4.1.1. Mean velocity

Figure 4 shows the mean fluid velocity profiles, with the inset showing the close-up of the velocities in the immediate vicinity of the wall. The mean velocity profile becomes flatter as $u_i/u_b$ increases from zero to 0.25. The mean velocity at $u_i/u_b=0.25$ becomes uniform for $y>0.3H$. As $u_i/u_b$ further increases, the mean velocity profile remains flat in the bulk region, but at $y$ around $0.2H$ the velocity decreases. We suspect that this decrease is related to the local peak of the particle concentration distribution, which could cause a larger local flow drag. Figure 6 shows a local peak of the particle concentration distribution at $y \approx 0.2H$ for $u_i/u_b=0.4$.

Figure 4. Mean fluid velocity profiles of turbulent channel flows for different $u_i/u_b$ at $a/H=0.05$, $\varPhi _0=2.36\,\%$, $Re_b=5746$ and $\rho _{r}=2.0$. The inset shows the close-up of the velocity profiles near the wall.

It was shown that the presence of neutrally buoyant finite-size particles generally leads to drag enhancement (Shao et al. Reference Shao, Wu and Yu2012; Picano et al. Reference Picano, Breugem and Brandt2015; Wang et al. Reference Wang, Peng, Guo and Yu2016). However, heavy particles in upward channel flow could cause a decrease in the wall friction at $u_i/u_b=0.159$ due to the suppression of large-scale vortices (Zhu et al. Reference Zhu, Yu, Pan and Shao2020a). From the inset of figure 4, the velocity gradient on the wall (i.e. flow friction) reaches a minimum at $u_i/u_b=0.2$, and it starts to increase as $u_i/u_b$ increases beyond 0.2. The reason for this increase will be discussed later.

Figure 5 depicts the profiles of the fluid relative mean velocity $\langle u _f\rangle$ with respect to the particle mean velocity $\langle u _p\rangle$ (i.e. mean slip velocity), where the angle brackets denote the phase averaging. It was observed experimentally and numerically that the mean velocity of heavy particles was smaller than the fluid counterpart in the bulk region, but larger than the fluid counterpart near the wall (Tsuji et al. Reference Tsuji, Morikawa and Shiomi1984; Shokri et al. Reference Shokri, Ghaemi, Nobes and Sanders2017; Zhu et al. Reference Zhu, Yu, Pan and Shao2020a). The reason for the larger particle mean velocity near the wall was attributed to the facts that the particles can slip on the wall and are preferentially located in the high-speed streaks (Tsuji et al. Reference Tsuji, Morikawa and Shiomi1984; Zhu et al. Reference Zhu, Yu, Pan and Shao2020a). Figure 5 shows negative slip velocities in the near-wall region at $u_i/u_b=0.159$ and 0.2. For $u_i /u_b \ge 0.25$, there are very few particles in the region of negative slip velocity due to the migration of the particles away from the wall (see figure 6); thus the expected negative slip velocities are not plotted. The position of zero mean slip velocity shifts towards the wall, as the settling coefficient increases. Two facts might be responsible for this observation. Firstly, the preferential distribution of the heavy particles in the high-speed streaks is caused by the entrainment of the large-scale streamwise vortices from the bulk region to the near-wall region under the condition of a higher particle concentration closer to the channel centre, as illustrated by Zhu et al. (Reference Zhu, Yu, Pan and Shao2020a) for $u_i/u_b=0.159$. At $u_i /u_b \ge 0.25$, shear-induced large-scale vortices are significantly suppressed (see figure 10) and the particles are distributed largely uniformly in the bulk region. Thus the preferential distribution of the heavy particles in the high-speed streaks only occurs in the near-wall region for high settling coefficients. Secondly, a larger settling coefficient itself means a larger positive slip velocity, and it also means a larger particle Reynolds number and thereby larger particle inertia. Consequently, at a larger $u_i/u_b$, particles with upward velocity smaller than the mean fluid velocity can approach closer to the wall.

Figure 5. Differences between the fluid and solid mean velocities at $a/H=0.05$, $\varPhi _0=2.36\,\%$, $Re_b=5746$ and $\rho _{r}=2.0$.

Figure 6. Profiles of the particle volume fraction at $a/H=0.05$, $\varPhi _0=2.36\,\%$, $Re_b=5746$ and $\rho _{r}=2.0$. The numerical result of Uhlmann (Reference Uhlmann2008) for $a/H=0.025$, $\varPhi _0= 0.42\,\%$ and $Re_p \approx 136$ is shown for comparison.

The mean particle volume fraction, $\varPhi _s/\varPhi _0$, as a function of the wall-normal distance $y/H$ is depicted in figure 6. The particle migration towards the channel (or pipe) centre was observed in numerical simulations (Zhu et al. Reference Zhu, Yu, Pan and Shao2020a,) and in experiments (Goes Oliveira, van der Geld & Kuerten Reference Goes Oliveira, van der Geld and Kuerten2017; Shokri et al. Reference Shokri, Ghaemi, Nobes and Sanders2017) when the particle settling effect was not strong, as a result of the Saffman effect (Saffman Reference Saffman1965; Auton Reference Auton1987; Auton, Hunt & Prud'Homme Reference Auton, Hunt and Prud'Homme1988). Costa et al. (Reference Costa, Brandt and Picano2020) showed that there exists a strong shear-induced lift force acting on the low-inertia particles in the near-wall region in turbulent channel flow and this lift force was well captured by the Saffman model, even though there was no gravity effect. As mentioned earlier, at $u_i /u_b \ge 0.25$ the particle concentration distribution is almost uniform in the bulk region, and a peak appears in the near-wall region. In fact, the peak position for $0.25 \le u_i /u_b \le 0.4$ is not so close to the wall, and it is around $y \approx 0.3H$ at $u_i/u_b = 0.25$. In the near-wall region, a high shear rate always exists, and a larger slip velocity causes a larger Saffman lift force, which leads to a very small particle volume fraction near the wall and the aforementioned peak. The peak position shifts closer to the wall as $u_i/u_b$ increases. The reason should be related to larger particle inertia at a higher settling coefficient.

It is interesting that the peak position is well correlated with the trough (local minimum) of the particle r.m.s. velocity component in the wall-normal direction, from the comparison between figure 6 and figure 9(b). Hence, the peak of the particle concentration distribution may also be related to turbophoresis, which drives the particles from the region of high turbulent intensity to that of low turbulent intensity (Caporaloni et al. Reference Caporaloni, Tampieri, Trombetti and Vittori1975; Reeks Reference Reeks1983). The mean momentum equation in the wall-normal (i.e. $y$) direction for the particle phase at the statistically steady state from the spatial averaging theorem (Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011) can be written as

(4.1)\begin{equation} {{\partial (\varPhi _s \langle {\sigma _s } \rangle_{yy} )} }/{{\partial y}} - {{\partial (\varPhi _s \rho _s \langle {v'_s v'_s } \rangle )} }/{{\partial y}} + \bar F_y = 0, \end{equation}

where $\sigma _s$ denotes the particle inner (Cauchy) stress and $\bar F_y$ is the volume-averaged hydrodynamic force on the particles in the $y$ direction. The particle r.m.s. velocity in the wall-normal direction (the second term in the above equation) plays a role in the force balance, similar to the fluid Reynolds stress for the fluid phase. It is also possible that the peak of the particle concentration results in the occurrence of the local minimum of the particle r.m.s. velocity.

The particle concentration distribution obtained by Uhlmann (Reference Uhlmann2008) for $a/H=0.025$ and $Re_p \approx 136$ is plotted in figure 6 for comparison. The two results for high slip velocities are qualitatively similar. The peak position of Uhlmann (Reference Uhlmann2008) is closer to the wall, because his particle size is smaller than ours. As shown in figure 19, the peak position is closer to the wall for smaller particles.

4.1.2. Fluctuating velocity and turbulent kinetic energy

The r.m.s. values of the fluid velocity fluctuations and the Reynolds shear stress for different settling coefficients are depicted in figure 7. The r.m.s. velocities do not change monotonically as $u_i/u_b$ increases from zero to 0.45. The addition of the particles decreases all r.m.s. velocity components across the channel at $u_i/u_b=0.159$. At $u_i/u_b \ge 0.2$, the streamwise r.m.s. velocity at the channel centre is enhanced, whereas the wall-normal and spanwise r.m.s. velocities are attenuated across the channel at $u_i/u_b \le 0.4$, and begin to exceed the single-phase values at the channel centre for $u_i/u_b \ge 0.45$. There is a local minimum in the profile of each r.m.s. velocity component for $0.25 \le u_i /u_b \le 0.4$ in figure 7, and we attribute the reason to the effect of the peak of the particle concentration distribution in figure 6. The Reynolds stress profiles for the particle-free and particle-laden cases are compared in figure 7($d$). The Reynolds stress decreases as $u_i/u_b$ increases up to around 0.25, and then it increases with increasing $u_i/u_b$, although it is still much smaller than the single-phase value at $u_i/u_b=0.45$.

Figure 7. Profiles of the r.m.s. velocity in single-phase and particle-laden turbulent channel flows: (a) streamwise, (b) wall-normal, (c) spanwise and (d) Reynolds shear stress $-\langle u'v'\rangle$ at $a/H=0.05$, $\varPhi _0=2.36\,\%$, $Re_b=5746$ and $\rho _{r}=2.0$.

Figure 8 shows the profiles of TKE and its relative change with respect to the single-phase flow. Turbulence suppression across the entire channel is observed for ${u_i}/{u_b} \le 0.25$. For ${u_i}/{u_b} > 0.25$, TKE augmentation in the bulk region and attenuation in the near-wall region can be observed in figure 8, as also observed in previous experiments (Tsuji et al. Reference Tsuji, Morikawa and Shiomi1984; Hosokawa & Tomiyama Reference Hosokawa and Tomiyama2004; Mena & Curtis Reference Mena and Curtis2020). As $u_i/u_b$ increases, the region for turbulence enhancement is expanded. At $u_i/u_b=0.5$, TKE is enhanced across the entire channel (not shown in figure 8 for convenience of presentation). From table 4, the relative increase in the total TKE at $u_i/u_b=0.5$ reaches $173.7\,\%$.

Figure 8. Profiles of (a) the TKE in single-phase and particle-laden turbulent channel flows and (b) the relative change of the TKE with respect to the single-phase flow, i.e. $(k-k_{{sp}})/k_{{sp}}$, at $a/H=0.05$, $\varPhi _0=2.36\,\%$, $Re_b=5746$ and $\rho _{r}=2.0$.

Table 4. The average fluid-phase r.m.s. velocity components, turbulence intensity $\tilde I$ and TKE $\tilde k$ in the entire channel. The values in parentheses represent the relative differences with respect to the single-phase flow. Here $Re_p$ is the particle Reynolds number based on mean inter-phase slip velocity at the channel centre at $a/H=0.05$, $\varPhi _0=2.36\,\%$, $Re_b=5746$ and $\rho _{r}=2.0$.

In table 4, the average fluid-phase r.m.s. velocity components, turbulence intensity $\tilde I$ and TKE $\tilde k$ in the entire channel for $u_i/u_b\le 0.5$ are presented. We see that the average r.m.s. velocities in all three directions, the mean turbulence intensity and the mean TKE are attenuated at $u_i/u_b \le 0.3$, with the maximum attenuation at $u_i/u_b=0.25$. The mean TKE is enhanced for $u_i/u_b \ge 0.4$, corresponding to $Re_p \ge 96.6$. At $u_i/u_b=0.4$ and 0.5, the increase in the mean TKE is caused by the increase in the streamwise r.m.s. velocity, since the wall-normal and spanwise r.m.s. velocities are reduced, compared to the single-phase flow. As $u_i/u_b$ exceeds 0.5, all components of the average r.m.s. velocities are augmented by the particles.

The solid-phase r.m.s. velocities and kinematic Reynolds shear stresses ($-\langle u'_pv'_p\rangle$ without the density) are plotted in figure 9. The streamwise solid-phase r.m.s. velocities are smaller than those of the fluid in the unladen case for $u_i/u_b < 0.35$, and the wall-normal r.m.s. velocities are larger than those of the fluid in the near-wall region due to the collision between the particles and the wall. The solid-phase r.m.s. velocities generally decrease with increasing $u_i/u_b$ for small $u_i/u_b$, and then increase with increasing $u_i/u_b$, consistent with the results of the fluid r.m.s. velocities in figure 7. As mentioned earlier, the trough (local minimum) in the profile of the wall-normal r.m.s. velocity is correlated with the peak of particle volume fraction in figure 6.

Figure 9. Profiles of the r.m.s. velocity and Reynolds shear stress for the particle phase: (a) streamwise, (b) wall-normal, (c) spanwise and (d) the particle kinematic Reynolds shear stress $-\langle u'_pv'_p\rangle$ at $a/H=0.05$, $\varPhi _0=2.36\,\%$, $Re_b=5746$ and $\rho _{r}=2.0$.

4.1.3. Vortex structure

Figure 10 shows the typical vortex structures in the half channel, identified with the $Q$ criterion ($Q=1.5$), which is defined as

(4.2)\begin{equation} Q=\tfrac{1}{2}(\varOmega_{i j} \varOmega_{i j}-\textit{\textsf{S}}_{i j} \textit{\textsf{S}}_{i j}), \end{equation}

where $\varOmega _{i j}=(u_{i,j}-u_{j,i}) / 2$ and $\textit {\textsf {S}}_{i j}=(u_{i,j} +u_{j,i}) / 2$ are the fluid rotation-rate tensor and the strain-rate tensor, respectively. Zhu et al. (Reference Zhu, Yu, Pan and Shao2020a) showed that the large-scale vortices were weakened more significantly at $u_i/u_b=0.159$ than at $u_i/u_b=0.092$. From figure 10(b), at $u_i/u_b=0.25$ with the particle Reynolds number being approximately 50, the large-scale vortices appear to be completely suppressed. The suppression of the large-scale vortices is clearly responsible for the pronounced attenuation of the turbulence intensity at relatively small settling coefficients, as observed in figures 7 and 8, and table 4.

Figure 10. Vortex structures of the (a) single-phase flow and (b,c) the particle-laden flows for (b) $u_i/u_b=0.25$ and (c) $u_i/u_b=0.45$ at $a/H=0.05$, $\varPhi _0=2.36\,\%$, $Re_b=5746$ and $\rho _{r}=2.0$. The colour of the vortices represents the fluid streamwise velocity, with the same scale for all three cases.

The reason for the suppression of the large-scale vortices is related to the effects of the drag force, the enhanced viscous dissipation, or the particle-induced wake structures. These three factors are interrelated, since a higher settling coefficient implies a larger drag force, a higher viscous dissipation and stronger particle wakes. It is not clear which effect is more important. Significant suppression of the large-scale vortices was observed in the numerical simulations based on the point-particle model (Vreman Reference Vreman2015; Muramulla et al. Reference Muramulla, Tyagi, Goswami and Kumaran2020). Vreman (Reference Vreman2015) demonstrated that the non-uniform part of the mean feedback drag force contributed significantly to particle-induced turbulence attenuation. Xia, Yu & Deng (Reference Xia, Yu and Deng2019) observed that, at the same settling coefficient, light (bubble) particles attenuated the turbulence in a downward channel flow, where the particles migrated towards the channel centre, whereas light particles enhanced the turbulence in an upward channel, where the particles migrated towards the channel wall. From these observations, it seems that the particle feedback force distributed in the bulk region suppresses the self-generation of large-scale vortices in the channel flow. On the other hand, the particle wakes tend to enhance the turbulence intensity when they are sufficiently strong. Figure 10 shows strong particle wakes at $u_i/u_b=0.45$, corresponding to $Re_p=113.5$, in which case the intensity of the particle-induced turbulence becomes stronger than that of large-scale-vortex-induced single-phase turbulence in the bulk region, as shown in figure 8. The stronger particle wakes at a higher $u_i/u_b$ in the near-wall region where the velocity gradient is large are expected to be responsible for the increase in the fluid Reynolds stress with increasing $u_i/u_b$ at $u_i/u_b>0.25$ in figure 7($d$).

4.2.1. Effects of other parameters on TKE

We have shown that there is a transition from turbulence attenuation to augmentation as the settling coefficient increases for the typical case of $a/H=0.05$, $\varPhi _0=2.36\,\%$, $Re_b=5746$ and $\rho _{r}=2.0$ (i.e. group 5). The TKE of the particle-laden flow starts to exceed that of the single-phase flow at the channel centre at $u_i/u_b=0.3$($Re_p=64$) (figure 8), while the augmentation in the total TKE starts at $u_i/u_b=0.45$ ($Re_p=96.6$) (table 4). In the following, the effects of the bulk Reynolds number, the particle size, the particle volume fraction and the particle–fluid density ratio on the TKE and the critical particle Reynolds number are examined. We consider two bulk Reynolds numbers $Re_b=5746$ and 12 000, four particle sizes $a/H=0.05$, 0.75, 0.1 and 0.15, two particle volume fractions $\varPhi _0=0.84\,\%$ and $2.36\,\%$, and three density ratios 2, 10 and 100, to inspect the effects of the specific parameter, while keeping the other parameters constant, as shown in table 2. The maximum density ratio is 100, because at this density ratio the particle inertia is already considerably strong for the particle size $a/H=0.1$ (Yu et al. Reference Yu, Lin, Shao and Wang2017). In addition, to investigate the effect of the particle size at the same particle number density, the case of $a/H=0.05$ and $\varPhi _0=0.3\,\%$ (group 7) is added for comparison to the case of $a/H=0.1$ and $\varPhi _0=2.36\,\%$ (Group 1) at the same particle number of 360 in the channel cell.

Figure 11 shows the relative changes of the mean turbulence intensity $\widetilde {{CT}}$ in the entire channel and the turbulence intensity $CT$ at the channel centre as a function of $Re_p$. When the other parameters are kept constant, both critical particle Reynolds numbers increase with increasing particle size $a/H$, bulk Reynolds number $Re_b$ and particle–fluid density ratio $\rho _r$. Although both critical particle Reynolds numbers generally increase with decreasing particle volume fraction $\varPhi _0$, the dependence of the critical $Re_p$ on the particle volume fraction for the turbulence intensity at the channel centre is more pronounced than for the total turbulence intensity in the channel. From figure 11(a), the critical $Re_p$ for the total turbulence intensity is relatively insensitive to the particle volume fraction: its value for $\varPhi _0=0.84\,\%$ is slightly larger than that for $\varPhi _0=2.36\,\%$ at $a/H=0.1$, and almost the same for $\varPhi _0=0.3\,\%$ and $\varPhi _0=2.36\,\%$ at $a/H=0.05$. Nevertheless, the effect of the particle volume fraction on the TKE is significant: a higher particle volume fraction causes more enhancement in the augmentation regime or more reduction in the attenuation regime. In other words, the magnitude of the turbulence enhancement or attenuation increases with increasing particle volume fraction.

Figure 11. The relative change in (a) the mean turbulent intensity $\widetilde {{CT}}$ and (b) the local turbulent intensity at the channel centre $CT$ as a function of $Re_p$.

At a critical particle Reynolds number for the total TKE in the channel, the TKE is enhanced in the bulk region and attenuated in the near-wall region. Figure 12($a$) compares the TKE profiles for different particle volume fractions at comparable $Re_p$ near the critical value. One can see that the increase in the particle volume fraction leads to turbulence enhancement in the bulk region, but it leads to turbulence attenuation in the near-wall region. The amounts of TKE enhancement and attenuation are comparable. Therefore, the total TKE is not sensitive to the variation of the particle volume fraction near the critical state, resulting in the critical $Re_p$ for the total TKE being insensitive to the particle volume fraction. Since the increase in the particle volume fraction leads to the increase in the TKE at the channel centre near the critical state, the critical $Re_p$ for the turbulence augmentation at the channel centre decreases with increasing particle volume fraction.

Figure 12. The TKE profiles of single-phase and particle-laden turbulent channel flows for (a) different particle volume fractions and (b) different particle–fluid density ratios.

It can be observed in figure 12($b$) that, at comparable values of $Re_p$, the flow laden with particles with a larger density (or inertia) has a lower turbulence intensity in the bulk region, which is responsible for the larger critical particle Reynolds number at a larger density ratio.

Figure 11 indicates that the effect of the particle size on the critical particle Reynolds number is most pronounced, compared to those of the other parameters. Some TKE profiles for different particle sizes at the same volume fraction of $\varPhi _0=2.36\,\%$ are shown in figure 13($a$). The TKE in the bulk region decreases significantly with increasing particle size at comparable $Re_p$. Some TKE profiles at the same particle number density (360 particles in the channel domain) for $a/H=0.05$ and $a/H=0.1$ are depicted in figure 13($b$). One can see that the TKEs for $a/H=0.05$ and $\varPhi _0=0.3\,\%$ are still larger than those for $a/H=0.1$ and $\varPhi _0=2.36\,\%$ at comparable $Re_p$. We will explain this particle size effect later.

Figure 13. The TKE profiles of single-phase and particle-laden turbulent channel flows for different particle sizes at (a) the same particle volume fraction of $2.36\,\%$ and (b) the same particle number of $N_p=360$.

The TKE profiles of single-phase and particle-laden turbulent channel flows for $Re_b=5746$ and 12 000 at $a/H=0.05$ and $\varPhi _0=2.36\,\%$ are compared in figure 14. For $Re_b=5746$, the TKE of the particle-laden flow in the near-wall region is considerably attenuated and the large-scale vortices existing in the single-phase flow are completely suppressed at $Re_p=64$ (i.e. $u_i/u_b=0.3$) (see figure 10), and there is appreciable TKE augmentation in the bulk region at $Re_p=96.6$. By contrast, for $Re_b=12\,000$, the TKE of the particle-laden flow in the near-wall region is close to that of the single-phase flow and TKE enhancement takes place only near the channel centre at $Re_p=102.1$. This difference indicates that stronger particle wakes (or drag force) are needed to prevail over the large-scale vortices and result in the turbulence augmentation for a higher channel (or bulk) Reynolds number. This is understandable, considering that the large-scale vortices are stronger (despite being smaller in size) at a higher channel Reynolds number. The direct effect of the particle wakes on the TKE should be related to the particle-induced TKE production rate, and we will show that the particle-induced TKE production rate normalized with the bulk velocity decreases significantly with increasing $Re_b$ at the same particle Reynolds number from (4.14). Our effect of $Re_b$ is consistent with the experimental results of Mena & Curtis (Reference Mena and Curtis2020) on an upward pipe flow, who observed the change from augmentation in the total TKE at $Re_b=2 \times 10^5$ to attenuation at $Re_b=3.5 \times 10^5$.

Figure 14. The TKE profiles of single-phase and particle-laden turbulent channel flows for different bulk Reynolds numbers at $a/H=0.05$ and $\varPhi _0=2.36\,\%$.

4.2.2. Criterion for turbulence modulation

Based on our simulation data in figure 11, we propose a new criterion to distinguish between turbulence augmentation and attenuation in terms of the total TKE and the TKE at the channel centre, respectively. The parameter for the total TKE is

(4.3)\begin{equation} {\chi_{1}} = \frac{{R{e_p}}}{{Re_b^{\textrm{{0}}\textrm{{.33}}}{{( {{d_p / H}} )}^{\textrm{{0.61}}}}\rho _r^{0.05}}}. \end{equation}

We have shown that the particle volume fraction has little effect on the critical Reynolds number for the total TKE, hence it is not included in the new parameter $\chi _1$. The relative change of the mean turbulence intensity $\widetilde {CT}$ as a function of $\chi _1$ is depicted in figure 15($a$). It is observed that the attenuation of the total TKE occurs at $\chi _1<20$, and the augmentation takes place at $\chi _1>20$. We collect experimental data on upward pipe flows in the literature and the numerical data of Uhlmann (Reference Uhlmann2008) on upward channel flow, as summarized in table 5, and compare them with our criterion of $\chi _1=20$ in figure 15($b$). One can see that our criterion agrees well with the experimental data and the numerical data of Uhlmann (Reference Uhlmann2008). There is some discrepancy between our criterion and the recent experiment of Mena & Curtis (Reference Mena and Curtis2020). At the moment, the reason is not clear. The pipe Reynolds number of Mena & Curtis (Reference Mena and Curtis2020) ranges from $Re_b=2 \times 10^5$ to $Re_b=3.5 \times 10^5$, and is much higher than ours. Further studies on the particle–turbulence interactions at high bulk Reynolds numbers are needed.

Figure 15. Relative changes of the mean turbulence intensity in the channel $\widetilde {CT}$ as a function of ${\chi _1}$ for (a) the present numerical simulations and (b) previous experimental and numerical data.

Table 5. Summary of the previous experiments on upward pipe flow and the simulation of Uhlmann (Reference Uhlmann2008) on upward channel flow. Here $D$ is the diameter of the pipe and $D$ equals $2H$ for the channel flow.

The parameter of the criterion for the turbulence intensity at the channel centre is

(4.4)\begin{equation} {\chi _\textrm{{2}}} = \frac{{\varPhi_0^{0.1}}}{{Re_b^{0.53} {{( {{d_p / H}} )}^{\textrm{{0}}\textrm{{.61}}}}\rho _r^{0.065}}}R{e_p}. \end{equation}

Figure 16($a$) shows that a value of $\chi _2$ of approximately 1.55 can classify the turbulence augmentation and attenuation in our simulations. From figure 16($b$), our criterion of $\chi _2 = 1.55$ for the turbulence intensity at the channel centre also agrees well with the previous experiments.

Figure 16. Relative changes of the turbulence intensity at the channel centre $CT$ as a function of ${\chi _2}$ for (a) the present numerical simulations and (b) previous experimental and numerical data.

4.3.1. Turbulent kinetic energy budget

The volume-averaged fluid TKE equation for particle-laden flows can be derived as (Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011; Ma et al. Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017; Vreman & Kuerten Reference Vreman and Kuerten2018; Peng et al. Reference Peng, Ayala and Wang2019)

(4.5)\begin{align} &\frac{\partial(\varPhi \rho k)}{\partial t}+ \frac{\partial(\varPhi \rho\langle u_{j}\rangle k)}{\partial x_{j}} \nonumber\\ &\quad=\underbrace{-\varPhi \rho\langle u_{i}^{\prime} u_{j}^{\prime}\rangle \frac{\partial\langle u_{i}\rangle}{\partial x_{j}}}_{\mathcal{P}} \underbrace{-\varPhi\left\langle\tau_{i j}^{\prime} \frac{\partial u_{i}^{\prime}}{\partial x_{j}}\right\rangle}_\epsilon \nonumber\\ & \qquad \underbrace{-\frac{\partial\left(\varPhi \rho\left\langle u_{j}^{\prime} \dfrac{u_{i}^{\prime} u_{i}^{\prime}}{2}\right\rangle\right.}{\partial x_{j}}- \frac{\partial(\varPhi\langle u_{j}^{\prime} p^{\prime}\rangle)}{\partial x_{j}}+\frac{\partial(\varPhi\langle u_{i}^{\prime} \tau_{i j}^{\prime} \rangle)}{\partial x_{j}}}_\mathcal{D} \nonumber\\ & \qquad +\underbrace{\varPhi(-\langle p\rangle \delta_{i j}+\langle\tau_{i j}\rangle)\left(\frac{\partial\langle u_{i}\rangle}{\partial x_{i}}-\left\langle\frac{\partial u_{i}}{\partial x_{j}}\right\rangle\right)}_{\mathcal{I} {\rm 1}} \nonumber\\ & \qquad \underbrace{-\frac{1}{V} \int_{S_{I}} n_{j} u_{i}({-}p \delta_{i j}+\tau_{i j}) \, \textrm{d} s+\frac{\langle u_{i}\rangle}{V} \int_{S_{I}}({-}p \delta_{i j}+\tau_{i j}) n_{j}\, \textrm{d} s}_{\mathcal{I} {\rm 2}}. \end{align}

Here, $\varPhi$, $\rho$ and $\tau _{i j}$ represent the fluid volume fraction, density and viscous stress, respectively. For convenience, the subscript ‘$f$’ denoting the fluid phase is omitted. The subscripts $i$ and $j$ denote the coordinate component, $V$ represents the volume for averaging, $S_I$ is the interface between the fluid and solid inside the volume $V$, and $n$ is the unit vector on the particle surface pointing to the fluid. The prime in (4.5) and other places below indicates the fluctuating part of the quantity.

For the channel flow at the statistically steady state, the statistics are homogeneous in the streamwise and spanwise directions and do not change with time. Thus the derivative is only non-zero in the wall-normal (i.e. $y$) direction, and the averaging volume $V$ is a thin band stretching over the streamwise and spanwise directions. In (4.5), ${\mathcal {P}}$ and $\epsilon$ represent the production due to the shear flow and the viscous dissipation of TKE; $\mathcal {D}$ denotes the diffusion terms including the turbulent and viscous diffusion terms. The terms ${\mathcal {P}}$, $\epsilon$ and $\mathcal {D}$ exist in the single-phase flow. The terms ${\mathcal {I}}1$ and ${\mathcal {I}}2$ exist only for the multiphase flows due to the presence of the interface between the two phases, and thus are called the interfacial terms or the particle source terms in the TKE equation (Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011; Ma et al. Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017). We will show that the interfacial term can be approximated with two terms: the product of the mean inter-phase drag force and the mean slip velocity, and the correlation between the particle fluctuating velocity and the drag force on the particle.

Let $\sigma _{ij}$ denote the total stress $\sigma _{ij} = - p\delta _{ij} + \tau _{ij}$. From the spatial averaging theorem (Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011), one can derive

(4.6)\begin{align} \varPhi \left\langle {\frac{{\partial u_i } }{{\partial x_j }}} \right\rangle & = \frac{{\partial (\varPhi \langle {u_i } \rangle )} }{{\partial x_j }} - \frac{1 }{V}\int_{S_I } {n_j } u_i\, \textrm{d} s \nonumber\\ & = \varPhi \frac{{\partial \langle {u_i } \rangle } }{{\partial x_j }} + \frac{1 }{V}\int_{S_I } {n_j } \langle {u_i } \rangle \, \textrm{d} s - \frac{1 }{V}\int_{S_I } {n_j } u_i \, \textrm{d} s \nonumber\\ & = \varPhi \frac{{\partial \langle {u_i } \rangle } }{{\partial x_j }} - \frac{1 }{V}\int_{S_I } {n_j } u'_i \, \textrm{d} s. \end{align}

Then, the interfacial term ${\mathcal {I}}1$ can be written as

(4.7)\begin{equation} {\mathcal{I}}1 = \varPhi ( { - \langle \,p\rangle \delta _{ij} + \langle {\tau _{ij} } \rangle } )\left( {\frac{{\partial \langle {u_i } \rangle } }{{\partial x_i }} - \left\langle {\frac{{\partial u_i } }{{\partial x_j }}} \right\rangle } \right) = \frac{1 }{V}\int_{S_I } {\langle {\sigma _{ij} } \rangle n_j } u'_i \, \textrm{d} s. \end{equation}

The interfacial term ${\mathcal {I}}2$ can be written as

(4.8)\begin{equation} {\mathcal{I}}2 ={-} \frac{1 }{V}\int_{S_I } {\sigma _{ij} n_j } u'_i \, \textrm{d} s. \end{equation}

Thus, the total interfacial term has the form (Ma et al. Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017; Vreman & Kuerten Reference Vreman and Kuerten2018)

(4.9)\begin{equation} {\mathcal{I}} = {\mathcal{I}}\textrm{{1}} + {\mathcal{I}}\textrm{{2}} ={-} \frac{1}{V}\int_{S_I } {\sigma '_{ij} n_j } u'_i \, \textrm{d} s. \end{equation}

Accurate computation of the interfacial term via the integration of the stress on the particle surface is difficult for the fictitious domain method, and here we obtain it from the TKE equation at the statistically steady state: ${\mathcal {I}}=-({\mathcal {P}}+\epsilon + \mathcal {D})$. Since the terms $\mathcal {P}$, $\epsilon$ and $\mathcal {D}$ are computed by averaging over time, and over the streamwise and spanwise directions (not over the wall-normal direction), the averaging volume can be regarded as a very thin band containing the plane $y=y_i$, where $y=y_i$ is the $y$ coordinate of the computational grid.

The profiles of the production, dissipation, diffusion and interfacial terms in the TKE equation for the case of $a/H=0.05$, $\varPhi _0=2.36\,\%$, $\rho _r=2.0$ and $Re_b=5746$ are depicted in figure 17. Similar to the behaviour of the Reynolds stress and the r.m.s. velocities in figure 7, the production rate of TKE ${\mathcal {P}}$ first decreases and then increases, as $u_i/u_b$ increases; ${\mathcal {P}}$ is attenuated by the particles for $u_i/u_b \le 0.45$, similar to the Reynolds stress in figure 7(b). In the bulk region, the particle wakes do not contribute significantly to the Reynolds stress and ${\mathcal {P}}$, since the mean fluid velocity there is nearly uniform at a high $u_i/u_b$. The reduction in ${\mathcal {P}}$ is clearly caused by the suppression of the large-scale vortices. For the present particle-laden flow where there is a mean interphase slip velocity, the interfacial term is always larger than zero and increases with increasing $u_i/u_b$, as shown in figure 17(c). Hence, the interfacial term can be regarded as the particle-induced TKE production. At a high $u_i/u_b$, the strong particle wakes cause a large TKE and accordingly a large viscous dissipation $\epsilon$. The large drag force at a high $u_i/u_b$ also leads to a large interfacial term ${\mathcal {I}}$. In the bulk region, the shear-induced production ${\mathcal {P}}$ and the diffusion $\mathcal {D}$ are small; therefore, the dissipation $\epsilon$ is nearly equal to the particle-induced production ${\mathcal {I}}$ in magnitude.

Figure 17. Profiles of the terms in the TKE equation for the typical case of $a/H=0.05$, $\varPhi _\textrm {{0}}=2.36\,\%$ and $Re_b=5746$: (a) the production term ${\mathcal {P}}$, (b) the dissipation term $\epsilon$, (c) the interfacial term ${\mathcal {I}}$ and (d) the diffusion term ${\mathcal {D}}$. All terms are normalized with ${{\rho _f}u_b^3}/H$.

The production, dissipation and interfacial terms for different particle sizes at the same particle number density and comparable particle Reynolds numbers are compared in figure 18($a$). It has been shown in figure 13($b$) that smaller particles can cause more augmentation in the TKE at the same particle number density and comparable particle Reynolds numbers. Figure 18($a$) shows the same behaviour for the interfacial and dissipation terms: they are larger for smaller particles at the same particle number density and Reynolds number. It should be noted that the interfacial term is larger for smaller particles ($a/H=0.05$) in the bulk region even if the local particle number density there is significantly smaller for smaller particles (because the particle concentration is higher in the near-wall region for smaller particles), as shown in figure 19. Figure 18($b$) compares the production, dissipation and interfacial terms for $\rho _r=2$ and 100 at $a/H=0.1$. Overall, the interfacial term is greater for the less dense particles with lower inertia. The reason should be related to the particle concentration distribution. For $\rho _r=2$ ($Re_p=163.8$, $a/H=0.1$), there is a peak at around $y=0.4H$ mainly as a result of the Saffman effect, which drives the particles away from the wall, whereas for $\rho _r=100$, the particle concentration distribution is more homogeneous across the channel. The higher particle concentration at around $y=0.4H$ for $\rho _r=2$ ($a/H=0.1$) in figure 19 should be responsible for its larger interfacial term at around $y=0.4H$ in figure 18($b$).

Figure 18. Production, dissipation and interfacial terms in the TKE equation, normalized with ${{\rho _f}u_b^3/H}$, for (a) different particle sizes $a/H$ at the same particle number density and (b) different particle–fluid density ratios $\rho _r$.

Figure 19. Profiles of the particle volume fraction for different particle sizes at the same particle number density and different particle–fluid density ratios.

The reason for the general decrease in the critical Reynolds number (particularly the one for the TKE at the channel centre) with decreasing particle size and density ratio and with increasing particle volume fraction should be related to the larger interfacial term (i.e. particle-induced TKE production) for a smaller particle size or density ratio, or a larger particle volume fraction at the same particle Reynolds number. However, it should be noted that, although the interfacial term ${\mathcal {I}}$ or the dissipation rate $\epsilon$ is generally correlated with the TKE, a higher ${\mathcal {I}}$ or $\epsilon$ does not always mean a higher TKE. For example, the TKE at the channel centre for $\rho _r=2$ and $Re_p=163.8$ is higher than that for $\rho _r=100$ and $Re_p=169.9$, as shown in figure 12($b$), but ${\mathcal {I}}$ and $\epsilon$ at the channel centre are almost the same for the two cases, as shown in 18($b$). Hence, the dependence of the interfacial term on the parameters can be regarded as the main reason for the dependence of the critical particle Reynolds number on the corresponding parameters, but it is not the only reason. As for the effect of the density ratio, another reason may be that less dense particles have stronger hydrodynamic interactions and consequently may have larger variations in the instantaneous slip velocity at the channel centre where the flow is dominated by the particle wakes, when the mean interphase slip velocity is the same.

4.3.2. Analysis of the interfacial term in the dilute limit

We here attempt to provide an explanation for the effects of the particle size and the bulk Reynolds number on the TKE at the same particle Reynolds number by simplifying the interfacial term in the dilute limit. The interfacial term ${\mathcal {I}}1$ is not important, since the average of the velocity gradient should be close to the gradient of the average velocity, as shown in the calculation of Peng et al. (Reference Peng, Ayala and Wang2019). Following the analysis of Crowe et al. (Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011), we neglect the effects of the particle rotation. We further assume that all particle surfaces are inside the volume $V$, which is reasonable for a homogeneous system such as particle sedimentation in the bulk region. The fluid velocity on the surface for a non-rotating particle $k$ is equal to the particle translational velocity, which is a constant for the integration on the particle surface. Then, the interfacial term ${\mathcal {I}}2$ can be simplified as

(4.10)\begin{equation} - \frac{1}{V}\int_{{S_I}} {{u_i}( - p{\delta _{ij}} + {\tau _{ij}})} \, \textrm{d} s + \frac{{\langle {{u_i}} \rangle }}{V}\int_{{S_I}} {( - p{\delta _{ij}} + {\tau _{ij}}){n_j}} \, \textrm{d} s = \frac{1}{V}\sum_k {(\langle {{u_i}} \rangle - {u_{pki}}){F_{ki}}}, \end{equation}

where ${F_{ki}}$ is the drag force on the particle $k$ and $u_{pki}$ is the translational velocity of the particle $k$.

For vertical channel flow, the drag force is mainly directed in the streamwise direction and thus we neglect the contributions from the other two directions:

(4.11)\begin{equation} \frac{1}{V}\sum_k {(\langle {{u_i}} \rangle - {u_{pki}}){F_{ki}}} \approx \frac{1}{V}\sum_k {(\langle u_f \rangle - {u_{pk}})F_k}, \end{equation}

where $u_f$, $u_{pk}$ and $F_k$ represent the fluid velocity, the particle translational velocity and the drag force on the particle $k$ in the streamwise direction, respectively. Considering that ${\langle u_f \rangle - {u_{pk}}}={\langle {{u_f}} \rangle - \langle {{u_p}} \rangle - u{'_{pk}}}$ (with $u{'_{pk}}$ being the particle velocity fluctuation), we have

(4.12)\begin{equation} {\mathcal{I}} \approx \frac{1}{V}\sum_k {(\langle {{u_f}} \rangle - \langle {{u_p}} \rangle - u{'_{pk}}){F_k}} = (\langle {{u_f}} \rangle - \langle {{u_p}} \rangle )\frac{1}{V}\sum_k {{F_k} - \frac{1}{V}\sum_k {u{'_{pk}}{F_k}}}.\end{equation}

The first term on the right-hand side is the product of the mean inter-phase drag force and the mean slip velocity, and was regarded as the main contribution to the interfacial term (Troshko & Hassan Reference Troshko and Hassan2001; Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011).

For particle-laden flows in the dilute limit, one may assume that each particle settles with the same slip velocity $u_r=\langle {{u_f}} \rangle - \langle {{u_p}} \rangle$, and then

(4.13)\begin{equation} (\langle {{u_f}} \rangle - \langle {{u_p}} \rangle ) \frac{1}{V}\sum_k {{F_k}}= n\frac{{u_r^3{C_D}{\rm \pi} {a^2}{\rho _f}}}{2} = n\frac{Re_p^3{\nu ^3}{C_D}{\rm \pi} {\rho _f}}{{16a}},\end{equation}

in which $n = \varPhi _s /(4{\rm \pi} {a^3}/3)$ is the particle number density and $C_D$ is the standard drag coefficient. The interfacial term in our simulations is normalized with ${{\rho _f}u_b^3/H}$. From (4.13), one has

(4.14)\begin{equation} \frac{(\langle u_f \rangle - \langle u_p \rangle)\displaystyle\frac{1}{V} \sum\nolimits_k F_k} { \rho _f u_b^3/H} \sim n^*\frac{{Re_p^3{C_D}}}{Re_b^3(a/H)} \sim \varPhi_s\frac{{Re_p^3{C_D}}}{Re_b^3(a/H)^4},\end{equation}

where $n^*$ denotes the dimensionless particle number density. The drag coefficient $C_D$ is a function of the particle Reynolds number alone.

Hence, the above dimensionless interfacial term is inversely proportional to the particle size at the same particle number density $n$, $Re_p$ and $Re_b$, which provides an explanation for the early observations that the interface term is larger for smaller particles at the same particle Reynolds number and particle number density in figure 18($a$) and the critical particle Reynolds number increases significantly with increasing particle size in figure 11. Equation (4.14) also explains why the dimensionless TKE in the bulk region is significantly smaller for a higher bulk Reynolds number at comparable particle Reynolds numbers in figure 14 and why the critical particle Reynolds number is higher at a higher bulk Reynolds number in figure 11.

4.3.3. Modelling of the interfacial term

Modelling of the interfacial term is crucial for the Reynolds-averaged Navier–Stokes (RANS) model and the Reynolds stress model (i.e. second-order moment model) for the simulations of industrial-scale multiphase turbulent flows (Ansys 2013; Ma et al. Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017, Reference Ma, Lucas, Jakirlić and Fröhlich2020). The Troshko–Hassan model (Troshko & Hassan Reference Troshko and Hassan2001) is one of the widely used models for the interfacial term:

(4.15)\begin{equation} {\mathcal{I}} = C_n \bar F{u_r}, \end{equation}

where $\bar F$ is the volume-averaged drag force, $\bar F = ({1 }/{V})\sum _k {F_k }$, and $C_n$ is a coefficient. In the Ansys FLUENT software, a value of 0.75 is recommended for $C_n$ (Ansys 2013). Ma et al. (Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017) determined $C_n$ as a function of $Re_p$, from the interface-resolved DNS data of Santarelli, Roussel & Fröhlich (Reference Santarelli, Roussel and Fröhlich2016) on the bubbly flow in a vertical channel:

(4.16)\begin{equation} C_n = 0.18Re_p^{0.23}.\end{equation}

In the following, we attempt to evaluate $C_n$ from our DNS data. Following Ma et al. (Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017), the drag force is determined directly from the simulations, rather than from the drag model. For the channel flow at the statistically stationary state, the fluid mean momentum equation in the streamwise direction can be written as follows (Yu et al. Reference Yu, Lin, Shao and Wang2017):

(4.17)\begin{equation} \frac{\textrm{{d}}}{{\textrm{{d}}y}}\left( {{\varPhi _f}\mu \frac{{\textrm{{d}}\langle {{u_f}} \rangle }}{{\textrm{d} y}}} \right) + {\varPhi _f}\left( { - \frac{{\textrm{{d}}{p_e}}}{{\textrm{{d}}x}}} \right) + \frac{\textrm{{d}}}{{\textrm{{d}}y}}( {{\varPhi _f}{\rho _f}\langle { - {{u'}_f}{{v'}_f}} \rangle } ) - \bar F = 0. \end{equation}

For upward flow, one can compute the pressure gradient from the force balance,

(4.18)\begin{equation} -\textrm{d} p_{e}/\textrm{d} x=\tau_{w}/H+\varPhi_{0} (\rho_{s}-\rho_{f})g, \end{equation}

where $\tau _{w}$ is the mean shear stress on the sidewalls. Then, the volume-averaged drag force can be calculated from (4.17). The calculation of the drag force at the channel centre can be simplified, if the flow is homogeneous there (i.e. $\textrm {d}(\,\cdot \,)/{\textrm {d} y}=0$). The solid-phase mean momentum equation in the streamwise direction is (Yu et al. Reference Yu, Lin, Shao and Wang2017)

(4.19)\begin{equation} \frac{\textrm{d}}{{\textrm{d} y}}(\varPhi_s \langle {\sigma _s } \rangle _{xy}) + \varPhi _s \left( - \frac{{\textrm{d}p_e } }{{\textrm{d} x}}\right) - \varPhi _s (\rho _s - \rho _f )g + \frac{\textrm{d}}{{\textrm{d} y}}(\varPhi _s \rho_s \langle { - {{u'}_s}{{v'}_s}} \rangle ) + \bar F = 0. \end{equation}

For a homogeneous flow, the fluid momentum equation (4.17) reduces to $\bar F = \varPhi _f ( - {{\textrm {d}p_e } }/{{\textrm {d} x}})$, and the solid-phase momentum equation (4.19) reduces to $\bar F = \varPhi _s (\rho _s - \rho _f )g - \varPhi _s ( - {{\textrm {d}p_e } }/{{\textrm {d} x}})$. From these two equations, one can obtain the pressure gradient and the drag force as

(4.20)\begin{equation} \bar F = \varPhi_f \left( - \frac{{\textrm{d}p_e } }{{\textrm{d} x}}\right) = \varPhi _f \varPhi _s (\rho _s - \rho _f )g, \end{equation}

which was used by Ma et al. (Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017) for the calculation of the drag force. Equation (4.17) is used in the present study, because, for $\rho _r=100$, the gradient of the particle Reynolds stress at the channel centre in (4.19) is found to be not negligibly small. For $\rho _r=2$, the results obtained with the two methods are close to each other for relatively high particle Reynolds numbers.

The mean $C_n$ values for the channel centre region $- 0.1 \le y/H \le 0.1$ (here $y=0$ corresponding to the channel centre) as a function of $Re_p$ for different particle sizes, different volume fractions and different density ratios are plotted in figure 20($a$). When we ran the cases for the study of the criterion on the TKE augmentation and attenuation, the data for the TKE budget were not computed. Owing to our limited computer resources, only some cases were selected for the generation of the data for the modelling of the interfacial term. Figure 20($a$) shows that the coefficient $C_n$ decreases with increasing particle Reynolds number for the same group of control parameters, and is higher for a lower particle volume fraction or a larger particle size at the same particle Reynolds number. The effect of the density ratio is small. The fitting of our data yields

(4.21)\begin{equation} {C_n}=2.1 Re_p^{ - 0.09}{\left( {\frac{a}{H}} \right)^{0.12}}\rho _r^{0.007}{({1 - {\varPhi_s}} )^{8.0}}. \end{equation}

Figure 20($b$) shows that the above model of $C_n$ fits our data well.

Figure 20. Plots of (a) $C_n$ versus $Re_p$ for different particle sizes, particle volume fractions and density ratios, and (b) ${C_n}=2.1 Re_p^{-0.09}{({{a}/{H}})^{0.12}}\rho _r^{0.007}{({1 - {\varPhi _s}})^{8.0}}$ as a function of $Re_p$.

The dependence of our $C_n$ on $Re_p$ appears to be inconsistent with that of Ma et al. (Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017). Our $C_n$ (4.21) decreases with increasing $Re_p$, whereas theirs (4.16) increases with increasing $Re_p$. However, our dependence of $C_n$ on $Re_p$ is examined at the same particle size and particle volume fraction, whereas theirs is examined at different particle sizes and particle volume fractions. Their results are consistent with ours in that $C_n$ is larger for larger particles and a lower particle volume fraction. Thus, we suspect that the dependence of their $C_n$ on $Re_p$ is caused by the effects of the particle size and the particle volume fraction, rather than $Re_p$. In addition, we note that they used the data on bubbly flow at $Re_p>230$, whereas we consider solid particles at $Re_p<230$.

The reasons for the dependence of $C_n$ on the particle Reynolds number, the particle size, the particle volume fraction and the density ratio are not clear, but should be related to the effects of the second term on the right-hand side in (4.12), i.e. the particle velocity fluctuation term: $({1}/{V})\sum _k {u{'_{pk}}{F_k}}$. If the effect of the particle rotation is not important and (4.12) is a good approximation to the interfacial term, then $C_n \approx 1 - ({1 }/{V})\sum _k {u'_{pk} F_k /} (u_r \bar F)$. Particles may be accelerated by the particle wakes, resulting in the positive correlation between the particle fluctuating velocity and the force on the particle. An increase in the particle Reynolds number or the particle volume fraction, or a decrease in the particle size or the density ratio, may cause a larger particle velocity fluctuation and the increase in $({1 }/{V})\sum _k {u'_{pk} F_k /} (u_r \bar F)$, i.e. the decrease in $C_n$. A more significant effect of the density ratio on $C_n$ than that shown in figure 20($a$) was expected from this reasoning. In fact, if the drag force is calculated with (4.20), instead of (4.17), the difference in $C_n$ for $\rho _r=100$ and $\rho _r=2$ would be more pronounced. Using (4.17), $C_n$ can be determined from:

(4.22)\begin{equation} \frac{{\mathcal{I}} }{{\rho _f u_b^3 /H}} = \frac{{C_n \varPhi _f \varPhi _s (\rho _s - \rho _f )gu_r H} }{{\rho _f u_b^3 }} = C_n \varPhi _f \varPhi _s \frac{3}{{8(a/H)}}\left(\frac{{u_i } }{{u_b }}\right)^2 \frac{{u_r}}{{u_b}}. \end{equation}

Now we compare two cases: $(a/H,\rho _r,u_i/u_b)=(0.1,2,0.3)$ in group 1 and $(a/H, \rho _r,u_i/u_b)=(0.1,100,0.25)$ in group 3. From table 6, the particle Reynolds numbers for these two cases are $Re_p=163.8$ and $Re_p=169.9$, respectively. From figures 18($b$) and 19, the interfacial term and the particle volume fraction at the channel centre for these two cases are nearly the same. The mean slip velocities are also close to each other, since the particle Reynolds numbers are close. The relative difference in $C_n$ is around $16\,\%$, since the settling coefficients $u_i/u_b$ are 0.3 and 0.275, respectively. As mentioned earlier, the effect of the particle Reynolds stress or the shear flow on the drag force even at the channel centre cannot be neglected for $\rho _r=100$ in the range of $Re_p$ studied. When gravity is not considered, i.e. $u_i/u_b=0$, there is still a pronounced mean drag force on the particles at the channel centre for $\rho _r=100$ (Yu et al. Reference Yu, Lin, Shao and Wang2017), which is mainly balanced with the gradient of the particle Reynolds stress.

Table 6. Results for the simulation cases: groups 1–5. Pressure gradients $(-\mathrm {d} p_{e} / \mathrm {d} x)$ are normalized by $\rho _{f} u_{b}^{2} / H$.

Only a preliminary work on the modelling of the interfacial term in the TKE equation is conducted here. We intend to do further studies on the models for the near-wall region and lower particle Reynolds numbers in the future. The drag and lift models, the RANS model and the second-order model for particle-laden turbulent flows based on the interface-resolved DNS data are good subjects for extensive future work.