2.1 Turbulent channel flow

The governing equations for particle-laden incompressible flow in Cartesian tensor notation are given by

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}u_{i}}{\text{D}t}=-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}u_{i}}{\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}x_{j}}+\frac{1}{\unicode[STIX]{x1D70C}}f_{i}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

in which $t$ is time; $x_{1},x_{2}$ and $x_{3}$ are the streamwise $(x)$ , wall-normal $(y)$ and spanwise coordinates $(z)$ in a channel, respectively; and $p,\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D708}$ are the pressure, fluid density and kinematic viscosity, respectively; $u_{i}$ and $f_{i}$ represent the fluid velocity and the reaction of particles with the fluid in the $x_{i}$ direction, respectively.

A DNS was performed for turbulent channel flow using the pseudo-spectral method. The dealiased Fourier expansion in the $x$ and $z$ directions and the Chebyshev expansion in the $y$ direction were used. For time advancing, the viscous and nonlinear terms were discretized using the Crank–Nicolson and third-order three-stage Runge–Kutta schemes, respectively.

A Reynolds number of $Re_{\unicode[STIX]{x1D70F}}=u_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D708}=180$ was considered in this study, where $u_{\unicode[STIX]{x1D70F}}$ and $\unicode[STIX]{x1D6FF}$ are the friction velocity and channel half-width, respectively. The same external pressure gradient $\text{d}P/\text{d}x$ was imposed on the flow using the wall shear stress of the particle-free flow $\unicode[STIX]{x1D70F}_{w}=\unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F}}^{2}$ as

(2.3) $$\begin{eqnarray}\frac{\text{d}P}{\text{d}x}=-\frac{\unicode[STIX]{x1D70F}_{w}}{\unicode[STIX]{x1D6FF}},\end{eqnarray}$$

in all the particle-free and particle-laden cases simulated. In the $(x,y,z)$ directions, the domain size was $(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF},2\unicode[STIX]{x1D6FF},2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF})$ and the mesh was $(128,129,128)$ . Then, the grid spacing was $\unicode[STIX]{x0394}x^{+}=17.7$ , $\unicode[STIX]{x0394}y^{+}=0.05$ –4.4 and $\unicode[STIX]{x0394}z^{+}=8.8$ in the streamwise, wall-normal and spanwise directions, respectively, in wall units. According to Kim, Moin & Moser (Reference Kim, Moin and Moser1987), this choice of resolution at $Re_{\unicode[STIX]{x1D70F}}=180$ is enough to predict a number of statistical correlations accurately. In the horizontal $(x,z)$ directions, periodic boundary conditions were employed.

2.2 Particles

For small (sub-Kolmogorov sized) spherical particles much heavier than the fluid in a horizontal turbulent channel flow, the particle equation of motion can be given in the Lagrangian frame as follows (Maxey & Riley Reference Maxey and Riley1983; Gerashchenko et al. Reference Gerashchenko, Sharp, Neuscamman and Warhaft2008):

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}v_{i}}{\text{d}t}=\frac{(1+0.15Re_{p}^{0.687})}{\unicode[STIX]{x1D70F}_{p}}(\tilde{u} _{i}-v_{i})-g\unicode[STIX]{x1D6FF}_{i2}, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}q_{i}}{\text{d}t}=v_{i}, & \displaystyle\end{eqnarray}$$

in which $v_{i}$ and $q_{i}$ are the particle velocity and the particle position in the $x_{i}$ direction, respectively, and $\tilde{u} _{i}$ indicates the fluid velocity at the particle position in the $x_{i}$ direction. The particle Reynolds number is defined as $Re_{p}=|\tilde{\boldsymbol{u}}-\boldsymbol{v}|d_{p}/\unicode[STIX]{x1D708}$ , where $d_{p}$ is the particle diameter and $\unicode[STIX]{x1D70F}_{p}$ is the particle response time; i.e. $\unicode[STIX]{x1D70F}_{p}=d_{p}^{2}\unicode[STIX]{x1D70C}_{p}/(18\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708})$ , where $\unicode[STIX]{x1D70C}_{p}$ is the particle density; $g$ is the magnitude of the gravitational acceleration and $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta.

Table 1. Parameters for comparison with the experiment by Gerashchenko et al. (Reference Gerashchenko, Sharp, Neuscamman and Warhaft2008) (cases A–C), for investigation of the two-way interaction between the settling particles and near-wall turbulence (cases D–H), and for comparison with the results without considering gravitational settling (cases I and J).

In this study, the two-way coupling effect at the $n$ th node position $\boldsymbol{x}^{n}$ in the $x_{i}$ direction was approximated by a point-force method, namely

(2.6) $$\begin{eqnarray}f_{i}(\boldsymbol{x}^{n})=-\mathop{\sum }_{k}w(\boldsymbol{x}^{n}-\boldsymbol{q}^{k})\frac{F_{i}^{k}}{\mathop{\vee }^{k}},\end{eqnarray}$$

where $F_{i}^{k}$ and $\vee ^{k}$ are the drag force of the $k$ th particle in the $x_{i}$ direction and the volume of the computational cell including the particle, respectively, and $(F_{i}^{k}/\vee ^{k})$ is exerted back on the fluid at the nearby eight grid points via the weight function $w$ based on linear interpolation from the $k$ th particle position $\boldsymbol{q}^{k}$ (Boivin, Simonin & Squires Reference Boivin, Simonin and Squires1998; Richter & Sullivan Reference Richter and Sullivan2013). Although more advanced methods are now available (see Maxey & Patel Reference Maxey and Patel2001; Capecelatro & Desjardins Reference Capecelatro and Desjardins2013; Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2015; Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2016), this classical point-force method has been extensively used for two-way coupled DNS of various flows such as homogeneous isotropic turbulence (Ferrante & Elghobashi Reference Ferrante and Elghobashi2003; Bosse et al. Reference Bosse, Kleiser and Meiburg2006; Abdelsamie & Lee Reference Abdelsamie and Lee2012), turbulent boundary layers (Li, Luo & Fan Reference Li, Luo and Fan2016a ; Li et al. Reference Li, Wei, Luo and Fan2016b ) and turbulent channel or Couette flows (Lee & Lee Reference Lee and Lee2015; Wang & Richter Reference Wang and Richter2019). However, the two-way force distribution $\boldsymbol{f}(\boldsymbol{x}^{n})$ obtained by this method is known to suffer from grid dependency as the number of particles per cell becomes $N_{p}/N_{c}<1$ , where $N_{p}$ is the number of total particles and $N_{c}$ is the number of computational cells (Eaton Reference Eaton2009; Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2013; Gualtieri, Battista & Casciola Reference Gualtieri, Battista and Casciola2017). As shown below in table 1, $N_{p}/N_{c}>1$ for the low Stokes number cases ( $St^{+}\leqslant 5.3$ ), whereas $N_{p}/N_{c}<1$ for the highest two Stokes numbers ( $St^{+}=10.6$ and 21.2) in our study. In order to ensure that our results for the highest Stokes number ( $St^{+}=21.1$ ) are reliable, we considered an extra case with particles of a very high density compared to the fluid ( $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}\approx 40\,000$ ) on the same resolution as in the DNS study by Nasr et al. (Reference Nasr, Ahmadi and McLaughlin2009) so that $N_{p}/N_{c}>1$ for the same $St^{+}=21.2$ . Upon comparison with case H ( $N_{p}/N_{c}<1$ , $St^{+}=21.2$ ), we found no qualitative difference between them (not shown here). Therefore, we believe that our results for the highest Stokes numbers are not significantly affected by numerical artefacts.

For time-efficient, high-accuracy calculation of $\tilde{u} _{i}$ , we employed the four-point Hermite interpolation scheme in the $x$ and $z$ directions and the fifth-order Lagrange polynomial interpolation scheme in the $y$ direction (Choi, Yeo & Lee Reference Choi, Yeo and Lee2004). Time advancement for (2.4) and (2.5) was performed using the third-order three-stage Runge–Kutta scheme, simultaneously with that for the flow field (Rouson & Eaton Reference Rouson and Eaton2001).

Figure 1. Particle vertical locations and streamwise velocities ensemble averaged over particles that have hit the bottom of the channel after injection at the upper wall as a function of time from the moment of the first particle–wall collision: (a) $St^{+}=0.81$ , (b) $St^{+}=2.65$ , (c) $St^{+}=5.3$ , (d) $St^{+}=10.6$ and (e) $St^{+}=21.2$ . The particle trajectory and velocity statistics were normalized by the particle radius $r_{p}(=d_{p}/2)$ and $r_{p}u_{\unicode[STIX]{x1D70F}}^{2}/\unicode[STIX]{x1D708}$ , respectively. For all the cases, $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}=833$ and $g=0.077u_{\unicode[STIX]{x1D70F}}^{3}/\unicode[STIX]{x1D708}$ .

2.3 Particle boundary condition

In order to determine the proper boundary condition for particles hitting the bottom wall, we examine the behaviour of heavy particles near the walls in horizontal turbulent channel flow in this subsection. We tracked the trajectories of 100 000 particles that were released simultaneously at random locations in the horizontal plane at the upper wall (i.e. a distance of a particle radius from the upper wall) with the interpolated horizontal fluid velocities and the vertical Stokes terminal velocity. Then, ensemble averages over the trajectories were performed without considering the feedback force of particles on the fluid. Figure 1 shows the average particle vertical locations and streamwise velocities for particles that have hit the bottom of the channel after injection at the upper wall as a function of time from the moment of the first particle–wall collision, for the several particle parameters which will be used in the two-way coupled simulations. Note that a perfectly elastic collision of particles with the walls was assumed in this study, as in most previous DNS studies for turbulent channel flow. In the figure, the conditionally ensemble-averaged trajectories eventually reach $y=r_{p}(=d_{p}/2)$ after a collision process between the particles and the channel bottom (although not shown here, the trajectories are never re-entrained into the flow even after several thousand wall time units), and the average particle streamwise velocities asymptotically approach the mean fluid velocity at $y=r_{p}$ , indicating that the particle forces on the fluid become negligible on average at some time after the occurrence of the first collision. Therefore, a particle that has ever collided with the lower wall is to be removed after its reaction with the fluid is no longer significant, i.e. after 6, 20, 40, 80, 160 wall time units following the first collision, respectively, for $St^{+}=0.81$ , 2.65, 5.3, 10.6, 21.2 (on the basis of the results of figure 1) in our two-way coupled simulations. Then, we release another particle at a distance of a particle radius from the upper wall (i.e. at $y=2\unicode[STIX]{x1D6FF}-r_{p}$ ), maintaining the total number of particles constant. The horizontal location of this new particle is chosen at random, and it begins settling with the interpolated horizontal fluid velocities and vertical Stokes terminal velocity.

In order to minimize the effect of this artificial initial condition of particle velocity, the feedback effect of the particles released at the upper wall on the fluid is taken into account only after they have reached on average a distance of $12\unicode[STIX]{x1D708}/u_{\unicode[STIX]{x1D70F}}$ away from the upper wall. By doing this, we can focus on the two-way interaction in regions where the coherent ejection motion of the fluid (away from the upper wall) becomes dominant, because its contribution to the Reynolds shear stress becomes significant at distances larger than 12 viscous lengths from the wall (Kim et al. Reference Kim, Moin and Moser1987). In figure 2, the horizontally averaged particle number density obtained when the ensemble-averaged trajectory of the particles released at the upper wall reaches $y=2\unicode[STIX]{x1D6FF}-12\unicode[STIX]{x1D708}/u_{\unicode[STIX]{x1D70F}}$ , normalized by the average number density over the entire channel, is shown along the vertical direction for each Stokes number. Here, the particle number density is defined as the number of particles per unit volume. In the figure, despite the same mean vertical location, the instantaneous concentration distribution of the particles depends strongly on the Stokes number. With decreasing Stokes number, the peak concentration near the upper wall is reduced, while the concentration distribution becomes more skewed toward the core region (away from the wall), consistent with Dorgan & Loth (Reference Dorgan and Loth2004) and Dorgan et al. (Reference Dorgan, Loth, Bocksell and Yeung2005).

Figure 2. Vertical profiles of the horizontally averaged particle number density when the ensemble-averaged trajectory of particles released at the upper wall reaches $y=2\unicode[STIX]{x1D6FF}-12\unicode[STIX]{x1D708}/u_{\unicode[STIX]{x1D70F}}$ , normalized by the average number density over the entire channel, for the same particle parameters as in figure 1.

In our two-way coupled DNS, a fully developed flow field for the $Re_{\unicode[STIX]{x1D70F}}=180$ channel simulation is used as the initial flow data. Then, the initial distribution for particles is a random uniform distribution over the entire channel domain, and the initial particle velocities are given by the interpolated fluid velocities at the particle locations. With the above-mentioned boundary conditions for the particles, the number of effective particles that do affect the fluid must be somewhat smaller than the total number of particles. After the initial uniform distribution, two-way coupling starts when the appearance of the effective particles near the upper wall and the particle removal at the lower wall are balanced at a statistically steady rate. Note that the injected particles at the upper wall always constituted less than 0.03 % of the total particles per time step for all the considered cases; the ratio was kept less than 0.002 %, 0.007 %, 0.013 %, 0.02 %, 0.03 % for $St^{+}=0.81$ , 2.65, 5.3, 10.6 and 21.2, respectively, throughout the entire simulation time. A time average for turbulence statistics is performed during 1500 wall time units. For particles moving outside the flow domain in the horizontal directions, periodic boundary conditions are applied, and inter-particle collisions are not included because a dilute suspension is assumed in our simulations.

Previously, a similar treatment was considered for particle sediments at the bottom. In their two-way coupled DNS study, Oresta & Prosperetti (Reference Oresta and Prosperetti2013) considered a liquid spray falling through a weakly turbulent Rayleigh–Bénard gas flow under gravity by removing and reinjecting the particles at the hot bottom and cold top plates, respectively, of a cylinder. Recently, Park, O’Keefe & Richter (Reference Park, O’Keefe and Richter2018) considered aeolian saltation by removing and reinserting non-isothermal particles at the lower boundary and in the lower 10 % of the domain, respectively, in a two-way coupled DNS study of Rayleigh–Bénard turbulence.

By using the point-particle approach without particle–particle collisions, which is adopted in this study, it is impossible to correctly simulate the presence of a sediment layer formed on the bottom wall over a long period of time (Soldati & Marchioli Reference Soldati and Marchioli2012). Instead, as described above, we eliminated deposited particles after several collisions with the lower wall and reinjected new particles into the flow near the upper wall to observe particle gravitational settling through the mean shear, which is the main subject of this paper. We believe that our investigation of this settling process will extend not only the Gerashchenko et al. (Reference Gerashchenko, Sharp, Neuscamman and Warhaft2008) experiment, but also the DNS studies by Dorgan & Loth (Reference Dorgan and Loth2004) and Dorgan et al. (Reference Dorgan, Loth, Bocksell and Yeung2005), in which small, heavy particles were released at the edge of the viscous sublayer near the upper wall in a horizontal boundary layer, to the two-way coupling regime. Furthermore, streamwise velocity streaks that are comparable with those of near-wall turbulence exist even in homogeneous shear turbulence (Lee, Kim & Moin Reference Lee, Kim and Moin1990) and in liquid shear layers formed by gas flows at a free surface (Rashidi & Banerjee Reference Rashidi and Banerjee1990). Since our focus is on the direct interaction between settling particles and those structures, the present study is also relevant to such turbulent shear flows in the absence of a wall.

Finally, it is worth mentioning that almost no momentum is added or eliminated by our particle treatment at both walls. When particles are released very near the top wall with the streamwise fluid velocity, the added momentum is negligible. Elimination of particles hitting the bottom wall was performed when the momentum of the particles was almost suppressed as explained above. Only when the particles which are released from the top wall reach on average 12 wall unit distances from the top wall, is a certain amount of momentum added to the fluid by the particles since the particle feedback force is turned on suddenly. However, an estimation of this momentum for the worst case of $St^{+}=21.1$ yields that it is at most 0.15 % of the mean pressure gradient force.

4.1 Modification of turbulence statistics

Figure 5. Mean streamwise fluid velocities for the particle-free case and cases D–H.

Figure 6. Mean streamwise fluid velocity profiles $\langle u_{1}^{\ast }\rangle$ normalized by $u_{\unicode[STIX]{x1D70F}}$ and $\unicode[STIX]{x1D708}$ of the respective particle-laden flow as a function of the distance (a) $2\unicode[STIX]{x1D6FF}^{\ast }-y^{\ast }$ and (b) $y^{\ast }$ from the upper and lower walls, respectively, in a semi-log plot for cases D–H. The thin solid line indicates the law of the wall, $\langle u_{1}^{\ast }\rangle =y^{\ast }$ and $\langle u_{1}^{\ast }\rangle =2.44\ln y^{\ast }+5.2$ .

Figure 5 shows the modifications of the mean streamwise fluid velocity $\langle u_{1}^{+}\rangle$ due to the presence of particles. With increasing Stokes number, the mean fluid velocities decrease in the upper part of the channel, while the maximum velocity is reduced and observed further below the geometric channel centreline ( $y/\unicode[STIX]{x1D6FF}=1$ ). Near the lower wall, the mean fluid velocities are little changed when $St^{+}=0.81$ and they are enhanced as the Stokes number increases further up to $St^{+}=10.6$ . The mean values between the two largest Stokes numbers ( $St^{+}=10.6$ and 21.2) are approximately the same in the region $0.2<y/\unicode[STIX]{x1D6FF}<0.6$ , but the enhancement at $y^{+}\approx 10$ is more pronounced for $St^{+}=5.3$ compared to that for $St^{+}=21.2$ (this is not visible in the figure, however, we will show this later on an increased scale in figure 18). The results of figure 5 indicate that the streamwise momentum of the fluid is carried in the direction of gravity by settling particles, as shown later.

Figure 6(a,b) displays the mean streamwise fluid velocity profiles $\langle u_{1}^{\ast }\rangle$ normalized by wall units of the respective particle-laden flow (represented by the superscript $\ast$ ) in the upper and lower parts of the channel, respectively, separated by the location of the peak $\langle u_{1}^{+}\rangle$ , as a function of the distance $2\unicode[STIX]{x1D6FF}^{\ast }-y^{\ast }$ and $y^{\ast }$ from the upper and lower walls, respectively. In the figure, the law of the wall with the von Kármán constant 0.41 and the additive constant 5.2 is also illustrated for comparison. Note that the mismatch of $\langle u_{1}^{\ast }\rangle$ at the location of the peak $\langle u_{1}^{+}\rangle$ indicates that $u_{\unicode[STIX]{x1D70F}}$ , used for normalization in figure 6(a,b), is differently modified at the upper and lower walls, respectively. In the upper part of the channel, the linear law still holds for the particle-laden cases, while as the Stokes number increases, $\langle u_{1}^{\ast }\rangle$ is slightly reduced in the buffer layer and lower part of the log layer, and its slope is enhanced around the core region compared to the log law. Note that, for $St^{+}=10.6$ and 21.2, the log law region where $\langle u_{1}^{\ast }\rangle \propto \ln y^{\ast }$ nearly disappears, while steeper slopes are instead observed in the profiles. On the other hand, $\langle u_{1}^{\ast }\rangle$ is significantly reduced with increasing Stokes number throughout the lower part of the channel.

Figure 7. Root-mean-square (a) streamwise, (b) wall-normal and (c) spanwise fluid velocity fluctuations, and (d) Reynolds shear stresses for the particle-free case and cases D–H.

In figure 7(a–d), the changes in the r.m.s. fluid velocity fluctuations $\langle (u_{i}^{\prime +})^{2}\rangle ^{1/2}$ (turbulence intensities) and the Reynolds shear stresses $\langle u_{1}^{+}u_{2}^{+}\rangle$ due to the presence of particles are illustrated. In the figure, a notable observation is that, with increasing Stokes number, all of the three components of turbulence intensities and the Reynolds shear stresses are suppressed near the lower wall while they are augmented in the upper part of the channel for $St^{+}\geqslant 2.65$ . Only the peak value of the streamwise fluctuation near the top wall slightly decreases with increasing Stokes number, which is probably not directly caused by the particles settling in that location because there are almost no effective particles that actually exert their feedback forces on the fluid there, as shown in figure 2. The results shown in figure 7 indicate that near-wall turbulence is altered in different ways depending on the sign of the mean shear rate $\text{d}\langle u_{1}\rangle /\text{d}y$ that the settling particles experience. In the next subsection, we will provide an explanation for this different modification behaviour. In particular, turbulence augmentation is observed even at $y/\unicode[STIX]{x1D6FF}=1$ , and hence the local minimum location in the turbulence intensity profiles, which is $y=\unicode[STIX]{x1D6FF}$ for the particle-free case, shifts toward the channel bottom. Accordingly, the Reynolds shear stress profiles cross zero below $y=\unicode[STIX]{x1D6FF}$ . The attenuation of turbulence near the lower wall is consistent with previous experimental observations (Tsuji & Morikawa Reference Tsuji and Morikawa1982; Kussin & Sommerfeld Reference Kussin and Sommerfeld2002). On the other hand, for the upper part of the channel, the turbulence augmentation is affected by the injection of particles at the upper wall. Interestingly, a similar trend was experimentally observed by Wu et al. (Reference Wu, Wang, Liu, Li, Zhang and Zheng2006), who measured the turbulence modification with polythene beads that settled in a horizontal turbulent channel flow, using a two-phase particle image velocimetry. The particle sizes were $60~\unicode[STIX]{x03BC}\text{m}$ and $110~\unicode[STIX]{x03BC}\text{m}$ , which were slightly smaller and larger than the Kolmogorov length scale of the flow, respectively, and the particle-to-fluid density ratio was $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}\approx 830$ . However, it is unclear whether the experimentally observed turbulence enhancement is physically associated with our numerical result.

Figure 8(a–c) shows the r.m.s. fluid vorticity fluctuations $\langle (\unicode[STIX]{x1D714}_{i}^{\prime +})^{2}\rangle ^{1/2}$ for the particle-free and particle-laden cases. In the outer layer (including the core region) in the upper part of the channel and near the lower wall, the overall trends for the modifications of the r.m.s. vorticity fluctuations are similar to those for the turbulence intensities and Reynolds shear stresses shown in figure 7. As will be shown in the next subsection by investigating instantaneous flow fields, the vorticity fluctuation enhancement and attenuation are consistent with an increased occurrence and suppression of small-scale vortical structures, respectively. A more intricate behaviour is observed very close to the upper wall, where no suppression is observed for the streamwise component of the r.m.s. vorticity fluctuations, while their wall-normal (including the local peak value) and spanwise components are reduced with increasing Stokes number. In the inset of figure 8(c), the spanwise r.m.s. vorticity fluctuations increase with increasing Stokes number up to $St^{+}=5.3$ , for which the Stokes number $St_{K}$ is approximately unity on the basis of the Kolmogorov time scale near the wall, and then decrease with further increasing $St^{+}$ in the viscous sublayer near the lower wall ( $y^{+}<5$ ). For the largest Stokes number, the spanwise vorticity fluctuations are reduced compared to the particle-free case. This spanwise r.m.s. vorticity fluctuation enhancement near the lower wall is due to enhanced shear fluctuations, because $\unicode[STIX]{x1D714}_{3}^{\prime }=-\text{d}u_{1}^{\prime }/\text{d}y$ at the wall. The $\langle (\text{d}u_{1}^{\prime }/\text{d}y)^{2}\rangle$ enhancement close to the lower wall is due to the fact that the settling particles carry the higher streamwise fluid momentum of the outer layer to the near-wall region, mainly via ‘sweeps’ by the near-wall streamwise vortices (Marchioli & Soldati Reference Marchioli and Soldati2002; Righetti & Romano Reference Righetti and Romano2004). This momentum transfer may be maximized when the particle time scale is comparable to the characteristic fluid time scale of the vortical structures. Among the Stokes numbers considered in this study, the particle response time for $St^{+}=5.3$ best matches the vortex time scale, defined as the inverse of the r.m.s. streamwise fluid vorticity $1/\langle (\unicode[STIX]{x1D714}_{x}^{\prime +})^{2}\rangle ^{1/2}$ , in the buffer layer, where near-wall streamwise vortices mostly occur (Soldati & Marchioli Reference Soldati and Marchioli2009; Lee & Lee Reference Lee and Lee2015). Therefore, the increase in $\langle (\text{d}u_{1}^{\prime }/\text{d}y)^{2}\rangle$ is most pronounced at $St^{+}=5.3$ , as shown in the inset of figure 8(c).

Figure 8. Root-mean-square (a) streamwise, (b) wall-normal and (c) spanwise fluid vorticity fluctuations for the particle-free case and cases D–H. In (c), the inset indicates $\langle (\unicode[STIX]{x1D714}_{3}^{\prime +})^{2}\rangle ^{1/2}$ near the lower wall on an increased scale.

4.2 Particle behaviour and turbulence modification mechanism

In near-wall turbulence, the spatial distribution of inertial particles under the influence of hydrodynamic drag forces exhibits two main patterns in terms of their local accumulation, one due to preferential concentration out of turbulent vortical regions and the other due to turbophoresis (preferential accumulation in lifted low-speed streaks) (Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012; Lee & Lee Reference Lee and Lee2015; Richter et al. Reference Richter, Garcia and Astephen2016). In order to gain further detailed information on the role of wall-normal gravity in the turbulence modification behaviour shown in § 4.1, in this subsection we investigate the gravitational effect on these local accumulation patterns and the consequent modification of coherent turbulent structures for two different cases, $St^{+}=5.3$ and $St^{+}=21.2$ , at which preferential concentration and turbophoresis are distinctly present, respectively (Lee & Lee Reference Lee and Lee2015).

Figure 9. Instantaneous distributions of particles near a streamwise fluid vorticity field in the $y^{+}$ – $z^{+}$ plane (left) and profiles of the horizontally averaged particle number density normalized by its initial value (right) for (a) $St^{+}=5.3$ (case F) and (b)  $St^{+}=21.2$ (case H). Red and blue colours represent, respectively, positive and negative values for the intense streamwise fluid vorticity, larger in magnitude than the maximum r.m.s. value. Black dots denote the locations of particles collected over a streamwise distance of approximately 35 viscous lengths.

Figure 10. Instantaneous velocity vectors of rapidly settling particles near a streamwise fluid velocity field in the $y^{+}$ – $z^{+}$ plane for (a) $St^{+}=5.3$ (case F) and (b)  $St^{+}=21.2$ (case H). Only particles that settle faster than (a) the maximum r.m.s. vertical fluid velocity fluctuation $\langle (u_{2}^{\prime })^{2}\rangle ^{1/2}\approx 1$ and (b) the Stokes terminal velocity $\unicode[STIX]{x1D70F}_{p}^{+}g^{+}=1.6324$ are shown. The results are for the particle distributions and flow fields shown in figure 9.

Figure 9(a,b) displays the instantaneous distributions of settling particles in comparison with a streamwise fluid vorticity field in the $y^{+}$ – $z^{+}$ plane and the profiles of the horizontally averaged particle number density, $n(y,t)$ , normalized by its initial value, $n_{0}(y)$ , for $St^{+}=5.3$ and 21.2, respectively. For the same particle distributions and flow fields as in figure 9(a,b), the velocity vectors for rapidly settling particles are shown, along with the streamwise fluid velocity contours, in figure 10(a,b), respectively. Note that in figure 10, only particles that settle faster than the maximum r.m.s. vertical fluid velocity fluctuation $\langle (u_{2}^{\prime })^{2}\rangle ^{1/2}\approx 1$ and the Stokes terminal velocity $\unicode[STIX]{x1D70F}_{p}^{+}g^{+}=1.6324$ were selected to distinguish the rapidly settling particles for $St^{+}=5.3$ and 21.2, respectively. In figure 9(a,b), the average particle number density remains almost constant at $n\approx n_{0}$ for $St^{+}=5.3$ and at $n\approx 0.6n_{0}$ for $St^{+}=21.2$ throughout the channel width, except in the viscous sublayer near the upper and lower walls. This indicates that a volume fraction of at least $\unicode[STIX]{x1D719}_{v}\approx 3.6\times 10^{-5}$ is maintained over the bulk of the flow by the particle injection at the upper wall. The extremely high peak of the average particle number density is observed at the lower wall (i.e. at $y=r_{p}$ ), in particular for $St^{+}=21.2$ , but this may not significantly affect the fluid because the influence of the no-slip condition is still dominant there. In figure 9(a), preferential concentration of particles with $St^{+}=5.3$ in regions of low streamwise fluid vorticity fluctuations is clearly observed. This is because the Stokes number $St^{+}=5.3$ matches the vortex time scale related to the intense streamwise vorticity fluctuations (Richter & Sullivan Reference Richter and Sullivan2013; Lee & Lee Reference Lee and Lee2015). Furthermore, it corresponds to a Stokes number of $St_{K}\approx 1$ in the buffer layer. In homogeneous isotropic turbulence under gravity, preferential concentration of heavy particles is responsible for their preferential sweeping in regions of descending fluid associated with local vortices and thus an increase in their average settling velocity over a similar Stokes number range ( $St_{K}\approx 1$ ) (Wang & Maxey Reference Wang and Maxey1993; Yang & Shy Reference Yang and Shy2005). This preferential sweeping effect also leads to a large downward particle force on the fluid in one side of a local vortical structure, and as a result, the vortical structures are stretched in the direction of gravity (Ferrante & Elghobashi Reference Ferrante and Elghobashi2003; Yang & Shy Reference Yang and Shy2005; Bosse et al. Reference Bosse, Kleiser and Meiburg2006). In near-wall turbulence, ‘ejections’ of low-speed fluid away from the upper wall and ‘sweeps’ of high-speed fluid toward the lower wall associated with streamwise vortical structures serve as such descending fluid for preferential sweeping. Therefore, particles rapidly settle in regions of ejections and sweeps in the upper and lower parts of the channel, respectively, as clearly seen in figure 10(a). These rapidly settling particles may enhance the surrounding fluid motions by exerting a large downward force on the fluid. In particular, lift-up of low-speed fluid (away from the wall) by ejections associated with streamwise vortices is responsible for a lifted low-speed streak that has been regarded as a fundamental element for the creation of small-scale vortical structures, such as hairpin vortices, in near-wall turbulence (Robinson Reference Robinson1991; Panton Reference Panton2001). As will be confirmed later, this indicates that more small-scale vortical structures are produced far away from the upper wall as low-speed fluid is pulled through preferential sweeping in ejections (see figures 13 and 14 a,b).

Figure 11. Instantaneous distributions of particles near coherent structures in the lower part of the channel in the $y^{+}$ – $z^{+}$ plane for $St^{+}=5.3$ (case F). Solid contour lines of streamwise fluid velocities with a contour level of $u_{1}^{+}=(3,5,7,9,11)$ and white arrow vectors of the vertical and spanwise components of the fluid velocity $(u_{2}^{+},u_{3}^{+})$ are shown. In (a), black arrow vectors and colours on dots represent the vertical and spanwise components $(v_{2}^{+},v_{3}^{+})$ and the streamwise component of the particle velocity $v_{1}^{+}$ , respectively. In (b), black arrow vectors and colours on dots indicate the vertical and spanwise components $(-(F_{2}/m_{p})^{+},-(F_{3}/m_{p})^{+})$ and the streamwise component of the feedback force of particles on the fluid per the particle mass $-(F_{1}/m_{p})^{+}$ , respectively.

Figure 12. Instantaneous distributions of particles near coherent structures in the lower part of the channel in the $y^{+}$ – $z^{+}$ plane for $St^{+}=21.2$ (case H). Solid contour lines, white and black arrow vectors, and colours on dots in (a,b) have the same meaning as in figure 11(a,b), respectively.

Figure 13. Instantaneous distributions of isosurfaces of (a) $\unicode[STIX]{x1D706}_{2}^{+}=-0.012$ and (b) $u_{1}^{+}=10$ in the lower half of the channel for the particle-free case. Colour contours of the vertical distance from the lower wall $y^{+}$ are plotted on the $\unicode[STIX]{x1D706}_{2}^{+}=-0.012$ and $u_{1}^{+}=10$ isosurfaces. The results are shown throughout the full domain in the horizontal directions.

Figure 14. Instantaneous distributions of isosurfaces of (a,c)  $\unicode[STIX]{x1D706}_{2}^{+}=-0.012$ and (b,d)  $u_{1}^{+}=10$ in the upper part of the channel based on the peak of $\langle u_{1}^{+}\rangle$ (upside-down view) for (a,b)  $St^{+}=5.3$ (case F) and for (c,d) $St^{+}=21.2$ (case H). Colour contours of the vertical distance from the upper wall $2\unicode[STIX]{x1D6FF}^{+}-y^{+}$ are plotted on the $\unicode[STIX]{x1D706}_{2}^{+}=-0.012$ and $u_{1}^{+}=10$ isosurfaces. The results are shown throughout the full domain in the horizontal directions.

Figure 15. Instantaneous distributions of isosurfaces of (a,c) $\unicode[STIX]{x1D706}_{2}^{+}=-0.012$ and (b,d)  $u_{1}^{+}=10$ in the lower part of the channel (based on the peak of $\langle u_{1}^{+}\rangle$ ) for (a,b) $St^{+}=5.3$ (case F) and for (c,d)  $St^{+}=21.2$ (case H). Colour contours of the vertical distance from the lower wall $y^{+}$ are plotted on the $\unicode[STIX]{x1D706}_{2}^{+}=-0.012$ and $u_{1}^{+}=10$ isosurfaces. The results are shown throughout the full domain in the horizontal directions.

Several previous investigators have shown that the local accumulation of heavy particles due to turbophoresis, i.e. the preferential accumulation in ejections associated with near-wall streamwise vortices (see figure 2 of Soldati & Marchioli Reference Soldati and Marchioli2009), becomes efficient for similar Stokes numbers to our $St^{+}=21.2$ case, assuming that the effect of gravity is negligible or is directed in the mean flow direction (Marchioli et al. Reference Marchioli, Giusti, Salvetti and Soldati2003, Reference Marchioli, Picciotto and Soldati2007; Lee & Lee Reference Lee and Lee2015). However, when the gravitational settling occurs at this Stokes number, it appears in figure 9(b) that this accumulation pattern disappears near the upper wall due to a large Stokes terminal fall velocity of $\unicode[STIX]{x1D70F}_{p}^{+}g^{+}=1.6324$ , which is even greater than the maximum $\langle (u_{i}^{\prime +})^{2}\rangle ^{1/2}$ . Furthermore, in figure 10(b), the rapidly settling particles exhibit very little spanwise motion until their arrival at the bottom, indicating the significance of the effect of crossing trajectories, which accounts for a decrease in the turbulent dispersion of heavy particles with a large terminal fall velocity in comparison with that of fluid or light particles (Yudine Reference Yudine1959; Csanady Reference Csanady1963; Reeks Reference Reeks1977; Wells & Stock Reference Wells and Stock1983; Elghobashi & Truesdell Reference Elghobashi and Truesdell1992). Dorgan & Loth (Reference Dorgan and Loth2004) also observed the effect of crossing trajectories for heavy particles settling in a horizontal turbulent boundary layer. A distinctive feature of flow modification due to the particles is that near-wall low-speed fluid is very efficiently transported away from the upper wall with rapidly settling particles, across almost the channel half-width $\unicode[STIX]{x1D6FF}$ in the vertical direction (in particular, at $z^{+}\approx 600$ ), and approximately $2\unicode[STIX]{x1D6FF}$ apart intermittently in the spanwise direction, as observed in figure 10(b). Furthermore, a comparison between figures 9(b) and 10(b) shows that streamwise vortex structures appear along the regions of low-speed fluid, vertically extended by rapidly settling particles in the upper part of the channel. In fact, as will be shown later, the number of small-scale vortical structures increases at quite large distances from the upper wall, as in the case of $St^{+}=5.3$ (see figure 14 c).

Figures 11 and 12 show the instantaneous particle velocities and feedback accelerations along with coherent near-wall streamwise vortices and a lifted low-speed streak near the lower wall in the $y^{+}$ – $z^{+}$ plane for $St^{+}=5.3$ and $St^{+}=21.2$ , respectively. It appears from figure 11(a) that when settling particles reach the wall region (the buffer layer and lower part of the log layer) where near-wall streamwise vortices are commonly observed, preferential sweeping occurs in coherent sweeps, consistent with the observation near the lower wall shown in figure 10(a). During the fast settling via sweeps, the particles maintain the high streamwise momentum of the outer layer (figure 11 a). Therefore, they experience a deceleration in the streamwise direction as the lower wall is approached, thereby exerting a positive streamwise feedback force on the fluid, as indicated by the red dots close to the lower wall in figure 11(b). These high-speed particles (i.e. swept to the wall by the preferential sweeping effect) contribute to the strong mean streamwise deceleration and negatively skewed streamwise acceleration distribution of the droplets near the lower wall found in the Gerashchenko et al.’s (Reference Gerashchenko, Sharp, Neuscamman and Warhaft2008) experiment and in Lavezzo et al.’s (Reference Lavezzo, Soldati, Gerashchenko, Warhaft and Collins2010) DNS and also shown in the present study (figures 3 and 4). Since the settling particles in sweeps can acquire enough vertical momentum to reach the wall without being affected by the ejection motion associated with the streamwise vortex from both gravity and the sweep motion of the fluid, it seems from figure 12 that only a few particles overcome the effect of gravity and move to the ejection region.

In figure 12, the effect of preferential concentration due to small-scale vortical structures is almost imperceptible, and the settling particles simply cross the regions of not only streamwise vortices but also a lifted low-speed streak with very little spanwise motion. Then, the feedback forces act to suppress the fluid motions associated with the coherent structures. In particular, the particles just passing vertically through a lifted low-speed streak exert large streamwise and vertical feedback forces against the turbulent velocities, as indicated by the red dots in the ejection region in figure 12(b). This eventually affects the formation of near-wall streamwise vortices, as will be shown later (see figure 15).

In figure 13, two typical coherent structures responsible for production of near-wall turbulence, lifted low-speed streamwise streaks and near-wall streamwise vortices, are visualized for an instantaneous flow field of the particle-free case, through an isosurface of the streamwise fluid velocity and the $\unicode[STIX]{x1D706}_{2}$ method, respectively (Jeong & Hussain Reference Jeong and Hussain1995; Jeong et al. Reference Jeong, Hussain, Schoppa and Kim1997). Figures 14 and 15 exhibit the modifications of these structures in the upper and lower parts of the channel, respectively, for two particle-laden cases of $St^{+}=5.3$ and 21.2. When gravitational settling occurs in the upper part of the channel, more small-scale vortical structures are generated at the log layer and even in the vicinity of the geometric channel centreline (at $y/\unicode[STIX]{x1D6FF}\approx 1$ ) (figure 14 a,c), as low-speed streaks are lifted farther in the direction of gravity (figure 14 b,d), particularly, for the largest Stokes number case. On the other hand, when particles settle in the lower part of the channel, the opposite is true; lifted low-speed streaks and near-wall vortical structures are damped compared to the particle-free case (figure 15). This trend is much stronger, again, for the largest Stokes number case (figure 15 c,d). The modification behaviour of coherent turbulent structures shown here is consistent with that of the turbulence statistics shown in figures 7 and 8.