1 Introduction
Turbulent flows laden with small inertial particles are found in many environmental and industrial contexts. These flows are inherently chaotic and multi-scale, with the interphase coupling categorized by the relevance of the dispersed phase to the overall dynamics (Elghobashi Reference Elghobashi1994). The so-called one-way-coupling regime corresponds to low volume and mass fractions of the solid phase, when particle–fluid interactions are negligible, and particle–particle interactions unlikely. Increasing the solid mass fraction while keeping the volume fraction at the same order of magnitude results in a regime where the overall particle load becomes high enough to modulate the turbulent flow, while particle–particle interactions remain negligible – the two-way-coupling regime. Finally, further increasing the volume fraction results in a regime where particle–particle interactions are also important – the four-way-coupling regime. From a modelling perspective, another important distinction concerns the particle size. When the ratio between the particle size and the Kolmogorov scale is smaller than one, the term ‘point particle’ is used, and particle–fluid coupling is considered to take place at a single point. Conversely, when the size ratio is large, the particles are termed ‘finite sized’ (Balachandar & Eaton Reference Balachandar and Eaton2010).
Particle-laden turbulence in the one-way-coupling and point-particle limit has been the subject of numerous studies throughout the last decades (for recent reviews, see e.g. Toschi & Bodenschatz (Reference Toschi and Bodenschatz2009), Balachandar & Eaton (Reference Balachandar and Eaton2010)). In these cases, it is assumed that the local properties of an undisturbed flow at the particle position drive the dispersed-phase dynamics (Gatignol Reference Gatignol1983; Maxey & Riley Reference Maxey and Riley1983). For relatively high particle-to-fluid density ratios, the particle dynamics is often simplified to a balance between particle inertial and drag forces, where the latter is a function of the so-called particle-to-fluid slip velocity. Under these conditions, particles display preferential clustering even in homogeneous and isotropic flows (Toschi & Bodenschatz Reference Toschi and Bodenschatz2009). When the flow is inhomogeneous, particles tend to migrate from regions of high to low turbulence intensity due to turbophoresis (Reeks Reference Reeks1983). In turbulent wall-bounded flows, in particular, the particle distribution is driven by the interplay between small-scale clustering, turbophoresis and the interaction between the particles and near-wall turbulence structures (Soldati & Marchioli Reference Soldati and Marchioli2009; Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012). When the system reaches a statistical equilibrium, particles tend to accumulate in the low-speed regions near the wall, resulting in a very inhomogeneous local particle concentration (see e.g. Fessler, Kulick & Eaton (Reference Fessler, Kulick and Eaton1994), Uijttewaal & Oliemans (Reference Uijttewaal and Oliemans1996), Marchioli et al. (Reference Marchioli, Giusti, Salvetti and Soldati2003, Reference Marchioli, Soldati, Kuerten, Arcen, Taniere, Goldensoph, Squires, Cargnelutti and Portela2008) and Kuerten (Reference Kuerten2006)).
In flows with locally higher mass loading, two-way-coupling effects may become important. In addition to solving the particle dynamics, a two-way-coupling point-particle algorithm must impose a localized momentum source/sink corresponding to the particle back-reaction to the flow. The challenge here is to determine the local characteristics of the undisturbed flow field, while the flow itself is subjected to a local disturbance due to the presence of the dispersed phase (see e.g. Gualtieri et al. (Reference Gualtieri, Picano, Sardina and Casciola2015)). The classical approach is the particle-in-cell method, developed by Crowe, Sharma & Stock (Reference Crowe, Sharma and Stock1977). Although widely used, the success of this method strongly depends on the number of particles per grid cell, i.e. it does not converge with grid refinement. Approaches for a consistent and more robust treatment are the object of active research, as shown by the number of recent studies (e.g. Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2015; Horwitz & Mani Reference Horwitz and Mani2016; Ireland & Desjardins Reference Ireland and Desjardins2017). We should also note the recent efforts to account for pairwise particle–particle hydrodynamic interactions in point-particle models (Akiki, Moore & Balachandar Reference Akiki, Moore and Balachandar2017). Investigations of particle-laden turbulent flows in the two-way-coupling regime are found in, for example, Vreman et al. (Reference Vreman, Geurts, Deen, Kuipers and Kuerten2009) and Capecelatro, Desjardins & Fox (Reference Capecelatro, Desjardins and Fox2018).
Despite the numerous studies involving direct numerical simulations (DNS) with point-particle methods, validations of the underlying assumptions remain scarce, even for simple flows. For instance, the parameter range where the Maxey–Riley–Gatignol equations are valid remains elusive (Bergougnoux et al. Reference Bergougnoux, Bouchet, Lopez and Guazzelli2014). Though experiments in particle-laden turbulence are insightful (Eaton & Fessler Reference Eaton and Fessler1994; Kaftori, Hetsroni & Banerjee Reference Kaftori, Hetsroni and Banerjee1995), parameter-matched numerical simulations remain challenging, because of the numerical limitations in terms of Reynolds number and the need for experiments in well-controlled, often idealized, configurations. Recent efforts in this direction have been undertaken in the recent study by Wang et al. (Reference Wang, Fong, Coletti, Capecelatro and Richter2019).
Another possible reference for new models is a particle-resolved (also denoted interface-resolved) DNS, i.e. a DNS that resolves the flow around the surface of each small particle. Despite the great computational challenge of such computations, the first direct comparisons between point-particle models and particle-resolved simulations have started to appear for forced homogeneous isotropic turbulence (HIT) with fixed particles (Vreman Reference Vreman2016) and decaying HIT with moving particles (Schneiders, Meinke & Schröder Reference Schneiders, Meinke and Schröder2017a,Reference Schneiders, Meinke and Schröderb; Fröhlich et al. Reference Fröhlich, Schneiders, Meinke and Schröder2018; Mehrabadi et al. Reference Mehrabadi, Horwitz, Subramaniam and Mani2018).
We consider particle-laden turbulent channel flow, a widely studied case with minimal governing parameters and particularly important to benchmark models for wall-bounded particle transport. Though interface-resolved DNS of these flows are in general quite demanding, recent studies have demonstrated that massively parallel simulations of wall-bounded flows with 
$O(10^{6})$ interface-resolved particles and 
$O(10^{9}{-}10^{10})$ grid points have become feasible (see Costa et al. Reference Costa, Picano, Brandt and Breugem2016; Kidanemariam & Uhlmann Reference Kidanemariam and Uhlmann2017; Horne & Mahesh Reference Horne and Mahesh2019).
We present interface-resolved DNS of gravity-free turbulent channel flow laden with small inertial particles (with a size of three viscous units, and 100 times denser than the fluid), in the dilute regime. Two cases are considered, with bulk volume solid fractions that approach the one-way-coupling regime: 0.003 % and 0.03 %. Both cases show a turbophoretic particle drift towards the wall, as expected. However, due to the inhomogeneous particle distribution, non-negligible two-way-coupling effects are found for the case with the larger volume fraction, but not in the more dilute case. The interface-resolved DNS are complemented with corresponding one-way-coupled point-particle DNS at the same Reynolds number using models for the drag and lift forces. Considering only the balance between inertia and nonlinear drag, the turbophoretic wall concentration is overestimated by a factor of 2–3 with respect to the more diluted case considering fully resolved particles. We show that the missing ingredient is a shear-induced lift force that reduces the turbophoretic drift. Finally, we demonstrate that these dynamics can be well captured using the model proposed in the seminal work of Saffman (Reference Saffman1965).
2 Methods and computational set-up
The numerical method solves the continuity and Navier–Stokes equations for an incompressible Newtonian fluid with density 
$\unicode[STIX]{x1D70C}_{f}$ and dynamic viscosity 
$\unicode[STIX]{x1D707}$ (kinematic viscosity 
$\unicode[STIX]{x1D708}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}_{f}$),
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{f}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}, & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{u}$ is the fluid velocity vector, 
$\unicode[STIX]{x1D748}=-(p+p_{e})\boldsymbol{I}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})$, with 
$p+p_{e}$ being the fluid pressure with respect to an arbitrary constant reference value; the resulting term 
$\unicode[STIX]{x1D735}p_{e}$ in (2.2) corresponds to a uniform pressure gradient that may serve as driving force for the flow.
2.1 Interface-resolved particle simulations
The fluid equations are coupled to the Newton–Euler equations governing the motion of a spherical particle with mass 
$m$, and moment of inertia 
$I$,
$$\begin{eqnarray}\displaystyle & \displaystyle m\dot{\boldsymbol{U}}=\oint _{\unicode[STIX]{x2202}V}\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}A+\boldsymbol{F}_{c}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle I\dot{\unicode[STIX]{x1D734}}=\oint _{\unicode[STIX]{x2202}V}\boldsymbol{r}\times (\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n})\,\text{d}A+\boldsymbol{T}_{c}, & \displaystyle\end{eqnarray}$$ where we use Newton’s dot notation for time differentiation. Here 
$\boldsymbol{U}$ and 
$\unicode[STIX]{x1D734}$ denote the particle linear and angular velocity vectors, 
$\boldsymbol{r}=\boldsymbol{x}-\boldsymbol{X}$ the position vector with respect to the particle centroid, 
$\boldsymbol{n}$ the outward-pointing unit vector normal to the particle surface 
$\unicode[STIX]{x2202}V$, 
$A$ the surface area of the particle, and 
$\boldsymbol{F}_{c}$ and 
$\boldsymbol{T}_{c}$ correspond to external forces and torques associated with short-range inter-particle or particle–wall interactions (such as solid–solid contact).
Equations (2.1)–(2.2) and (2.3)–(2.4) are coupled through the imposition of no-slip and no-penetration boundary conditions at the particle surface:
$$\begin{eqnarray}\boldsymbol{U}+\unicode[STIX]{x1D734}\times \boldsymbol{r}\stackrel{!}{=}\boldsymbol{u},\quad \forall ~\boldsymbol{x}\in \unicode[STIX]{x2202}V.\end{eqnarray}$$The Navier–Stokes equations are solved with a 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). The solver was extended with a direct forcing immersed-boundary method (IBM) for particle-laden flows developed by Breugem (Reference Breugem2012) (see also Uhlmann (Reference Uhlmann2005)) and the lubrication/soft-sphere collision model for short-range particle–particle and particle–wall interactions in Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015). Several recent studies describe the method, present validations and assess its computational performance (Picano, Breugem & Brandt Reference Picano, Breugem and Brandt2015; Costa Reference Costa2018; de Motta et al. Reference de Motta, Costa, Derksen, Peng, Wang, Breugem, Estivalezes, Vincent, Climent and Fede2019).
Turbulent channel flow is simulated in a domain periodic in the streamwise (
$x$) and spanwise (
$z$) directions, with no-slip/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 physical and computational parameters are reported in table 1. Since the IBM requires a fixed, regular Eulerian grid, resolving the spherical particles with 
$O(10)$ grid points over the diameter is, by far, what dictates the grid resolution. Hence, the spatial resolution for the interface-resolved simulations is, in each direction, approximately one order of magnitude larger than what is required for single-phase simulation, leading to 
${\sim}10^{10}$ grid points.
Computational parameters of the DNS dataset. Here 
$\unicode[STIX]{x1D6F7}$ denotes the bulk solid volume fraction, 
$\unicode[STIX]{x1D6F9}$ the bulk solid mass fraction and 
$N_{p}$ the total number of particles. For all cases, the bulk Reynolds number 
$Re_{b}=5600$ (i.e. unladen friction Reynolds number 
$Re_{\unicode[STIX]{x1D70F}}^{sph}\approx 180$); particle size ratio 
$D/(2h)=1/120$; particle-to-fluid mass density ratio 
$\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D70C}}=100$, corresponding to a particle diameter (in viscous units) 
$D^{+}=3$ and a Stokes number 
$St=50$. For the interface-resolved cases, 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/\unicode[STIX]{x0394}x=12$ grid points over the particle diameter (420 Lagrangian grid points in total). For case PP, 500 000 point particles have been simulated in a grid that is five times coarser in each direction. Note that the last column reports the friction Reynolds number 
$Re_{\unicode[STIX]{x1D70F}}$ extracted from the simulations.
The solid–solid contact between the rigid particles is frictionless, with a normal dry coefficient of restitution, i.e. the ratio of rebound to impact velocity in the absence of lubrication effects 
$e_{n,d}=0.97$. As suggested in Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015), the particle–particle/wall collision time is set to 
$N_{c}\approx 8$ times the time step of the Navier–Stokes solver. Lubrication corrections for the normal lubrication force when the distance between two solid surfaces is small are used as described in Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015), with the same model parameters as those used in Costa et al. (Reference Costa, Picano, Brandt and Breugem2018).
To mimic a flow close to the one-way-coupling regime, we consider two low values of solid volume fraction, 
$\unicode[STIX]{x1D6F7}\simeq 3\times 10^{-5}$, denoted very dilute (VD), and 
$\unicode[STIX]{x1D6F7}\simeq 3\times 10^{-4}$, dilute (D). The bulk Reynolds number 
$Re_{b}=U_{b}(2h)/\unicode[STIX]{x1D708}=5600$, which corresponds to an unladen friction Reynolds number 
$Re_{\unicode[STIX]{x1D70F}}^{sph}=u_{\unicode[STIX]{x1D70F}}h/\unicode[STIX]{x1D708}\approx 180$; where 
$U_{b}$ is the bulk flow velocity and 
$u_{\unicode[STIX]{x1D70F}}$ the wall friction velocity. The particle properties are chosen close to those used in point-particle simulations at the same 
$Re_{\unicode[STIX]{x1D70F}}$, which result in strong turbophoresis and large-scale clustering (Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012). This choice corresponds to a particle Reynolds number 
$Re_{p}=Du_{\unicode[STIX]{x1D70F}}/\unicode[STIX]{x1D708}=D^{+}=3$, and Stokes number 
$St_{p}=\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D70C}}Re_{p}^{2}/18=50$, where 
$D$ is the particle diameter and 
$\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D70C}}$ the particle-to-fluid mass density ratio.
Two strategies have been considered for achieving an initial condition close to the fully developed turbulent state with little computational cost: (1) using a fully developed one-way-coupled point-particle simulation at a much coarser grid, interpolated onto the grid of the interface-resolved case; and (2) using the fully developed state of a coarse interface-resolved DNS, where the fluid is well resolved, but the particles are under-resolved (four grid points over the particle diameter). The latter approach was preferred, as both the particle distribution and two-way-coupling effects are closer to the fully developed turbulent state of the interface-resolved DNS.
Finally, the mesoscale-averaged profiles reported in this paper correspond to mean intrinsic averages. The intrinsic average of an observable 
$o$ is given by
$$\begin{eqnarray}\langle o\rangle (y_{j})=\frac{\displaystyle \mathop{\sum }_{ik,t}o_{ijk,t}C_{ijk,t}}{\displaystyle \mathop{\sum }_{ik,t}C_{ijk,t}},\end{eqnarray}$$ where 
$C_{ijk,t}$ is the volume fraction of the specific phase at the grid cell 
$ijk$ and instant 
$t$, and 
$y_{j}$ the wall-normal location of the averaging bin, which extends over the entire domain in the two homogeneous directions. These quantities have been obtained from an ensemble average on approximately 200 instances sampled over a period of 
$100h/U_{b}$ once a fully developed state has been reached. The velocity statistics of the solid phase consider the rigid-body motion throughout the volume of the particles, i.e. account for both translation and rotation. We should note that higher-order particle statistics pertaining to case VD still show some high-frequency fluctuation away from the wall, due to the extremely low local volume fraction.
Each interface-resolved simulation cost approximately 10 million core hours on the supercomputer Marconi based in Italy at CINECA.
2.2 One-way-coupled point-particle simulations
The interface-resolved simulations are complemented with one-way-coupled point-particle simulations. Given the large density ratio, we first assume that the particle dynamics simplifies to a balance between inertial and nonlinear Schiller–Naumann (Schiller & Nauman Reference Schiller and Nauman1933) drag forces, as often done in the literature. Test simulations confirmed that the results are not sensitive to Faxén corrections in the particle dynamics. In addition to this standard case, denoted PP (point particle), we also investigate the role of shear-induced lift forces, which may be important close to the wall.
The point-particle dynamics is driven by the conservation of linear momentum of the particle,
$$\begin{eqnarray}\displaystyle & m\dot{\boldsymbol{U}}=\boldsymbol{F}_{d}+\boldsymbol{F}_{l}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \dot{\boldsymbol{X}}=\boldsymbol{U}, & \displaystyle\end{eqnarray}$$where
$$\begin{eqnarray}\boldsymbol{F}_{d}=-3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}D\boldsymbol{U}_{s}(1+0.15Re_{p}^{0.687})\end{eqnarray}$$and
$$\begin{eqnarray}\boldsymbol{F}_{l}=1.615J\unicode[STIX]{x1D707}D|\boldsymbol{U}_{s}|\sqrt{{\displaystyle \frac{D^{2}|\unicode[STIX]{x1D74E}|}{\unicode[STIX]{x1D708}}}}\,{\displaystyle \frac{\unicode[STIX]{x1D74E}\times \boldsymbol{U}_{s}}{|\unicode[STIX]{x1D74E}||\boldsymbol{U}_{s}|}},\end{eqnarray}$$ with 
$\boldsymbol{U}_{s}=\boldsymbol{U}-\boldsymbol{u}|_{\boldsymbol{x}=\boldsymbol{X}}$ the particle-to-undisturbed-flow slip velocity, and 
$\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D74E}|_{\boldsymbol{x}=\boldsymbol{X}}$ the undisturbed flow vorticity evaluated at the particle position. Two widely used approaches are considered for modelling the lift force 
$\boldsymbol{F}_{l}$. First, 
$\boldsymbol{F}_{l}$ is modelled simply by the Saffman lift force (Saffman Reference Saffman1965), i.e. 
$J=1$ in (2.10). Second, we consider a correction for finite-Reynolds-number effects that fits the tabulated results in McLaughlin (Reference McLaughlin1991) proposed by Mei (Reference Mei1992), where 
$J$ in (2.10) is a function of a parameter 
$\unicode[STIX]{x1D700}=\sqrt{|\unicode[STIX]{x1D74E}|\unicode[STIX]{x1D708}}/|\boldsymbol{U}_{s}|$ (note that the formula in the original reference has a typo (cf. Loth & Dorgan Reference Loth and Dorgan2009)):
$$\begin{eqnarray}J=0.3\left(1+\tanh [{\textstyle \frac{5}{2}}(\log _{10}\unicode[STIX]{x1D700}+0.191)]\right)\left({\textstyle \frac{2}{3}}+\tanh (6\unicode[STIX]{x1D700}-1.92)\right).\end{eqnarray}$$We should note that this type of correction for finite-Reynolds-number effects has been employed in numerous studies, and often as basis for more elaborate models (see e.g. Uijttewaal & Oliemans (Reference Uijttewaal and Oliemans1996), Wang et al. (Reference Wang, Squires, Chen and McLaughlin1997), Marchioli & Soldati (Reference Marchioli and Soldati2002), Marchioli, Picciotto & Soldati (Reference Marchioli, Picciotto and Soldati2007), Loth & Dorgan (Reference Loth and Dorgan2009)). The results obtained using the simple Saffman lift term are here denoted PP-Saffman, and those obtained with the correction as PP-McLaughlin.
In all point-particle simulations, the fluid observables evaluated at the particle position are computed using trilinear interpolation, and the particle positions integrated in time with the same third-order low-storage Runge–Kutta (RK) scheme used for the interface-resolved simulations. As regards particle–wall collisions, a perfectly elastic hard-sphere rebound is adopted: if a particle is about to overlap with the wall during an RK substep with size 
$\unicode[STIX]{x0394}t_{s}$, the particle final wall-normal position 
$Y$ and velocity 
$V$ are computed analytically from the two-body kinematics assuming a piecewise-constant particle velocity. For instance, for a particle colliding onto the lower wall (
$y=0$) with centroid position 
$Y_{in}$ and impact wall-normal velocity 
$V_{in}<0$, we have, if 
$Y_{in}+V_{in}\unicode[STIX]{x0394}t_{s}<D/2$,
$$\begin{eqnarray}\displaystyle & V=-V_{in}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & Y=V\unicode[STIX]{x0394}t_{s}-(Y_{in}-D/2), & \displaystyle\end{eqnarray}$$3 Results
3.1 Interface-resolved simulations and the role of lift force
Visualization of case D. Isocontours of the 
$Q$-criterion, 
$Q=20(U_{b}/h)^{2}$ (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988) coloured by the local wall-normal velocity 
$v$ (increasing from blue to white with minimum/maximum values of 
$\mp 0.1U_{b}$). The interface-resolved particles are depicted in orange. The isosurfaces have been rendered from a field in a grid four times coarser in each direction than that used in the DNS.
Contours of streamwise fluid velocity (flow from left to right) in the plane 
$y/h=0.056$, i.e. 
$y^{+}\approx 10$. Particles with wall distance 
$y_{p}/h<0.06$ are depicted in white colour. Panels (a) and (b) correspond to cases D and PP with matched bulk number density.
Figure 1 shows a visualization of the particle-laden flow for case D. Similarly to case VD (not shown), the near-wall large-scale structures resemble those of the unladen flow. Nonetheless, the localized effect of the particles is evident, as depicted by the high-vorticity trail due to their wakes. The overall drag is a measure of the cumulative effect of the localized disturbances: this is reported in the last column of table 1 in terms of a friction Reynolds number 
$Re_{\unicode[STIX]{x1D70F}}=u_{\unicode[STIX]{x1D70F}}h/\unicode[STIX]{x1D708}$. While the very dilute (VD) case shows, as expected, approximately the same drag as the unladen flow, the dilute case (D) shows approximately 5 % higher friction Reynolds number (i.e. approximately 10 % increase in pressure drop). This is a remarkable drag increase, better comparable to what has been observed for finite-sized neutrally buoyant particle suspensions with volume fraction of 5 % (Picano et al. Reference Picano, Breugem and Brandt2015). The modulation of the near-wall structures for case D is noticeable in figure 2, where we display the streamwise velocity contours close to the wall. In addition to showing a microscopic footprint, the resolved particles disrupt streamwise-correlated near-wall structures, as can be noticed in the spanwise autocorrelation of the wall-normal fluid velocity, given by 
$R_{vv}^{z}(\unicode[STIX]{x1D6FF}_{z})=\langle v^{\prime }(z)v^{\prime }(z+\unicode[STIX]{x1D6FF}_{z})\rangle /\langle v^{\prime 2}(z)\rangle$ shown in figure 3(a). Clearly, the magnitude of the negative peak in correlation – footprint of the near-wall low- and high-speed streaks – is dampened for case D, while remaining very close to that of single-phase flow in case VD. We will see that these observations are closely connected to the inhomogeneous particle distribution near the wall.
(a) Spanwise autocorrelations of the wall-normal velocity 
$R_{vv}^{z}$, at 
$y/h=0.056$, i.e. 
$y^{+}\approx 10$. (b) Local solid mass fraction as a function of the outer-scaled wall-normal distance (same legend as the former panel). The inset shows the corresponding bulk-normalized profile, versus the wall-distance in particle diameters.
(a) Inner- and (b) outer-scaled mean velocity profiles for the different cases (lines, fluid velocity; symbols (colour-matched), particle velocity). The inset in panel (a) shows the inner-scaled difference between fluid and particle velocity profiles 
$\unicode[STIX]{x0394}\langle u\rangle$.
Outer-scaled second-order moments of particle velocity: (a) streamwise velocity r.m.s., (b) wall-normal velocity r.m.s. and (c) ditto for spanwise velocity. (d) Reynolds stress profile with inner-scaled inset. Lines, fluid; symbols (colour-matched), particles.
Figure 3(b) depicts the wall-normal profiles of the local solid mass fraction as a function of the wall-normal distance. The inset shows the same quantity divided by the bulk values versus the wall-normal distance in particle diameters. The profiles show a near-wall peak at 
$y\approx D/2$ approximately one order of magnitude larger than the bulk value, which explains the drag increase observed in case D. Despite the low bulk mass loading, the near-wall mass fraction becomes high enough to modulate the flow in this critical region (
$\unicode[STIX]{x1D6F9}\simeq 0.2$). We should note that the corresponding local volume fraction is still too low (
$\unicode[STIX]{x1D6F7}\simeq 0.002$) for particle–particle interactions to be significant: the collision frequency in the viscous sublayer is estimated to be less than 
$10^{-5}$ collisions per particle pair per bulk time unit 
$h/U_{b}$ for the less dilute case D. The interface-resolved simulations show pronounced differences when compared to the point-particle results (PP), which consider only the nonlinear drag in the dynamics. In particular, the latter overpredicts the concentration peak by a factor of 2–3. Moreover, the concentration profile of case PP shows a more gentle decrease away from the wall, while in the other cases the concentration peak corresponds to a single particle layer. We attribute this localized peak to the combination of turbophoresis, the kinematic constraint that the wall imposes on the particles, and the stabilizing effect of lubrication forces (Picano et al. Reference Picano, Breugem and Brandt2015). Clearly the mechanisms causing preferential concentration in the standard point-particle simulation are greatly weakened. At least for the very dilute case VD, these differences cannot be explained in terms of turbulence modulation nor particle–particle interactions (i.e. two/four-way-coupling effects). Hence, the cause must be the different near-wall dynamics of (isolated) particles.
Figure 4(a) shows the inner-scaled profiles of mean streamwise fluid and particle velocity for the different cases. While the fluid velocity profile of case VD matches that of the unladen flow, case D shows significant deviations. Since the inner-scaled profiles of all cases agree in the outer region (see figure 4b), these deviations are attributed to the increase in wall shear. We thus confirm that, somewhat unexpectedly, only case VD satisfies the one-way-coupling assumption, i.e. the turbulence modulation by the dispersed phase is negligible. This case serves therefore as benchmark to assess the validity of the point-particle models to predict the dispersed-phase dynamics in the one-way-coupling regime.
As for the particles, they tend to flow slower than the fluid in the buffer layer (
$10\lesssim y/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}\lesssim 40{-}50$, with 
$\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}=\unicode[STIX]{x1D708}/u_{\unicode[STIX]{x1D70F}}$; see also the inset in figure 4a, showing the inner-scaled difference between the profiles of each phase). This velocity difference is more pronounced in case PP, where the particle velocity reduction is clear also very close to the wall. This has been observed in previous studies using point-particle DNS (Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012), and is attributed to the preferential sampling of the fluid low-speed streaks near the wall. The particle-resolved cases show a much weaker reduction of the mean particle velocity, which suggests that particles reside in the low-speed regions for shorter times before resuspending into the bulk; see also the discussion below about the particle dynamics. In the viscous sublayer, particles show a slightly higher mean velocity in the interface-resolved simulations. This higher slip velocity causes hot-spots of higher wall shear stress, which favour an increase in overall drag (Costa et al. Reference Costa, Picano, Brandt and Breugem2016, Reference Costa, Picano, Brandt and Breugem2018). Clearly this effect is significant in case D, where the near-wall number density is high enough, but not in case VD.
Figure 5 shows the second-order statistics of fluid and particle velocity. Focusing first on the fluid phase, we see once more that the data for VD tend to those of the single-phase flow, with the minor differences attributed to a better statistical convergence of case PP, and possibly a slight two-way-coupling effect. Conversely, turbulence modulation is evident for case D. Here the Reynolds stresses are higher, consistently with the overall drag increase – see inset in figure 5(d). Moreover, small differences are found for all the velocity root mean square (r.m.s.) values of case D near the wall, where the velocity fluctuations become less anisotropic, i.e. 
$u_{r}$ decreases and 
$v_{r}$ and 
$w_{r}$ increase. This is attributed to the enhanced mixing due to the near-wall particles, whose local mass fraction is high enough for two-way-coupling effects to be significant.
Interestingly, the second-order moments of the particle velocity for the fully resolved one-way-coupling case, VD, strongly differ from those of the point-particle simulations near the wall, while in the bulk the two cases display a similar behaviour. In the bulk, where the local shear is relatively low, the point-particle model succeeds in predicting the particle dynamics. We should note that the same closure for the point-particle dynamics was used in Mehrabadi et al. (Reference Mehrabadi, Horwitz, Subramaniam and Mani2018) for decaying HIT, and the results also compared well to the corresponding interface-resolved case. Closer to the wall (
$y/h\lesssim 0.1$), however, the interface-resolved simulations show higher fluctuation levels than the point-particle reference. One can depict a clear change in trend in case VD for the profiles pertaining to the quantities in the plane at 
$y/h\approx 0.1$. There is a clear local minimum of 
$v_{r}$, and a sudden change in slope for 
$\langle u^{\prime }v^{\prime }\rangle$. The exception is the spanwise velocity r.m.s. 
$w_{r}$, which attains similar values also close to the wall for cases VD and PP. We should note that similar trends for the streamwise and wall-normal particle velocity r.m.s. have been observed in recent experiments of particle-laden turbulent downward flow in a vertical channel (see Fong, Amili & Coletti Reference Fong, Amili and Coletti2019). All these observations suggest differences in the single-particle dynamics, in particular in the way particles approach and depart from the wall in the two models. In the one-way-coupled point-particle DNS, the particle dynamics is modelled by a simple drag law without considering shear-induced lift forces. In this case, particles are driven towards the wall with high velocity by turbophoresis. Their inertia prevents resuspension, resulting in long periods of wall accumulation in low-speed regions, while only few of them drift back into the bulk (Soldati & Marchioli Reference Soldati and Marchioli2009; Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012).
Since at equilibrium the net wall-normal particle flux is zero, a large number of particles accumulate at the wall. Conversely, when the flow around particles is resolved, these tend to reside for much shorter times near the wall before resuspending. Hence, point particles tend to skim along the wall in low-speed streaks for long periods (Soldati & Marchioli Reference Soldati and Marchioli2009; Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012), whereas resolved particles show shorter residence times at the wall, quickly take off and are not preferentially localized in low-speed streaks (see figure 2b).
This faster cycle explains the larger value of 
$v_{r}$ near the wall, and consequently the larger values of 
$u_{r}$ and 
$\langle u^{\prime }v^{\prime }\rangle$ since the fluctuations are correlated through the mean shear. To better quantify this effect, figure 6 shows the average time that a particle close to the wall (i.e. located at 
$y\approx D/2$) needs to exit the viscous sublayer (i.e. to reach a wall-normal position 
$y>5\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$), 
$\unicode[STIX]{x0394}t^{up}$. Figure 6 shows that near-wall particles in the fully resolved cases take approximately the same time to exit the viscous sublayer, which is approximately one order of magnitude shorter than that of the point-particle DNS PP.
The particle dynamics just described suggests that a shear-induced lift force is the missing key ingredient absent in the point-particle model. Such a force plays a very important role in the particle dynamics near the wall, where the mean shear is high (Wang et al. Reference Wang, Squires, Chen and McLaughlin1997; Soldati & Marchioli Reference Soldati and Marchioli2009). It is known that a particle flowing near a wall in a shear flow experiences a strong lift force (Saffman Reference Saffman1965; McLaughlin Reference McLaughlin1991; Cherukat & McLaughlin Reference Cherukat and McLaughlin1994; Bagchi & Balachandar Reference Bagchi and Balachandar2002; Magnaudet Reference Magnaudet2003). A similar mechanism should work for small particles in the viscous sublayer of a turbulent flow.
Inner-scaled average time that a wall-skimming particle takes to reach a wall-normal distance 
$y>5\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ (i.e. to exit the viscous sublayer), 
$\unicode[STIX]{x0394}t^{up}$. The dashed yellow line corresponds to the time that a wall-skimming particle takes, in a DNS of a model laminar Couette flow at equivalent Reynolds number, to reach the same inner-scaled wall-normal distance.
To confirm the nature of the mechanism for particle detachment, we performed an auxiliary DNS of laminar Couette flow at the same particle Reynolds number. The computational domain has size 
$L_{x}\times L_{y}\times L_{z}=20D\times 10D\times 10D$ with a regular grid where 
$D/\unicode[STIX]{x0394}x=16$. The boundary conditions are the same as for the turbulent channel flow, except that the flow is now driven by a non-zero streamwise velocity 
$U_{w}$ at 
$y=L_{y}$. The Reynolds number based on the local shear rate 
$\dot{\unicode[STIX]{x1D6FE}}=U_{w}/L_{y}$ and particle size is set to match that of the particle in the viscous sublayer, i.e. 
$\dot{\unicode[STIX]{x1D6FE}}D^{2}/\unicode[STIX]{x1D708}=(u_{\unicode[STIX]{x1D70F}}/\unicode[STIX]{x1D6FF}_{v})D^{2}/\unicode[STIX]{x1D708}$. A single particle with the same physical properties is placed at the bottom wall, with initial linear and angular velocity conforming to the local flow velocity and vorticity. The flight time 
$\unicode[STIX]{x0394}t^{up}$ for the particle to detach from the wall and travel five viscous units in 
$y$ is reported by the dashed line of figure 6. The measured time is remarkably close to the average value measured in the interface-resolved DNS for the two turbulent cases under consideration. This strongly suggests that the mechanism for particle detachment from the wall is, to first approximation, purely shear-driven. Moreover, the following scaling considerations suggest that shear-induced lift forces can be important close to the wall.
Let us neglect short-range particle–wall interactions and consider the particle dynamics in the viscous sublayer to be modelled by an unbounded linear shear flow. In the limit of vanishing particle Reynolds number, the ratio between the drag and shear-induced lift force scales with 
$D^{+}$ (cf. (2.9) and (2.10) with 
$J=1$ for Saffman lift, and note that in this case 
$|\unicode[STIX]{x1D74E}\times \boldsymbol{U}_{s}|=|\unicode[STIX]{x1D74E}||\boldsymbol{U}_{s}|=(u_{\unicode[STIX]{x1D70F}}^{2}/\unicode[STIX]{x1D708})|\boldsymbol{U}_{s}|$):
$$\begin{eqnarray}\left|\frac{\boldsymbol{F}_{l}}{\boldsymbol{F}_{d}}\right|_{y^{+}<5}=\frac{1.615}{3\unicode[STIX]{x03C0}}\sqrt{\frac{D^{2}|\unicode[STIX]{x1D74E}|}{\unicode[STIX]{x1D708}}}=0.171D^{+},\end{eqnarray}$$ with a proportionality coefficient 
$O(1)$ in the present set-up. This suggests that shear-induced lift forces cannot be neglected near the wall when 
$D^{+}\gtrsim 1$, as these can be as high as the streamwise drag forces. It should be noted, however, that the lift force might remain important for lower values of 
$D^{+}$, since particles in the viscous sublayer would tend to flow for a longer time while subjected to this wall-normal force. Away from the wall, instead, the order of magnitude of the lift force should become negligible with respect to that of the viscous drag, just like the other terms in the Maxey–Riley–Gatignol equations that are typically neglected for small inertial particles in wall-bounded turbulent flows (see e.g. Wang et al. (Reference Wang, Squires, Chen and McLaughlin1997)).
We have seen that the dynamics of resolved and point particles is similar in the bulk, and thus the drift towards the wall is well described by point-particle methods. When moving close to the wall, the high shear induces strong lift forces that quickly dislodge particles. The lower value of the near-wall volume fraction peak in figure 3(b) for case D, and its shorter average take-off time in figure 6 are consistent with this picture, as the mean wall shear is higher in this case. This wall detachment mechanism is absent in models considering only the drag, as in the PP case, which instead display slow resuspension and particles spending a long time in the near-wall low-speed streaks. In the next section we assess the validity of simple shear-induced lift models to predict the particle dynamics.
(a) Local solid mass fraction as a function of the outer-scaled wall-normal distance for the different point-particle cases. (b) Corresponding inner-scaled mean particle velocity profiles. The inset in panel (b) shows the inner-scaled difference between fluid and particle velocity profiles 
$\unicode[STIX]{x0394}\langle u\rangle$.
Outer-scaled second-order moments of particle velocity for the different point-particle cases: (a) streamwise velocity r.m.s., (b) wall-normal velocity r.m.s. and (c) ditto for spanwise velocity. (d) Reynolds stress profile.
3.2 Assessment of lift models for point-particle simulations
The most well-known models proposed to describe a shear-induced lift force in the one-way-coupling point-particle regime introduced at the end of § 2 are here compared with fully resolved simulations. In all cases, the standard Eulerian particle statistics for the point-particle models with lift force are compared to those of the model without lift, and the interface-resolved case VD.
Figure 7(a) shows the normalized local volume fraction profiles for the different cases considered. The point-particle models accounting for lift forces predict well the wall concentration peak. The peak location is predicted slightly away from the wall, probably due to a minor stabilizing effect of short-range particle–wall interactions and softer contact, absent in the point-particle model. These secondary mechanisms can possibly be modelled with near-wall closures such as an effective (wet) restitution coefficient that is a function of the particle impact Stokes number (see e.g. Joseph et al. (Reference Joseph, Zenit, Hunt and Rosenwinkel2001), Legendre et al. (Reference Legendre, Zenit, Daniel and Guiraud2006), Fong et al. (Reference Fong, Amili and Coletti2019)), or by more sophisticated near-wall corrections for the particle dynamics (Gondret, Lance & Petit Reference Gondret, Lance and Petit2002; Lee & Balachandar Reference Lee and Balachandar2010). We should note that the former option may be too simple for realistically predicting the particle dynamics, as the drag force already accounts for part of the fluid effects (Gondret et al. Reference Gondret, Lance and Petit2002). Indeed, we tested a point-particle simulation with the same dynamics as PP-Saffman but with a lower coefficient of restitution of 0.9; the results showed a shift towards the wall of the concentration peak that, however, overpredicts the maximum concentration by a factor of 1.3.
Away from the wall, the case using the seminal lift model, PP-Saffman, optimally predicts the concentration profile. The mean particle velocity pertaining to the cases where a lift force is accounted for is much closer to the interface-resolved case (see figure 7b). This is particularly evident in the apparent mean particle-to-fluid velocity difference displayed in the figure inset, where the results with lift force show a much better agreement with case VD, in particular in the region where the negative slip is highest. Overall, both lift models predict relatively well the first-order statistics pertaining to case VD shown in figure 7. This is not the case for the second-order moments of the particle velocity.
Figure 8 shows the second-order moments of particle velocity for the different cases considered. As seen in the previous section, the profiles pertaining to the streamwise and wall-normal velocity fluctuations (figure 8a,b,d) of the interface-resolved simulations, case VD, show significant differences with case PP close to the wall, with a peculiar non-monotonic trend as compared to the smooth decrease obtained from the point-particle simulations. All panels point to a somewhat surprising result: the low-Reynolds-number shear-induced lift model of Saffman (Reference Saffman1965) predicts almost perfectly the second-order moments of particle velocity, while the model accounting for finite-Reynolds-number effects shows major qualitative differences. Actually, the somewhat better agreement of the Saffman model compared to that of McLaughlin (Reference McLaughlin1991) for interface-resolved simulations of a particle in a shear flow has also been noticed in the work of Mei (Reference Mei1992). Still, we should stress here the better performance of Saffman’s model in this more demanding benchmark, in view of the number of references that have implemented/extended lift models based on the work of McLaughlin (Reference McLaughlin1991) for point-particle simulations. We also tested other approaches that have been suggested in the literature based on the work of McLaughlin (Reference McLaughlin1991) and observed similar discrepancies (not shown here).
The very good prediction of Saffman’s lift model also sheds light on the cause of the non-monotonic trend of the fluctuation intensities and shear stresses in figure 8(a,b,d) for 
$y/h<0.1$. In this region, the magnitude of the shear-induced lift force is highest. Hence, a particle approaching the wall is forced to suddenly depart from it, enhancing the velocity fluctuations in the wall-normal and streamwise directions. This mechanism, however, does not increase the spanwise velocity fluctuation.
4 Conclusions
We have presented two particle-resolved DNS of turbulent channel flow laden with small inertial particles with volume fractions of 
$3\times 10^{-4}$ and 
$3\times 10^{-5}$. Since the latter case hardly shows flow modulation, it can be considered to fall into the one-way-coupling assumption, i.e. the particles do not influence the statistics of the turbulence. The less dilute case, instead, shows approximately 10 % drag increase compared to the unladen case. This striking increase is attributed to a significant particle mass fraction near the wall caused by turbophoresis. These particles flow with high particle-to-fluid (apparent) slip velocity, producing hot-spots of large wall shear.
We also examine the differences between the most dilute case, which satisfies the one-way-coupling assumption, and a classic one-way-coupled point-particle simulation where the particle acceleration is balanced by the nonlinear (Schiller–Naumann) drag. While, in the bulk of the channel, concentration profiles and moments of the particle velocity agree well, most of these quantities show clear differences close to the wall. This disagreement is attributed to a missing key ingredient in this point-particle simulation: a shear-induced lift force for the near-wall particle dynamics. All differences between resolved and point-particle statistics can be related to this near-wall effect. First, the high concentration at the wall in the interface-resolved simulations is limited to a single particle layer, since particles tend to promptly depart from the wall. Second, particles are resuspended before they can accumulate for a long time in low-speed streaks, which lowers the apparent slip velocity in the near-wall region. Finally, because of this lift force, particles show higher near-wall velocity fluctuations.
To test this hypothesis, the average particle residence time in the viscous sublayer has been measured and compared to that of a single resolved sphere in a laminar Couette flow with matched (particle shear) Reynolds number. The average flight time of the resolved particles in a turbulent channel is extremely close to that of the laminar simulation, strongly supporting our conclusion. Conversely, point particles reside in the viscous sublayer for a time approximately one order of magnitude longer. Properly accounting for this reduced residence time is therefore fundamental to accurately predict the particle statistics without resolving the flow conforming to the particles. Moreover, a scaling analysis shows that the lift force in the viscous sublayer can be as high as the streamwise drag force: in the limit of vanishing particle Reynolds number, the ratio between the two forces in the viscous sublayer scales with the inner-scaled particle diameter, 
$D^{+}$, which is 
$O(1)$ in the cases addressed here.
Accordingly, we tested the validity of simple but widely used shear-induced lift models: that from the seminal work of Saffman (Reference Saffman1965), and a model based on the work of McLaughlin (Reference McLaughlin1991) that considers finite-Reynolds-number effects (see Mei Reference Mei1992). Our results show that both models predict reasonably well the wall peak in particle concentration and mean particle velocity, where the minor differences can be attributed to short-range particle–wall interactions that are not incorporated in the model. As regards the second-order moments of the particle velocity, the interface-resolved results from the simulations show a more complex, non-monotonic trend near the wall for the streamwise and wall-normal r.m.s. and the Reynolds shear stress in a region where the point-particle simulation without lift predicts a monotonic decrease. Quite remarkably, the oldest lift model available in the literature, that due to Saffman, quantitatively predicts all these trends, while the other model shows differences with respect to the most dilute interface-resolved simulation. We should remark that we have tested other approaches that have been suggested in the literature based on the work of McLaughlin (Reference McLaughlin1991), and observed similar discrepancies. The present work opens new perspectives towards a revision of the point-particle models. In particular, trajectories and force time-histories of the point particles could be compared to those of fully resolved particles in order to assess possible local discrepancies and potential improvements. Nevertheless, the good agreement between statistics of interface-resolved particle simulations and of the point-particle model including the Saffman lift is an indirect proof that the dominant dynamics is well captured. The success of the Saffman lift force model also suggests that short-range particle–wall interactions play a minor role in the near-wall particle dynamics when compared to that of the local shear near the wall.
We hope that the present results can be exploited for the development of improved point-particle models for one- and two-way-coupling regimes in wall-bounded turbulent flows.
Acknowledgements
We acknowledge PRACE for awarding us access to the supercomputer Marconi, based in Italy at CINECA under project 2017174185–DILPART, and the computing time provided by SNIC (Swedish National Infrastructure for Computing). This work was supported by the European Research Council grant no. ERC-2013-CoG-616186, TRITOS, the Swedish Research Council grant no. VR 2014-5001, and the grant BIRD192032/19 from the University of Padova. The authors acknowledge the anonymous referees for useful comments on an earlier version of the manuscript. B. Vreman, H. Kuerten and P. Pakseresht are also thanked for their feedback on an earlier version of the manuscript.
Declaration of interests
The authors report no conflict of interest.
Open access
