1 Introduction
The modification of fluid turbulence by the particles in a gas–particle suspension is an active area of research, due to the importance of these suspensions in geophysical phenomena such as snowstorms and sand avalanches, as well as in industrial applications involving fluidisation and pneumatic conveying. It is well established, since the early work of Gore & Crowe (Reference Gore and Crowe1989) and Elghobashi & Truesdell (Reference Elghobashi and Truesdell1993), that there is significant turbulence modification even when the particle volume fraction is in the range 
$10^{-4}$–
$10^{-3}$. The initial studies were primarily experimental, and relied on direct measurements of average and local velocity statistics. Advances in high-speed computation subsequently facilitated direct simulations of turbulent flows with large numbers of suspended particles. There is now growing evidence that large particles enhance turbulent fluctuations due to the formation of a wake behind the particles, whereas small particles attenuate turbulence due to increased dissipation by particle drag. However, sharp criteria delineating regimes of turbulence intensification or attenuation, as well as the extent of turbulence modification, are not yet available. Some of the earlier studies are first summarised, followed by a discussion of the objective of the present study.
In an early review, Gore & Crowe (Reference Gore and Crowe1989) analysed the experiments done on turbulence modification in multi-phase flows, and proposed a classification using a single critical parameter 
$d_{p}/L$ where 
$d_{p}$ is the diameter of the particle and 
$L$ is the integral length scale of turbulence. They reported that turbulence is attenuated for 
$(d_{p}/L)<0.1$, and augmented for 
$(d_{p}/L)>0.1$. Hetsroni (Reference Hetsroni1989) was perhaps the first to suggest that particles with low Reynolds number suppress the turbulence, while particle with high Reynolds number increase turbulence due to the wake formation. In a review of the experimental and numerical results, Tanaka & Eaton (Reference Tanaka and Eaton2008) have introduced the particle momentum number (
$Pa$),
$$\begin{eqnarray}Pa=\frac{1}{18}\frac{Re_{L}^{2}}{Re_{p}}\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}_{f}}\left(\frac{d_{p}}{L}\right)^{3},\end{eqnarray}$$ to classify and understand different regimes of turbulence modification in the particle laden flows. Here, 
$\unicode[STIX]{x1D70C}_{p}$ and 
$\unicode[STIX]{x1D70C}_{f}$ are the mass densities of the particle and fluid, 
$L$ is the characteristic flow scale, 
$d_{p}$ is the particle diameter, 
$Re_{L}$ is the flow Reynolds number and 
$Re_{p}$ is the particle Reynolds number. They proposed that the particle momentum number, and not the Reynolds number, determines the turbulence modification. For particles with extremely small 
$Pa$, there is no turbulence modification. Fluid turbulence is augmented at small 
$Pa$, attenuated at moderate 
$Pa$ and again augmented at high 
$Pa$.
In numerical simulations, Li & McLaughlin (Reference Li and McLaughlin2001) studied the effect of particle feedback on fluid turbulence and found that at small mass loading, the particles tend to increase the turbulent intensities, while they suppress the turbulence at high mass loadings. They also found that particle–particle collisions tend to reduce the accumulation of the particles at the wall, even for inelastic collisions. Vreman et al. (Reference Vreman, Geurts, Deen, Kuipers and Kuerten2009) performed large eddy simulations (LES) of the particle-laden turbulent flow in a vertical riser with a volume fraction of 1.5 %. The presence of a large number of interacting particles through collisions led to a strong modulation of the turbulence in the channel. Capecelatro, Desjardins & Fox (Reference Capecelatro, Desjardins and Fox2015, Reference Capecelatro, Desjardins and Fox2018) studied the effect of particle clustering on the turbulence modification in a channel, and suggested a transition from a dilute limit where mean-shear production is the dominant mechanism of turbulence generation, to a dense limit where particle-induced fluctuations are the dominant mechanism. Vreman (Reference Vreman2015) studied the effect of wall roughness on turbulence modification, and concluded that turbulence attenuation is enhanced due to wall roughness. The turbulence attenuation due to an increase in particle loading has also been reported in more recent simulations (Vreman Reference Vreman2015; Capecelatro et al. Reference Capecelatro, Desjardins and Fox2018). A gradual decrease in the turbulence intensities with increasing particle loading is considered to be caused by the increasing dissipation of energy due to the particles compensating for the transfer of energy to turbulent fluctuations. Reynolds-averaged equations have been developed by Capecelatro, Desjardins & Fox (Reference Capecelatro, Desjardins and Fox2016) for particle–gas suspension flows in a vertical channel at high mass loading which include the particle-phase turbulent kinetic energy. Richter (Reference Richter2015) found that the presence of particles reduces the production of kinetic energy at all wavenumbers, and that the turbulence production is much larger than the exchange of energy between the fluid and the particles. Wang & Richter (Reference Wang and Richter2019) reported that the high inertia particles dampen the large scale vortices, whereas low inertia particles amplify the large scale vortices and reduce the transition Reynolds number.
The broad consensus is that small particles moving through the fluid at low particle Reynolds number tend to attenuate turbulence due to the drag force, whereas large particles moving at high Reynolds number increase turbulence due to the wake effect. The turbulence attenuation is considered to be a gradual process as the particle loading increases, due to the increase in the dissipation of energy due to the particle drag. Here, we attempt to examine the turbulence attenuation as a function of particle loading in a vertical channel. The force exerted by a particle on the fluid is modelled as a point force, and the particle Reynolds number is sufficiently small that there could a steady separation bubble but no vortex shedding at the rear.
In order to detect precisely the changes in the turbulence attenuation due to the particles, a large number of simulations are carried out for small increments in the particle loading at different particle Stokes numbers and for two different forms of the drag law. The channel Reynolds number based on the gas density, average flow velocity and the channel width is maintained at a two relatively low values of approximately 3300 and 5600 to enable the necessary number of computations in reasonable time. The particle Reynolds number based on the particle mean velocity and diameter has values of 60 and 103. Taneda (Reference Taneda1956) has shown that, at these Reynolds numbers, there is a steady separation bubble at the rear (which appears at a Reynolds number of 24), but no unsteady vortex shedding (which starts at a Reynolds number of 130). The particle Reynolds number based on the channel-averaged root mean square particle–fluid velocity difference is in the range 4–15 for all the Reynolds and Stokes numbers reported here. The particle Stokes number based on the particle relaxation time, the bulk velocity and channel width is varied in the range 1–420. The particle size is comparable to the Kolmogorov scale, so that the point-particle approximation is used for the particle force on the fluid. The particle volume fraction is varied in the range 
$0{-}3.5\times 10^{-3}$, and the mass loading is in the range 0–13.6. The parameters used here are compared with those in earlier studies in table 1. The channel Reynolds number and the friction Reynolds number studied here are comparable to or lower than those in earlier studies, and the particle volume fraction is comparable to those used in earlier studies. However, the range of particle Stokes numbers and mass loading is significantly higher than those in previous studies.
Table 1. The channel Reynolds number 
$Re_{b}=(\unicode[STIX]{x1D70C}_{f}\bar{u}h/\unicode[STIX]{x1D702}_{f})$, friction Reynolds number 
$Re_{\ast }=(\unicode[STIX]{x1D70C}_{f}u_{\ast }h/2\unicode[STIX]{x1D702}_{f})$, the Reynolds number based on the Taylor microscale 
$Re_{\unicode[STIX]{x1D706}}=(\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D706}u^{\prime }/\unicode[STIX]{x1D702}_{f})$, the particle Reynolds number 
$Re_{p}=(\unicode[STIX]{x1D70C}_{f}d_{p}\unicode[STIX]{x0394}u/\unicode[STIX]{x1D702}_{f})$, the particle Stokes number 
$St=(\unicode[STIX]{x1D71A}_{p}d_{p}^{2}\bar{u}/18h\unicode[STIX]{x1D702}_{f})$, the mass loading 
$\unicode[STIX]{x1D719}_{m}$ (mass of particles per unit mass of fluid) and the volume fraction 
$\unicode[STIX]{x1D719}$ in the present study compared with select earlier studies. Here, 
$\unicode[STIX]{x1D70C}_{f}$ and 
$\unicode[STIX]{x1D702}_{f}$ are the fluid density and viscosity, 
$\bar{u}$ is the gas mean velocity, 
$\unicode[STIX]{x0394}u$ is the characteristic velocity difference between the particles and gas, 
$\unicode[STIX]{x1D71A}_{p}$ is the particle material density, 
$h$ is the channel width and 
$d_{p}$ is the particle diameter.
In contrast to earlier studies, the focus of the present study is to re-examine the conventional wisdom of a continuous decrease in the turbulence intensities with the particle volume fraction, and whether the decrease in the turbulence intensities is a result of an increase in the dissipation due to the drag force on the particles. Direct numerical simulation (DNS) of the Navier–Stokes equations is used for the fluid turbulence, and discrete particles are simulated using Newton’s laws of motion. The fluid–particle interaction is modelled using drag and lift forces based on the difference between the particle and fluid velocities at the particle location.
The particle size is comparable to the Kolmogorov scale for the fluid turbulence, and the point-particle approximation is used for the force exerted by the particles on the fluid. Instead of the ratio of the particle diameter and Kolmogorov scale, recent studies (Balachandar Reference Balachandar2009; Balachandar & Eaton Reference Balachandar and Eaton2010) suggest that the point-particle approximation can be used when the ratio of the Kolmogorov time and the particle relaxation is small. The ratio of the particle relaxation time and the Kolmogorov time, which is a Stokes number based on the Kolmogorov time, is 
$(\unicode[STIX]{x1D70C}_{p}/18\unicode[STIX]{x1D70C}_{f})(d_{p}/\unicode[STIX]{x1D702})^{2}$ when the mass density of the fluid is much larger than that of the particle, where 
$(d_{p}/\unicode[STIX]{x1D702})$ is the ratio of the particle diameter and the Kolmogorov scale 
$\unicode[STIX]{x1D702}$. Mehrabadi et al. (Reference Mehrabadi, Horwitz, Subramaniam and Mani2018) compared fully resolved particle simulations with point-particle simulations, and found that, even for 
$(d/\unicode[STIX]{x1D702})=1$, the results for the energy dissipation rate and the kinetic energy are accurately predicted by the point-particle simulation provided the Stokes number based on the Kolmogorov time is 
$O(100)$. They also found that it is essential to correct for the undisturbed velocity at the particle position while calculating the particle drag. The turbulent Stokes number based on the Kolmogorov scale in our study is 
$O(100)$; therefore, we have used the point-particle approximation along with the algorithm of Esmaily & Horwitz (Reference Esmaily and Horwitz2018) for the undisturbed fluid velocity at the particle position while calculating the drag.
While it is desirable to use the most accurate force model for the particle force, it should be noted that particle force models are still evolving. The force models have become relatively advanced in incorporating the correction for the inertial effects at high Reynolds number (Naumann & Schiller Reference Naumann and Schiller1935), the correction for the ‘undisturbed’ velocity at the particle centre (Esmaily & Horwitz Reference Esmaily and Horwitz2018; Mehrabadi et al. Reference Mehrabadi, Horwitz, Subramaniam and Mani2018), lift due to the fluid strain rate and the relative velocity between particle and fluid (Saffman Reference Saffman1965; Wang et al. Reference Wang, Squires, Chen and McLaughlin1997; Zhang & Ahmadi Reference Zhang and Ahmadi2000) and the correction to the drag and lift forces due to the presence of a nearby wall (Zeng et al. Reference Zeng, Najjar, Balachandar and Fischer2009). The exact regularised point particle (ERPP), a rigorous method for the regularisation of the singularities in a point-particle description, has been developed (Gualtieri, Battista & Casciola Reference Gualtieri, Battista and Casciola2017; Battista et al. Reference Battista, Gualtieri, Mollicone and Casciola2018) and extended to incorporate the effect of a nearby wall (Battista et al. Reference Battista, Mollicone, Gualtieri, Messina and Casciola2019). A comparison of simulations using the point-particle approximation and experimental results by Wang et al. (Reference Wang, Fong, Coletti, Capecelatro and Richter2019) showed that the results are in good agreement when the volume loading is 
$3\times 10^{-6}$, but the particle clustering is under-predicted at a higher volume loading of 
$5\times 10^{-5}$. More recently, Costa, Brandt & Picano (Reference Costa, Brandt and Picano2020) carried out a comparison of point-particle simulations with fully resolved particle simulations for inertial particles in a turbulent channel flow. They found that the turbulence statistics are well captured by the point-particle simulations in the bulk, but it is necessary to include the effects of the Saffman lift force in order to capture the near-wall dynamics. Moreover, the Saffman expression for the lift force provides a better prediction in comparison to more sophisticated lift models.
Despite the large amount of work done, further advances are expected in the force model, which currently do not incorporate the corrections to the torque on a particle due to inertia, the undisturbed vorticity at the particle centre or due to a nearby wall. Advances in computation will result in further refinement of the models that are currently available. The objective is to uncover phenomena that are independent of the force model used, rather than do a computation with a force model which is currently the most sophisticated but which may be further refined in the future. The approach here is to study turbulence modification at the large scale using different models for the particle force, and examine whether the collective effect of particles on turbulence modification is qualitatively the same and independent of the force model. Therefore, multiple proposed models for the drag and lift force on the particles have been used in the computations; these are described in § 2.2. The qualitative and quantitative natures of the turbulence attenuation due to these models are examined. The simulation procedures for the fluid and particles are discussed in the next section, followed by a discussion of the results in § 3. The important conclusions are summarised in § 4.
2 Formulation
2.1 Fluid equations
The fluid is an incompressible Newtonian fluid described using the Navier–Stokes mass and momentum equations,
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\frac{1}{\unicode[STIX]{x1D70C}_{f}}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\frac{1}{\unicode[STIX]{x1D70C}_{f}}\boldsymbol{f}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{u}(\boldsymbol{x},t)$ is the velocity, 
$\unicode[STIX]{x1D70C}_{f}$ is the fluid density, 
$p(\boldsymbol{x},t)$ is the pressure field, 
$\unicode[STIX]{x1D708}$ is the kinematic viscosity and 
$\boldsymbol{f}(\boldsymbol{x},t)$ is the reverse force on the fluid exerted by the particles, which is discussed in § 2.2. In all equations, bold type is used to represent vectors. The volume fraction of particles is of 
$O(10^{-4})$, so the volume occupied by the particle is neglected in the fluid mass conservation equation. Hydrodynamic interactions between particles are not explicitly included in the formulation. However, the point force due to a particle in the fluid momentum equations will result in an effect on neighbouring particles, mediated by the fluid flow, which is turbulent. The force density due to the particles, 
$\boldsymbol{f}(\boldsymbol{x},t)$, is the negative of the sum of the drag and lift forces on the particles,
$$\begin{eqnarray}\boldsymbol{f}(\boldsymbol{x},t)=-\mathop{\sum }_{I}(\boldsymbol{F}_{I}^{D}+\boldsymbol{F}_{I}^{L})\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{I}),\end{eqnarray}$$ where 
$\boldsymbol{x}_{I}$ is the position of particle 
$I$, 
$\boldsymbol{F}_{I}^{D}$ and 
$\boldsymbol{F}_{I}^{L}$ are the drag and lift forces on the particle 
$I$, specified later in § 2.2, 
$m_{p}$ is the particle mass and 
$\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{I})$ is the Dirac delta function in three dimensions.
Brief details of the simulation procedure are provided in appendix A, the interpolation procedure for calculating the fluid velocity at the particle location is discussed in § A.1 and the projection procedure for the particle force on to the fluid grid points is explained in § A.2. The standard methods for the interpolation of the particle force on the fluid grid points are the particle-in-cell method, where the force due to all the particles in a (usually cubic) control volume around a grid point is imposed at the grid point (Crowe Reference Crowe1982; Boivin, Simonin & Squires Reference Boivin, Simonin and Squires1998) and the projection on nearest neighbours (PNN) method, where the particle force is projected onto the nearest neighbouring grid points (Squires & Eaton Reference Squires and Eaton1990; Elghobashi & Truesdell Reference Elghobashi and Truesdell1993). Here, a variant of the PNN method is used, where three coordinate planes with origin at the particle location are used to divide the control volume into eight cuboids, and the fraction of the force projected onto a grid point is the ratio of the volume of the cuboid opposite to the grid point and the total volume. This has the advantage that the force on a fluid grid point varies continuously as the particle crosses the surface of a control volume. The details of the implementation are provided in § A.2, and the results of the projection procedure are validated against the results obtained using a delta function at the particle location.
2.2 Particle equations
Discrete spherical particles are simulated using Newton’s law for the evolution of the particle velocity,
$$\begin{eqnarray}m_{p}\frac{\text{d}\boldsymbol{v}_{I}}{\text{d}t}=\boldsymbol{F}_{I}^{D}+\boldsymbol{F}_{I}^{L}+\mathop{\sum }_{I\neq J}\boldsymbol{F}_{IJ}+\boldsymbol{F}_{Iw}+m_{p}\boldsymbol{g},\end{eqnarray}$$ where 
$\boldsymbol{v}_{I}$ is the particle velocity, 
$m_{p}$ is the particle mass, 
$\boldsymbol{F}_{IJ}$ is the interaction force exerted on particle 
$I$ by particle 
$J$, 
$\boldsymbol{F}_{Iw}$ is the force exerted on the particle 
$I$ due to the wall in a particle–wall interaction, 
$\boldsymbol{g}$ is the gravitational acceleration and 
$\boldsymbol{F}_{I}^{D}$ and 
$\boldsymbol{F}_{I}^{L}$ are the drag and lift forces on the particle. The drag force is given by the Stokes drag law at small Reynolds number,
$$\begin{eqnarray}\boldsymbol{F}_{I}^{D}=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{f}d_{p}(\boldsymbol{u}(\boldsymbol{x}_{I},t)-\boldsymbol{v}_{I}).\end{eqnarray}$$Several nonlinear equations proposed earlier have been used to study the drag and lift forces on a particle.
(i) The most commonly used drag law that incorporates the correction due to inertia at relatively low Reynolds number, called the Schiller–Naumann correlation (Naumann & Schiller Reference Naumann and Schiller1935), is of the form
$$\begin{eqnarray}\boldsymbol{F}_{I}^{D}=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{f}d_{p}(\boldsymbol{u}_{f}(\boldsymbol{x}_{I},t)-\boldsymbol{v}_{I})(1+0.15Re_{p}^{0.687}),\end{eqnarray}$$ where 
$d_{p}$ is the particle diameter and the particle Reynolds number 
$Re_{p}=(\unicode[STIX]{x1D70C}_{f}|\boldsymbol{u}_{f}-\boldsymbol{v}_{I}|d_{p}/\unicode[STIX]{x1D702}_{f})$.
(ii) In earlier studies, the velocity 
$\boldsymbol{u}_{f}(\boldsymbol{x}_{I},t)$ in (2.6) was considered to be the fluid velocity at the particle location. It has now been realised that the appropriate velocity is the ‘undisturbed’ fluid velocity at the particle position (Esmaily & Horwitz Reference Esmaily and Horwitz2018; Mehrabadi et al. Reference Mehrabadi, Horwitz, Subramaniam and Mani2018), which is the far-field velocity extrapolated to the particle centre in the absence of the particle. The undisturbed fluid velocity at the particle location is determined here using the scheme of Esmaily & Horwitz (Reference Esmaily and Horwitz2018). The undisturbed fluid velocity 
$\boldsymbol{u}_{f}$ is defined as the difference between the actual interpolated fluid velocity 
$\boldsymbol{u}(\boldsymbol{x}_{I},t)$ and a ‘cell’ velocity 
$\boldsymbol{u}_{c}^{(i)}$, which is the velocity of a volume of fluid of size comparable to the computational cell 
$i$ subjected to a point force. The cell velocity 
$\boldsymbol{u}_{c}^{(i)}$ is determined using the evolution equation
$$\begin{eqnarray}\frac{3m_{c}}{2}\frac{\text{d}\boldsymbol{u}_{c}^{(i)}}{\text{d}t}=-3\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{f}d_{c}K_{t}^{(i)}\boldsymbol{u}_{c}^{(i)}-\boldsymbol{F}^{(i)},\end{eqnarray}$$ where 
$m_{c}$ is the mass of fluid in the cell, 
$d_{c}$ is the effective cell diameter, 
$\boldsymbol{F}^{(i)}$ is the sum of the fluid drag and lift forces in cell 
$i$ due to the particle and 
$K_{t}^{(i)}$ is a correction factor for the drag force, which is expressed as 
$K_{t}^{(i)}=(K_{p}^{(i)}C_{r}/K_{p}^{(i)}C_{t}^{(i)})$. Here, the correction factor 
$K_{c}^{(i)}$ accounts for the non-spherical shape of the computational cell, 
$C_{r}$ incorporates the correction to the Stokes drag law for the cell due to finite Reynolds number, 
$K_{p}^{(i)}$ accounts for the distribution of the drag force among different computational cells and 
$C_{t}^{(i)}$ is a correction for the finite residence time of the particle in the cell. The algorithm provided in Esmaily & Horwitz (Reference Esmaily and Horwitz2018) is used here for determining the undisturbed fluid velocity.
(iii) When the particles are close to solid walls, there is a correction to the drag force on the particle due to the presence of the wall. Two kinds of corrections were developed for the drag force calculated in Zeng et al. (Reference Zeng, Najjar, Balachandar and Fischer2009), that due to a stationary particle in a shear flow, and due to a particle moving parallel to a wall in a stationary fluid. The two expressions are nearly identical when the distance between the particle surface and the wall is greater than approximately 
$0.6d_{p}$, but the latter is much larger than the former when the distance from the wall is small. In order to use the larger estimate for the effect of a solid wall on the drag force, the following expression for the force due to a particle moving parallel to a stationary wall from Zeng et al. (Reference Zeng, Najjar, Balachandar and Fischer2009) is used:
$$\begin{eqnarray}\displaystyle \boldsymbol{F}_{I}^{D} & = & \displaystyle 3\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{f}d_{p}(\boldsymbol{u}_{f}(\boldsymbol{x}_{I},t)-\boldsymbol{v}_{I})(1+\unicode[STIX]{x1D6FC}Re_{p}^{\unicode[STIX]{x1D6FD}})\nonumber\\ \displaystyle & & \displaystyle \times \left(1.028-\frac{0.07}{1+4(\unicode[STIX]{x1D6FF}/d_{p})^{2}}-\frac{8}{15}\ln \left(\frac{270(\unicode[STIX]{x1D6FF}/d_{p})}{135+256(\unicode[STIX]{x1D6FF}/d_{p})}\right)\right),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}$ is the distance between the wall and the location on the particle surface closest to the wall, and 
$\unicode[STIX]{x1D6FC}$ and 
$\unicode[STIX]{x1D6FD}$ are parameters which depend on 
$\unicode[STIX]{x1D6FF}$ and reduce to 0.15 and 0.687 for 
$\unicode[STIX]{x1D6FF}\gg 1$,
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D6FC}=0.15(1-\exp (-\sqrt{\unicode[STIX]{x1D6FF}/d_{p}})),\\ \unicode[STIX]{x1D6FD}=0.687+0.313\exp (-2\sqrt{\unicode[STIX]{x1D6FF}/d_{p}}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$(iv) The Saffman lift force (Saffman Reference Saffman1965; Wang et al. Reference Wang, Squires, Chen and McLaughlin1997; Zhang & Ahmadi Reference Zhang and Ahmadi2000) is
$$\begin{eqnarray}F_{I}^{L}=\frac{1.61\unicode[STIX]{x1D702}_{f}d_{p}^{2}(\text{d}u_{x}/\text{d}y)|\boldsymbol{v}_{I}-\boldsymbol{u}_{f}|}{|\text{d}u_{x}/\text{d}y|^{1/2}\unicode[STIX]{x1D708}^{1/2}},\end{eqnarray}$$ where 
$u_{x}$ is the streamwise velocity, 
$y$ is the cross-stream direction and 
$\unicode[STIX]{x1D708}$ is the fluid kinematic viscosity. There are two important correction factors that have been proposed for the Saffman lift force. The first correction factor is due to the Reynolds number 
$Re_{G}$ based on the strain rate and particle diameter, in contrast to the particle Reynolds number 
$Re_{p}$ which is proportional to the relative velocity between the particle and fluid (Wang et al. Reference Wang, Squires, Chen and McLaughlin1997). The lift force depends on the parameter (
$Re_{G}^{1/2}/Re_{p}$); it decreases to zero for 
$(Re_{G}^{1/2}/Re_{p})\ll 1$ and has the maximum value expressed in (2.10) for 
$(Re_{G}^{1/2}/Re_{p})\gg 1$. In the present analysis, the Reynolds number 
$Re_{G}$ based on the maximum strain rate at the wall is approximately 17, and so the ratio 
$(Re_{G}^{1/2}/Re_{p})$ is less than 0.1. In this case, the correction factor is proportional to 
$-32\unicode[STIX]{x03C0}^{2}(Re_{G}^{1/2}/Re_{p})^{5}\log (Re_{p}^{2}/Re_{G})$ (McLaughlin Reference McLaughlin1993; Cherukat, McLaughlin & Graham Reference Cherukat, McLaughlin and Graham1994), which is approximately 1.5 % of the Saffman lift. Therefore, the lift force is primarily due to the translation of the particle relative to the fluid.
(v) There is a second correction to the lift force due to the presence of a nearby wall. As discussed above, for 
$(Re_{G}^{1/2}/Re_{p})\ll 1$, the lift due to particle translation is dominant. Therefore, the expression for the lift force for a particle translating parallel to the wall (Zeng et al. Reference Zeng, Najjar, Balachandar and Fischer2009) has been used,
$$\begin{eqnarray}\displaystyle F_{I}^{L} & = & \displaystyle \frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}|\boldsymbol{v}_{p}|^{2}d_{p}^{2}}{8} (f(y/d_{p},Re_{p})\nonumber\\ \displaystyle & & \displaystyle +\,(0.313+0.812\exp (-0.125Re_{p}^{0.77})-f(0.5,Re_{p}))\text{e}^{(-11(\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{c})^{1.2})}\! ),\qquad\end{eqnarray}$$where
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}f(y/d_{p},Re_{p})=(1+0.329Re_{p}+0.00485Re_{p}^{2})C_{0}(Re_{p}L/d_{p})(L/d_{p})^{[-0.9\tanh (0.0022Re_{p})]},\\ \unicode[STIX]{x1D6FF}_{c}=3\exp (-0.17Re_{p}^{0.7}),\\ \begin{array}{@{}rcl@{}}C_{0}(Re_{p}L/d_{p})\ & =\ & (9/8+5.79\times 10^{-6}(Re_{p}L/d_{p})^{4.58})\\ \ & \ & \times \exp (-0.292Re_{p}L/d_{p})\quad \text{for}~0\leqslant Re_{p}L/d_{p}\leqslant 10,\\ \ & =\ & 8.94(Re_{p}L/d_{p})^{-2.09}\quad \text{for}~Re_{P}L/d_{p}>10.\end{array}\end{array}\right\}\end{eqnarray}$$The particles are considered as smooth rigid hard spheres, and the particle interactions are modelled as instantaneous elastic collisions. In inter-particle collisions, the velocity of the centre of mass is unaltered. The component of the relative velocity (difference in the velocities of the two colliding particles) in the direction of the line joining the centres of the particles at impact is reversed, and the component of the relative velocity perpendicular to the line joining centres is unchanged. It can be shown that these relations between the pre- and post-collisional velocities conserve momentum and energy. Inelastic collision models, where the magnitude of the post-collisional velocity is less than that of the pre-collisional velocity, have also been widely used. Here, the focus is on the turbulence modification due to the particle drag, and so the simplest elastic collision model has been used. The specular reflection model is used for particle–wall collisions, where the component of the velocity perpendicular to the wall is reversed in a collision, while the component of the velocity parallel to the wall is unchanged.
The particle equations are solved using an explicit time-stepping method with the same time increment as that for the fluid velocity field. An event-driven molecular dynamics algorithm (Allen & Tildesley Reference Allen and Tildesley2017) is used for the particle collisions, where the location of an impending collision is predicted based on the current positions and velocities of all the particles, and the simulation is advanced in time to execute the impending collision. The total number of particles in the simulation is proportional to the volume fraction of particles, and 8000 particles are used when the volume fraction is 
$10^{-3}$.
2.3 Flow configuration
The configuration used for the simulations, shown in figure 1, consists of a rectangular box of length 
$4\unicode[STIX]{x03C0}h$ in the flow (
$x$) direction, 
$h$ in the cross-stream (
$y$) direction and 
$(2\unicode[STIX]{x03C0}h/3)$ in the spanwise (
$z$) direction. Periodic boundary conditions are used in the flow and the spanwise direction, and there are solid walls in the 
$y$ direction. The gas flow is driven by a pressure gradient in the 
$x$ direction, and the pressure gradient is adjusted so that the average gas flow velocity 
$\bar{u}$ is the same in all the simulations. The Reynolds number 
$Re_{b}$ for the unladen flow is defined as 
$(\unicode[STIX]{x1D70C}_{f}\bar{u}h/\unicode[STIX]{x1D702}_{f})$, where 
$\unicode[STIX]{x1D70C}_{f}$ and 
$\unicode[STIX]{x1D702}_{f}$ are the fluid density and viscosity. The Reynolds number for the unladen flow is maintained at a fixed value of 
$Re_{b}=3300$ and 
$Re_{b}=5600$ for the two sets of simulations. When there is a change in the wall shear stress due to the addition of particles, the pressure gradient is adjusted so that the channel Reynolds number unchanged.
Figure 1. Configuration and coordinate system used for the simulations.
In order to resolve the smallest scales at 
$Re_{b}=3300$, 128 Fourier modes are used in the flow direction, 
$64$ Fourier modes in the spanwise direction and 65 Chebyshev modes in the cross-stream (
$y$) direction. For 
$Re_{b}=5600$, 192 Fourier modes are used in the flow direction, 160 Fourier modes in the spanwise direction and 129 Chebyshev modes in the cross-stream (
$y$) direction have been used. The time step in the simulations is 
$(0.0033h/\bar{u})$, where 
$h$ is the channel width and 
$\bar{u}$ the average velocity. An explicit scheme (for nonlinear term) is used for the simulations; the time step has been decided on the basis of the Courant–Friedrichs–Lewy condition for convergence. The simulation procedure is briefly explained in appendix A, and further details of the numerical scheme are provided in Goswami & Kumaran (Reference Goswami and Kumaran2011a). The results for the unladen flow are validated against the results of Goswami & Kumaran (Reference Goswami and Kumaran2011a) for a Reynolds number of 3300, and against the results of Kim, Moin & Moser (Reference Kim, Moin and Moser1987) for the Reynolds number 5600.
Although the maximum particle volume fraction considered in the simulations is 
$3.5\times 10^{-3}$, there is a significant increase in the mass density because the material density of the particles is three orders of magnitude greater than that of the fluid. For the particle-laden flow, it is necessary to define the Reynolds number based on the suspension density and viscosity. Based on the Einstein correction, the difference between the suspension and fluid viscosity is less than 0.5 % of the fluid viscosity for the range of volume fractions considered here, and so the fluid viscosity is used in the definition of the Reynolds number. However, the suspension mass density does increase significantly, and it is necessary to incorporate this correction in the definition of the Reynolds number. The corrected Reynolds number is then 
$(1+(\unicode[STIX]{x1D719}_{av}\unicode[STIX]{x1D71A}_{p}/\unicode[STIX]{x1D70C}_{f}))Re_{b}$, where 
$\unicode[STIX]{x1D71A}_{p}$ is the material density of the material comprising the particles, and the notation 
$\unicode[STIX]{x1D719}_{av}$ is used for the average volume fraction, to distinguish it from the local volume fraction denoted by 
$\unicode[STIX]{x1D719}$.
The ratio of the particle diameter and the channel width is set equal to 
$1.845\times 10^{-2}$ in all the simulations. The particle terminal velocity depends on the particle material density and the drag law. For the Stokes drag law, the ratio of the particle terminal velocity and the bulk gas velocity is 
$2.06\times 10^{-4}$ for the lowest particle material density 
$\unicode[STIX]{x1D71A}_{p}=100~\text{kg}~\text{m}^{-3}$, and 
$1.65\times 10^{-2}$ for the highest particle material density of 8000 kg m-3. Thus, the particle terminal velocity is always much smaller than the bulk gas velocity. The particle Reynolds number based on the particle diameter, terminal velocity and gas kinematic viscosity varies between 
$1.3\times 10^{-2}$ for particle material density 100 kg m-3 and 1.01 for particle material density 8000 kg m-3. Thus, the particle Reynolds number based on the terminal velocity of the particles is small. The results of the simulation show that the particle mean velocity is very nearly a constant across the channel, whereas the unladen fluid mean velocity for the turbulent flow is plug like at the centre, with sharp gradients at the walls. The particle Reynolds number 
$(\unicode[STIX]{x1D70C}_{f}\bar{u}d_{p}/\unicode[STIX]{x1D702}_{f})$ is approximately 60 when the channel Reynolds number is 3300, and approximately 103 when the channel Reynolds number is 5600. However, the particle Reynolds number based on the channel-averaged root mean square velocity difference between the particles and fluid is in the range 4–15. The Stokes drag law (2.5) is in error by approximately 50 % at this 
$Re_{p}=60$, but the Schiller–Naumann correlation (2.6) is in error by approximately 1.2 % at 
$Re_{p}=60$ and approximately 2 % at 
$Re_{p}=103$. The fluid time scale is estimated as 
$(h/\bar{u})$, the ratio of the channel width and the average flow velocity. The particle viscous relaxation time is
$$\begin{eqnarray}t_{v}=(\unicode[STIX]{x1D71A}_{p}d_{p}^{2}/18\unicode[STIX]{x1D702}_{f}),\end{eqnarray}$$for the Stokes drag law and
$$\begin{eqnarray}t_{v}=(\unicode[STIX]{x1D71A}_{p}d_{p}^{2}/18\unicode[STIX]{x1D702}_{f}(1+0.15Re_{p})^{0.687}),\end{eqnarray}$$ for the Schiller–Naumann correlation (2.6). The viscous relaxation time is provided as a function of the particle material density in table 2. The particle Stokes number, 
$(t_{v}/t_{f})$ varies over nearly two orders of magnitude, from approximately 5 to approximately 420 for the Stokes drag law, and from approximately 1.5 to approximately 150 for the Schiller–Naumann correlation (2.6). The other time scales of importance are the average time between particle–particle collisions, 
$t_{pp}$, and the average time between particle–wall collisions, 
$t_{pw}$. These depend on the particle loading, and the times determined from the simulations are shown in figure 21 in § 3. Here, it suffices to note that we have covered all parameter regimes, where the viscous relaxation time is smaller and larger than the collision time, and where the particle–particle collision time is smaller and larger than the particle–wall collision time.
Table 2. The different values of the Stokes number calculated for the Stokes drag law (2.5) with the viscous relaxation time given by (2.13), and the Schiller–Naumann correlation (2.6) using the viscous relaxation time in (2.14). The particle Reynolds number 
$Re_{p}$ is set equal to 60 for channel Reynolds number 3300, and 103 for channel Reynolds number 5600 in the calculation of the nonlinear drag coefficient in (2.14).
The Kolmogorov scale for the unladen flow is estimated as 
$(\unicode[STIX]{x1D708}^{3}/(\unicode[STIX]{x1D716}/\unicode[STIX]{x1D70C}))^{1/4}$, where 
$\unicode[STIX]{x1D708}=(\unicode[STIX]{x1D702}_{f}/\unicode[STIX]{x1D70C}_{f})$ is the kinematic velocity, and 
$\unicode[STIX]{x1D716}$ is the rate of dissipation of energy per unit volume due to the turbulent velocity fluctuations. For the unladen flow, the average value of 
$\unicode[STIX]{x1D716}$, determined by averaging across the channel width, is approximately 
$4\times 10^{-3}(\unicode[STIX]{x1D70C}_{f}\bar{u}^{3}/h)$. Therefore, the Kolmogorov scale is 
$(4\times 10^{-3}Re_{b}^{3})^{1/4}$, which is approximately 
$1.6\times 10^{-2}$ times the channel width. This is comparable to the particle diameter, and so the point-force approximation is used for the force on the fluid due to the particles.
It is important to note that, in order to generate a turbulent flow in DNS simulations for the fluid, it is necessary to impose a relatively large perturbation on the fluid velocity profile, otherwise the fluctuations decay and the velocity profile reverts to a laminar profile. In the present simulations, perturbations are introduced in the Fourier–Chebyshev coefficients of the velocities following the standard procedure of Gibson, Halcrow & Cvitanovi (Reference Gibson, Halcrow and Cvitanovi2008). For each point in Fourier–Chebyshev space, in the two Fourier directions, the perturbation is introduced in the form of a random number in the interval -0.1 to 0.1 at the starting time. The perturbation to the Chebyshev coefficient is determined from the incompressibility condition. This amplitude of the perturbation repeatably generates a turbulent steady state after a development time of approximately 
$50(h/\bar{u})$, where 
$h$ is the channel width and 
$\bar{u}$ is the average velocity.
For the simulation of a particle-laden flow, the instantaneous velocity for a steady unladen flow is taken as the initial condition, and the particles are added at random locations with the same velocity as the interpolated fluid velocity at that location. The flow is evolved for approximately 
$1000(h/\bar{u})$ for the particle statistics to reach steady state for a particle-laden flow. Averages are then calculated for a further 
$300(h/\bar{u})$. If necessary, small adjustments are made to the pressure gradient so that the average gas flow velocity is equal to that for the unladen flow. In the primary results in figures 2, 3 and 20, the mean values and the standard deviations in the profiles of the mean and root mean square fluctuating velocities are shown. The standard deviations are estimated by carrying out three independent simulation runs for the same particle volume fraction but with different random initial particle locations and fluid perturbations. The error bars show one standard deviation below and above the mean value.
3 Results
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
$$\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}$,
$$\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.
(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.
(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.(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.
(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.
(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
$$\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}$$
$$\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,
$$\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
$$\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
$$\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
$$\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
$$\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
$$\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}$$
$$\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.
4 Conclusions
In the present study, we have investigated the phenomenon of turbulence collapse for two different fluid Reynolds numbers and over a range of particle Stokes numbers. The study reveals a hitherto unknown facet of turbulent particle–gas suspensions, which is the discontinuous decrease in the turbulence intensities at a critical volume fraction. There is a reduction of 1–2 orders of magnitude in all components of the second moment of the turbulent velocity fluctuations when the volume fraction exceeds a critical value. This transition is robust in many ways. For a steady flow with volume fraction above the critical value, upon imposition of velocity perturbations of relatively large magnitude (50 % larger than that used to generate turbulence in DNS in an unladen flow), the perturbation decays and the velocity profile spontaneously returns to its initial value. When the volume fraction is below the critical value, the fluid velocity statistics revert to those for an unladen turbulent flow when the particles are deleted. In contrast, when the volume fraction is above the critical value, the flow reverts to the laminar parabolic profile when the particles are deleted. This suggests that, above the critical volume fraction, the flow is essentially a laminar flow with superposed fluid velocity fluctuations generated due to the force exerted by the particles.
The phenomenon of turbulence collapse is remarkably robust and has been found to occur irrespective of the drag law. This has been tested using both the linear Stokes drag law and the Schiller–Naumann correlation (2.6) and the turbulence collapse phenomenon is observed in both cases, although the critical volume fraction does depend on the drag law. It is observed that there is very little change in the fluid statistics, and no significant change in the critical volume fraction, due to the inclusion of the correction for the undisturbed velocity at the particle centre, the lift force and the effect of a nearby wall. Thus, this phenomenon is not dependent on the detailed models for the drag and lift forces on the particles, which are likely to become more refined as time progresses, but is a basic characteristic of turbulent particle–gas suspensions.
There is no comparable change in the particle fluctuation intensities at the critical volume fraction. The particle concentration and mean velocity show no change despite the large decrease in the fluid turbulence. The particle mean square velocities show an increase, in the range of 10 %–50 %, at the critical volume fraction, in contrast to the decrease of 1–2 orders of magnitude in the fluid fluctuating velocities. There is also a small but discernible decrease in the particle–particle and particle–wall collision times, commensurate with the increase in the particle fluctuating velocities. There is virtually no change in the distribution functions for the particle fluctuating velocities. Thus, the particle dynamics does not undergo a significant change at the critical volume fraction.
An important observation is that the magnitudes of the particle velocity fluctuations are smaller than those of the fluid velocity fluctuations below the critical volume fraction, whereas the former are larger than the latter above the critical volume fraction. This indicates that the fluid velocity fluctuations are causing the particle velocity fluctuations via the drag coupling below the critical volume fraction, whereas the particle fluctuations are driving the fluid fluctuations above the critical volume fraction.
One of the important limitations of the present study is that it has been conducted at two Reynolds numbers, 3300 and 5600. Even though these are relatively low, they are of importance in practical applications like pneumatic conveying (Goswami & Kumaran Reference Goswami and Kumaran2011b). Due to the relatively low Reynolds number, it was feasible to carry out a sufficiently large number of simulations in order to detect the discontinuous changes in the turbulence intensities – there have been five drag laws used here, for each drag law there have been six different particle Stokes numbers used, at least five volume fractions have been examined for each Stokes number to detect the critical volume fraction and three different production runs have been carried out for each volume fraction to ensure repeatability, resulting in approximately 500 different simulation runs. Moreover, the simulation time required for a particle-laden flow is significantly higher than that for an unladen flow. Whereas an unladen flow repeatably reaches the turbulent steady state within approximately 
$(50h/\bar{u})$, it is necessary to run the simulation for approximately 
$(1000h/\bar{u})$ for the particle statistics to reach steady state for a particle-laden flow, followed by averaging for another 
$(300h/\bar{u})$.
The second important limitation is the assumption of a steady drag force on the particles. Although the particle Reynolds numbers based on particle diameter and average fluid mean velocity are approximately 60 and 103, the particle Reynolds number based on the root mean square difference between the fluid and particle velocity is in the range 4–15. For particle Reynolds numbers greater than 24, there is a separation bubble at the rear of the sphere but no vortex shedding, and inertial effects cannot be neglected. One important objective of the present analysis was to examine whether the turbulence collapse phenomenon does depend on the specific form of the drag law used. For this reason, simulations have been carried out using two different forms of the drag law, the Stokes drag law which is applicable in the absence of inertia, and the Schiller–Naumann correlation (2.6). The force on the fluid due to the particle has been modelled as a point force, and the particle shape has not been resolved in the simulations. In order to examine the accuracy of the point-force approximation, an attempt has been made to include the symmetric and antisymmetric point dipoles due to the particles in a related study (Tyagi Reference Tyagi2017). The results of that study are qualitatively similar to those from the present simulations, although small quantitative differences less than 10 % have been observed. Another important limitation is that the channel width is sufficiently small that the particle volume fraction and the mean velocity are nearly constant across the channel. The small channel width, approximately 60 times the particle diameter, was necessary to limit the number of particles to approximately 36 000 for the largest volume fraction, so that hard-particle simulations that include the effect of particle collisions could be carried out in reasonable time. However, as shown in figure 21, the range of parameters considered includes cases where the collision time is smaller and larger than the viscous relaxation time of a particle, as well as cases where the particle–particle collision time is smaller and larger than the particle–wall collision time. The particle Stokes number has been varied over two orders of magnitude.
An analysis of the fluid momentum equation indicates that there is a sharp decline, by a factor of 10 or more, in the divergence of the Reynolds stress. This is accompanied by a small increase in the drag force exerted on the fluid by the particles. Similarly, in the energy equation, there is a dramatic decrease in the turbulent energy production rate at the transition volume fraction. There is a small decrease in the energy dissipation rate due to the mean flow, and a small increase in the dissipation rate due to the particle drag. It is significant to note that the increase in the energy dissipation rate due to the particles is only approximately one half of the decrease in the turbulent energy production. This results in a decrease in the total energy dissipation rate in the fluid at the critical volume fraction. Moreover, there is no increase in the energy dissipation rate due to particle drag in the zone of maximum turbulence production, as shown in figure 14(b); rather, there is an increase in the dissipation rate in the centre of the channel due to the increase in the difference between the particle and fluid mean velocities. This clearly indicates that the turbulence collapse is not due to the increased dissipation by the particles, but rather due to the disruption of the turbulent production by the particles. This is at variance with the conventional wisdom in turbulent gas–particle flows, that the turbulence attenuation is due to the increase in the energy dissipation due to the particles. Rather, the present study shows that the turbulence attenuation is due to the collapse in the Reynolds stress and the turbulent production rate in the fluid. This collapse occurs discontinuously at a critical volume fraction, in contrast to the current belief that there is a steady decline in the turbulence intensity. The exact mechanism of this disruption is the subject of further study.
Acknowledgements
The authors would like to thank the Department of Science and Technology, Government of India for financial support. V.K. would like to thank the JRD Tata Trust for financial support.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Direct numerical simulations
The primitive variable formulation along with the coupled method has been used to integrate Navier–Stokes equation using the Kleiser–Schumann algorithm (Kleiser & Schumann Reference Kleiser and Schumann1980; Canuto, Hussaini & Zang Reference Canuto, Hussaini and Zang1988). The numerical formulation for the fluid turbulence is adapted from chapter 7 of Canuto et al. (Reference Canuto, Hussaini, Quarteroni and Zhang2007). The procedure has been validated earlier in Goswami (Reference Goswami2008) for one-way coupling for wall-bounded two phase flow based on the open source incompressible Navier–Stokes solver Channelflow (Gibson et al. Reference Gibson, Halcrow and Cvitanovi2008). This has been extended to include two-way coupling between the particles and fluid. The simulation procedure is the same as that used in Goswami & Kumaran (Reference Goswami and Kumaran2011a). Here, a brief summary of the simulation procedure and validation is provided, followed by the interpolation procedure for calculating the fluid velocity at the particle location and the particle force on the fluid.
The instantaneous velocity field can be decomposed into the mean velocity and the fluctuating velocity,
$$\begin{eqnarray}\boldsymbol{u}(x,t)=U(y)\boldsymbol{e}_{x}+\boldsymbol{u}^{\prime }(x,t),\end{eqnarray}$$ where 
$x$ and 
$y$ are the streamwise and cross-stream directions, 
$U(y)$ is the mean (time-averaged) velocity and 
$\boldsymbol{e}_{x}$ is the unit vector in the streamwise (
$x$) direction. Similarly, the pressure field can be decomposed into the linearly varying mean pressure and the fluctuating pressure as
$$\begin{eqnarray}p(x,t)=\frac{\text{d}P}{\text{d}x}x+p^{\prime }(x,t).\end{eqnarray}$$The pressure gradient term can be written as
$$\begin{eqnarray}\unicode[STIX]{x1D735}p(x,t)=\frac{\text{d}P}{\text{d}x}\boldsymbol{e}_{x}+\unicode[STIX]{x1D735}p^{\prime }(x,t).\end{eqnarray}$$Substituting (A 1) and (A 2) in (2.2), we get
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}^{\prime }}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D70C}_{f}}\unicode[STIX]{x1D735}p^{\prime }=\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{\prime }-\left[\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }+U\frac{\unicode[STIX]{x2202}\boldsymbol{u}^{\prime }}{\unicode[STIX]{x2202}x}+u_{y}^{\prime }\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}\boldsymbol{e}_{x}\right]+\left[\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}y^{2}}-\frac{1}{\unicode[STIX]{x1D70C}_{f}}\frac{\text{d}P}{\text{d}x}\right]\boldsymbol{e}_{x}+\frac{\boldsymbol{f}}{\unicode[STIX]{x1D70C}_{f}}.\end{eqnarray}$$Different forms have been used for the nonlinear term in the momentum conservation equation (A 3). Here, we have used the rotational form of Kleiser & Schumann (Reference Kleiser and Schumann1980). Since this introduces errors at higher wavenumbers, the dealiasing procedure of Zang (Reference Zang1991) has been used. Crank–Nicholson time discretisation has been used for the linear terms in (A 3), and a second-order Adams–Bashforth scheme has been used for the nonlinear term.
The zero-velocity boundary conditions at the walls are enforced using the influence matrix method. After carrying out the Fourier transforms in the streamwise and spanwise directions, the momentum equations are reduced to two independent one-dimensional Helmholtz equations for the divergence of the velocity and the wall-normal component of the velocity. The inhomogeneous pressure boundary condition in the equation for the divergence of the velocity is adjusted in order to ensure that the tangential velocity boundary condition is satisfied. The procedure is explained in Kleiser & Schumann (Reference Kleiser and Schumann1980).
A.1 Fluid velocity interpolation
The fluid velocity is calculated at finite number of fluid grid points, but the calculation of fluid drag on a particle requires the interpolation of fluid velocities at the particle location. The most accurate method to calculate the interpolated velocities is direct summation of the Fourier and Chebyshev series, known as spectral interpolation, used, for example, in McLaughlin (Reference McLaughlin1989). Spectral interpolation requires a significant computational cost 
$O(N^{3}N_{p})$, where 
$N$ is the number of grid points in each direction and 
$N_{p}$ is the number of particles, so it not feasible to simulate large numbers of particles. Yeung & Pope (Reference Yeung and Pope1988) concluded that linear interpolation is not accurate to calculate the interpolated fluid velocity because there could be high velocity gradients near the wall. Yeung & Pope (Reference Yeung and Pope1988) and Balachandar & Maxey (Reference Balachandar and Maxey1989) studied different interpolation schemes such as linear interpolation, Lagrangian interpolation, Hermite and cubic interpolation for accuracy and found that the interpolation scheme should at least be third-order accurate in the spatial directions. In the case of wall-bounded flows, Chebyshev polynomials are used in the wall-normal direction to capture the near-wall physics. Since we are using a pseudo-spectral method, it is easy to implement Chebyshev polynomial interpolation in the wall-normal direction. For this reason, Kontomaris, Hanratty & McLaughlin (Reference Kontomaris, Hanratty and McLaughlin1992) studied mixed interpolation schemes like linear–Chebyshev, third-order Lagrangian–Chebyshev, third-order Lagrangian–Chebyshev, Hermite–Chebyshev etc. They found that Hermite–Chebyshev is the most accurate interpolation scheme but it requires the storage of the first and second derivatives of the three-dimensional velocity field in the periodic directions. This requires a large amount of computational memory. Garg et al. (Reference Garg, Narayanan, Lakehal and Subramaniam2007) have compared different interpolation schemes, and have concluded that a higher-order polynomial interpolation scheme could result in larger errors than a lower-order scheme in some cases, if there is no correction for the undisturbed velocity at the particle location. Since we do include the correction for the undisturbed velocity (Esmaily & Horwitz Reference Esmaily and Horwitz2018), and considering both accuracy and memory optimisation, we have used a fifth-order Lagrangian Chebyshev interpolation scheme in our study, as used by Goswami & Kumaran (Reference Goswami and Kumaran2010).
In the pseudo-spectral method, the fluid velocity can be written as,
$$\begin{eqnarray}\displaystyle \boldsymbol{u}^{\prime }(x,y,z,t)=\mathop{\sum }_{n_{x}=-N_{x/2}}^{N_{x/2}-1}\mathop{\sum }_{n_{z}=-N_{z/2}}^{N_{z/2}-1}\mathop{\sum }_{n_{y}=0}^{N_{y}-1}\hat{\tilde{\boldsymbol{u}}}(n_{x},n_{y},n_{z},t)\text{e}^{2\unicode[STIX]{x03C0}\imath [(n_{x}x/L_{x})+(n_{z}z/L_{z})]}T_{n}(y), & & \displaystyle\end{eqnarray}$$ where 
$N_{x}$, 
$N_{y}$ and 
$N_{z}$ are the total number of grid points in the streamwise, wall-normal and vorticity directions, all lengths are scaled by the half-width of the channel, 
$L_{x}$ and 
$L_{z}$ are the box lengths in the 
$x$ and 
$z$ directions. The summations over 
$n_{x}$ and 
$n_{z}$ are computationally expensive if the particles are not at the fluid grid points, since fast Fourier transforms cannot be used. In the present simulations, these have been carried out by the less computationally expensive Lagrangian interpolation,
$$\begin{eqnarray}\boldsymbol{u}(x,y,z,t)=\mathop{\sum }_{j=0}^{k}\boldsymbol{u}_{j}l_{j}(x)l_{j}(z),\end{eqnarray}$$ where 
$l_{j}$ are the Lagrange basis functions. Here, we have used fifth-order Lagrangian basis sets, 
$k=5$,
$$\begin{eqnarray}\displaystyle & \displaystyle l_{j}(x)=\mathop{\prod }_{0\leqslant m\leqslant k}\frac{x-x_{m}}{x_{j}-x_{m}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle l_{j}(z)=\mathop{\prod }_{0\leqslant m\leqslant k}\frac{z-z_{m}}{z_{j}-z_{m}}. & \displaystyle\end{eqnarray}$$ After summation over 
$x$ and 
$z$, we get the velocity in the 
$x$–
$y$ plane,
$$\begin{eqnarray}\displaystyle \boldsymbol{u}(x,y,z,t)=\mathop{\sum }_{n_{y}=0}^{N_{y}-1}\hat{\boldsymbol{u}}(x,n_{y},z,t)T_{n}(y). & & \displaystyle\end{eqnarray}$$The summation in (A 8) is carried out to obtain velocities at 36 wall-normal grid points surrounding the particle location. Then, fifth-order Lagrangian interpolation is used to obtain the velocity at the particle location.
A.2 Particle reverse force
In spectral method, the reverse force calculated at the particle centre is represented by delta functions, that is, the force is non-zero only at the particle centre, and is zero otherwise. The force is represented as the sum of forces exerted by individual particles,
$$\begin{eqnarray}\boldsymbol{f}=-\mathop{\sum }_{I=1}^{N_{p}}(\boldsymbol{F}_{I}^{D}+\boldsymbol{F}_{I}^{L})\unicode[STIX]{x1D6FF}(x-x_{I})\unicode[STIX]{x1D6FF}(y-y_{I})\unicode[STIX]{x1D6FF}(z-z_{I}),\end{eqnarray}$$ where 
$F_{I}^{D}$ and 
$F_{I}^{L}$ are the drag and lift forces on the particles (see (2.4)). The force in spectral space after Fourier–Chebyshev transformation is
$$\begin{eqnarray}\displaystyle \hat{\boldsymbol{f}}(n_{x},n_{y},n_{z},t)=-\mathop{\sum }_{I=1}^{N_{p}}(\boldsymbol{F}_{I}^{D}+\boldsymbol{F}_{I}^{L})\text{e}^{-2\unicode[STIX]{x03C0}\imath [((n_{x}x_{I}/L_{x})+(n_{z}z_{I}/L_{z}))]}\frac{T_{n_{y}}(y_{I})}{\sqrt{1-y_{I}^{2}}}. & & \displaystyle\end{eqnarray}$$ Here, 
$N_{p}$ is the number of particles in the simulation box. This method involves a very high computational cost, since complex exponential functions are evaluated at the particle centre, which is off grid in general.
To avoid the high computational cost of the spectral representation of the reverse force, equation (A 10), the reverse force in the simulations is determined using the projection on nearest neighbours (PNN) method, which is computationally faster. A projection technique is used where the particle force is projected onto the eight surrounding fluid grid points by a volume weighted method, as illustrated in figure 15,
$$\begin{eqnarray}\boldsymbol{f}=-\frac{1}{\unicode[STIX]{x0394}V}\mathop{\sum }_{j=1}^{8}(\boldsymbol{F}_{I}^{D}+\boldsymbol{F}_{I}^{L})(\unicode[STIX]{x0394}V)_{j}^{\prime }\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{j}).\end{eqnarray}$$ Here, the summation is over the eight vertices 
$j$ of the volume 
$\unicode[STIX]{x0394}V$ within which the particle is located, and 
$(\unicode[STIX]{x0394}V)_{j}^{\prime }$ is the volume of the cuboid opposite to the vertex 
$j$.
Figure 15. The PNN method for projecting the force on a particle located within the volume 
$\unicode[STIX]{x0394}V$ onto the vertices of the volume.
Figure 16 shows a comparison of the fluid mean square fluctuating velocities obtained using the PNN projection and the direct Fourier–Chebyshev transform for a suspension with channel Reynolds number 
$Re_{b}=3300$ and particle Stokes number 105.47 in which the particle drag is given by Stokes law. The results of the PNN method are in excellent agreement with those of the spectral interpolation for the range of volume fractions considered here, and so the volume interpolation has been used to obtain all of our results.
Figure 16. The scaled mean square velocities 
$(\overline{u_{x}^{\prime 2}}/\bar{u}^{2})$ (a) and 
$(-\overline{u_{x}^{\prime }u_{y}^{\prime }}/\bar{u}^{2})$ (b) calculated using the direct Fourier–Chebyshev transform (solid lines) and the PNN projection method (dashed lines) for a suspension with particle volume fractions 
$10^{-3}$ (○) and 
$2\times 10^{-3}$ (▵) when the channel Reynolds number is 3300, the particle drag is given by Stokes law (2.5) and the particle Stokes number is 105.47.
Figure 17. The fluid mean velocity (a) and root mean square velocity fluctuation (b) scaled by the friction velocity 
$u_{\ast }$ as a function of the scaled coordinate 
$(yu_{\ast }/\unicode[STIX]{x1D708})$ for an unladen channel flow at 
$Re_{b}=5600$ and 
$Re_{\ast }=180$. In (a), current (——) DNS simulation are compared with (○) Kim et al. (Reference Kim, Moin and Moser1987) for 
$Re_{\ast }=180$. In (b), (——) 
$\sqrt{\overline{u_{x}^{\prime 2}}}/u_{\ast }$, (
$\cdots$)
$\sqrt{\overline{u_{y}^{\prime 2}}}/u_{\ast }$ and (– – –) 
$\sqrt{\overline{u_{z}^{\prime 2}}}/u_{\ast }$, obtained here are compared with (○) 
$\sqrt{\overline{u_{x}^{\prime 2}}}/u_{\ast }$, (▵) 
$\sqrt{\overline{u_{y}^{\prime 2}}}/u_{\ast }$ and (♢) 
$\sqrt{\overline{u_{z}^{\prime 2}}}/u_{\ast }$ from Kim et al. (Reference Kim, Moin and Moser1987).
Figure 18. A comparison between the current simulations (lines) and the experimental data from PIV of Goswami & Kumaran (Reference Goswami and Kumaran2011b) of the fluid mean velocity 
$\bar{u}_{x}$ (a), root mean square velocity fluctuation in the streamwise 
$\sqrt{\overline{u_{x}^{\prime 2}}}$ and wall-normal 
$\sqrt{\overline{u_{y}^{\prime 2}}}$ directions (b) and the correlation 
$\overline{u_{x}^{\prime }u_{y}^{\prime }}$ (c), scaled by suitable powers of the turbulent centre line velocity 
$\bar{u}_{c}$, as a function of the scaled coordinate 
$(y/h)$. The channel Reynolds number 
$Re_{b}=3300$, the particle Stokes number is 2.5 and the particle volume fractions is 
$10^{-4}$. In (b) (——) 
$\sqrt{\overline{u_{x}^{\prime 2}}}/\bar{u}_{c}$ and (– – –) 
$\sqrt{\overline{u_{y}^{\prime 2}}}/\bar{u}_{c}$, from simulations are compared with (○) 
$\sqrt{\overline{u_{x}^{\prime 2}}}/\bar{u}_{c}$, and (●) 
$\sqrt{\overline{u_{y}^{\prime 2}}}/\bar{u}_{c}$ of Goswami & Kumaran (Reference Goswami and Kumaran2011b).
The present results for the mean fluid velocity and root mean square fluctuations, scaled by friction velocity, for an unladen flow for Reynolds number (
$Re_{b}$) of 5600 (
$Re_{\unicode[STIX]{x1D70F}}=180$) are compared with the results of Kim et al. (Reference Kim, Moin and Moser1987) in figure 17. This figure shows that the results of the present simulation are in quantitative agreement with previous simulations for the gas turbulence. The results have also been validated against those of Goswami & Kumaran (Reference Goswami and Kumaran2011a) for an unladen flow at a 
$Re_{b}=3300$, and quantitative agreement was found. The fluid turbulence intensities for 
$Re_{b}=3300$ for a particle-laden flow with Stokes number 2.5 and volume fraction 
$\unicode[STIX]{x1D719}_{av}=10^{-4}$ are compared with experimental results obtained by Goswami & Kumaran (Reference Goswami and Kumaran2011b) using particle image velocimetry (PIV) in figure 18. This figure shows that the fluid mean velocity and the root means square of the fluid velocity fluctuations in the streamwise and wall-normal directions obtained from the present simulations are in good quantitative agreement with experimental data of Goswami & Kumaran (Reference Goswami and Kumaran2011b). Figure 18(c) shows that the error bars in the experimental results for the correlation 
$\overline{u_{x}^{\prime }u_{y}^{\prime }}$ are relatively large compared to the mean values due to the limited resolution in the experiments. Subject to these error bars, there is good quantitative agreement between the simulation and experimental results.
Appendix B. Particle statistics
The cross-stream variations of the particle concentration and velocity fields are shown in figure 19. There is virtually no variation in either the particle concentration or the average velocity with cross-stream distance. It should be noted that the particle–wall collisions are considered to be smooth and elastic; consequently, there is no momentum exchange between the particles and the wall in the flow direction, resulting in a ‘zero-stress’ condition for the particle phase at the wall. There is also virtually no change in the particle mean velocity as the particle loading is increased, and the particle mean velocity is virtually equal to the average fluid velocity. This is because the terminal velocity of the particles is, at maximum, approximately 0.01 times the flow velocity, and so the average slip between the particle and fluid phase is negligible.
Figure 19. The volume fraction 
$\unicode[STIX]{x1D719}$ (a) and the scaled mean particle velocity 
$(\bar{v}_{x}/\bar{u})$ (b) as a function of 
$(y/h)$ for cross-sectionally averaged particle volume fractions 
$\unicode[STIX]{x1D719}_{av}=5\times 10^{-4}$ (○), 
$\unicode[STIX]{x1D719}_{av}=10^{-3}$ (▵), 
$\unicode[STIX]{x1D719}_{av}=1.3\times 10^{-3}$ (▿), 
$\unicode[STIX]{x1D719}_{av}=1.4\times 10^{-3}$ (◃), 
$\unicode[STIX]{x1D719}_{av}=1.6\times 10^{-3}$ (▹) and 
$\unicode[STIX]{x1D719}_{av}=2\times 10^{-3}$ (♢). The Stokes law is used for the drag force on the particle, the channel Reynolds number 
$Re_{b}=3300$ and the particle Stokes number is 105.47. The dashed line in (b) is the fluid mean velocity for the unladen flow.
Figure 19 shows that the particle mean velocity profiles are also independent of the fluid turbulence intensities. In contrast, the particle fluctuating velocities in figure 20 show a strong dependence on the turbulence intensities. The mean square of the streamwise fluctuating velocity 
$\overline{v_{x}^{\prime 2}}$, shown in figure 20(a), decreases with distance from the wall of the channel. When the particle loading is increased from 
$5\times 10^{-4}$ to 
$1.3\times 10^{-3}$, there is a decrease in 
$\overline{v_{x}^{\prime 2}}$. However, when the particle loading is increased from 
$1.3\times 10^{-3}$ to 
$1.4\times 10^{-3}$, there is a sharp increase of approximately 20 % in 
$\overline{v_{x}^{\prime 2}}$. This is followed by a decrease upon further increase in the volume fraction from 
$1.3\times 10^{-3}$ to 
$2\times 10^{-3}$. A comparison of figures 2(a) and 20(a) indicates that 
$\overline{v_{x}^{\prime 2}}$ is significantly smaller than 
$\overline{u_{x}^{\prime 2}}$ for 
$\unicode[STIX]{x1D719}_{av}\leqslant 1.3\times 10^{-3}$, but the particle fluctuating velocity 
$\overline{v_{x}^{\prime 2}}$ becomes larger than 
$\overline{u_{x}^{\prime 2}}$ for 
$\unicode[STIX]{x1D719}_{av}\geqslant 1.4\times 10^{-3}$. This suggests that the fluid velocity fluctuations drive the particle velocity fluctuations for the turbulent state, whereas the particle velocity fluctuations drive the fluid velocity fluctuations after turbulence collapse. A similar inference can be drawn for the fluctuations in the other two directions as discussed below.
Figure 20. The mean square of the particle fluctuating velocities in the streamwise (a), cross-stream (b) and spanwise (c) directions, and correlation 
$\overline{v_{x}^{\prime }v_{y}^{\prime }}$ (d), non-dimensionalised 
$\bar{u}^{2}$, as a function of the scaled coordinate 
$(y/h)$, for average particle volume fractions 
$\unicode[STIX]{x1D719}_{av}=5\times 10^{-4}$ (○), 
$\unicode[STIX]{x1D719}_{av}=10^{-3}$ (▵), 
$\unicode[STIX]{x1D719}_{av}=1.3\times 10^{-3}$ (▿), 
$\unicode[STIX]{x1D719}_{av}=1.4\times 10^{-3}$ (◃), 
$\unicode[STIX]{x1D719}_{av}=1.6\times 10^{-3}$ (▹), 
$\unicode[STIX]{x1D719}_{av}=2\times 10^{-3}$ (♢). The Stokes law is used for the drag force on the particles, the channel Reynolds number is 
$Re_{b}=3300$ and the particle Stokes number is 105.47.
The mean square of the fluctuating velocities in the cross-streamwise and spanwise directions, 
$\overline{v_{y}^{\prime 2}}$ and 
$\overline{v_{z}^{\prime 2}}$, in figure 20(b,c), exhibit remarkably little variation with cross-stream distance. This lack of variation in the cross-stream direction is understandable, because momentum conservation requires that the product of the concentration and 
$\overline{v_{y}^{\prime 2}}$ has to be a constant. Figure 20(b,c) shows that there is a monotonic increase in 
$\overline{v_{y}^{\prime 2}}$ and 
$\overline{v_{z}^{\prime 2}}$ when the average volume fraction is increased from 
$5\times 10^{-4}$ to 
$1.3\times 10^{-3}$, and then a discontinuous increase of approximately 30 % when the average volume fraction is increased from 
$1.3\times 10^{-3}$ to 
$1.4\times 10^{-3}$. A further increase in the volume fraction results in a slight decrease in the mean square of the fluctuating velocity. A comparison of figures 2(b,c) and 20(b,c) shows that 
$\overline{v_{y}^{\prime 2}}$ and 
$\overline{v_{z}^{\prime 2}}$ are comparable to the fluid mean square velocities 
$\overline{u_{y}^{\prime 2}}$ and 
$\overline{u_{z}^{\prime 2}}$ for 
$\unicode[STIX]{x1D719}_{av}\leqslant 1.3\times 10^{-3}$ when the fluid flow is turbulent. For 
$\unicode[STIX]{x1D719}_{av}\geqslant 1.4\times 10^{-3}$, the moments of the particle velocity fluctuations are significantly higher than the moments of the fluid velocity fluctuations.
Figure 20(d) shows that there is a discontinuous increase in the magnitude of the velocity correlation 
$\overline{v_{x}^{\prime }v_{y}^{\prime }}$ when the average volume fraction is increased from 
$1.3\times 10^{-4}$ to 
$1.4\times 10^{-4}$. The magnitude of the particle velocity moment 
$\overline{v_{x}^{\prime }v_{y}^{\prime }}$ is significantly smaller than the fluid velocity moment 
$\overline{u_{x}^{\prime }u_{y}^{\prime }}$ (figure 2) for the turbulent flow for 
$\unicode[STIX]{x1D719}_{av}\leqslant 1.3\times 10^{-3}$, but it becomes significantly larger after turbulence collapse for 
$\unicode[STIX]{x1D719}_{av}\geqslant 1.4\times 10^{-3}$. Since the mass loading in our simulations is 
$O(1)$, this implies that the stress transmitted due to the particle velocity fluctuations is significantly larger than that due to the fluid velocity fluctuations for 
$\unicode[STIX]{x1D719}_{av}\geqslant 1.4\times 10^{-3}$. Thus, the discontinuous decrease in the turbulence intensity at the critical volume loading also causes a discontinuous increase in all the second moments of the particle fluctuating velocities. Whilst the second moments of the particle velocity fluctuations are lower than those for the fluid velocity fluctuations for 
$\unicode[STIX]{x1D719}_{av}\leqslant 1.3\times 10^{-3}$, the reverse is true for 
$\unicode[STIX]{x1D719}_{av}\geqslant 1.4\times 10^{-3}$.
Figure 21. The scaled time between particle–particle collisions 
$(t_{pp}\bar{u}/h)$ (solid lines), time between particle–wall collisions 
$(t_{pw}\bar{u}/h)$ (dashed lines) and the viscous relaxation time of the particles 
$(t_{v}\bar{u}/h)$ (horizontal dotted lines) as a function of the volume fraction for channel Reynolds number 
$Re_{b}=3300$ and Stokes numbers of 21.09 (○), 105.47 (▵) and 316.40 (▿), when the Stokes law (2.5) is used for the drag force on the particles (a); and for Stokes numbers of 6.49 (○), 32.45 (▵) and 97.34 (▿) when the Schiller–Naumann correlation (2.6) is used for the drag force on the particles (b). The vertical dotted lines show the range of 
$\unicode[STIX]{x1D719}_{av}$ over which there is turbulence collapse for the Stokes number indicated by the symbols on the abscissa.
Figure 22. The measures of the mean square of the fluctuating velocities (3.2) scaled by the square of the average flow velocity, 
$(\langle \overline{v_{x}^{\prime 2}}\rangle _{s}/\bar{u}^{2})$ (a), 
$(\langle \overline{v_{y}^{\prime 2}}\rangle _{s}/\bar{u}^{2})$ (b), 
$(\langle \overline{v_{z}^{\prime 2}}\rangle _{s}/\bar{u}^{2})$ (c) and 
$(-\langle \overline{v_{x}^{\prime }v_{y}^{\prime }}\rangle _{s}/\bar{u}^{2})$ (d) as a function of the average volume fraction for channel Reynolds number 
$Re_{b}=3300$ and particle Stokes numbers of 21.09 (○), 105.47 (▵) and 316.40 (▿), when the Stokes law (2.5) is used for the drag force on the particles (solid lines); and for Stokes numbers of 6.49 (●), 32.45 (▴) and 97.34 (▾) when the Schiller–Naumann correlation (2.6) is used for the drag force on the particles (dashed lines). The vertical dotted lines show the range of 
$\unicode[STIX]{x1D719}_{av}$ over which there is turbulence collapse for the particle drag law and Stokes number indicated by the symbols on the abscissa.
The time between particle–particle collisions averaged across the channel width, 
$t_{pp}$, and the time between particle–wall collisions, 
$t_{pw}$, all scaled by the fluid integral time 
$(h/\bar{u})$, are shown as a function of the volume fraction in figure 21. The viscous relaxation time, equations (2.13) and (2.14), which is independent of the particle volume fraction, is shown by the horizontal dotted lines in figure 21. The results in figure 21(a) are obtained when the Stokes drag law is used for the particle drag. Here, the viscous relaxation time is larger than the particle–particle and particle–wall collision times for Stokes numbers of 105.47 and 316.40, implying that the particle does not relax to the local fluid velocity between successive collisions. For the Stokes number 21.09, the viscous relaxation time is smaller than the collision time for low volume fraction, and larger than the collision time for high volume fraction. When the Schiller–Naumann correlation (2.6) is used for the drag force, the viscous relaxation time is smaller than the collision time for the Stokes number 6.49, comparable to the collision time for the Stokes number 32.45 and larger for the Stokes number 97.34. Thus, figure 21 shows that our simulations span the range 
$t_{pw}<t_{pp}$ where the particle travels across the channel before colliding with another particle, and 
$t_{pp}<t_{pw}$, where the particle collides with other particles before crossing the width of the channel. The collapse in turbulence intensities at a critical volume loading is observed for this entire range of ratios of time scales.
Figure 23. The probability distributions for the particle velocity fluctuations in the streamwise (a), cross-stream (b) and spanwise (c) directions in the zone near the wall (○) and at the centre (▵) for particle volume fractions of 
$1.3\times 10^{-3}$ (solid line) and 
$1.4\times 10^{-3}$ (dashed line) when the channel Reynolds number is 3300, the Stokes drag law is used for the drag force and the particle Stokes number is 105.47. The dotted line is the Gaussian distribution.
The spatially averaged mean square of the particle velocity fluctuations, scaled by the square of the average flow velocity, is shown as a function of the average particle volume fraction in figure 22. The volume fractions bracketing the turbulence collapse transition are shown by the vertical lines, with the symbols on the abscissa indicating the drag law and Stokes number. There is a decrease in the streamwise fluctuating velocity and an increase in the cross-stream and spanwise fluctuating velocities as the volume fraction is increased, due to the enhanced transfer of fluctuating energy from the streamwise to the cross-stream directions by collisions. At the transition volume fraction, there is a significant increase in the mean square velocity in all three directions. This increase is relatively small, approximately 10 %, in the streamwise fluctuating velocity, but it is quite large, up to 50 %, in the cross-steam and spanwise directions. There is also a moderate increase of approximately 20 % in magnitude of the correlation 
$\langle \overline{v_{x}^{\prime }v_{y}^{\prime }}\rangle _{s}$. However, the increase in the particle velocity fluctuations is significantly smaller than the decrease by 1–2 orders of magnitude in the fluid velocity fluctuations in figure 6.
The probability distributions for the particle fluctuating velocities are shown in figure 23 for the turbulent flow at volume fraction 
$1.3\times 10^{-3}$, and after turbulence collapse at volume fraction 
$1.4\times 10^{-3}$. Figure 23 shows that there is virtually no change in the distribution functions for the particle fluctuating velocities, in contrast to the qualitative change in the distribution function for the fluid fluctuating velocities observed in figure 8. Thus, the particle statistics are unchanged even when there is a significant change in the fluid statistics at turbulence collapse.
