3.1 Turbulence attenuation

Figure 2.

The mean velocity scaled by the average velocity ($\bar{u}_{x}/\bar{u}$) (a), the mean square of the fluid fluctuating velocities in the streamwise (b) and cross-stream (c) directions and the correlation $\overline{u_{x}^{\prime }u_{y}^{\prime }}$ (d), all non-dimensionalised $\bar{u}^{2}$, as a function of the scaled coordinate $(y/h)$, for average particle volume fractions $\unicode[STIX]{x1D719}_{av}=0$ (○), $5\times 10^{-4}$ (▵), $10^{-3}$ (▿), $1.3\times 10^{-3}$ (◃), $1.4\times 10^{-3}$ (▹), $1.6\times 10^{-3}$ (♢) and $2\times 10^{-3}$ (▫). The channel Reynolds number is 3300, Stokes law (2.5) is used for the drag force on the particles, and the Stokes number is 105.47.

The mean velocity scaled by the average flow velocity is shown as a function of the cross-stream distance for different particle loadings in figure 2(a). The Stokes law (2.5) is used for the particle drag and the particle Stokes number is 105.5 in all cases. The velocity profiles are shown for the left half of the channel, because the profiles are symmetric about the centre line. Two distinct velocity profiles are observed for $\unicode[STIX]{x1D719}_{av}\leqslant 1.3\times 10^{-3}$ and $\unicode[STIX]{x1D719}_{av}\geqslant 1.4\times 10^{-3}$, with a discontinuous transition between these two profiles. For $\unicode[STIX]{x1D719}_{av}\leqslant 1.3\times 10^{-3}$, the fluid mean velocity profile is close to that for an unladen flow. When the average volume fraction is $1.4\times 10^{-3}$ or higher, the velocity profile is noticeably different from that for an unladen flow, with a smaller slope at the wall and a larger curvature at the centre. The mean square of the velocity fluctuations in the streamwise direction, $\overline{u_{x}^{\prime 2}}$, is shown in figure 2(b). The characteristic near-wall maximum in the streamwise mean square velocity is observed for the unladen flow, and for the particle laden flow for $\unicode[STIX]{x1D719}_{av}\leqslant 1.3\times 10^{-3}$. The height of this maximum decreases by approximately 30 % when the particle volume fraction is increased from 0 to $1.3\times 10^{-3}$. When the particle volume fraction is increased from $1.3\times 10^{-3}$ to $1.4\times 10^{-3}$, there is a dramatic collapse by two orders of magnitude in the mean square velocity. There is very little change in $\overline{u_{x}^{\prime 2}}$ when the average volume fraction is further increased from $1.4\times 10^{-3}$ to $2\times 10^{-3}$. This dramatic collapse is observed in the mean square velocities in the spanwise direction, $\overline{u_{y}^{\prime 2}}$, shown in figure 2(c). There is a larger decrease of approximately 75 % in the values of $\overline{u_{y}^{\prime 2}}$ when the particle volume fraction is increased from 0 to $1.3\times 10^{-3}$. However, when the volume fraction is increased from $1.3\times 10^{-3}$ to $1.4\times 10^{-3}$, there is a decrease of two orders of magnitude in the cross-stream mean square velocity. The variation in the correlation $\overline{u_{x}^{\prime }u_{y}^{\prime }}$ with volume fraction, shown in figure 2(d), is similar to that for the mean square velocities. There is a gradual decrease of approximately 50 % in $\overline{u_{x}^{\prime }u_{y}^{\prime }}$ when the particle volume fraction is increased from 0 to $1.3\times 10^{-3}$, and a discontinuous decrease by over two orders of magnitude when the volume fraction is increased from $1.3\times 10^{-3}$ to $1.4\times 10^{-3}$. The variation in $\overline{u_{x}^{\prime }u_{y}^{\prime }}$ is of significance, since it implies that the Reynolds stress is virtually zero for $\unicode[STIX]{x1D719}_{av}\geqslant 1.4\times 10^{-3}$, and momentum transport is entirely due to the viscous stress or due to transport by the particles.

The discontinuous decrease in the turbulence intensities at a critical particle volume fraction is also observed for other values of the Stokes number, although the critical volume fraction does vary with Stokes number, as discussed a little later. More importantly, the discontinuous decrease is also observed for the Schiller–Naumann correlation (2.6) which has an inertial correction, as shown in figure 3. There is a collapse of the turbulence intensities when the volume fraction is increased from $9\times 10^{-4}$ to $1\times 10^{-3}$; in contrast, there is a continuous variation in the velocity correlations for $\unicode[STIX]{x1D719}_{av}\leqslant 9\times 10^{-4}$ and $\unicode[STIX]{x1D719}_{av}\geqslant 1\times 10^{-3}$. All the features of the second moments of velocity fluctuations and correlation $\overline{u_{x}^{\prime }u_{y}^{\prime }}$ are qualitatively similar to those for the linear drag law, but the critical volume fraction is different, and it does depend on the drag law. However, the phenomenon of turbulence collapse itself appears to be independent of the form of the drag law.

Figure 3.

The mean velocity scaled by the average velocity ($\bar{u}_{x}/\bar{u}$) (a), the mean square of the fluid fluctuating velocities in the streamwise (b) and cross-stream (c) directions, scaled by $\bar{u}^{2}$, and the correlation $(\overline{u_{x}^{\prime }u_{y}^{\prime }}/\bar{u}^{2})$ (d) as a function of the scaled coordinate $(y/h)$, for average particle volume fractions $\unicode[STIX]{x1D719}_{av}=0$ (○), $5\times 10^{-4}$ (▵), $9\times 10^{-4}$ (▿), $10^{-3}$ (◃), $1.1\times 10^{-3}$ (▹) and $1.4\times 10^{-3}$ (♢). The particle force is given by ($\cdots$) Schiller–Naumann (2.6); (– – –) (2.6) with correction for undisturbed velocity at the particle centre (item (ii) in § 2.2); (——) (2.6) with correction for undisturbed velocity at the particle centre (item (ii) in § 2.2), wall correction for drag force (item (iii) in § 2.2) and translational lift force with wall correction (item (v) in § 2.2). The channel Reynolds number is 3300 and the particle Stokes number is 32.45.

Figure 3 also shows the variation in the fluid fluctuating velocities when different particle force models are used. The results of different force models for the mean velocity are not shown in figure 3(a) in order to enhance clarity. In the other panels, the dotted lines are the results when the drag force is given by (2.6), the dashed lines are the results when the correction for the undisturbed velocity profile is included as described in item (ii) in § 2.2 and the solid lines are the results when the lift force and the wall corrections to the drag and lift forces (items (iii) and (v) in § 2.2) are included. There is a discernible reduction of approximately 10 % in the mean square velocities when the correction for the undisturbed velocity profile is included. However, there is very little change in the mean and mean square velocities when the lift force and the wall corrections are included. The modification of the force law has the same qualitative effect on the velocity fluctuations at other particle Stokes numbers as well. More importantly, the critical volume fraction for the turbulence collapse does not change when corrections are included in the force law. For other Stokes numbers, it is found that the critical volume fraction is either unchanged or is altered by a maximum of $10^{-4}$ when the force laws are modified.

The variation of the mean and mean square of the fluctuating velocity at a higher channel Reynolds number of 5600 (friction Reynolds number 180) is shown in figure 4. In this case, turbulence collapse happens at a significantly higher volume fraction of approximately $2.9\times 10^{-3}$. In comparison to the unladen flow, there is a gradual attenuation of the mean square fluctuating velocities of approximately 35 % in the streamwise direction and approximately 50 % in the cross-stream direction and in the Reynolds stress when the particle volume fraction is increased from 0 to $2.9\times 10^{-3}$. When the volume fraction is further increased from $2.9\times 10^{-3}$ to $3\times 10^{-3}$, there is a discontinuous decrease by an order of magnitude. A further increase in the particle volume fraction does not appreciably change the turbulence intensities.

Figure 4.

The mean fluid velocity (a), the mean square of the fluid fluctuating velocities in the streamwise (b) and cross-stream (c) directions and the correlation $\overline{u_{x}^{\prime }u_{y}^{\prime }}$ (d), scaled by suitable powers of $\bar{u}$, as a function of the scaled coordinate $(y/h)$, for average particle volume fractions 0 (○), $2\times 10^{-3}$ (▵), $2.5\times 10^{-3}$ (▿), $2.9\times 10^{-3}$ (◃), $3\times 10^{-3}$ (▹) and $3.5\times 10^{-3}$ (♢). The Schiller–Naumann correlation (2.6) is used for the drag force on the particles. The channel Reynolds number is 5600 and the particle Stokes number is 25.11.

Two other tests have been done to verify that the turbulence collapse observed in the simulations is a fluid-dynamical phenomenon rather than a computational artefact. In the first test, an attempt has been made to determine whether the flow for particle loading just above the critical volume fraction becomes turbulent if it is perturbed in a manner similar to that used for an unladen flow. For an accurate interpretation of the simulation results, it is important to note that, in DNS simulations of an unladen flow, if the initial state is a laminar velocity profile, the flow continues to be in the laminar state even when the Reynolds number exceeds the transition Reynolds number. It is necessary to impose a reasonably large perturbation in order to cause a transition to the turbulent steady state. In the simulations, the perturbation is imposed on the Fourier–Chebyshev transforms of the fluid velocity field, as indicated in § 2.3. The perturbations to the Fourier–Chebyshev coefficients of the streamwise and spanwise velocities, scaled by the maximum of the parabolic velocity profile for a laminar flow, are set to an amplitude $A$ times a random number in the interval $-1$ to $+1$. The cross-stream component of the velocity is determined from the incompressibility condition. The perturbation amplitude $A=0.1$ is used here for the unladen flow. If the Reynolds number exceeds the transition Reynolds number, the turbulent steady state is sustained after the perturbation is imposed. If the Reynolds number is less than the transition Reynolds number, the laminar state is recovered even after the perturbation is imposed.

We have examined the effect of perturbations on a particle laden flow at volume fraction $1.4\times 10^{-3}$ or higher for the Stokes drag law (2.5), and $10^{-3}$ and higher for the nonlinear drag law (2.6). The procedure for generating the perturbations is identical to that used to generate turbulence in an unladen flow. Even when the amplitude is $A=0.15$ (in comparison to the amplitude 0.1 that is adequate for an unladen flow), we observe that the perturbations decay, and the final steady state is as shown in figures 2 and 3.

In the second test, simulations have been carried out for decreasing volume fraction, where the initial state is the steady particle-laden flow with volume fraction $1.4\times 10^{-3}$ using the Stokes drag law at particle Stokes number 105.47. Particles are deleted at random so that the final volume fraction is $1.3\times 10^{-3}$. The state with low fluctuation intensities persists at final volume fraction $1.3\times 10^{-3}$ in the absence of perturbations. However, when we impose perturbations of the same form as those used to trigger turbulence in an unladen flow with the same amplitude $A=0.1$, it is found that the turbulence is sustained, and the fluctuation amplitudes are identical to those for a flow with volume fraction $1.3\times 10^{-3}$ in figure 2. A similar result is obtained when the Schiller–Naumann correlation (2.6) is used for the drag law, and the corrections for the drag and lift forces are included. These tests indicate that the turbulence collapse transition is a robust steady-state phenomenon, independent of numerical details or the path used to attain the final state.

In another series of tests, the steady fluid velocity field for the particle-laden flow is used as the initial condition, and all the particles are deleted. In this case, the flow is unladen, but the initial condition is that for a particle-laden flow at steady state. The fluid velocity field is then evolved in time until a steady state is reached. The results of the tests for the linear drag law with particle Stokes number 105.47 are as follows. When the particle volume fraction is $1.3\times 10^{-3}$ or less, the velocity profile evolves to the unladen turbulent state, as shown in figure 5. When the particle volume fraction is $1.4\times 10^{-3}$ or more, the velocity profile evolves to a parabolic profile for the laminar state. This indicates that when the volume fraction exceeds the critical loading for the particle-laden flow, the velocity fluctuations (of very low magnitude) are due to the presence of particles in a laminar flow and not due to the fluid turbulence; the dramatic decrease in the turbulence intensities is because the flow has completely relaminarised due to the presence of the particles.

Figure 5.

The solid lines show the scaled mean fluid velocity $(\bar{u}_{x}/\bar{u})$ as a function of scaled cross-stream distance for average particle volume fractions $\unicode[STIX]{x1D719}_{av}=1.3\times 10^{-3}$ (○) and $1.4\times 10^{-3}$ (▵). The dashed lines show the steady mean velocity profiles when the solid lines are used as initial conditions and all the particles are deleted from the simulation domain. The channel Reynolds number is 3300, Stokes law is used for the particle drag and the Stokes number is 105.47.

The turbulence collapse is observed for particles with different particle Stokes numbers, as shown in figure 6. A measure of the departure from the unladen mean velocity profile, $(\unicode[STIX]{x0394}\bar{u}_{x})^{2}$, is defined as

(3.1)$$\begin{eqnarray}(\unicode[STIX]{x0394}\bar{u}_{x})^{2}(\unicode[STIX]{x1D719}_{av})=\frac{\displaystyle \int _{0}^{h}\,\text{d}y(\bar{u}_{x}-\bar{u}_{x0})^{2}}{\displaystyle \int _{0}^{h}\,\text{d}y(\bar{u}_{x0})^{2}},\end{eqnarray}$$

where $\bar{u}_{x}$ is the mean velocity for average volume fraction $\unicode[STIX]{x1D719}_{av}$ and $\bar{u}_{x0}$ is the mean velocity for the unladen turbulent flow. The averages of the mean square of the fluctuating velocities over the channel width are defined as measures of the mean square velocities denoted by the angle brackets $\langle \,\rangle _{s}$,

(3.2)$$\begin{eqnarray}\langle \star \rangle _{s}=\frac{2}{h}\int _{0}^{h/2}\,\text{d}y\star .\end{eqnarray}$$

These measures, scaled by the square of the average velocity, are shown as a function of volume fraction in figure 6 for different particle Stokes numbers and for two different forms of the particle drag law. The discontinuous change in the velocity profile and the discontinuous decrease in the turbulence intensities are observed at a critical volume fraction for all particle Stokes numbers, and for two different drag laws. The measures of the turbulence intensity clearly decrease by two orders of magnitude at the transition volume fraction when the Stokes drag law is used to calculate particle drag. Even when the Schiller–Naumann correlation (2.6), with corrections for the undisturbed velocity at the particle centre, lift and the presence of the wall, is used, there is a decrease of approximately 1.5 orders of magnitude in the turbulence intensities.

Figure 6.

The measure $\unicode[STIX]{x0394}\bar{u}_{x}$ (3.1) for the difference between the mean velocity profile and the unladen turbulent velocity profile (a), the measures of the fluid fluctuating velocity (3.2) $(\langle \overline{u_{x}^{\prime 2}}\rangle _{s}/\bar{u}^{2})$ (b), $(\langle \overline{u_{y}^{\prime 2}}\rangle _{s}/\bar{u}^{2})$ (c) and $(\langle |\overline{u_{x}^{\prime }u_{y}^{\prime }}|\rangle _{s}/\bar{u}^{2})$ (d) as a function of the average particle volume fraction $\unicode[STIX]{x1D719}_{av}$. The solid lines are the results for channel Reynolds number 3300 when the Stokes law (2.5) is used for the drag force on the particles and for Stokes numbers of 21.09 (○), 105.47 (▵) and 316.40 (▿). The dashed lines are the results for channel Reynolds number 3300 when the Schiller–Naumann correlation (2.6) is used for the drag force, the correction for the undisturbed velocity and the wall effect on the particle drag and lift forces are included and for Stokes numbers of 6.49 (●), 32.45 (▴) and 97.34 (▾). The dotted lines are the results for channel Reynolds number 5600 when the particle force is given by the Schiller–Naumann correlation (2.6), for Stokes numbers of 25.11 (♢), 50.22 (◃) and 75.33 (▹).

There is also a large reduction in the wall shear stress at the critical volume fraction. It should be recalled that, in our protocol, the average fluid velocity is maintained constant as the particle loading is increased by varying the pressure gradient. The variation of the scaled wall shear stress, $(\unicode[STIX]{x1D70F}_{w}/(\unicode[STIX]{x1D70C}_{f}\bar{u}^{2}))$ is shown as a function of the average volume fraction in figure 7(a). There is a step decrease in the scaled wall shear stress at the critical volume fraction, and then a gradual increase as the volume fraction is further increased. The scaled wall shear stress is also equal to $(u_{\ast }/\bar{u})^{2}$, where $u_{\ast }=(\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D70C}_{f})^{1/2}$ is the friction velocity based on the fluid density. The friction Reynolds number $Re_{\ast }=(\unicode[STIX]{x1D70C}_{f}u_{\ast }h/2\unicode[STIX]{x1D702}_{f})$, decreases sharply at the transition, due to the decrease in the wall shear stress, even though the Reynolds number based on the average flow velocity and the channel width is unchanged.

Figure 7.

The wall shear stress $\unicode[STIX]{x1D70F}_{w}$ scaled by $\unicode[STIX]{x1D70C}_{f}\bar{u}^{2}$ and the Reynolds number $Re_{\ast }=(\unicode[STIX]{x1D70C}_{f}u_{\ast }h/2\unicode[STIX]{x1D702}_{f})$ based on the friction velocity (b) as a function of the average particle volume fraction $\unicode[STIX]{x1D719}_{av}$. The solid lines, referenced to the left $y$ axis in (b), are the results for channel Reynolds number 3300 when the Stokes law (2.5) is used for the drag force on the particles for Stokes numbers 21.09 (○), 105.47 (▵) and 316.40 (▿). The dashed lines, referenced to the left $y$ axis in (b), are the results for channel Reynolds number 3300 when the Schiller–Naumann correlation (2.6) is used for the drag force, the correction for the undisturbed velocity and the wall effect on particle drag and lift forces are included for Stokes numbers of 6.49 (●), 32.45 (▴) and 97.34 (▾). The dotted lines, referenced to the right $y$ axis in (b), are the results for channel Reynolds number 5600 when the particle force is given by the Schiller–Naumann correlation (2.6), for Stokes numbers of 25.11 (♢), 50.22 (◃) and 75.33 (▹). The vertical dotted lines show the range of $\unicode[STIX]{x1D719}_{av}$ over which there is turbulence collapse for the particle Stokes number and drag law indicated by the symbols on the abscissa.

The fluid velocity distributions are examined both at the centre and near the wall of the channel. Since there is a variation in the flow properties in the cross-stream $y$ direction, the distributions are calculated in two zones, one of extent 20 % of the channel width at the wall, and the second spanning 20 % of the channel width at the centre. In these two zones, the distributions for the fluctuating velocities in the three directions, $P(u_{x}^{\prime })$, $P(u_{y}^{\prime })$ and $P(u_{z}^{\prime })$, are shown as a function of the velocity fluctuations scaled by the standard deviation in figure 8. In figure 8(a) for the streamwise velocity distribution, there is not much change in the distribution function close to the wall. The distributions at $\unicode[STIX]{x1D719}_{av}=1.3\times 10^{-3}$ and $1.4\times 10^{-3}$ are both close to a Gaussian distribution, but with a small positive skewness. The distribution function for the streamwise fluctuations at the centre does show some change at transition; the distribution has a pronounced negative skewness at $\unicode[STIX]{x1D719}_{av}=1.3\times 10^{-3}$, but it becomes close to a symmetric Gaussian distribution after turbulence collapse at $\unicode[STIX]{x1D719}_{av}=1.4\times 10^{-3}$. There is a striking change in the probability distribution for the cross-stream velocity fluctuations in figure 8(b). The peaked distribution function with high-velocity tails for the turbulent flow at $\unicode[STIX]{x1D719}_{av}=1.3\times 10^{-3}$ transitions to a distribution function that is closer to a Gaussian after turbulence collapse at $\unicode[STIX]{x1D719}_{av}=1.4\times 10^{-3}$ in both zones. A similar, but less pronounced, transition is observed for the spanwise velocity fluctuations in figure 8(c). Thus, the large attenuation of the velocity fluctuations is also accompanied by a change in the forms of the velocity distribution functions.

Figure 8.

The probability distributions for the fluid velocity fluctuations in the streamwise (a), cross-stream (b) and spanwise (c) directions in the zone near the wall (○) and at the centre (▵) for channel Reynolds number $Re_{b}=3300$, particle volume fractions $1.3\times 10^{-3}$ (red) and $1.4\times 10^{-3}$ (blue). The Stokes drag law is used for the particles, and the Stokes number is 105.47. The dotted line is the Gaussian distribution.

The variation in the critical volume fraction with particle Stokes number is shown in figure 9(a). For the Stokes drag law, there is an initial decrease in the critical volume fraction as the Stokes number is increased for $St\lesssim 100$, and then the critical volume fraction tends to a constant value for $St\gtrsim 100$. Similarly, for the Schiller–Naumann correlation (2.6), there is a decrease in the critical volume fraction with Stokes number for $St\lesssim 40$, and the critical volume fraction is independent of the Stokes number for $St\gtrsim 40$.

An examination of (2.2) and (2.4) reveals the reason for the lack of dependence of the critical volume fraction on the particle Stokes number at high Stokes number. The force $\boldsymbol{f}$ in (2.2) contains two parts, the drag force and the gravitational force. The gravitational force $m_{p}\boldsymbol{g}$ in (2.3) does depend on the particle material density, but it turns out that the gravitational force is much smaller than the drag force. For the largest Stokes number considered here, the ratio $(|\boldsymbol{v}_{t}|/\bar{u})$ is $1.3\times 10^{-2}$, and therefore the gravitational force is approximately two orders of magnitude smaller than the drag due to the gas flow. The drag force, given by (2.5) or (2.6), is independent of the particle material density. The particle acceleration in (2.4) does depend on the particle material density. However, in a steady flow, the average particle acceleration is zero. Since the Stokes number is varied here by changing the particle material density and keeping the particle diameter and fluid viscosity unchanged, the critical volume fraction is independent of Stokes number for high Stokes number.

Figure 9.

The critical average volume fraction $\unicode[STIX]{x1D719}_{cr}$ for turbulence collapse as a function of the Stokes number for (○) the Stokes drag law (2.5) and $Re_{b}=3300$, (▵) the Schiller–Naumann correlation (2.6) and $Re_{b}=3300$, (▿) the Schiller–Naumann correlation (2.6) with the correction to undisturbed velocity at the particle centre and the wall effect on drag and lift force and $Re_{b}=3300$ and (♢) for the Schiller–Naumann correlation (2.6) and $Re_{b}=5600$.

The increase in the critical volume fraction at low Stokes number is less easily understood. It is expected that, when the Stokes number becomes smaller, there is less force exerted by the particles and fluid, and the difference in the local particle and fluid velocities is smaller. Due to this, a higher number of particles are required for turbulence collapse. This expectation is consistent with the increase in the critical volume fraction as the Stokes number is decreased. However, the data in figure 9 are insufficient to deduce a clear scaling law for the dependence of the critical volume fraction on Stokes number in this limit.

3.2 Particle fluctuations

The details of the mean and mean square velocities and the particle velocity distributions are summarised in appendix B. The important conclusions that are as follows.

1. (i) The particle volume fraction is approximately constant across the channel, and there is virtually no change in the cross-stream variation of the volume fraction at turbulence collapse.

2. (ii) The particle mean velocity is nearly constant across the channel, below and above the critical volume fraction. The particle mean velocity is approximately equal to the average flow velocity $\bar{u}$, because the terminal velocity of the particles is small compared to the flow velocity.

3. (iii) There is a small increase of approximately 15 %–30 % in the mean and mean square of the particle fluctuating velocities at the critical volume fraction. This is much smaller than the decrease, by approximately an order of magnitude, in the fluid turbulence intensities.

4. (iv) Below the critical volume fraction, the particle fluctuating velocities are smaller than the fluid fluctuating velocities, suggesting that the fluid turbulence is driving the particle fluctuations. Above the critical volume fraction, the particle fluctuating velocities are higher than the fluid fluctuating velocities, indicating that the particle fluctuations are driving the fluid fluctuations.

5. (v) There is no qualitative change in the distribution functions for the particle velocities at the critical volume fraction.

3.3 Momentum balance

The fluid mean momentum equations are

(3.3)$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}\bar{p}}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D702}_{f}\frac{\text{d}^{2}\bar{u}_{x}}{\text{d}y^{2}}+\frac{\text{d}\unicode[STIX]{x1D70F}_{xy}^{Rf}}{\text{d}y}-\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}, & \displaystyle\end{eqnarray}$$

(3.4)$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}\bar{p}}{\unicode[STIX]{x2202}y}+\frac{\text{d}\unicode[STIX]{x1D70F}_{yy}^{Rf}}{\text{d}y}-\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{y}-v_{y})}, & \displaystyle\end{eqnarray}$$

where $n_{p}$ is the number density of the particles, $u_{i}$ and $v_{i}$ are the local undisturbed fluid velocity and the particle velocity,

(3.5)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{xy}^{Rf}=-\unicode[STIX]{x1D70C}_{f}\overline{u_{x}^{\prime }u_{y}^{\prime }}\end{eqnarray}$$

is the Reynolds stress, and $\unicode[STIX]{x1D707}_{f}$ is the drag coefficient, which is a function of the difference between the fluid and particle velocities for the nonlinear drag law (2.6). The fourth term on the right in (3.3) and the third term on the right in (3.4) are due to the drag force exerted by the particles on the fluid.

The terms in the momentum balance equation (3.3), are shown in figure 10 for the Schiller–Naumann correlation (2.6), with correction to the undisturbed velocity at particle centre and the wall effect on drag and lift force. The force density due to the particles, $-\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}$, exhibits a sharp decrease to zero very close to the wall, because the particles are excluded within a region of width one particle radius from the wall. Close to the wall, the separation between the fluid grid points used for Chebyshev collocation becomes smaller than the particle diameter; the minimum grid spacing is approximately 0.7 times the particle diameter. In order to account for the finite size of the particles, the results in figure 10 have been calculated by relaxing the point-force approximation, and calculating the force within a Chebyshev grid interval based on the fraction of the particle surface within the cell, as prescribed by the Faxen theorem. This procedure results in the continuous decrease in the force density due to the particles close to the wall. The procedure of Esmaily & Horwitz (Reference Esmaily and Horwitz2018) used here was validated against fully resolved particle simulations, and was found to accurately predict the turbulence statistics even when the particle size is comparable to or larger than the grid size. It is also important to note that the maximum in the divergence of the Reynolds stress is at a distance greater than the particle diameter from the wall.

When the volume fraction is increased from $0.9\times 10^{-3}$ to $1\times 10^{-3}$, there is little change in the magnitude and shape of the terms $(\unicode[STIX]{x1D702}_{f}\text{d}^{2}\bar{u}_{x}/\text{d}y^{2})$ (divergence of the viscous stress) and $-\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}$ (force density due to the particles). There is a significant change in the divergence of the Reynolds stress, which is larger than the force density due to the particles near the wall at $\unicode[STIX]{x1D719}_{av}=0.9\times 10^{-3}$ (figure 10a), but decreases to virtually zero across the entire channel for $\unicode[STIX]{x1D719}_{av}=1\times 10^{-3}$ (figure 10b). For $\unicode[STIX]{x1D719}_{av}=1\times 10^{-3}$, there is a balance between the divergence of the viscous stress due to the mean flow, the drag force due to the particles and the pressure gradient, as shown in figure 10(b).

Figure 10.

The terms in the fluid momentum equation (3.3) scaled by $(\unicode[STIX]{x1D70C}_{f}\bar{u}^{2}/h)$ as a function of scaled cross-stream distance for volume fractions $9\times 10^{-4}$ (a) and $10^{-3}$ (b). The channel Reynolds number is 3300, the particle Stokes number is 32.45, the particle force is given by the Schiller–Naumann correlation (2.6), with the correction to the undisturbed velocity at the particle centre and the wall effect on drag and lift force are included. The different symbols are ○ — $-(\unicode[STIX]{x2202}\bar{p}/\unicode[STIX]{x2202}x)$, ▵ — $(\unicode[STIX]{x1D702}_{f}\text{d}^{2}\bar{u}_{x}/\text{d}y^{2})$, ▿ — $(\text{d}\unicode[STIX]{x1D70F}_{xy}^{Rf}/\text{d}y)$, ♢ — $-\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}$.

3.4 Energy balance

The equation for the fluid mean kinetic energy is

(3.6)$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle -\bar{u}_{x}\frac{\unicode[STIX]{x2202}\bar{p}}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D702}_{f}\frac{\text{d}}{\text{d}y}\left(\bar{u}_{x}\frac{\text{d}\bar{u}_{x}}{\text{d}y}\right)-\unicode[STIX]{x1D702}_{f}\left(\frac{\text{d}\bar{u}_{x}}{\text{d}y}\right)^{2}+\frac{\text{d}(\unicode[STIX]{x1D70F}_{xy}^{Rf}\bar{u}_{x})}{\text{d}y}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D70F}_{xy}^{Rf}\frac{\text{d}\bar{u}_{x}}{\text{d}y}-\bar{u}_{x}\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}.\end{eqnarray}$$

In (3.6), the first term on the right is the rate of increase of energy due to the pressure gradient, the third term on the right is the viscous dissipation of energy due to the mean shear, the fifth term on the right is the turbulent production of energy due to the Reynolds stress, which results in a transfer of energy from the mean flow to the fluctuations, and the sixth term is the rate of dissipation of energy due to the mean particle drag. The second and fourth terms on the right of (3.6) are the divergences of energy fluxes due to the fluid viscous stress and the Reynolds stress, respectively. The fluctuating kinetic energy balance equation is

(3.7)$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle -\frac{\text{d}q_{y}}{\text{d}y}-\frac{\text{d}\overline{p^{\prime }u_{y}^{\prime }}}{\text{d}y}+\unicode[STIX]{x1D702}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\overline{\boldsymbol{u}^{\prime }\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime })^{\text{T}}})+\unicode[STIX]{x1D70F}_{xy}^{Rf}\frac{\text{d}\bar{u}_{x}}{\text{d}y}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D702}_{f}\overline{(\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime })\boldsymbol{ : }(\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }+(\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime })^{\text{T}})}-D_{f},\end{eqnarray}$$

where

(3.8)$$\begin{eqnarray}q_{y}={\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}_{f}\overline{u_{y}^{\prime }u^{\prime 2}}\end{eqnarray}$$

is the energy flux due to the fluid fluctuating velocity, and

(3.9)$$\begin{eqnarray}\displaystyle D_{f}(\boldsymbol{x},t) & = & \displaystyle \overline{\unicode[STIX]{x1D707}_{f}n_{p}u_{x}(u_{x}-v_{x})}+\overline{\unicode[STIX]{x1D707}_{f}n_{p}u_{y}(u_{y}-v_{y})}+\overline{\unicode[STIX]{x1D707}_{f}n_{p}u_{z}(u_{z}-v_{z})}\nonumber\\ \displaystyle & & \displaystyle -\,\bar{u}_{x}\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}\end{eqnarray}$$

is the difference between the total energy dissipation rate due to the particles and the mean energy dissipation rate in (3.6).

When (3.6) and (3.7) are averaged across the channel cross-section using (3.2), divergences of the energy fluxes average to zero due to no-slip conditions at the solid walls, and we obtain

(3.10)$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\bar{u}\frac{\unicode[STIX]{x2202}\bar{p}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D702}_{f}\left\langle \left(\frac{\unicode[STIX]{x2202}\bar{u}_{x}}{\unicode[STIX]{x2202}y}\right)^{2}\right\rangle _{s}-\left\langle \unicode[STIX]{x1D70F}_{xy}^{Rf}\frac{\unicode[STIX]{x2202}\bar{u}_{x}}{\unicode[STIX]{x2202}y}\right\rangle _{s}-\langle \bar{u}_{x}\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}\rangle _{s}, & \displaystyle\end{eqnarray}$$

(3.11)$$\begin{eqnarray}\displaystyle & \displaystyle 0=\left\langle \unicode[STIX]{x1D70F}_{xy}^{Rf}\frac{\unicode[STIX]{x2202}\bar{u}_{x}}{\unicode[STIX]{x2202}y}\right\rangle _{s}-\unicode[STIX]{x1D702}_{f}\langle \overline{(\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime })\boldsymbol{ : }((\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }+\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime })^{\text{T}})}\rangle _{s}-\langle D_{f}\rangle _{s}. & \displaystyle\end{eqnarray}$$

The terms in the energy balance equations, equations (3.6) and (3.7), are shown in figure 11. The divergences of the energy fluxes, which are the second and fourth terms on the right in (3.6) and the first, second and third terms on the right in (3.7), are not shown, since these are the transport terms which do not alter the total fluid energy dissipation rate. For volume fraction $\unicode[STIX]{x1D719}=0.9\times 10^{-3}$, figure 11(a) shows that the largest terms are the mean pressure work $\bar{u}_{x}(\text{d}\bar{p}/\text{d}x)$, the mean viscous dissipation $(\unicode[STIX]{x1D702}(\text{d}\bar{u}_{x}/\text{d}y)^{2})$ and the dissipation due to the particle drag $\bar{u}_{x}\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}$. The transport of energy from the mean flow to the fluctuations due to the Reynolds stress $-\unicode[STIX]{x1D70F}_{xy}^{Rf}(\text{d}\bar{u}_{x}/\text{d}y)$, and the viscous dissipation due to the fluid fluctuating velocity $\unicode[STIX]{x1D702}_{f}(\overline{\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }\boldsymbol{ : }(\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }+(\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime })^{\text{T}})})$ are comparable in magnitude close to the wall, while the dissipation $D_{f}$ (3.9) due to the drag force is smaller. When the volume fraction is increased to $1\times 10^{-3}$, the transfer of energy from mean flow to the fluctuations $-\unicode[STIX]{x1D70F}_{xy}^{Rf}(\text{d}\bar{u}_{x}/\text{d}y)$ decreases virtually to zero throughout the channel, and consequently, the energy dissipation due to the fluid fluctuations also decreases to zero. There is a balance between the work done due to the mean pressure gradient, the viscous energy dissipation and the dissipation due to the mean particle drag.

Figure 11.

The terms in the fluid energy equations (3.6) and (3.7), scaled by $(\unicode[STIX]{x1D70C}_{f}\bar{u}^{3}/h)$ as a function of scaled cross-stream distance for volume fractions $9\times 10^{-4}$ (a) and $10^{-3}$ (b). The channel Reynolds number is 3300, the particle force is given by the Schiller–Naumann correlation (2.6), with correction to the undisturbed velocity at the particle centre and the wall effect on drag and lift force, and the particle Stokes number is 32.45. The different symbols are ○ $-\bar{u}_{x}(\text{d}p/\text{d}x)$, ▵ $(\unicode[STIX]{x1D702}_{f}(\text{d}\bar{u}_{x}/\text{d}y)^{2})$, ▿ $-\unicode[STIX]{x1D70F}_{xy}^{Rf}(\text{d}\bar{u}_{x}/\text{d}y)$, ♢ $-\bar{u}_{c}\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}$, ◃ $\unicode[STIX]{x1D702}_{f}(\overline{\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }\boldsymbol{ : }(\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }+(\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime })^{\text{T}})})$, ▹ $D_{f}$.

The terms in the spatially averaged fluid energy conservation equation, equation (3.10), averaged across the channel cross-section, $\bar{u}(\text{d}\bar{p}/\text{d}x)$, the turbulent production $\langle \unicode[STIX]{x1D70F}_{xy}^{Rf}(\unicode[STIX]{x2202}\bar{u}_{x}/\unicode[STIX]{x2202}y)\rangle _{s}$ and the dissipation due to the particles $\langle \bar{u}_{x}\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}\rangle _{s}$, all scaled by $(\unicode[STIX]{x1D70C}_{f}\bar{u}^{3}/h)$, are shown in figure 12. The turbulent energy production rate, shown by the ○ symbol, shows a large decrease at the critical volume fraction, consistent with the results in figure 11. The total energy dissipation rate due to the particle drag, shown by the ▵ symbol, does increase at the critical volume fraction. However, this increase is not more than one half of the decrease in the turbulent production rate. There is also a decrease in the energy dissipation rate due to the mean shear, with the consequence that the total energy dissipation rate in the fluid actually decreases significantly at the transition volume fraction. This indicates that the decrease in the turbulent energy production is not due to a compensating increase in the rate of dissipation of energy due to the particles. Consequently, turbulence attenuation is not due to the increase in dissipation due to particle drag, but rather due to the disruption of the turbulence production mechanism.

Figure 12.

The cross-sectionally averaged turbulent energy production rate $\langle \unicode[STIX]{x1D70F}_{xy}^{Rf}(\text{d}\bar{u}_{x}/\text{d}y)\rangle _{s}$ (○), the rate of dissipation of energy due to particle drag $\langle \bar{u}_{x}\overline{\unicode[STIX]{x1D707}_{f}n_{p}(u_{x}-v_{x})}\rangle _{s}$ (▵), the total rate of input of energy $\bar{u}(\unicode[STIX]{x2202}\bar{p}/\unicode[STIX]{x2202}x)$ (▿) and the rate of dissipation of energy due to mean shear $\unicode[STIX]{x1D702}_{f}\langle (\text{d}\bar{u}_{x}/\text{d}y)^{2}\rangle _{s}$ (♢), all scaled by $(\unicode[STIX]{x1D70C}_{f}\bar{u}^{3}/h)$, as a function of volume fraction. In (a), the Stokes drag law (2.5) is used for the drag force, the channel Reynolds number is 3300, with particle Stokes numbers of 21.09 (dotted lines), 105.47 (solid lines) and 316.40 (dashed lines). In (b), the Schiller–Naumann correlation (2.6) is used for the drag force, the channel Reynolds number is 5600 and for particle Stokes numbers of 25.11 (dotted lines), 50.22 (dashed lines) and 75.33 (solid lines).

3.5 Turbulence decay

The time evolution of the terms in the energy balance equation, averaged over the half-width of the channel, due to the addition of particles from an initially unladen flow is examined in figure 13. Here, the initial condition is a fully developed unladen turbulent flow at steady state. The particles are instantaneously added at $t=0$ at random locations with a velocity equal to the interpolated fluid velocity at the particle location, so that the final volume fraction is $1\times 10^{-3}$. The terms in the kinetic energy equation are shown as a function of the scaled time, $(t\bar{u}/h)$. There is a sharp decrease in the turbulent energy production rate upon addition of particles, and the production rate virtually decreases to zero within approximately 200 integral times. There is a smaller and less sharp increase in the dissipation rate due to the particle drag. The mean viscous dissipation rate also decreases, resulting in a decrease in the total energy dissipation rate of energy. This clearly shows that the addition of particles disrupts the turbulent energy production rate, and the collapse of the energy production is not accompanied by a commensurate increase in the energy dissipation due to the drag force exerted by the particles. This disruption of the energy production results in a decrease in the pressure drop at constant flow rate. This is consistent with the relaminarisation of the flow – because the wall shear stress decreases when the flow is relaminarised, a lower pressure drop is required to drive the flow.

Figure 13.

The spatially averaged turbulent energy production rate $\langle \unicode[STIX]{x1D70F}_{xy}^{Rf}(\text{d}\bar{u}_{x}/\text{d}y)\rangle _{s}$ (○), the rate of dissipation of energy due to particle drag $\langle \bar{u}_{x}\overline{\unicode[STIX]{x1D707}_{f}n_{p}(v_{x}-u_{x})}\rangle _{s}$ (▵), the total rate of input of energy $\bar{u}(\unicode[STIX]{x2202}\bar{p}/\unicode[STIX]{x2202}x)$ (♢) and the rate of dissipation of energy due to mean shear $\unicode[STIX]{x1D702}_{f}\langle (\text{d}\bar{u}_{x}/\text{d}y)^{2}\rangle _{s}$ (▿), all scaled by $(\unicode[STIX]{x1D70C}_{f}\bar{u}^{3}/h)$, as a function of scaled time $(t\bar{u}/h)$. The particle force is given by the Schiller–Naumann correlation (2.6), the channel Reynolds number $Re_{b}$ is 3300 and the particle Stokes number is 32.45.

Figure 14.

The turbulent energy production rate $\unicode[STIX]{x1D70F}_{xy}^{Rf}(\text{d}\bar{u}_{x}/\text{d}y)$ (a), the rate of dissipation of energy due to particle drag $\bar{u}_{x}\overline{\unicode[STIX]{x1D707}_{f}n_{p}(v_{x}-u_{x})}$ (b) and the rate of dissipation of energy due to mean shear $\unicode[STIX]{x1D702}_{f}(\text{d}\bar{u}_{x}/\text{d}y)^{2}$, all scaled by $(\unicode[STIX]{x1D70C}_{f}\bar{u}^{3}/h)$, averaged over scaled time periods $(t\bar{u}/h)=0{-}100$ (○), 100–200 (▵), 200–300 (▿), 300–400 (◃), 400–500 (▹) and 500–600 (♢). The channel Reynolds number is 3300, the particle force is given by the Schiller–Naumann correlation (2.6), with the correction to the undisturbed velocity at the particle centre and the wall effect on drag and lift force, and the particle Stokes number is 32.45.

The time evolution of the profiles for the turbulent energy production rate, the dissipation rate due to the mean shear and the dissipation rate due to the particle drag are shown in figure 14. The profiles are averaged over six time intervals of extent $(100h/\bar{u})$, in order to obtain profiles that are relatively smooth in comparison to the instantaneous profiles. The sharp near-wall maximum in the turbulent production rate is evident in the first time interval in figure 14(a). However, there is a rapid decrease in the energy production rate in the second time interval, and the energy production rate is virtually zero within $(300h/\bar{u})$. In figure 14(b) there is an increase in the magnitude of the energy dissipation due to the drag near the centre of the channel. This is because there is an increase in the fluid mean velocity at the centre of the channel, shown in figure 2(a), when the particles are added to an unladen flow, whereas there is virtually no change in the particle mean velocity, shown in figure 19. This results in an increase in $\bar{u}_{x}-\bar{v}_{x}$ at the centre of the channel, and therefore an increase in the magnitude of the energy dissipation due to particle drag. However, in the near-wall region, where there is a significant decrease in the turbulent energy production, there is virtually no change in the energy dissipation due to the particles. There is, however, a significant decrease in the viscous dissipation of energy near the wall, as shown in figure 14(c), due to a decrease in the fluid strain rate at the wall when compared to an unladen flow (see figure 2a). The profiles of the production and dissipation rates clearly show a collapse in the turbulent energy production, with no comparable increase in the dissipation due to the particles, and an overall decrease in the energy dissipation rate. This indicates that the turbulence is not attenuated due to an increase in the energy dissipation rate due to particle drag in the zone of maximum turbulent production, but rather the particles seem to have a direct role in disrupting the turbulence production mechanism.