2.1. Simulation cases

To investigate the topics mentioned in the introduction, numerical simulations are performed. The configuration and parameters of the simulation cases are based on two experimental cases in KFE1994: the unladen case and the particle-laden case for which the strongest reduction of turbulence was measured. Although the basic presentation will be in dimensional variables, to retain the advantage of dimensional checking, the relevant non-dimensional numbers will be specified. Table 1 shows an overview of four simulations, A0–A3. The unladen flow (A0) is a turbulent channel flow at Reynolds number $Re_{{\it\tau}}=642$ , based on channel half-width $H=0.02~\text{m}$ , wall friction velocity $u_{{\it\tau},0}=0.4896~\text{m}~\text{s}^{-1}$ and kinematic viscosity ${\it\nu}=1.525\times 10^{-5}~\text{m}^{2}~\text{s}^{-1}$ . Thus, the viscous length scale of the reference flow is equal to ${\it\delta}_{{\it\nu}}=H/Re_{{\it\tau}}=31.15~{\rm\mu}\text{m}$ . The centreline velocity equals $u_{cl,0}=10.42~\text{m}~\text{s}^{-1}$ and the bulk velocity is $9.212~\text{m}~\text{s}^{-1}$ . In addition, the gas is air with a constant density of ${\it\rho}=1.2~\text{kg}~\text{m}^{-3}$ , and the flow in the vertically positioned channel is downward (gravity acceleration $g=9.8~\text{m}~\text{s}^{-2}$ ). The diameter of the (free) particles in the laden cases is $d_{p}=70~{\rm\mu}\text{m}$ , and the density of the particles is the density of copper, ${\it\rho}_{p}=8800~\text{kg}~\text{m}^{-3}$ . The Stokes response time, ${\it\tau}_{p}={\it\rho}_{p}d_{p}^{2}/(18{\it\rho}{\it\nu})$ , is equal to 0.131 s, which corresponds to a Stokes number ${\it\tau}_{p}u_{{\it\tau},0}/{\it\delta}_{{\it\nu}}\approx 2060$ . The particle diameter is small compared with the length scale of most turbulent eddies, $d_{p}\approx 0.45{\it\eta}$ at the centre and $d_{p}\approx 1.55{\it\eta}$ at the wall, where ${\it\eta}=({\it\nu}^{3}/{\it\epsilon})^{1/4}$ is the Kolmogorov length scale, based on the unladen turbulence dissipation rate ( ${\it\epsilon}$ ), which is a function of the distance to the nearest wall. The domain-averaged solid volume fraction of the particles is the same in the three laden cases, namely $1.092\times 10^{-4}$ . By definition, the mass loading ratio parameter ${\it\phi}$ is equal to the overall solid volume fraction multiplied by ${\it\rho}_{p}/{\it\rho}$ , thus ${\it\phi}=0.801$ (80 %). The walls are smooth in case A1 and rough in cases A2 and A3. The wall in case A3 is rougher than in case A2. The wall roughness will be defined in a separate section.

2.2. The Eulerian–Lagrangian approach

The unladen case (A0) is simulated with DNS, while the laden cases (A1–A3) are simulated with an Eulerian–Lagrangian method, PP-DNS with two-way coupling and inter-particle collisions. For an introduction to the Eulerian–Lagrangian method, the reader is referred to Crowe, Sommerfeld & Tsuji (Reference Crowe, Sommerfeld and Tsuji1998), Deen et al. (Reference Deen, van Sint Annaland, van der Hoef and Kuipers2007) and van der Hoef et al. (Reference van der Hoef, van Sint Annaland, Deen and Kuipers2008). The method can be one-way coupled (the particle–fluid interaction appears in the particle equations only), two-way coupled (the particle–fluid interaction appears in the fluid and particle equations) or two-way coupled with inter-particle collisions included. In PP-LES studies of dilute vertical downward particle-laden channel flow at $Re_{{\it\tau}}\approx 640$ in the literature, all three types have been used: one-way coupling (Wang & Squires Reference Wang and Squires1996), two-way coupling (Segura Reference Segura2004) and two-way coupling with inter-particle collisions (Yamamoto et al. Reference Yamamoto, Potthoff, Tanaka, Kajishima and Tsuji2001).

PP-DNS studies of vertical downward particle-laden plane channel flow in the literature were hitherto limited to lower Reynolds number: $Re_{{\it\tau}}=180$ with one-way coupling (Rouson & Eaton Reference Rouson and Eaton2001), $Re_{{\it\tau}}=210$ with two-way coupling (Kubik & Kleiser Reference Kubik and Kleiser2006), and $Re_{{\it\tau}}=125$ (Li et al. Reference Li, McLaughlin, Kontomaris and Portela2001) and $Re_{{\it\tau}}=150$ (Mito & Hanratty Reference Mito and Hanratty2006) with two-way coupling and inter-particle collisions. If we include PP-DNS of horizontal or non-gravitational particle-laden plane channel flows in our discussion, one-way coupled simulations of such flows were (for example) performed at $Re_{{\it\tau}}=150$ (Pedinotti, Mariotti & Banerjee Reference Pedinotti, Mariotti and Banerjee1992; Marchioli & Soldati Reference Marchioli and Soldati2002), $Re_{{\it\tau}}=300$ (Lavezzo et al. Reference Lavezzo, Soldati, Gerashchenko, Warhaft and Collins2010) and $Re_{{\it\tau}}=950$ (Geurts & Kuerten Reference Geurts and Kuerten2012; Kuerten & Brouwers Reference Kuerten and Brouwers2013), two-way coupled simulations up to $Re_{{\it\tau}}=395$ (Kuerten, van der Geld & Geurts Reference Kuerten, van der Geld and Geurts2011; Zhao, Andersson & Gillissen Reference Zhao, Andersson and Gillissen2013) and two-way coupled simulations with inter-particle collisions at $Re_{{\it\tau}}=125$ (Nasr, Ahmadi & McLaughlin Reference Nasr, Ahmadi and McLaughlin2009).

In view of the above review of the literature, PP-DNS with two-way coupling and inter-particle collisions is a state-of-the-art method to simulate particle-laden channel flow at $Re_{{\it\tau}}=642$ . For most particle-laden flows, including the present flow, it is not yet feasible to perform a DNS that resolves the boundary layers around the particles. In the class of feasible methods, two-way coupled PP-DNS with inter-particle collisions is probably the most advanced one. In addition, the present PP-DNS cases have several new features. First, although $Re_{{\it\tau}}=642$ is not high, the Reynolds number is approximately four times higher than in previous studies of turbulent channel flows using the same method (PP-DNS with two-way coupling and inter-particle collisions). Although a detailed study of the effect of Reynolds number falls outside the scope of this paper, the effect will be briefly addressed at the end of § 3.3. Second, it is the first time that wall roughness effects have been included in PP-DNS plane channel flow. It is not expected that the results would be very different if PP-LES at the same Reynolds number were used. Point-particle DNS has been chosen because it is feasible for this Reynolds number and it does reduce uncertainties caused by the simulation method to some extent, simply because the range of scales taken into account is larger than in PP-LES. Without correction, the turbulence intensity and turbulence kinetic energy in LES are generally lower than in DNS (how much lower depends on the resolution of the LES). Third, it is the first time that all interactions between the particles and the rough walls have been simulated with a fully deterministic approach. A comparison between the deterministic model and existing stochastic models falls outside the scope of this paper.

2.3. Governing equations

The governing equations of the gas phase are $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0$ and

where $\boldsymbol{u}$ is the gas velocity, $P$ is the part of the pressure that is periodic in the streamwise direction and $\boldsymbol{f}$ the feedback force, the force exerted by the particles on the gas per unit mass of gas. The gravity term points in the direction of the mean flow; the unit vector in the streamwise direction is denoted with $\boldsymbol{e}_{\mathbf{1}}$ . In addition, ${\it\rho}a(t)$ is the (time-dependent) domain-averaged streamwise component of the pressure gradient, discussed in a later section. Time is indicated by $t$ , and the streamwise, normal and spanwise directions are indicated by $x_{1}$ , $x_{2}$ and $x_{3}$ respectively. The computational domain is given by $[0,L_{1}]\times [0,L_{2}]\times [0,L_{3}]$ , with $L_{1}=6H$ and $L_{2}=L_{3}=2H$ . Periodic boundary conditions are imposed in the streamwise and spanwise directions, while in the normal direction no-slip boundary conditions are applied to the gas at the walls.

For the discretization of these equations a staggered central differencing method is used, second order in the normal and fourth order in the homogeneous directions (Vreman & Kuerten Reference Vreman and Kuerten2014a ). The implementation of the no-slip boundary conditions of the tangential velocities at the wall involves third-order extrapolations of these velocities to ghost cells, based on the boundary condition and the velocities in the first three grid cells off the wall (see Vreman (Reference Vreman2014) for a discussion of different implementations of the no-slip condition on a staggered grid). The pressure Poisson equation is solved by a direct method (fast Fourier transforms in the homogeneous directions). The temporal discretization is an explicit compact storage Runge–Kutta method with Runge–Kutta coefficients $1/3$ , $1/2$ and 1, in fact the three-stage variant of the four-stage method of Jameson & Baker (Reference Jameson and Baker1983), see also Vreman (Reference Vreman2014). The time step ${\rm\Delta}t$ is equal to $10^{-5}~\text{s}$ . The spatial grid contains $512\times 320\times 256$ cells and is stretched in the normal direction by the tangent hyperbolical function (Vreman Reference Vreman2014). The grid points $x_{2,j}$ of the cell centres are located at

The uniform streamwise and spanwise grid spacings are equal to $7.5{\it\delta}_{{\it\nu}}$ and $5.0{\it\delta}_{{\it\nu}}$ respectively. The grid spacing in the normal direction varies from $0.99{\it\delta}_{{\it\nu}}$ at the wall to $7{\it\delta}_{{\it\nu}}$ at the centre of the channel. These are standard grid sizes for DNS of single-phase channel flow (Moser, Kim & Mansour Reference Moser, Kim and Mansour1999) and are sufficiently small to compute mean velocity and Reynolds stress profiles accurately (Vreman & Kuerten Reference Vreman and Kuerten2014a ,Reference Vreman and Kuerten b ).

The particles are tracked in a Lagrangian framework. The following equations are solved for each particle $p$ :

In scalar particle variables, $p$ is written as a subscript, while in vector particle variables $p$ is written as superscript.

The particle position, velocity and angular velocity are respectively denoted by $\boldsymbol{x}^{p}$ , $\boldsymbol{v}^{p}$ and ${\bf\omega}^{p}$ , while $m_{p}={\rm\pi}{\it\rho}_{p}d_{p}^{3}/6$ is the mass of the particle and $I_{p}=(m_{p}d_{p}^{2})/10$ is the moment of inertia of the particle. The force and torque caused by the inter-particle collisions and particle–wall collisions on particle $p$ at time $t$ are represented by $\boldsymbol{F}_{coll}^{p}$ and $\boldsymbol{T}_{coll}^{p}$ , while the force and torque exerted by the gas on the particle are represented by $\boldsymbol{F}^{p}$ and $\boldsymbol{T}^{p}$ . These forces and torques will be defined in later subsections.

The time integration of the particle position and velocity is performed with forward Euler with (at least) three time steps per time step ${\rm\Delta}t$ of the Runge–Kutta method used for the gas phase. The forward Euler method has been selected because the time step is not the same for each particle, due to collisions. Since there are many collisions per time step, it is not efficient to update all particles after each collision. The time levels after stages 1, 2 and 3 of the Runge–Kutta method are $t+{\rm\Delta}t/3$ , $t+{\rm\Delta}t/2$ and $t+{\rm\Delta}t$ respectively. The basic particle time step, ${\rm\Delta}t_{p}$ , is set to ${\rm\Delta}t/3$ in stage 1, ${\rm\Delta}t/6$ in stage 2 and ${\rm\Delta}t/2$ in stage 3. For particles that experience one or more collisions within ${\rm\Delta}t_{p}$ , the time integration over ${\rm\Delta}t_{p}$ is divided into two or more substeps.

The feedback force at location $\boldsymbol{x}$ and time $t$ is defined by

where the ${\it\delta}_{x}$ is the delta function with unit $\text{m}^{-3}$ . The discretization of this equation is

where $|V_{\boldsymbol{x}}|$ is the volume of the Eulerian (gas phase) grid cell $V_{\boldsymbol{x}}$ around $\boldsymbol{x}$ , and the coefficients $w_{p}(\boldsymbol{x},t)$ represent the force distribution coefficients. The simplest way is to define $w_{p}(\boldsymbol{x},t)=1$ if the particle location $\boldsymbol{x}^{p}(t)$ is inside $V_{\boldsymbol{x}}$ and $w_{p}(\boldsymbol{x},t)=0$ otherwise. However, we use a slightly different approach, which is more accurate and also more convenient on a staggered grid. In this approach the three gas velocity components are transferred from staggered to cell-centre positions (the cell centre is the position where also the pressure and the velocity divergence of the grid are defined). Then the gas velocity in $\boldsymbol{x}^{p}$ is defined by linear interpolation of the cell-centre gas velocity to particle locations (particle centres). Then the particle forces $\boldsymbol{F}^{p}$ are computed. The weights of the linear interpolation just mentioned are reused to distribute the particle force, which is defined at the particle centre, to the eight surrounding cell-centre positions. Finally, the force $\boldsymbol{f}$ is transferred to the staggered positions at cell faces by taking the average over the appropriate two cells for each component. No special treatment is used in the wall proximity. There, the cell size spacing in the normal direction is approximately two times smaller than the particle size, but the effective volume of fluid over which a particle force is distributed is still considerably larger than the particle volume. At the wall, the Eulerian cell volume $|V_{x}|$ is equal to 6.3 particle volumes. Due to the averaging of $\boldsymbol{f}$ to the staggered positions, the particle force is distributed over at least $2|V_{x}|$ . At least, because due to the use of weight coefficients $w_{p}$ in (2.7), the effective fluid volume over which each particle force is distributed is generally larger than $2|V_{x}|$ .

The unladen case A0 is initialized by a parabolic mean profile with a suitable divergence-free perturbation. The particle-laden cases A1–A3 are initialized by a snapshot of the turbulent velocity field of case A0. Statistical averaging is started only after transient effects have disappeared. In each case the statistical averaging time is at least 0.5 s, which corresponds to 38 flow-through times, with the flow-through time defined by the streamwise channel length divided by the gas bulk velocity. Each particle-laden case contains 116 756 particles.

2.4. Particle–gas forces

The force $\boldsymbol{F}^{p}$ , exerted by the gas on a single free particle, is parametrized by the well-known Schiller–Naumann standard drag correlation:

where $Re_{p}$ is the Reynolds number of the particle, defined by $d_{p}|\boldsymbol{u}-\boldsymbol{v}^{p}|/{\it\nu}$ . Based upon recent literature on the relevance of the lift force (Zeng et al. Reference Zeng, Balachandar, Fischer and Najjar2008), the inertial, added mass and history force (Armenio & Fiorotto Reference Armenio and Fiorotto2001; Bagchi & Balachandar Reference Bagchi and Balachandar2003; Burton & Eaton Reference Burton and Eaton2005), and hydrodynamic particle–particle interactions (Vreman Reference Vreman2007; Nasr et al. Reference Nasr, Ahmadi and McLaughlin2009), explicit models for these forces are not included. In one-way coupled simulations of similar case, Rouson & Eaton (Reference Rouson and Eaton2001) tested the lift force expression derived by McLaughlin (Reference McLaughlin1993) and found that the term is negligible compared with the particle drag in the same direction. Fully resolved simulations of a single small particle in isotropic turbulence (Bagchi & Balachandar Reference Bagchi and Balachandar2003; Burton & Eaton Reference Burton and Eaton2005) and in wall-bounded turbulence (Zeng et al. Reference Zeng, Balachandar, Fischer and Najjar2008) indicate that the standard drag law provides a good description for the mean force on the particle. However, these references also show that the error in the fluctuating force may be large and that incorporation of other terms of the Maxey Riley equation is no remedy for this limitation. A simulation method that fully resolves the boundary layers around the particles would of course be a more accurate method to simulate the turbulence modification. However, such a simulation is not possible yet for the present case, in which the flow domain contains more than $10^{5}$ small moving particles.

The particle angular velocity ${\bf\omega}^{p}$ is included in the PP-DNS, since it is required to describe non-ideal collisions, in particular collisions with non-zero tangential friction. The particle angular velocity is both generated and dissipated by non-ideal collisions via $\boldsymbol{T}^{p,c}$ . In addition, the particle angular velocity is influenced by the particle drag torque $\boldsymbol{T}^{p}$ . The drag torque is included, because it was found to dampen the particle angular velocity fluctuation considerably. The drag torque is modelled by (Takagi Reference Takagi1977; Dennis, Singh & Ingham Reference Dennis, Singh and Ingham1980; Yamamoto et al. Reference Yamamoto, Potthoff, Tanaka, Kajishima and Tsuji2001)

In this equation ${\bf\zeta}^{p}={\bf\omega}^{p}-(\boldsymbol{{\rm\nabla}}\times \boldsymbol{u})/2$ and $Re_{p,a}=d_{p}^{2}|{\bf\zeta}^{p}|/4{\it\nu}$ is the particle Reynolds number based on the angular velocity difference. The torque coefficients $(C_{1},C_{2},C_{3})$ are given by (0, 50.27, 0) if $Re_{p,a}<1$ , (0, 50.27, 0.0418) if $1<Re_{p,a}<10$ , (5.32, 37.2, 5.32) if $10<Re_{p,a}<20$ , (6.44, 32.2, 6.44) if $20<Re_{p,a}<50$ and (6.45, 32.1, 6.45) if $Re_{p,a}>50$ .

2.5. Particle collisions

Although the volume fraction of the particle-laden simulations is only of the order of $10^{-4}$ , all collisions between particles are taken into account. It is often mentioned that inter-particle collisions are negligible if the volume fraction is smaller than $10^{-3}$ . This may be the case in isotropic turbulence; for wall-bounded turbulence it has repeatedly been shown that collisions do significantly influence the results if the volume fraction is of the order of $10^{-4}$ or even lower (Li et al. Reference Li, McLaughlin, Kontomaris and Portela2001; Yamamoto et al. Reference Yamamoto, Potthoff, Tanaka, Kajishima and Tsuji2001; Vreman Reference Vreman2007; Nasr et al. Reference Nasr, Ahmadi and McLaughlin2009). These references show that collisions reduce the tendency of preferential concentration of particles. Without collisions all particles slowly accumulate either at the walls or at the centre (Rouson & Eaton Reference Rouson and Eaton2001), while with collisions a statistically stationary steady state can be achieved.

The particle collisions are modelled with a hard-sphere collision model. The normal restitution coefficient is given by $e_{r}=0.95$ and the tangential friction coefficient by ${\it\mu}=0.3$ , which are the same values as in Yamamoto et al. (Reference Yamamoto, Potthoff, Tanaka, Kajishima and Tsuji2001). The tangential restitution coefficient is given by ${\it\beta}_{0}=0.33$ (Hoomans Reference Hoomans1999). Two particles $p$ and $q$ collide when the distance between their point-particle locations is equal to $(d_{p}+d_{q})/2$ . Conservation of momentum and angular momentum during collision between two particles $p$ and $q$ implies the following modification of the particle velocity and angular velocity (Hoomans et al. Reference Hoomans, Kuipers, Briels and van Swaaij1996; Hoomans Reference Hoomans1999):

where the superscript zero denotes the values before collision. The normal unit vector $\boldsymbol{n}$ is defined by $(\boldsymbol{x}^{p}-\boldsymbol{x}^{q})/|\boldsymbol{x}^{p}-\boldsymbol{x}^{q}|$ . The relative velocity vector before collision is defined by

with normal component $v_{n}=\boldsymbol{v}^{pq,0}\boldsymbol{\cdot }\boldsymbol{n}$ . The tangential unit vector $\boldsymbol{t}$ is defined by $(\boldsymbol{v}^{pq,0}-v_{n}\boldsymbol{n})/|\boldsymbol{v}^{pq,0}-v_{n}\boldsymbol{n}|$ . The impulse vector $\boldsymbol{J}$ is then defined as $\boldsymbol{J}=J_{n}\boldsymbol{n}+J_{t}\boldsymbol{t}$ , where the normal component $J_{n}$ is given by $-(1+e_{r})v_{n}/B_{2}$ and $B_{2}=1/m_{p}+1/m_{q}$ . The tangential component $J_{t}$ depends on the maximum Coulomb friction ${\it\mu}J_{n}$ and $J_{{\it\beta}}=(1+{\it\beta}_{0})(\boldsymbol{v}^{pq,0}\boldsymbol{\cdot }\boldsymbol{t})/B_{1}$ with $B_{1}=7B_{2}/2$ . It can be proven that at contact ${\it\mu}J_{n}\geqslant 0$ and $J_{{\it\beta}}\geqslant 0$ . If $J_{{\it\beta}}\leqslant {\it\mu}J_{n}$ , the collision is of the sticking type and $J_{t}=-J_{{\it\beta}}$ . However, if $J_{{\it\beta}}>{\it\mu}J_{n}$ , the collision is of the sliding type and $J_{t}=-{\it\mu}J_{n}$ (Hoomans et al. Reference Hoomans, Kuipers, Briels and van Swaaij1996; Hoomans Reference Hoomans1999). Since (2.13) with a plus instead of a minus sign was also encountered in the literature, it is stressed that both minus signs in (2.12) and (2.13) are correct. This directly follows from the conservation of angular momentum about the centre of mass of the two particles during a collision.

If particle $p$ is involved in collision $k$ , we define $\boldsymbol{J}^{p,k}=m_{p}(\boldsymbol{v}^{p}-\boldsymbol{v}^{p,0})$ , where $\boldsymbol{v}^{p}$ is the post- and $\boldsymbol{v}^{p,0}$ the pre-collision velocity of particle $p$ at collision $k$ . Thus, for the particles $p$ and $q$ involved in collision $k$ , we have $\boldsymbol{J}^{q,k}=-\boldsymbol{J}^{p,k}=-\boldsymbol{J}$ . For any particle $p$ not involved in collision $k$ , $\boldsymbol{J}^{p,k}=0$ . If particle $p$ collides with a plane wall in the $x_{2}$ direction, then the relative velocity vector is aligned with the $x_{2}$ direction. The equations for a particle–wall collision can be derived from those of a particle–particle collision in the limit $d_{q}\rightarrow \infty$ and $m_{q}\rightarrow \infty$ .

The collision force and torque vectors, $\boldsymbol{F}_{coll}^{p}$ and $\boldsymbol{T}_{coll}^{p}$ , can be expressed as a sum of temporal delta functions, one delta function for each collision,

where $k$ is the index that runs over all collisions (note that $\boldsymbol{J}^{p,k}=0$ for any particle $p$ not involved in collision $k$ ). The time of collision $k$ is denoted by $t_{coll,k}$ . The unit of the delta function ${\it\delta}_{t}$ is $\text{s}^{-1}$ . If the collisions were ideal, $e_{r}$ would be 1 and ${\it\beta}_{0}$ would be $-1$ , then $J_{t}$ , $\boldsymbol{n}\times \boldsymbol{J}$ and $\boldsymbol{T}_{coll}^{p}$ would be zero. If in addition the rotational drag force were neglected, the angular particle velocities would remain constant in time and there would be no need to include them in the simulation.

An efficient search for collision partners is facilitated by a uniform secondary grid ( $120\times 40\times 40$ cells). Each cell of this secondary grid has a list attached that contains the particle indices for the particles located in that cell. Collisions that involve a free particle in a certain cell of the particle grid are detected by checking all particles in the neighbouring cells. The lists are set up before each particle time step and updated after each collision. To ensure that no collision is missed, the absolute velocity of each particle is verified to remain smaller than the size of the particle grid cells divided by twice the basic particle time step. In addition, for each collision it is verified that the particle involved touches the other particle (or the wall) without any overlap.

2.6. Forcing of the flow

Before we describe the forcing technique, the two averaging operators used in this paper are defined. The statistical mean of a quantity $Q$ is defined by

where $t_{1}$ and $t_{2}$ are two times, such that (theoretically) $t_{2}-t_{1}\rightarrow \infty$ . To increase the statistical sample size in practice, both halves of the channel are included in the average, using a parity coefficient $c_{Q}$ , equal to 1 for even and $-1$ for odd statistical quantities. Thus, $c_{Q}=1$ in the computation of, for example, $\overline{u}_{1}$ , $\overline{u}_{3}$ , $\overline{u_{1}^{2}}$ , $\overline{u_{2}^{2}}$ , $\overline{u_{3}^{2}}$ , $\overline{\partial u_{2}/\partial x_{2}}$ , $\overline{f_{1}}$ , $\overline{u_{2}f_{2}}$ , while $c_{Q}=-1$ in the computation of, for example, $\overline{u}_{2}$ , $\overline{u_{1}u_{2}}$ , $\overline{\partial u_{1}/\partial x_{2}}$ , $\overline{f_{2}}$ . The ensemble average of particle properties, also denoted with an overbar, is similarly defined, except that the integrals over the $x_{1}{-}x_{3}$ planes are replaced by averages over all particles in (or near) these planes. The fluctuation of a quantity is defined by $Q^{\prime }=Q-\overline{Q}$ . In addition, we introduce the domain average of a quantity $Q$ :

The domain average applied to a mean quantity, $\langle \overline{Q}\rangle$ , is the same as the mean quantity averaged over the $x_{2}$ direction and does not depend on $t$ , $x_{1}$ , $x_{2}$ , $x_{3}$ .

The overall force balance of the flow is given by the sum of the domain-averaged mean gas and domain-averaged mean particle momentum equation in the streamwise direction:

where $u_{{\it\tau}}=({\it\nu}\,\text{d}\overline{u}/\text{d}x_{2})^{1/2}$ is the friction velocity. The term $F_{wall}$ represents the streamwise force exerted by the walls on the particles per unit of volume,

where the index $p$ runs over all free particles. If this index ran over the free particles that collide with the wall (or collide with fixed wall particles), then the force would be equal and opposite, since the sum over the forces of all collisions between free particles is zero ( $\boldsymbol{J}^{q,k}=-\boldsymbol{J}^{p,k}$ ).

The overall balance, (2.19), shows that it is not easy to decide what should remain constant in a study where laden cases are compared with an unladen case. Suppose that the streamwise pressure gradient, ${\it\rho}\overline{a}$ , is held constant. Then an increase of ${\it\phi}$ necessarily implies an increase of the wall friction forces, ${\it\rho}u_{{\it\tau}}^{2}/H-F_{wall}$ , which means that the fluid friction velocity, $u_{{\it\tau}}$ , or the particle–wall friction force, $-F_{wall}$ , increases. If the latter were negligible, $u_{{\it\tau}}$ at ${\it\phi}=0.8$ would become $1.8^{1/2}=1.34$ times larger than in the unladen case, i.e. the friction Reynolds number and probably also the bulk Reynolds number would become much larger, such that the grid would have to be finer than in the unladen case.

The literature on vertical particle-laden channel or pipe flows shows that, instead of the streamwise pressure gradient, another quantity is usually fixed in a comparison between laden and unladen cases. In the experiments of KFE1994, the mean centreline velocity of the gas was fixed, since the authors state on pp. 116–117: ‘for each set of particle conditions, the tunnel speed was adjusted so that the laden-flow velocity matched the unladen-flow velocity at the centreline’. In the simulations of Yamamoto et al. (Reference Yamamoto, Potthoff, Tanaka, Kajishima and Tsuji2001) the quantity $(1+{\it\phi})g+\overline{a}$ was fixed. In the simulations of Segura (Reference Segura2004) the bulk flow was apparently fixed, since the authors mention that the simulations were run at the same bulk Reynolds number. Mito & Hanratty (Reference Mito and Hanratty2006) explicitly mention that the bulk velocity was fixed. In the pipe flow experiments of Caraman et al. (Reference Caraman, Borée and Simonin2003) and Borée & Caraman (Reference Borée and Caraman2005) and corresponding simulations of Vreman (Reference Vreman2007), the bulk flow was fixed as well.

In view of the above discussion of the literature, the present simulations are performed for a prescribed bulk velocity, which means that $a(t)$ is adapted at each time step to keep the domain-average velocity constant. The same bulk velocity is prescribed in each simulation case. Since the numerical scheme is fully explicit and conserves momentum, there is no need for a control algorithm; $L_{1}L_{2}L_{3}(g+a(t))$ is set equal to the instantaneous viscous shear stress integrated over the walls minus the instantaneous streamwise feedback force integrated over the spatial domain. Forcing to maintain a constant volumetric flow rate is also a standard approach in DNS of unladen turbulent channel flows. Since the effects of the particles on the mean gas velocity are relatively small (see the next section), a comparison for constant bulk velocity can be expected to lead to the same conclusions as a comparison for constant centreline velocity.

2.7. Wall roughness

A rough wall is defined as a plane coated with tiny wall particles, half-spheres of diameter $d_{p,w}$ , with the circular side of each half-sphere fixed on the plane. The wall particles are arranged in square pitch of size $d_{p,w}$ ; each wall particle touches four neighbouring wall particles, two of these neighbours are aligned in the $x_{1}$ direction, the other two are aligned in the $x_{3}$ direction. For this pattern, the standard deviation of the wall roughness $k_{rms}$ is approximated by $0.172d_{p,w}$ , while the peak-to-trough wall roughness $k_{l}$ is equal to $0.5d_{p,w}$ . The diameter $d_{p,w}$ of the wall particles is smaller than the wall unit ${\it\delta}_{{\it\nu}}=H/Re_{{\it\tau}}=31~{\rm\mu}\text{m}$ : $d_{p,w}=10~{\rm\mu}\text{m}$ in case A2 and $d_{p,w}=20~{\rm\mu}\text{m}$ in case A3. Thus, $k_{rms}^{+}=0.055$ and $k_{l}^{+}=0.16$ in case A2, while $k_{rms}^{+}=0.11$ and $k_{l}^{+}=0.32$ in case A3. As discussed in the introduction, these values are sufficiently small to justify the application of smooth boundary conditions for the gas phase.

The equations for a collision between a free particle $p$ and a wall particle $q$ can be derived from the equations for a collision between two free particles by taking the limit $m_{q}\rightarrow \infty$ , such that $\boldsymbol{v}^{q}$ and ${\bf\omega}^{q}$ remain zero. There are a large number of wall particles in the rough cases (96 million in case A2, 24 million in case A3). Fortunately, due to the regular structure of the locations of these particles, it is not necessary to store the positions of all wall particles. In the simulations, wall-particle positions were stored for a square of $9~\text{mm}^{2}$ ; from these positions the positions of the other wall particles were recomputed whenever they were required.

In the present paper, a rough wall denotes a wall with a bumpy surface, while a smooth wall denotes a wall with a flat surface. In both the smooth and rough cases, the coefficient of friction ${\it\mu}$ is 0.3 for particle–particle and particle–wall collisions. Thus, the smooth wall is in fact a flat frictional wall (Jenkins & Louge Reference Jenkins and Louge1997). In addition, it is remarked that bumpiness in the form of spheres fixed to a plane is a well-known concept in granular flow experiments (Pouliquen Reference Pouliquen1999). In the granular flow literature, a so-called specularity coefficient has been used to describe phenomenologically the diffusive scatter of particles caused by bumpiness of the surface of the boundary (Johnson & Jackson Reference Johnson and Jackson1987).

3.1. Velocity and volume fraction statistics

The basic characteristics of the mean velocity profiles and particle collision frequencies are shown in table 2. Columns two to six contain the gas friction velocity $u_{{\it\tau}}$ , gas bulk velocity $u_{b}=\langle \overline{u}_{1}\rangle$ , gas mean centreline velocity $u_{cl}$ , particle bulk velocity $v_{b}=\langle \overline{v}_{1}\rangle$ and particle mean centreline velocity $v_{cl}$ . The last two columns contain the overall collision frequencies. A p–p collision is a collision between two free particles somewhere in the domain. A p–w collision is a collision between a free particle and a wall or, in the rough-wall cases, a particle fixed on a wall. We will refer to the values in the table during the discussion of figures 1–4. Like in KFE1994 and in the work of Benson et al. (Reference Benson, Tanaka and Eaton2005), the mean and fluctuation velocity profiles in the figures in this section have been normalized with the mean centreline velocity of the unladen case, $u_{cl,0}=10.42~\text{m}~\text{s}^{-1}$ . A typical run of a particle-laden case with rough walls, which was performed on an eight-core (16-thread) processor of a Linux cluster, took approximately two weeks (wall-clock time); roughly 50 % of that time was used by the particle collision module.

Figure 1. Mean gas (a) and mean particle streamwise velocity profile (b), normalized with the unladen centreline velocity $u_{cl,0}$ . Unladen DNS A0 (dashed) and the laden simulations A1 (solid, no markers), A2 (solid, circles) and A3 (solid, triangles). The filled circles and filled triangles represent experiments of Kulick et al. (Reference Kulick, Fessler and Eaton1994) for copper particles at ${\it\phi}=0.4$ and ${\it\phi}=0.2$ respectively; the stars represent a fully rough experiment of Benson et al. (Reference Benson, Tanaka and Eaton2005) at ${\it\phi}=0.15$ . The unladen mean velocity profile is shown in both (a) and (b).

The mean velocity profiles of the two phases are shown in figure 1. According to KFE1994 the mean gas velocity profile in the experiments was not significantly modified in the laden cases up to a mass loading ratio of 0.4. The experimental uncertainty in the mean gas-phase velocity profiles was reported to be approximately 2 %. The near-wall flow, including $u_{{\it\tau}}$ , could not be measured in the laden cases. In the present simulations, the effect of the particles on the mean gas velocity profile appears also to be small, in particular in case A2 (rough wall, $10~{\rm\mu}\text{m}$ ); the mean gas velocity profiles of case A2 and the unladen case apparently coincide in figure 1(a). Quantitatively, the difference between the mean gas velocity profiles of A2 and A0 is less than 1 % for $u_{{\it\tau}}$ and $u_{cl}$ , while the maximum deviation, 2.6 %, is attained in the near-wall region (at $y^{+}\approx 7$ ). Figure 1(a) indicates that, compared with case A2, the particles in cases A1 and A3 have a slightly larger effect on the mean gas velocity profile, but quantitatively the effects are also small in these cases (see table 2).

The mean particle velocity profile is flatter than the mean gas velocity profile in simulations A1–A3. An explanation of this feature of wall-bounded flows laden with heavy particles was already given in the introduction. The wall roughness is observed to lower and to flatten the mean particle velocity, like in experiments (Benson et al. Reference Benson, Tanaka and Eaton2005). The mean particle velocity becomes lower because the rough wall causes more friction. The particle centreline velocity $v_{cl}$ is $0.86u_{cl,0}$ in case A2 and $0.70u_{cl,0}$ in case A3. The mean particle velocity becomes flatter due to stronger transverse mixing; the interaction of the particles with the rough wall enhances the wall-normal (transverse) particle velocity fluctuation in particular (Kussin & Sommerfeld Reference Kussin and Sommerfeld2002; Benson et al. Reference Benson, Tanaka and Eaton2005; Squires & Simonin Reference Squires and Simonin2006; Vreman Reference Vreman2007; Breuer et al. Reference Breuer, Alletto and Langfeldt2012), as will be shown in figure 3. The mean particle profiles in cases A2 and A3 are nearly horizontal, like the mean particle velocity profile measured in a fully rough channel (Benson et al. Reference Benson, Tanaka and Eaton2005), shown in figure 1(b). In addition, two mean particle velocity profiles from Kulick et al. (Reference Kulick, Fessler and Eaton1994) are shown in figure 1(b). Since no particle velocity statistics for ${\it\phi}=0.8$ are shown in KFE1994, we include experimental data for ${\it\phi}=0.4$ and ${\it\phi}=0.2$ . These data probably give indications for the behaviour at ${\it\phi}=0.8$ , since the data in KFE1994 indicate that particle velocity statistics are not as sensitive to the mass loading ratio as the gas velocity fluctuations. The KFE1994 mean particle velocity profile for ${\it\phi}=0.4$ is almost as low as the one of simulation A3, but apparently not as flat. According to Benson et al. (Reference Benson, Tanaka and Eaton2005), the mean particle velocity profiles in the experiments of KFE1994 were lowered due to (poorly defined) wall roughness in the development section of the flow. The flow had 0.30 m to recover from these effects, which corresponds to approximately 0.04 s (based on a particle bulk velocity of approximately $7.5~\text{m}~\text{s}^{-1}$ ). The time for recovery was significantly smaller than the Stokes response time ( ${\it\tau}_{p}=0.131~\text{s}$ ).

Figure 2 shows that wall roughness of sufficiently small size is able to enhance particle-induced turbulence attenuation in vertical channels. The root-mean-square gas velocity fluctuations, $\text{RMS}(u_{i})$ , and the Reynolds shear stress, $R_{12}$ , are shown; by definition $\text{RMS}(u_{i})=R_{ii}^{1/2}$ , where $R_{ij}=\overline{u_{i}^{\prime }u_{j}^{\prime }}$ is the Reynolds stress tensor. Figure 2 clearly shows attenuation of the turbulence in all three laden cases. In each of the laden cases, the three fluctuations and the Reynolds shear stress are attenuated compared with the unladen counterparts. The attenuation in the rough cases (A2 and A3) is clearly stronger than in the smooth case (A1). The strongest turbulence attenuation is found at the channel centre, where the root mean square of the streamwise velocity fluctuation in the rough cases is reduced by approximately 75 % of the unladen value. As indicated in the introduction, enhancement of turbulence attenuation by wall roughness has been observed in experiments on horizontal particle-laden channel flow (Kussin & Sommerfeld Reference Kussin and Sommerfeld2002). In these experiments the size of the wall roughness was also smaller than the viscous length scale, the maximum peak-to-trough roughness was $k_{l}^{+}\approx 0.5$ . In the experiments on vertical particle-laden channel flow by Benson et al. (Reference Benson, Tanaka and Eaton2005), however, the turbulence was amplified in the fully rough case. The mass loading ratio in this experiment was relatively low (0.15), while the walls were roughened by a wire mesh with $250~{\rm\mu}\text{m}$ thick wires in a square pattern with squares of $1~\text{mm}^{2}$ . This corresponds to a peak-to-trough roughness of $k_{l}^{+}\approx 8$ and a roughness standard deviation of $k_{rms}^{+}\approx 3$ . Roughness of this size might have modified the turbulence in the unladen-flow experiment, see the criteria of Flack & Schultz (Reference Flack and Schultz2014), discussed in the introduction.

Figure 2. Root mean square of gas velocity fluctuations (a–c) and gas Reynolds shear stress (d), normalized with the unladen centreline velocity $u_{cl,0}$ . Unladen DNS A0 (dashed) and the laden simulations A1 (solid, no markers), A2 (solid, circles) and A3 (solid, triangles). The filled triangles represent the experiment of Kulick et al. (Reference Kulick, Fessler and Eaton1994) for copper particles at ${\it\phi}=0.8$ .

Figure 3. Root mean square of particle velocity fluctuations (normalized with the unladen centreline velocity $u_{cl,0}$ ) (a–c) and mean particle volume fraction (d). Laden simulations A1 (solid, no markers), A2 (solid, circles) and A3 (solid, triangles). The filled circles represent the experiment of Kulick et al. (Reference Kulick, Fessler and Eaton1994) for copper particles at ${\it\phi}=0.4$ ; the stars represent a fully rough experiment of Benson et al. (Reference Benson, Tanaka and Eaton2005) at ${\it\phi}=0.15$ . The filled squares denote representative values of the very flat particle fluctuation profiles in the high-roughness experiment of Kussin & Sommerfeld (Reference Kussin and Sommerfeld2002) for $d_{p}=100~{\rm\mu}\text{m}$ and ${\it\phi}=0.5$ .

Figure 4. Mean (a) and RMS (b) of ${\it\omega}_{3}$ , the spanwise component of the particle angular velocity, normalized with $u_{cl,0}/(0.316d_{p})$ . Laden simulations A1 (solid, no markers), A2 (solid, circles) and A3 (solid, triangles).

Compared with the unladen simulation, the centreline value of the streamwise (wall-normal) gas velocity fluctuation in case A1 is reduced by 34 % (44 %). The turbulence attenuation in this case is stronger than predicted by the two-way coupled LES (Segura Reference Segura2004), which did not include inter-particle collisions. In that work the centreline streamwise (wall-normal) velocity fluctuation was reduced by 15 % (25 %) at ${\it\phi}=0.8$ . Thus, two-way coupled PP-DNS with inter-particle collisions predicts stronger turbulence attenuation than two-way coupled PP-LES without inter-particle collisions. Nevertheless, the discrepancy between the simulation data for smooth channels and the available experimental data for turbulence attenuation in vertical particle-laden channel flow remains large. It is stressed that we do not conclude that the turbulence attenuation found in KFE1994 was enhanced by wall roughness. In these experiments, only the development section was probably rough. For a mass loading ratio of ${\it\phi}=0.15$ , Benson et al. (Reference Benson, Tanaka and Eaton2005) showed that, unlike the mean particle velocity, the gas-phase statistics measured in a relatively smooth measurement section were hardly influenced by a rough development section. In addition, the roughness size of the walls of the development section in the experiments of Benson et al. (Reference Benson, Tanaka and Eaton2005), and probably in the experiments of KFE1994, was in a different regime from the roughness size used in the present simulations. However, we can conclude, after detailed reading of Benson et al. (Reference Benson, Tanaka and Eaton2005), that there is probably a lack of experimental data on turbulence attenuation in particle-laden flows in vertical channels for which both the development and measurement section are known to be smooth or have the same specified wall roughness. The only experiments on such vertical particle-laden channel flows seem to be the experiments at ${\it\phi}=0.15$ of Benson et al. (Reference Benson, Tanaka and Eaton2005). However, attenuation of turbulence is not reported there; no (clear) attenuation of turbulence is observed when the centreline streamwise velocity fluctuations of the laden smooth case and the unladen smooth case are compared (the filled right-pointing triangles in figures 3 and 7 in the reference).

It is possible that the quantitative discrepancy between the turbulence attenuation in simulation and experiment is caused by a limitation of the simulation method: PP-DNS does not resolve the boundary layers around the particles. As the wall is approached, the ratio of particle diameter and Kolmogorov length scale $d_{p}/{\it\eta}$ , which is approximately 0.45 at the centre, increases. It equals 1 at $x_{2}\approx 0.09H$ and attains a maximum of approximately 1.55 at the wall. This applies to the unladen case; the ratio is slightly reduced in the laden cases. With increasing $d_{p}/{\it\eta}$ , the error in the fluctuating part of the particle drag force model increases. This error was shown to be between 15 % and 30 % at $d_{p}/{\it\eta}=2$ (Burton & Eaton Reference Burton and Eaton2005). Due to this error, the dissipation via the particle-induced source term in the turbulence kinetic energy equation, introduced in § 3.3, is often underpredicted in point-particle simulations (Hwang & Eaton Reference Hwang and Eaton2006; Balachandar & Eaton Reference Balachandar and Eaton2010). However, as mentioned in § 2, PP-DNS is the most advanced method to date applicable to a simulation at this Reynolds number and particle size; a particle-resolved DNS of the statistically stationary state of a channel flow at $Re_{{\it\tau}}\approx 640$ with more than $10^{5}$ particles sized around the Kolmogorov scale is not feasible yet. In addition, a particle-resolved DNS will still have to deal with the issue of how an inelastic collision between two solid particles in a gas or between a (rough) wall and a solid particle in gas should be modelled.

The discussion above shows that more research is needed to draw definitive conclusions about the reason for the large quantitative discrepancy between simulation results for smooth vertical channels and available experimental results. On the one hand, it cannot be ruled out that the rough development section in KFE1994 had an effect on the turbulence attenuation statistics. On the other hand, the simulated turbulence attenuation probably underpredicts the turbulence attenuation by an unknown amount. It is not within the scope of the present paper to settle this issue. The purpose of this paper is to show the qualitative effect of wall roughness of relatively small size on turbulence attenuation in vertical particle-laden channel flow and to show that uniformity of the mean feedback force can contribute to turbulence attenuation. Since at least the qualitative effect of particles on the turbulence is reproduced by the simulations, the fundamental limitation of PP-DNS is not expected to alter the conclusions of the present paper.

Before we proceed to analyse the simulated turbulence attenuation in more detail in the next sections, we discuss a few more statistics of the particles. Experiments have shown that the streamwise and wall-normal spanwise particle velocity fluctuations are strongly increased by wall roughness (Kussin & Sommerfeld Reference Kussin and Sommerfeld2002; Benson et al. Reference Benson, Tanaka and Eaton2005). The same trend is found in the simulations, see figure 3(a–c), in which experimental results are included for comparison. It is observed that the spanwise particle velocity fluctuation, for which no experimental results could be found, is only approximately 50 % of the wall-normal particle velocity fluctuation, in the cases with rough walls. Figure 3(b) shows that the wall-normal fluctuation in the experiment of KFE1994 quickly recovered from the wall roughness of the development section that, according to Benson et al. (Reference Benson, Tanaka and Eaton2005), caused the low mean particle velocity profile. An explanation could be that if particles entered the measurement section with a wall-normal velocity of approximately $1~\text{m}~\text{s}^{-1}$ (the level of the fluctuation in the fully rough experiment), they had a high probability to collide with the smooth wall of the measurement section before the velocity was measured (the channel width was 0.04 m and the residence time in the measurement section before measurement was approximately 0.04 s).

The simulated mean particle volume fraction profiles, shown in figure 3(d), are fairly uniform. A weak tendency of particles to concentrate at the channel centre is observed, most clearly in the case with smooth walls (A1).

The increase of the kinetic energy in the particle velocity fluctuations with increasing wall roughness is accompanied by an increase of the overall collision frequencies, which are shown in the last columns of table 2. The first observation is that the collision frequencies strongly increase with increasing wall roughness. The second observation is that a significant proportion of all particle collisions are p–p collisions, 58 %, 34 % and 33 % in cases A1, A2 and A3 respectively. Thus, if p–w collisions are included, it seems appropriate to include also p–p collisions in simulations of this type, although the overall particle volume fraction is low. Others have shown that p–p collisions significantly flatten the mean particle velocity profile and the mean particle volume fraction velocity profile at low volume fractions (Li et al. Reference Li, McLaughlin, Kontomaris and Portela2001; Yamamoto et al. Reference Yamamoto, Potthoff, Tanaka, Kajishima and Tsuji2001; Vreman Reference Vreman2007; Nasr et al. Reference Nasr, Ahmadi and McLaughlin2009). At larger volume fractions p–p collisions were shown to enhance the transfer from the kinetic energy in the streamwise particle velocity fluctuations to the kinetic energy in both the normal and spanwise particle velocity fluctuations (Vreman et al. Reference Vreman, Geurts, Deen, Kuipers and Kuerten2009).

The mean and fluctuation of the spanwise component of the angular velocity of the particles are shown in figure 4. The mean streamwise and mean normal components are zero, while the fluctuation of each of these components is roughly 50 % of the spanwise fluctuation. The magnitude of the angular velocity increases with wall roughness. This observation and the observation that the angular velocity profiles peak at the wall indicate that particle–wall collisions play an important role in the production of angular velocity. It is remarked that the total kinetic energy of a spherical particle is given by $(m_{p}(|\boldsymbol{v}^{p}|^{2}+(d_{p}^{2}|{\bf\omega}|^{2})/10))/2$ . For this reason the particular angular velocity profiles were multiplied by $(1/10)^{1/2}d_{p}/u_{cl,0}$ to allow comparison with the particle translative velocity fluctuations shown in figure 3. The angular kinetic energy appears to be roughly 10 % of the kinetic energy in the translative velocity fluctuations.

3.2. The feedback force and the streamwise mean momentum equation

The attenuation of the gas-phase turbulence observed in the simulations must be caused by the feedback force, the term $\boldsymbol{f}$ in the equations that govern the continuous phase. To investigate this term in more detail, it is decomposed into three parts,

The three parts are the uniform part of the mean feedback force, the non-uniform part of the mean feedback force,

and the fluctuating part of the feedback force, $\boldsymbol{f}^{\prime }$ . The normal and spanwise components of the first two parts are negligible. The decomposition of the streamwise component of the feedback force is shown in the second column of table 3 (the uniform part of the mean force), in figure 5(a) (the non-uniform part of the mean) and in figure 5(b) (the root-mean-square profile of the fluctuation).

Figure 5. Non-uniform part (a) and RMS (b) of the particle drag force in the streamwise component of the gas momentum equation, $\overline{f}_{1}^{nu}$ and $\text{RMS}(f_{1})$ , normalized with ${\it\rho}u_{{\it\tau},0}^{2}/H$ . Laden simulations A1 (solid, no markers), A2 (solid, circles) and A3 (solid, triangles). In (a) an estimate based on experimental results (KFE1994) is also shown, as described in the text.

Table 3 shows an overview of all terms in the domain-averaged streamwise momentum equations. The sum of the terms is zero in each case, for both the gas and particle momentum equations. The force exerted by the wall on the gas strongly increases with wall roughness. As a result the mean particle velocity profile becomes lower, and the particles exert a negative domain-averaged feedback force on the gas in the rough cases. Therefore, the streamwise pressure gradient ${\it\rho}\overline{a}$ needs to increase to maintain the constant gas bulk velocity in the rough cases. The effect of a larger forcing term $\overline{a}$ in single-phase flow is a larger friction velocity $u_{{\it\tau}}$ and also a larger bulk velocity $u_{b}$ and a larger mean centreline velocity $u_{cl}$ or, in other words, a larger Reynolds number of the flow.

By definition, the uniform part of the mean feedback force does not depend on time and space, therefore it affects the flow in a similar way to $\overline{a}$ ; in fact $g+\overline{a}+\langle \,\overline{f}_{1}\rangle$ could be lumped into a single term. The magnitude of that single term sets $u_{{\it\tau}}$ and thereby the Reynolds number $Re_{{\it\tau}}$ of the flow. Thus, the effect of the mean feedback force is a change of the Reynolds number if the mean streamwise pressure gradient ( ${\it\rho}\overline{a}$ ) is not changed. As discussed in § 2.6, $a$ is modified to maintain a constant volume flow and to have the same bulk velocity in each case. In this way the mean feedback force is not able to change the bulk Reynolds number, while the variation of $u_{{\it\tau}}$ (and $Re_{{\it\tau}}$ ) among the different cases remains small (table 2). Comparison of the values of $\overline{a}$ in table 3 with the unladen value indicates that if we had fixed the pressure gradient, the (bulk) Reynolds number would have increased in the smooth laden case but decreased in the rough laden cases.

The non-uniform part of the feedback force, shown in figure 5(a), appears to be a monotonically decreasing function in the channel half $(0,H)$ for simulations A1–A3. The slope of the profile is steep in the near-wall region and zero at the centre. Thus, the non-uniformity is much stronger near the wall than in the centre region. The non-uniformity of the feedback force is apparent in each laden simulation and increases with increasing wall roughness. Provided that the particle volume fraction is uniform, the non-uniformity of the profile is directly related to the uniformity of the mean relative velocity profile, $\overline{v}_{1}-\overline{u}_{1}$ . As explained in the introduction, the mean particle velocity tends to be flatter than the mean gas velocity, such that the mean relative velocity tends to be non-uniform. It is remarked that the definition of $\overline{f}_{1}^{nu}$ implies that the volume average or the integral over $x_{2}$ of $\overline{f}_{1}^{nu}$ is zero. This means that even for a profile $\overline{f}_{1}$ that is non-uniform in the near-wall region only, the non-uniform part $\overline{f}_{1}^{nu}$ can be non-zero in the centre region. The effect of the non-uniformity on the turbulence is not only local. Instead of the slope (derivative) of $\overline{f}_{1}^{nu}$ , it is the primitive of $\overline{f}_{1}^{nu}$ that will appear to be tied to the Reynolds shear stress. As a global measure of the non-uniformity of the feedback force, we introduce the parameter ${\it\beta}$ ,

which is in fact the normalized $L_{1}$ -norm of $\overline{f}_{1}^{nu}$ . For simulations A1, A2 and A3, we find ${\it\beta}=0.43$ , ${\it\beta}=0.73$ and ${\it\beta}=0.87$ respectively.

Very flat mean particle velocity profiles have been observed in experiments in rough channels (Benson et al. Reference Benson, Tanaka and Eaton2005). Kulick et al. (Reference Kulick, Fessler and Eaton1994) also reported that the mean particle velocity profile was flatter than the mean gas velocity profile, in particular in the near-wall region (see figure 1 b). It is appropriate to consider figure 1(b) again, since according to the experimental data for ${\it\phi}=0.4$ the slope of the mean particle velocity profile is not smaller than the slope of the unladen mean gas velocity profile for $0.7H<x_{2}<H$ . To show that nonetheless the experiments in KFE1994 do support the notion of a non-uniform part of the mean feedback force, we provide an estimate of $\overline{f}_{1}^{nu}$ based on experimental data. We define

and substitute ${\it\phi}=0.8$ and ${\it\tau}_{p}=0.131~\text{s}$ . Since for the mean particle velocity profile at ${\it\phi}=0.8$ no experimental data are available, we approximate $\overline{v}_{1}$ by the data for ${\it\phi}=0.4$ and ${\it\phi}=0.2$ . The two approximated $\overline{v}_{1}$ profiles are linearly interpolated from the experimental data. Between the wall and the measurement point nearest to the wall, the profile is horizontally extrapolated from the measurement point. This is also reasonable for the ${\it\phi}=0.4$ data, since the ${\it\phi}=0.2$ data show that the mean particle velocity profile is flat in the near-wall region, except for a cusp very close to the wall. For $\overline{v}_{1}$ taken from the ${\it\phi}=0.4$ experiment, $\overline{f}_{exp}-\langle \overline{f}_{exp}\rangle$ is shown in figure 5(a) (the domain-averaged factor in (3.4) is then 1.69). The shape of the profile is more or less similar to the shapes of $\overline{f}_{1}^{nu}$ obtained from simulations A1–A3, although it is not monotonically decreasing on the entire channel half. The corresponding global non-uniformity parameter ${\it\beta}$ is equal to 0.59, in between those of simulations A1 and A2. Alternatively, if $\overline{v}_{1}$ is taken from the ${\it\phi}=0.2$ experiment, we find ${\it\beta}=0.47$ . Since the profile of the non-uniform part of the feedback force obtained with the latter estimate is very similar to the profile of simulation A1, we do not show the curve in figure 5(a), to keep the figure clearer. It is concluded that the existence of a non-uniform part of the mean feedback force is supported by experimental data.

Since the Reynolds shear stress is strongly modified by particles (figure 2 d), the balance of the mean streamwise momentum of the gas should also be modified. Owen (Reference Owen1969) noted the importance of the mean streamwise momentum equation, since in the context of the reduction of the gas turbulence in wall-bounded particle-laden flows he wrote in a footnote ‘in a statistically steady flow, an increasing proportion of the shear stress developed in the gas phase is transferred to the particles, in the form of a momentum flux, as the wall is approached’. In the following we consider the mean streamwise momentum equation of the gas in detail. For lower $Re_{{\it\tau}}$ this has been done before by Mito & Hanratty (Reference Mito and Hanratty2006) in PP-DNS of channel flow ( $Re_{{\it\tau}}=150$ ) and by Vreman (Reference Vreman2007) in PP-DNS of pipe flow ( $Re_{{\it\tau}}=140$ ). The derivation below shows similarities with the derivations in these papers.

The streamwise component of the mean momentum equation of the gas phase can be written as

The domain-averaged mean momentum equation of the gas phase is then derived as

If this equation is subtracted from (3.5), we obtain the following equation for the Reynolds shear stress:

This is an interesting equation, as it shows the dependence of the Reynolds shear stress on the non-uniform part of the mean feedback force $\overline{f}_{1}^{nu}$ . This term is zero for a uniform mean feedback force, and then $R_{12}-{\it\nu}\,\text{d}\overline{u}_{1}/\text{d}x_{1}$ is a straight line with slope $u_{{\it\tau}}^{2}/H$ , like in the unladen flow. After non-dimensionalization with ${\it\rho}u_{{\it\tau}}^{2}/H$ , (3.7) has the form

All terms in (3.7) have been computed for cases A1–A3: the first term on the right-hand side is shown in table 3 (after multiplication by the density), the last term on the right-hand side is shown in figure 5(a) and the left-hand side and second term on the right-hand side are shown in figure 6. The mean viscous term appears to be hardly influenced by particles and is negligible at the centre (figure 6 b); the normalized centre value is roughly $-0.02$ . Provided that the mean viscous term at the centre remains negligible, the slope of the Reynolds shear stress normalized with $u_{{\it\tau}}^{2}/H$ must decrease at the centre of the channel due to the negativity of the non-uniform part of the mean feedback force at the centre. In the present cases the effect of the particles on $u_{{\it\tau}}^{2}$ is not very large; in particular, in case A2 the effect is small, approximately 1 %.

Figure 6. Differentiated Reynolds shear stress (a) and the mean of the viscous term in the streamwise component of the gas momentum equation (b), normalized with ${\it\rho}u_{{\it\tau},0}^{2}/H$ . Laden simulations A1 (solid, no markers), A2 (solid, circles) and A3 (solid, triangles).

The integration of (3.8) between position $x_{2}$ and $H$ yields

In the derivation the properties $R_{12}=0$ and $\text{d}\overline{u}_{1}/\text{d}x_{2}=0$ at $x_{2}^{\prime }=H$ have been used. Since $\overline{f}_{1}^{nu}$ is a monotonically decreasing function of between 0 and $H$ and its integral between 0 and $H$ is by definition zero, we have

Thus, in the case where the mean velocity profile does not change and the non-uniform part of the mean feedback force is monotonically decreasing in the interval $(0,H)$ , the effect of particles is a reduction of minus the Reynolds shear stress ( $-R_{12}$ ) in the interval $(0,H)$ . As long as the reduction of $-R_{12}$ is such that $-R_{12}$ remains positive in the interval $(0,H)$ , this reduction is a reduction of the absolute Reynolds shear stress and thus an attenuation of an important statistical quantity of the gas turbulence. Conversely, if the Reynolds shear stress is attenuated but the mean gas velocity profile remains the same, then (3.9) shows that the mean feedback force must be non-uniform. If we define the attenuation of the Reynolds shear stress by ${\rm\Delta}R_{12}=R_{12}-R_{12}^{0}$ ( $R_{12}^{0}$ is the unladen profile), then the assumption of unchanged mean velocity profile implies

which is a direct relation between the attenuation of the Reynolds shear stress and the non-uniform part of the mean feedback force. The mean gas velocity profiles in KFE1994 are not modified by particles for mass loading ratios up to ${\it\phi}=0.4$ . Equation (3.11) shows that in cases where the mean gas velocity remains the same, a relatively large attenuation of $R_{12}$ is equivalent to a relatively large non-uniform part of the mean feedback force. Since the turbulence attenuation of streamwise and wall-normal gas velocity fluctuations was very strong in the KFE1994 experiments, the attenuation of the Reynolds shear stress was probably also strong (the Cauchy–Schwarz inequality implies $|R_{12}|\leqslant \text{RMS}(u_{1})\text{RMS}(u_{2})$ ). It seems therefore that also the non-uniform part of the mean feedback force in these experiments must have been quite large, despite the fact that the measured mean particle velocity profile was not yet as flat as in simulations and experiments with (fully) rough walls.

Although the derivation above shows that the Reynolds shear stress and the non-uniform part of the mean feedback force are related, it does not prove that this non-uniformity causes the turbulence attenuation or part of the turbulence attenuation. The combination of turbulence attenuation and an unchanged mean velocity profile could also produce the non-uniformity of the drag force. Intuitively this may not be very likely, but at least there is another cause of turbulence attenuation that should be investigated, which is the fluctuating drag force. Figure 5 shows that the magnitude of the fluctuation of the drag force $\boldsymbol{f}^{\prime }$ is larger than the magnitude of the non-uniform part of the drag force. To investigate the role of the fluctuating drag force, we proceed with an investigation of the turbulence kinetic energy budget.

3.3. The feedback force and the turbulence kinetic energy budget

Unlike the mean feedback force, the fluctuating feedback force occurs in the turbulence kinetic energy equation of the gas phase. The corresponding term has in general a dissipative effect. It has hitherto been regarded as the mechanism by which the particle-induced turbulence attenuation should be explained.

The turbulence kinetic energy of the gas phase is defined by $K=(\overline{\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}^{\prime }})/2=(R_{11}+R_{22}+R_{33})/2$ . The budget or balance equation of $K$ is given by

where

Here, $P_{K}$ is the turbulence production, $D_{K}$ is the sum of turbulent diffusion, pressure diffusion and viscous diffusion, ${\it\epsilon}_{K}$ is minus the turbulence dissipation and $S_{K}$ is a particle-induced energy source term (often negative). The term $A_{K}$ is the correlation between the streamwise velocity fluctuation and the fluctuation of the domain-averaged pressure gradient. Since $A_{K}$ is expected to be small, the term is usually omitted, but its negligibility should be verified if $a$ is a function of time.

Figure 7. Turbulence kinetic energy budget in the near-wall region: turbulence production $P_{K}$ (a), total diffusion $D_{K}$ (b), turbulence dissipation ${\it\epsilon}_{K}$ (c) and particle-induced source term $S_{K}$ (d), normalized with $u_{{\it\tau},0}^{3}/H$ . Unladen DNS A0 (dashed) and the laden simulations A1 (solid, no markers), A2 (solid, circles) and A3 (solid, triangles).

The budget of turbulence kinetic energy is shown in figure 7 (for the near-wall region) and in figure 8 (for the core region). The term $A_{K}$ is not shown, since its absolute value normalized with $u_{{\it\tau},0}^{3}/H$ was less than 0.003, negligible indeed. It is observed that all other terms in the budget are strongly affected by the particles. In particular, the turbulence production is reduced across the entire channel, the turbulence dissipation is reduced, while the reduction of $D_{K}$ in the centre region of the channel shows that also turbulence transport from the wall to the centre is strongly inhibited by the particle motion. The term $S_{K}$ is interesting in particular; it has since long been recognized that this term plays an important role in turbulence attenuation (Elghobashi & Truesdell Reference Elghobashi and Truesdell1993). It is the correlation between the fluctuating feedback force and the fluctuating gas velocity. The correlation tends to be negative if the standard drag law is used and the particle fluctuation and fluid fluctuation are decorrelated. If the second condition is not fulfilled, the term may be positive. In regions with weak gas turbulence but strong particle velocity fluctuation, the instantaneous oscillations of gas regions produced by particles that cross the region may be larger than the turbulence generated by (nonlinear) channel flow instabilities. Gas fluctuations created by particle fluctuations are probably positively correlated, which tends to make $S_{K}$ positive. Very close to the wall $S_{K}$ is positive for A1–A3, while at the channel centre $S_{K}$ is positive in cases A2 and A3. In cases A2 and A3, the turbulence attenuation is strong, and the particle velocity fluctuations are large. Thus, the particles produce some turbulence kinetic energy, which is immediately dissipated by viscosity. The latter contributes to the viscous dissipation, which possibly explains why the centreline dissipation of case A3 is hardly reduced compared with the (small) unladen dissipation.

Figure 8. Turbulence kinetic energy budget in the core region: turbulence production $P_{K}$ (a), total diffusion $D_{K}$ (b), turbulence dissipation ${\it\epsilon}_{K}$ (c) and particle-induced source term $S_{K}$ (d), normalized with $u_{{\it\tau},0}^{3}/H$ . Unladen DNS A0 (dashed) and the laden simulations A1 (solid, no markers), A2 (solid, circles) and A3 (solid, triangles).

Perhaps the most surprising feature of the energy budget equation is that strong turbulence attenuation is possible, while the quantity $S_{K}$ remains relatively small. In case A3, the domain average of $S_{K}$ is approximately zero, such that the term has no net dissipative effect, although the turbulence attenuation is strong. As discussed in § 3.1, Hwang & Eaton (Reference Hwang and Eaton2006) found that a point-particle model of $S_{K}$ underpredicts the dissipation by particles in isotropic turbulence if $d_{p}\gtrapprox {\it\eta}$ , see also Balachandar & Eaton (Reference Balachandar and Eaton2010). Nevertheless, the observation that a simulation of the Navier–Stokes equations shows strong turbulence attenuation in a case where the feedback term has no net dissipative effect remains surprising.

We finish this section with a brief discussion of the influence of the Reynolds number on turbulence attenuation. In the present PP-DNS of gas–solid channel flow with smooth walls at mass loading ratio 0.8 and $Re_{{\it\tau}}=642$ , the turbulence is significantly reduced but still relevant. However, in a PP-DNS of gas–solid pipe flow with smooth walls at a similar mass loading ratio, 0.77, but lower Reynolds number, $Re_{{\it\tau}}=140$ , the gas-phase turbulence production, the Reynolds shear stress and also the radial diagonal component of the Reynolds stress tensor were negligibly small, less than 10 % of the unladen values (Vreman Reference Vreman2007). Although pipe flow and channel flow are not equivalent, this comparison seems to suggest that turbulence attenuation becomes less pronounced with increasing Reynolds number.

4.1. Direct numerical simulation of turbulent channel flow with non-uniform forcing

Turbulence attenuation in particle-laden channel flows is often attributed (solely) to the fluctuating feedback force, since that is the only part of the feedback force that appears in the particle-induced source term $S_{K}$ in the turbulence kinetic energy equation. However, (4.7) shows that an important statistical quantity of the turbulence, the Reynolds shear stress, is related to the primitive of the non-uniform part of the mean feedback force, in particular if the mean gas velocity profile is not much affected. This leads to the hypothesis that a non-uniform part of the mean feedback force acts as one of the causes of turbulence attenuation in particle-laden channel flow. So far, this is only an hypothesis, because the possibility that the non-uniform part of the mean feedback force is not a cause but just a consequence of turbulence attenuation has not been ruled out yet.

To demonstrate that a non-uniform feedback force is able to attenuate the turbulence by itself, three additional direct numerical simulations are performed, B1–B3, in which the effect of the mean feedback force is isolated. Unlike A1–A3, these simulations do not use the point-particle model, at least not explicitly. Simulation B1 is a DNS of the Navier–Stokes equations with forcing $\boldsymbol{f}=f_{1}(x_{2})\boldsymbol{e}_{\mathbf{1}}$ . Thus, in these cases $f_{2}=f_{3}=0$ and $f_{1}$ is a function of $x_{2}$ only. The function $f_{1}$ in simulation B1 is precisely the mean feedback force of simulation A1. Similarly, simulations B2 and B3 use the mean feedback forces of simulations A2 and A3 respectively. The fluctuating feedback force is not used in simulations B1–B3. Like in simulations A1–A3, the domain average of the streamwise pressure gradient ( ${\it\rho}a$ ) is adapted to fix the bulk velocity in B1–B3 to $9.212~\text{m}~\text{s}^{-1}$ . Thus, the turbulence modification in B1 (B2, B3) compared with case A0 can only be caused by the profile of the non-uniform part of the mean feedback force A1 (A2, A3). That profile is shown in figure 5(a).

The results of simulations B1–B3 are shown in figure 9. It appears that the turbulence is also significantly attenuated in these cases, although the attenuation is not as strong as in cases A1–A3. For example, in case B1 the centreline value of the streamwise (wall-normal) fluctuation is reduced by 33 % (26 %); in case A1 it was reduced by 34 % (44 %). The spanwise velocity fluctuation, not shown for these cases, is attenuated in a similar way to the wall-normal velocity fluctuation. The attenuation observed in cases B1–B3 cannot be caused by $S_{K}$ , the particle-induced source term in the turbulence kinetic energy equation, since this term is zero in these cases ( $\boldsymbol{f}^{\prime }=0$ ). This means that the attenuation in B1–B3 must be caused by the mean streamwise feedback force. Wall roughness makes the mean particle force more non-uniform (figure 5 a), and simulations B1–3 show that stronger non-uniformity leads to more turbulence attenuation.

Figure 9. Mean streamwise velocity $\overline{u}_{1}$ (a), $\text{RMS}(u_{1})$ (b), $\text{RMS}(u_{2})$ (c) and Reynolds shear stress $R_{12}$ (d), normalized with the unladen centreline velocity $u_{cl,0}$ . Results from four direct numerical simulations: A0 (dashed), B1 (solid, no markers), B2 (solid, circles) and B3 (solid, triangles). The forcing in B1–B3 is $\overline{f}_{1}\boldsymbol{e}_{\mathbf{1}}$ with the mean force $\overline{f}_{1}$ taken from A1, A2 and A3 respectively.

While the difference between case B1 (B2, B3) and the unladen case A0 gives an estimate of the contribution of the mean feedback force to the turbulence modification in case A1 (A2, A3), the difference between case A1 (A2, A3) and case B1 (B2, B3) gives an estimate of the contribution of the fluctuation of the feedback force to the turbulence modification in case A1 (A2, A3). Both contributions turn out to be significant.

Compared with the modification of the velocity fluctuations and the Reynolds shear stress, the modification of the mean velocity profile in cases B1–B3 is small, in particular in case B1; in cases B2 and B3 the effects on the linear profile are somewhat larger than in cases A2 and A3. The friction velocity $u_{{\it\tau}}$ is equal to $0.490~\text{m}~\text{s}^{-1}$ , $0.496~\text{m}~\text{s}^{-1}$ , $0.501~\text{m}~\text{s}^{-1}$ and $0.504~\text{m}~\text{s}^{-1}$ in cases A0, B1, B2 and B3 respectively. The bulk velocity is by definition the same in each case. The centreline velocity $u_{cl}$ is equal to $10.42~\text{m}~\text{s}^{-1}$ , $10.40~\text{m}~\text{s}^{-1}$ , $10.08~\text{m}~\text{s}^{-1}$ and $9.97~\text{m}~\text{s}^{-1}$ in cases A0, B1, B2 and B3 respectively. Thus, with increasing non-uniformity of the forcing, $u_{{\it\tau}}$ increases slightly and $u_{cl}$ decreases slightly; in other words the mean profile becomes less rounded. How can this be reconciled with case A2, in which both $u_{{\it\tau}}$ and $u_{cl}$ were almost identical to $u_{{\it\tau},0}$ and $u_{cl,0}$ ? The smaller change of the mean gas velocity profile in case A2 compared with B2 can only be explained by an interplay of the non-uniform part of the mean feedback force and the fluctuating feedback force. Both parts attenuate the turbulence. However, the non-uniform mean part tends to make the mean profile less rounded, while the fluctuating part tends to make the mean profile more rounded. These two effects on the mean gas velocity profile cancel to certain extent; this is most clearly visible in case A2.

4.2. Linear stability analysis

As another illustration of the response of the Navier–Stokes equations to a feedback force with a non-uniform mean part, a linear stability analysis is performed. We consider the linear instability of a laminar velocity profile to two-dimensional perturbations, as described by the Orr–Sommerfeld equations (Criminale, Jackson & Joslin Reference Criminale, Jackson and Joslin2003). In linear stability theory the flow is decomposed into a base profile $U(x_{2})$ and a velocity perturbation fluctuation proportional to the Fourier mode $\exp (\text{i}({\it\alpha}_{1}x_{1}-{\it\gamma}t))$ , where $\text{i}=(-1)^{1/2}$ , ${\it\alpha}_{1}$ is a prescribed streamwise spatial wavenumber and the frequency ${\it\gamma}$ is the complex eigenvalue.

For simplicity of the analysis we consider a linear version of the feedback force, $\boldsymbol{f}=c(\boldsymbol{v}-\boldsymbol{u})$ , with $c=({\it\phi}/{\it\tau}_{p})(1+0.15R^{0.687})$ , where the constant $R$ is a representative particle Reynolds number which does not depend on time and space. As a further simplification, we assume that the particles move with a uniform velocity, $\boldsymbol{v}=V\boldsymbol{e}_{\mathbf{1}}$ , where $V$ is a constant. Since $c$ is uniform, the value of the uniform particle velocity $V$ is irrelevant for the flow dynamics, since $cV$ , $g$ and $a$ can be lumped into a single parameter $b=g+a+cV$ .

The unperturbed momentum equation of the linearized problem becomes

and $U(0)=U(2H)=0$ . The analytical solution of this equation is

where $U_{0}$ denotes the centreline velocity and ${\it\lambda}^{2}=c/{\it\nu}$ . The constants $b$ and $U_{0}$ are related by

Three examples of base profiles $U(x_{2})$ are shown in figure 10(a). For ${\it\lambda}\rightarrow 0$ , the profile converges to the standard parabolic profile $U_{0}x_{2}(2H-x_{2})/H^{2}$ , while for ${\it\lambda}\rightarrow \infty$ , the profile converges to $U_{0}$ for $0<x_{2}<2H$ .

The Orr–Sommerfeld equation for the eigenfunction of the normal velocity component, denoted by $\hat{u} _{2}(x_{2})$ , reads

where the prime denotes differentiation with respect to $x_{2}$ . The non-conventional term proportional to $c$ represents the fluctuating feedback force. The equation above can be derived from the evolution equation of the Laplacian of the normal component of the velocity perturbation, following the procedure described in Criminale et al. (Reference Criminale, Jackson and Joslin2003).

The linear growth rate Im( ${\it\gamma}$ ), the imaginary part of ${\it\gamma}$ , is shown in figure 10(b) for several cases, all for $Re=U_{0}H/{\it\nu}=14\,000$ , roughly the same as the Reynolds number based on the mean centreline velocity of the unladen DNS. It appears that a very modest flattening of the base flow ( ${\it\lambda}=1$ ) leads to a strong reduction of the most unstable growth rate. Comparison with a case with the mean feedback force, but without the fluctuating feedback force ( ${\it\lambda}=1$ in (4.2) and $c=0$ in (4.4)), shows that the linear stabilization is primarily caused by the mean feedback component which slightly changes the parabolic profile. The critical value for ${\it\lambda}$ is approximately 1.3; for ${\it\lambda}>1.3$ the system is linearly stable for all ${\it\alpha}_{1}$ .

Figure 10. (a) Base velocity profile in the linear stability theory (4.2) for ${\it\lambda}\rightarrow 0$ (red dashed), ${\it\lambda}=1$ (blue solid) and ${\it\lambda}=10$ (black squares). (b) Linear growth rate (normalized by $U_{0}/H$ ) as a function of the streamwise wavenumber ${\it\alpha}$ (normalized by $H^{-1}$ ) for ${\it\lambda}\rightarrow 0$ (red dashed), for ${\it\lambda}=1$ and $c={\it\lambda}^{2}{\it\nu}$ (blue solid), and for ${\it\lambda}=1$ and $c=0$ (black dash–dotted).

4.3. The equation of the velocity fluctuation

Since the mean feedback force does not occur in the turbulence kinetic energy equation or in any of the equations of the Reynolds stresses, the question may arise of how a non-uniform mean feedback force can contribute to significant turbulence attenuation. It is tempting to assume that the mean feedback force can only modify the turbulence via a modification of the mean gas velocity profile, which does explicitly appear in the Reynolds stress equations. However, the mean gas velocity profile is not necessarily significantly modified, consider for example case A2. According to this argument, the mean feedback force cannot play an important role in the turbulence attenuation, at least not in case A2. In the following we will indicate where this intuitively logical argument fails.

The presumption in the argument is that the Reynolds stress equations and the mean flow are all that is required to describe the turbulence. This is not correct. Even if the mean velocity profile is prescribed, the Reynolds stress equations do not form a closed set. The argument would be valid if it were possible to close the unclosed terms (for example turbulence transport, pressure strain and turbulence dissipation) by algebraic models in terms of the Reynolds stresses and the mean velocity profile. However, an accurate closure of this type is unknown and may not exist. Of course, additional transport equations can be derived for the unclosed expressions, but then these additional equations will suffer from similar closure problems.

These closure problems are avoided if we consider the system of equations for the mean and instantaneous fluctuation of the gas velocity. For simplicity, we assume that the forcing term $\boldsymbol{f}$ is known as input as a vector field that varies in time and space. The following system is then a closed system for the mean and fluctuation of the gas velocity:

where $u_{b,0}$ is the prescribed gas bulk velocity and the expression for $a^{\prime }$ follows from the requirement that the domain-averaged streamwise velocity remains constant, $\text{d}\langle u_{1}\rangle /\text{d}t=0$ . This is a complete system of integro-differential equations for the unknowns $\boldsymbol{u}^{\prime }$ , $P^{\prime }$ , $\overline{u}_{1}$ and $a^{\prime }$ , for given $\boldsymbol{f}$ , $u_{b,0}$ , initial condition for $\boldsymbol{u}^{\prime }$ and no-slip boundary conditions for $\boldsymbol{u}^{\prime }$ and $\overline{u}_{1}$ . The unknown $a^{\prime }$ can easily be eliminated by substitution of the last equation into (4.5).

In the above system only (4.5) contains a time derivative; only this equation is an evolution equation, the other equations are additional constraints. One constraint is the incompressibility constraint, another one is (4.7). In other words, the field $\boldsymbol{u}^{\prime }$ can only evolve such that the statistical average of $u_{1}^{\prime }u_{2}^{\prime }$ satisfies a constraint in which the non-uniform part of the mean feedback force explicitly appears. Thus, although the mean feedback force does not explicitly appear in (4.5), this does not mean that $\boldsymbol{u}^{\prime }$ would not depend on it. To clarify this further, we substitute (4.7) into (4.5), assume $a^{\prime }=0$ (then $\text{d}\langle u_{1}\rangle /\text{d}t$ is not necessarily zero), skip (4.8) and prescribe $u_{{\it\tau}}=u_{{\it\tau},0}$ . We obtain

which is a closed system for $\boldsymbol{u}^{\prime }$ , $P^{\prime }$ and $\overline{u}_{1}$ , for given $\overline{f}_{1}^{nu}$ , $\boldsymbol{f}^{\prime }$ , $u_{{\it\tau},0}$ , initial condition for $\boldsymbol{u}^{\prime }$ and no-slip boundary conditions for $\boldsymbol{u}^{\prime }$ and $\overline{u}_{1}$ . In this formulation the fluctuating velocity explicitly depends on the non-uniform part of the mean feedback force. If we contract (4.10) with $\boldsymbol{u}^{\prime }$ , the product $u_{1}^{\prime }\overline{f}_{1}^{nu}$ appears. The mean of this product is zero, and that is why it does not appear in the statistical turbulence kinetic energy equation. Although $\overline{f}_{1}^{nu}$ has no mean dissipative effect, it can modify the dynamics of $\boldsymbol{u}^{\prime }$ , and the consequence may be that statistical correlations that do appear in the turbulence kinetic energy equation (production, transport and dissipation) are reduced. The appearance of $\overline{f}_{1}^{nu}$ in the equation of $\boldsymbol{u}^{\prime }$ is a nonlinear effect. In linear stability theory all terms that are quadratic in $\boldsymbol{u}^{\prime }$ vanish; in particular, the left-hand side of (4.12) is then replaced by zero, such that the sum of the three terms between brackets in (4.10) also vanishes.