1. Introduction
The turbulent channel flow laden with small, gravity-free inertial spheres is a paradigmatic multiphase flow, yet not fully understood and correctly modelled. This system has been widely studied using direct numerical simulations (DNS) of the Navier–Stokes equations using the point-particle approximation (Balachandar & Eaton Reference Balachandar and Eaton2010), which assumes that interphase coupling is localized at a single point and the particle dynamics driven by an undisturbed flow field, evaluated at the particle position. This approximation allows us to study the flow dynamics beyond the particle scale, without solving the actual flow around each particle. Despite the many studies, the validity of these particle–fluid coupling models has not been fully assessed for canonical turbulent wall flows, in particular when the underlying turbulence is altered by the dispersed phase.
When the particle feedback on the flow becomes important – the two-way coupling regime – care should be taken in the estimation of the undisturbed velocity. The challenge is conciliating the estimation of an undisturbed velocity sampled by the particle with the need of forcing the local velocity field. Indeed, approaches for accurate and realistic two-way coupling methods have been pursued since the work of Crowe, Sharma & Stock (Reference Crowe, Sharma and Stock1977), and are still object of active research (see, e.g. Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2015; Horwitz & Mani Reference Horwitz and Mani2016; Ireland & Desjardins Reference Ireland and Desjardins2017; Battista et al. Reference Battista, Mollicone, Gualtieri, Messina and Casciola2019; Pakseresht, Esmaily & Apte Reference Pakseresht, Esmaily and Apte2020). Moreover, results in the literature for integral observables show qualitatively different trends. For instance, while Vreman (Reference Vreman2007) and Zhao, Andersson & Gillissen (Reference Zhao, Andersson and Gillissen2010) have observed a drag-reducing behaviour in two-way coupled particle-laden turbulent channel flow, other studies have measured a drag increase (Battista et al. Reference Battista, Mollicone, Gualtieri, Messina and Casciola2019). These opposite trends suggest the presence of different regimes of momentum transfer in wall-bounded particle-laden turbulence. Indeed, Capecelatro, Desjardins & Fox (Reference Capecelatro, Desjardins and Fox2018) used Euler–Lagrange volume-filtered DNS to investigate in detail the turbulent kinetic energy (TKE) budget in vertical particle-laden channel flow at different volume fractions. Different regimes were observed with increasing mass fraction – first, a turbulent regime at low mass fractions; second, increasing the mass fraction, the turbulence activity decreases until the flow laminarizes, due to the growing importance of a ‘drag dissipation-and-exchange term’ associated with fluid–particle correlated velocity fluctuations; finally, increasing the mass fraction beyond this limit triggers a second mechanism for TKE increase, due to the average velocity difference between phases, which grows important and re-energizes fluid turbulence. Despite the exciting developments exploiting particle-modelled two-way coupling interactions, there is a clear need for high-fidelity data to support the body of work using these modelled simulations, and to possibly reconcile seemingly conflicting observations of drag increase and drag reduction in two-way coupling point-particle numerical simulations.
Indeed, the tremendous developments of approaches for DNS with modelled particle–fluid coupling have not been accompanied by the high-fidelity data which are essential for validating the underlying assumptions. These data can be obtained either from interface-resolved simulations of the Navier–Stokes equations, or well-controlled experiments matching the computational set-up. Nowadays, the first direct comparisons between point-particle DNS data and experiments or interface-resolved DNS are possible. Examples are the work by Wang et al. (Reference Wang, Fong, Coletti, Capecelatro and Richter2019), which compares one-to-one two-way coupled point-particle DNS with experimental measurements for vertical turbulent channel flow at different volume loads, and the interface-resolved (also denoted particle-resolved) DNS of spherical particles near the point-particle limit in homogeneous isotropic turbulence (Schneiders, Meinke & Schröder Reference Schneiders, Meinke and Schröder2017; Fröhlich et al. Reference Fröhlich, Schneiders, Meinke and Schröder2018; Mehrabadi et al. Reference Mehrabadi, Horwitz, Subramaniam and Mani2018) and in turbulent channel flow (Horne & Mahesh Reference Horne and Mahesh2019; Costa, Brandt & Picano Reference Costa, Brandt and Picano2020).
The present work also exploits interface-resolved simulations and focuses on near-wall turbulence modulation by small inertial particles. We adopt the configuration in Costa et al. (Reference Costa, Brandt and Picano2020) and add results at higher volume fraction. Three bulk mass fractions 
$\varPsi =0.34\,\%$, 
$3.37\,\%$ and 
$33.7\,\%$ at fixed particle to fluid density ratio (
$\varPi _\rho =100$) are hence considered, where the latter two show non-negligible two-way coupling effects. Our results reveal two distinct mechanisms for turbulence modulation. At lower volume fractions, the flow is only dense very close to the wall, and the higher flow inertia in this region, with particles travelling much faster than the fluid, results in a drag-increased flow resembling single-phase turbulence at slightly higher Reynolds number. At higher volume fractions, the dispersed phase is dynamically important over the entire channel. Here, the drag increasing effect is amplified by correlated particle velocity fluctuations, but counterbalanced by a substantial fluid turbulence drag reduction. This results in a milder increase in drag for the denser case. These observations may be explained in light of the streamwise momentum balance for vanishing volume fraction, but with non-negligible mass fraction.
2. Methods and computational set-up
Since we use here the tools and set-up in Costa et al. (Reference Costa, Brandt and Picano2020), with one additional case at the largest volume fraction, we briefly summarize the numerical method and refer to the previous work for more details. We solve the continuity and Navier–Stokes equations for an incompressible Newtonian fluid, together with the Newton–Euler equations driving the motion of the solid spherical particles. These two sets of equations are coupled directly using the immersed-boundary method developed by Breugem (Reference Breugem2012) (see also Uhlmann Reference Uhlmann2005), built on a standard second-order finite-difference method on a three-dimensional, staggered Cartesian grid, using a fast-Fourier-transform-based pressure-projection method (Kim & Moin Reference Kim and Moin1985; Costa Reference Costa2018). Short-range particle–particle/particle–wall interactions (lubrication and solid–solid contact) are modelled using the method of Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015), as in Costa et al. (Reference Costa, Brandt and Picano2020). More specifically, as regards solid–solid contact, the particles are frictionless, and with a normal solid–solid coefficient of restitution of 
$0.97$.
Turbulent channel flow is simulated in a domain periodic in the streamwise (
$x$) and spanwise (
$z$) directions, with no-slip and no-penetration boundary conditions imposed at the walls (
$y=h\mp h$), where 
$h$ is the channel half-height. The flow is driven by a uniform pressure gradient that ensures a constant bulk velocity. The bulk Reynolds number is equal to 
${Re}_b = U_b(2h)/\nu =5600$, which corresponds to an unladen friction Reynolds number 
${Re}_\tau ^{sph} = u_\tau h/\nu \approx 180$; where 
$U_b$ is the bulk flow velocity and 
$u_\tau$ the wall friction velocity. The particle properties are chosen to yield a particle Reynolds number 
${Re}_p = D u_\tau /\nu =D^{+}= 3$, and Stokes number 
$St_p = \varPi _\rho {Re}_p^{2}/18=50$, where 
$\nu$ is the fluid kinematic viscosity, 
$D$ is the particle diameter and 
$\varPi _\rho$ the particle-to-fluid mass density ratio; since the particle size is restricted by resolution requirements, 
$\varPi _\rho$ is used to achieve the target particle Stokes number, which was demonstrated to feature highly inhomogeneous particle distributions in wall turbulence (see, e.g. Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012). Three values of solid volume fraction are considered, increasing by factors of 
$10$: 
$\varPhi \simeq 3\times 10^{-5}$, denoted very dilute (VD), 
$\varPhi \simeq 3\times 10^{-4}$, dilute (D), and 
$\varPhi \simeq 3\times 10^{-3}$, semi-dilute (SD), with the data pertaining to the first two cases also used in the recent study by Costa et al. (Reference Costa, Brandt and Picano2020). Table 1 shows all relevant physical and computational parameters.
Table 1. Computational parameters; 
$\varPhi /\varPsi$ denote the bulk solid volume/mass fraction, and 
$N_p$ the total number of particles. Common to all cases: bulk Reynolds number 
${Re}_b = 5600$ (i.e. friction Reynolds number in the single-phase limit 
${Re}_\tau ^{sph}\approx 180$); particle size ratio 
$D/(2h)=1/120$; particle-to-fluid mass density ratio 
$\varPi _\rho =100$. These correspond to a particle Reynolds number based on the unladen reference values of 
$D^{+}=3$ and Stokes number 
$St=50$. The fluid domain is discretized on a regular Cartesian grid with 
$(L_x/N_x)\times (L_y/N_y)\times (L_z/N_z) = (6h/4320)\times (2h/1440)\times (3h/2160)$, while the particles are resolved with 
$D/\Delta x=12$ grid points over the particle diameter (in total 
$420$ Lagrangian grid points) (same as Costa et al. Reference Costa, Brandt and Picano2020).
Finally, unless otherwise stated, the mesoscale-averaged profiles reported in this manuscript correspond to intrinsic averages (i.e. averaged only over the corresponding phase) in time and along the two homogeneous directions for each phase (Costa et al. Reference Costa, Brandt and Picano2020), obtained with wall-parallel bins with thickness equal to the grid spacing.
3. Results
A three-dimensional visualization of the flow pertaining to the different cases is reported in figure 1, where iso-surfaces of the positive second invariant of the velocity gradient tensor Q (coloured by the local wall-normal velocity) and the particles are displayed. High-vorticity spots, footprints of the presence of the particles, are noticeable in all particle-laden cases, especially closer to the wall, where, as we will see, the mean slip velocity between the phases is highest. As expected, no significant qualitative differences can be seen between the single-phase and the very dilute case VD, apart from the very small number of dispersed particles in the latter case. Qualitative differences between these cases and case D are also small; however, an increased number of high-vorticity spots due to the particles is obvious over the entire domain. Finally, a strong modulation of the flow dynamics by the particles is seen when inspecting case SD, where the disruption of the flow coherent structures is evident. One of the main messages of the present work is that, although cases D and SD have macroscopic flow variations and, therefore, are formally in a two-way coupling regime, the mechanism of turbulence modulation is definitively different.
Figure 1. Three-dimensional flow visualizations of the different cases under consideration. We display surfaces of constant second invariant of the velocity gradient tensor 
$Q=20(U_b/h)^{2}$, coloured by the local wall-normal velocity (white – high and blue – low). The particles are shown to scale, in orange colour. See table 1 for a description of the different cases.
Figure 2(a) presents the friction Reynolds number as a function of the bulk volume fraction, while the change in drag relative to the unladen flow is quantified in panel (b). The friction Reynolds number increases significantly (
$5\,\%$) from case VD to case D, despite the relatively low mass fraction 
$\varPsi = \mathcal {O}(10^{-2})$. As discussed in Costa et al. (Reference Costa, Brandt and Picano2020), this increase is attributed to the higher inertia induced by the large local mass fraction near the wall; the particles are driven by turbophoresis towards the wall, where they experience a large apparent slip velocity. Further increasing the mass fraction by an order of magnitude (case SD) results in a mild increase in drag, only approximately 
$3\,\%$ with respect to case D. This milder increase indicates that a competing drag-reducing mechanism comes into play at higher mass loading. The blue symbols in figure 2 display quantities computed using a wall friction velocity determined from the centreline slope of the Reynolds stress, so as to quantify the decrease in turbulent momentum transfer, as discussed in detail later.
Figure 2. (a) Friction Reynolds number 
${Re}_\tau$ vs the bulk volume fraction 
$\varPsi$; 
${Re}_\tau ^{t}$ is based on a velocity scale computed from the slope of the outer-scaled Reynolds shear stresses evaluated at the channel centreline 
$u_\tau ^{t} = \sqrt {\partial _{y/h} \langle u^{\prime } v^{\prime } \rangle |_{y/h=1}}$, as a rough measure of the contribution of fluid velocity correlations to the mean wall shear (see Picano, Breugem & Brandt Reference Picano, Breugem and Brandt2015). (b) Drag modulation 
$\mathrm {DC}={Re}_\tau ^{2}/{Re}_\tau ^{2,sph}$, computed with 
${Re}_\tau$ and 
${Re}_\tau ^{t}$ (
$\mathrm {DC}^{t}$).
Figure 3(a) shows the mass fraction normalized by the corresponding bulk value vs the wall-normal distance in particle diameters. All the cases show a peak at the wall, corresponding to a particle layer. Note also that these peaks occur slightly beyond a wall distance of one particle radius, due to the bouncing dynamics caused by particle–wall collisions. Clearly, the fraction of particles near the wall decreases with increasing mass loading. This decrease in wall accumulation is attributed to the increased shear rate near the boundary, which enhances lift forces (Costa et al. Reference Costa, Brandt and Picano2020). However, case SD shows a disproportionally strong decrease of the relative particle concentration, considering the relatively low increase in wall shear with respect to case D. This is partly caused by two-way coupling effects, which dampen the intensity of the wall-normal fluid velocity fluctuations responsible – to first approximation – for the turbophoretic drift (Marchioli & Soldati Reference Marchioli and Soldati2002). In addition, particle–particle interactions may become significant. The inset of figure 3(a) shows the mass fraction profile as a function of the outer-scaled wall-normal distance. While the mass fraction profiles become uniform far from the wall for cases VD and D, the distribution shows a mild monotonic increase with the wall distance in case SD. These observations indicate that, for this densest case, non-negligible particle–particle interactions may be driving particles towards regions of low shear. The mechanism underlying this drift is most likely similar to that reported in Fornari et al. (Reference Fornari, Formenti, Picano and Brandt2016) – particle inertia and local high shear promote inter-particle interactions, causing a net migration towards low shear regions. We attribute the non-zero slope at the channel centreline of the profile pertaining to case SD to a lower statistical convergence in the channel core in combination with the amplification induced by the logarithmic scale, as explained in the figure caption.
Figure 3. (a) Local solid mass fraction 
$\psi$, normalized by the bulk value 
$\varPsi$ as a function of the wall-normal distance in particle diameters (note that the profile is clamped beyond 
$y/D=6$, to highlight the near-wall peaks); the inset shows the nominal (i.e. not normalized) local mass fraction 
$\psi$ as a function of the outer-scaled wall-normal distance. The non-zero slope at the centreline, for case SD, is attributed to insufficient statistical convergence in this region where the flow is more dilute and the volume fraction fluctuations higher, further highlighted by the logarithmic scaling. (b) Outer-scaled spanwise autocorrelation of streamwise velocity at 
$y/h=0.067$ and 
$y/h=0.167$, (
$y/\delta _{\nu }^{sph}\approx 12$ and 
$30$); ‘sph’ corresponds to the unladen flow.
The inset of figure 3(a) also illustrates where two-way coupling effects are expected to be important. These can be envisaged near the wall for the less dilute cases, as the solid mass fraction increases to 
$20\,\%$ for case D, and to approximately 
$80\,\%$ for case SD. Away from the wall, instead, the solid mass fraction retains high values only for case SD, 25 %–40 %, and therefore turbulence modulation is also expected in the bulk, which we will relate to the milder overall drag increase from case D to case SD. In the latter case, we observe a significant turbulence modulation by the solid particles: the spacing between low- and high-speed streaks is larger, the spanwise velocity variation smoother, and the maximum streak amplitude occurs at a much larger wall-normal distance, which is reflected in the spanwise autocorrelation of the streamwise velocity in wall-parallel planes shown in figure 3(b). This trend can be clearly seen in figure 4, showing the contours of the near-wall streamwise fluid velocity. Case SD also presents non-negligible four-way coupling effects. In particular, we estimated the particle collision frequency near the wall following Sundaram & Collins (Reference Sundaram and Collins1997)
\begin{equation} N_c = \tfrac{1}{2}{\rm \pi} R^{2}n^{2} g(R)\left\langle \Delta U^{n-}\right\rangle , \end{equation}
where 
$R$ is the particle radius, 
$n$ the local number density, 
$g(R)$ the particle radial distribution function at the imminence of contact and 
$\langle \Delta U^{n-}\rangle$ the mean particle relative velocity projected along the line-of-centres, conditioned to negative values (i.e. promoting a collision) at the imminence of contact (i.e. interparticle distance equal to one particle diameter). We found the collision frequency for case SD to be 
$5.5$ near the wall and to reduce to 
$0.53$ in the bulk, a quantity defined per unit volume 
$L_x L_z D$ and bulk eddy turnover time 
$(O(h/u_\tau ))$, more than one order of magnitude larger than in case D (
$0.29$ near the wall and 
$0.002$ in the bulk). We should note that, owing to the computational cost and the low volume fractions considered, we were unable to collect enough data to directly measure particle collision statistics. Also, a direct assessment of collision frequency can yield deceiving results due to the possibility of sustained contacts; see Kuerten & Vreman (Reference Kuerten and Vreman2016).
Figure 4. Contours of outer-scaled streamwise fluid velocity 
$u$ at a plane 
$y/h = 0.07$, together with the particles located at 
$y/h<0.15$. Perspective of an observer aligned with the negative wall-normal direction, with flow from left to right (the contour is shown with slight transparency).
The inner-scaled mean velocity profiles, see figure 5(a), resemble those of the single-phase flow, with a downward shift that indicates drag increase (the difference is less apparent in the corresponding outer-scaled profiles shown in panel b). The inset of the figure shows the particle-to-fluid apparent slip velocity, defined as the difference between the mean velocity profiles of the fluid and solid phase. The negative minimum in the near-wall region is due to the higher particle velocity where the fluid velocity is vanishing. The local maximum of positive slip velocity occurs in the buffer region, which is an indirect indication of the well-known tendency of near-wall inertial particles to oversample regions of low streamwise fluid velocity observed in numerous studies (e.g. Rouson & Eaton Reference Rouson and Eaton2001; Marchioli & Soldati Reference Marchioli and Soldati2002, among others). The higher the mass loading, the weaker this slip velocity is. Case SD, in particular, shows a deeper minimum of the slip velocity, and virtually no slip between the two phases beyond 
$y>10\delta _\nu$. This deeper minimum is reflected in more pronounced wakes for particles moving faster than the fluid, see figure 4. Interestingly, the fluid-to-particle apparent slip velocity does not attain positive values for case SD, meaning that particles are less prone to sampling low-speed regions. Since the low (and high) speed regions are located further away from the wall, particles departing from the wall may not reside for long enough in these regions. Finally, we should note that wider streak spacing with larger wall-normal extent is often a feature of turbulent drag reduction (see e.g. Tiederman, Luchik & Bogard Reference Tiederman, Luchik and Bogard1985).
Figure 5. Inner- (a) and outer-scaled (b) mean streamwise velocity profiles. Lines – fluid velocity; symbols – particle velocity. The insets show the apparent fluid-to-particle slip velocity, defined as the difference between the intrinsic mean velocity of each phase.
Figure 6 presents the second-order moments of the fluid and particle velocities. Particle fluctuations are usually higher in the near-wall region and smaller or similar to that of the fluid in the bulk. The modifications with respect to the unladen case are much smaller for cases VD and SD, compared with case SD. As discussed in Costa et al. (Reference Costa, Brandt and Picano2020), the minor deviations of VD with respect to the unladen case are attributed to a lower statistical convergence of the massive particle-resolved simulation, and a slight two-way coupling effect, while the more pronounced differences in case D are caused by the increase in mean wall shear. Compared with these two cases, case SD shows three significant differences: a strong reduction of the streamwise velocity fluctuations near the wall and enhancement in the bulk, a large decrease of the wall-normal and spanwise velocity variances across the channel, and a reduction in Reynolds shear stresses. This observation is consistent with particle-modelled DNS and experiments of particle-laden wall-bounded turbulence (see, e.g. Kulick, Fessler & Eaton Reference Kulick, Fessler and Eaton1994; Zhao et al. Reference Zhao, Andersson and Gillissen2010; Richter Reference Richter2015; Capecelatro et al. Reference Capecelatro, Desjardins and Fox2018; Wang et al. Reference Wang, Fong, Coletti, Capecelatro and Richter2019). Focusing on the fluid Reynolds shear stress, case D features a slight increase of the peak value and slope in the outer region, while the opposite applies to case SD. Note that the particle Reynolds stresses are always larger than the fluid ones.
Figure 6. Second-order moments of mean fluid (lines) and particle (symbols) velocity. (a) Inner-scaled streamwise, (b) wall-normal and (c) spanwise velocity variances, and (d) outer-scaled Reynolds stresses profile. The legend in (a) holds for all panels.
Despite the overall increase in drag, the significant decrease of the Reynolds shear stresses for case SD is consistent with the reduced contribution to the overall drag from the fluid turbulence shown in figure 2 by the friction Reynolds number based on the velocity based on the slope of the Reynolds stresses, 
$u_\tau ^{t} = \sqrt {\partial _{y/h} \langle u^{\prime } v^{\prime } \rangle |_{y/h=1}}$. This quantity initially increases with the mass loading (case VD to D), denoting an enhancement of turbulence and its induced drag. Further increasing the solid mass fraction, (case SD), 
$Re_{\tau }^{t}$ decreases – even below the single-phase flow – while the conventional friction Reynolds number increases. In other words, while turbulence is attenuated, a fundamentally different mechanism acts to anyway increase the drag. Also, despite showing drag increase, case D features statistics that highly resemble those of the single-phase flow. Hence, only case SD shows a more intricate two-way coupling mechanism, at play over the whole channel region.
To gain further insight, we examine the streamwise momentum balance of the fluid,
\begin{equation} {\tau}= {\rho}\left(1-\left\langle \phi\right\rangle \right) \left(\nu \frac{\mathrm{d}\left\langle u\right\rangle }{\mathrm{d}y} - \left\langle u_{f}^{\prime} v_{f}^{\prime}\right\rangle \right) + \left\langle \phi\right\rangle \tau_p = \rho u_\tau^2\left(1-\frac{y}{h}\right)\mathrm{,} \end{equation}
where the first and second terms denote the fluid viscous and Reynolds stresses, while the last term the total particle stress 
$\tau _p = \tau _{p,\nu }-\rho _p\langle u_p^{\prime } v_p^{\prime }\rangle$, with 
$\tau _{p,\nu }$ including the viscous and collisional particle stresses. The data are shown in figure 7(a–c), where the insets depict the local relative contribution of each term to the total momentum transfer; panel (d) reports the share of each momentum transfer mechanism to the overall drag, normalized by the overall single-phase drag. These contributions are related to the stress profiles via the Fukagata, Iwamoto and Kasagi (FIK) identity (see Fukagata, Iwamoto & Kasagi Reference Fukagata, Iwamoto and Kasagi2002; Yu et al. Reference Yu, Xia, Guo and Lin2021), with the contribution of each stress mechanism 
$\tau _i$ given by the following weighted integral: 
$\smallint _0^{1} 3(1-y)\tau _i\,\mathrm {d}\kern 0.05em y$. Expectedly, the stress budget pertaining to case VD resembles that of the single-phase flow, with negligible contribution from the particles. Case D displays a noticeable, though small, contribution of the particle stresses which have a peak of approximately 
$5\,\%$ close to the wall and are of the order of the viscous stresses in the bulk. Finally, the total particle stresses show a significant relative contribution to the total in case SD, of approximately 
$20\,\%$. Examining panel (d) of the same figure, we note that the sum of the viscous and turbulent contributions corresponds to the single-phase drag in case VD, with the particle stress showing an almost vanishing two-way effect. In case D, the sum of the viscous and turbulent stress contributions is higher than the overall single-phase drag, which is further increased by the particle stress. This indicates that the overall drag is increased by two different mechanisms: a direct one – the particles induce an extra stress, localized near the wall where the average slip velocity is not negligible; and an indirect effect – the solid phase enhances the Reynolds stress and hence the turbulence. Notably, the sum of the viscous and turbulent stress contributions is smaller than the overall single-phase drag in case SD; the particle stress, however, counterbalances this reduction and the overall drag is still higher than the corresponding unladen case. This confirms that we have turbulent drag reduction and a strong contribution of the total particle stress, which combine to a net drag increase.
Figure 7. Budget of streamwise momentum as a function of the outer-scaled wall-normal distance, for the different cases. The inset shows the same budget, but normalized by the local value of the total stresses. Panel (d) shows the relative contribution of each term in the stress budget to the mean wall friction, normalized by that of the unladen case, i.e. 
$\smallint _0^{1} 3(1-y)\tau _i\,\mathrm {d}y/(\rho u_\tau ^{2\,sph})$ (Fukagata et al. Reference Fukagata, Iwamoto and Kasagi2002; Yu et al. Reference Yu, Xia, Guo and Lin2021), with 
$\tau _i$ the different contributions, and 
$u_\tau ^{2\,sph}$ the wall friction velocity of the unladen case.
Let us take the limit of vanishing volume fraction, 
$\phi \to 0$, but finite mass fraction 
$\psi = \rho _p\phi$, of (3.2),
\begin{equation} \tau_{two\text{-}way} = \rho\left(\nu \frac{\mathrm{d}\left\langle u\right\rangle }{\mathrm{d}y} - \left\langle u_f^{\prime} v_f^{\prime}\right\rangle \right) - \left\langle \psi\right\rangle \left\langle u_p^{\prime} v_p^{\prime}\right\rangle\approx\rho u_\tau^{2}\left(1-\frac{y}{h}\right). \end{equation}This limit corresponds to negligible particle excluded volume, but finite effects of the particle mass, i.e. two-way coupling conditions. In this limit, the contribution from the correlated particle fluctuations corresponds to a particle direct contribution to the drag. When plotting the terms in (3.3), see figure 8, we note that the difference between the total stress and the sum of the different terms is very small, confirming that the contribution of the particle inertial shear stress to the budget is dynamically significant and increases with the mass fraction. Moreover, the ‘two-way coupling’ budget explains the true nature of the turbulence attenuation for case SD. As the mass loading is increased, the fluid Reynolds stresses progressively give in to particle Reynolds stresses, which alter the nature of the flow. A significant portion of the total stress is spent to accelerate the inertial particles, which in turn lowers the fluid Reynolds stresses. Consequently, the flow shows fluid turbulence attenuation. Unlike correlated fluid–fluid or fluid–particle velocity fluctuations, particle–particle correlated motions do not sustain near-wall turbulence, as demonstrated in Capecelatro et al. (Reference Capecelatro, Desjardins and Fox2018) for the fluid turbulence kinetic energy budget.
Figure 8. Budget of streamwise momentum in the two-way coupling limit, for cases D and SD, (see (3.3)). The dashed grey line corresponds to the difference between the sum of the terms in (3.3), and the total stress, to assess the validity of the approximation. Case VD was omitted from this figure, since given the negligible two-way coupling effect the budget is trivial, resembling that of the single-phase flow.
Since the net effect of the particle Reynolds stresses is an increased drag for case SD, we expect that at higher mass fractions the drag may eventually decrease – the contribution of Reynolds stresses to the momentum transport would decrease and the flow becomes smoother and eventually laminarizes. Consequently, the correlated particle velocity fluctuations would also decrease, drastically reducing the overall flow drag. We have actually observed such a laminarizing and drag-reducing trend in a preliminary simulation of a case with the same governing parameters as SD, but 
$10$ times the number of particles (note that here particle excluded volume effects may grow important). Note that, reductions of the skin friction coefficient have also been observed in the two-way coupling point-particle DNS by Vreman (Reference Vreman2007) and Zhao et al. (Reference Zhao, Andersson and Gillissen2010) and in the volume-filtered Euler–Lagrange simulations of Capecelatro et al. (Reference Capecelatro, Desjardins and Fox2018), with the latter study observing laminarization at sufficiently high mass fraction. In particular, Capecelatro et al. (Reference Capecelatro, Desjardins and Fox2018) reported that the fluid velocity fluctuations were re-energized at even higher mass fractions, due to the growing importance of turbulence-producing particle–fluid velocity correlations. New experiments and fully resolved simulations are hence needed to further support and explain these observations at higher mass fractions.
4. Conclusions
We have used particle-resolved DNS to study the near-wall turbulence modulation by small inertial particles. Three cases have been considered with volume fractions progressively increased by one order of magnitude, chosen in the one-way and two-way coupling regimes. The two densest cases show a non-negligible turbulence modulation, however, of a fundamentally different nature. The case with 
$\varPsi =3.4\,\%$ (
$\varPhi =0.034\,\%$) features drag increase due to increased inertia near the wall, with increased Reynolds stresses due to particles travelling at high particle-to-fluid slip velocity, and dynamics close to that of the single-phase flow at slightly higher Reynolds number. Increasing the volume and mass fraction by a factor of 
$10$, the particle presence is felt over the entire channel, however, the relative drag increase is significantly reduced with respect to that of the previous configuration. The streamwise momentum balance in the two-way coupling limit of vanishing volume fraction, but finite mass fraction, reveals that correlated particle fluctuations seize a larger share of the total stresses as the mass fraction increases. This results in a reduction of the fluid Reynolds stresses, and consequently of the turbulent drag. For the parameter setting considered here, this turbulent drag reduction is counterbalanced by the contribution of correlated particle velocity fluctuations; however, a further increase in mass fraction may lead to a net drag reduction, something which has been noticed in previous particle-modelled DNS, and should inspire future experimental and particle-resolved numerical studies.
Acknowledgements
We acknowledge the computing time provided by SNIC (Swedish National Infrastructure for Computing), and PRACE for awarding us access to the supercomputer Marconi, based in Italy at CINECA under project 2017174185 – DILPART.
Funding
This work was supported by the European Research Council grant no. ERC-2013-CoG-616186, TRITOS, the grant BIRD192032/19 from the University of Padova and the University of Iceland Recruitment Fund grant no. 1515-151341, TURBBLY.
Declaration of interests
The authors report no conflict of interest.
