2.1 Turbulent channel flow

The fluid flow is described by velocity $\boldsymbol{u}(\boldsymbol{x},t)$ and pressure $p(\boldsymbol{x},t)$ fields that evolve according to the incompressible Navier–Stokes equations,

(2.1a,b)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}(p/\unicode[STIX]{x1D70C}_{f})+\unicode[STIX]{x1D708}_{f}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{f},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{f}$ is the fluid mass density, $\unicode[STIX]{x1D708}_{f}$ is the kinematic viscosity and $\boldsymbol{f}$ signifies forcing by particles. The effect of including or neglecting the two-way coupling represented by $\boldsymbol{f}$ is documented in appendix D. For the mass fractions simulated in this paper, the two-way coupling is of secondary importance for the particle concentration profile compared with particle–particle collisions and so is neglected in the numerical calculations presented in the main body of the paper. Although the presence or absence of the two-way coupling does not change the form of the conservation expressions derived below for the particle phase, it is expected that it would quantitatively affect the drag force terms in those expressions when the mass fraction is sufficiently large.

For studying the physics of wall-bounded flows in the simplest case possible, this paper considers a turbulent channel flow. An imposed pressure gradient in the $x$ (streamwise) direction, $\text{d}p/\text{d}x$, drives the flow and no-slip, no-penetration conditions are imposed at the (smooth) walls separated by a distance $2h$ in the $y$ (wall-normal) direction. Thus, the flow and particle statistics are homogeneous in the $x$ and $z$ (spanwise) directions in addition to the mirror symmetry about the centreline. A statistically stationary channel flow is characterized by the friction Reynolds number, $Re_{\ast }=u_{\ast }h/\unicode[STIX]{x1D708}_{f}$, where $u_{\ast }=\sqrt{-(\text{d}p/\text{d}x)h/\unicode[STIX]{x1D70C}_{f}}$ is the friction velocity. The smallest active scales in the flow are characterized by the friction length scale, $\unicode[STIX]{x1D6FF}_{\ast }=\unicode[STIX]{x1D708}_{f}/u_{\ast }$, and the friction time scale, $\unicode[STIX]{x1D70F}_{\ast }=\unicode[STIX]{x1D708}_{f}/u_{\ast }^{2}$.

2.2 Lagrangian particle tracking

The particle phase is described by $N_{p}$ discrete particles with centre of mass $\boldsymbol{x}(t)$ and velocity $\boldsymbol{v}(t)$, evolving according to the trajectory equations,

(2.2a,b)$$\begin{eqnarray}\displaystyle \dot{\boldsymbol{x}}^{(i)}=\boldsymbol{v}^{(i)},\quad \dot{\boldsymbol{v}}^{(i)}=\boldsymbol{a}^{(i)}(\boldsymbol{x},\boldsymbol{v}),\quad i=1,2,\ldots ,N_{p}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{a}(\boldsymbol{x},\boldsymbol{v})$ is the acceleration of the particle due to fluid forces. Each spherical particle is characterized by its diameter ($d_{p}$) and mass density ($\unicode[STIX]{x1D70C}_{p}$). For small particles ($d_{p}<\unicode[STIX]{x1D6FF}_{\ast }$) at low particle slip Reynolds numbers ($Re_{p}=|\boldsymbol{u}-\boldsymbol{v}|d_{p}/\unicode[STIX]{x1D708}_{f}\ll 1$), Stokes drag can be used, $\boldsymbol{a}_{St}=(\boldsymbol{u}(\boldsymbol{x})-\boldsymbol{v})/\unicode[STIX]{x1D70F}_{p}$, with relaxation time $\unicode[STIX]{x1D70F}_{p}=\unicode[STIX]{x1D70C}_{p}d_{p}^{2}/(18\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D708}_{f})$ in the limit $\unicode[STIX]{x1D70C}_{p}\gg \unicode[STIX]{x1D70C}_{f}$. In DNS, the fluid velocity seen by the particle, $\boldsymbol{u}(\boldsymbol{x})$, must be interpolated from the Eulerian grid to particle locations for the drag law. In LES, the interpolated velocity only represents the filtered velocity seen by the particle. The remaining SGS velocity seen by the particles requires additional modelling (Fede et al. Reference Fede, Simonin, Villedieu and Squires2006; Kuerten Reference Kuerten2006; Marchioli et al. Reference Marchioli, Salvetti and Soldati2008b; Ray & Collins Reference Ray and Collins2014; Minier Reference Minier2015, Reference Minier2016; Marchioli Reference Marchioli2017; Park et al. Reference Park, Bassenne, Urzay and Moin2017). The simulations in this paper use the Schiller–Naumann drag, $\boldsymbol{a}_{SN}=\boldsymbol{a}_{St}(1+0.15Re_{p}^{0.687})$, which gives a finite Reynolds number correction (Schiller & Naumann Reference Schiller and Naumann1933; Balachandar & Eaton Reference Balachandar and Eaton2010). Note that the present work does not include a lift force, which could potentially play a role in the near-wall region where velocity gradients are largest. Simple lift models, such as the Saffman lift force, are not applicable here. Given the complexity and variety of lift models (Wang et al. Reference Wang, Squires, Chen and McLaughlin1997; Marchioli, Picciotto & Soldati Reference Marchioli, Picciotto and Soldati2007), we defer the effect of lift to future work.

At times within § 3, the Stokes drag form will be used to simplify expressions. The same procedures as shown below can also be performed using Schiller–Naumann drag to include finite Reynolds number correction terms. These expressions are used to generate the plots in this paper, to be consistent with the simulation, and are given in appendix A.

3.1 Single-particle probability density function evolution

The particle ensemble has statistical homogeneity in the $x$ and $z$ directions, so the particle statistics of interest include only wall-normal position and velocity. To this end, we study the single-particle position–velocity probability density function (PDF) defined by,

(3.1)$$\begin{eqnarray}\displaystyle f(y,v_{y};t)=\langle \unicode[STIX]{x1D6FF}(y-{\hat{y}}(t))\unicode[STIX]{x1D6FF}(v_{y}-\hat{v}_{y}(t))\rangle , & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}(x)$ is the Dirac delta function representing the fine-grained PDF. The angled brackets, $\langle \cdot \rangle$, are employed to denote ensemble averaging over a monodisperse system of particles. Practically speaking, these averages may be obtained by averaging over all identical particles in a given realization, as well as averaging in time. Such temporal averaging is more convenient for obtaining converged statistics from numerical simulations and is therefore used in subsequent sections.

For use with (3.1), the dynamics of each individual particle, equation (2.2), is projected in the wall-normal direction,

(3.2a,b)$$\begin{eqnarray}\displaystyle \dot{{\hat{y}}}=\hat{v}_{y}\quad \dot{\hat{v}}_{y}=\hat{a}_{y}=\frac{u_{y}(\hat{\boldsymbol{x}},t)-\hat{v}_{y}}{\unicode[STIX]{x1D70F}_{p}}. & & \displaystyle\end{eqnarray}$$

Differentiating (3.1) in time and substituting (3.2), one can obtain an evolution equation for $f$,

(3.3)$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(v_{y}f)}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}(\langle a_{y}|y,v_{y}\rangle f)}{\unicode[STIX]{x2202}v_{y}}={\dot{f}}_{coll}, & & \displaystyle\end{eqnarray}$$

where the conditional average $\langle a_{y}|y,v_{y}\rangle$ is shorthand for $\langle \hat{a}_{y}|{\hat{y}}=y,\hat{v}_{y}=v_{y}\rangle$, i.e. the average particle acceleration conditioned on wall-normal position and velocity. The right-hand side term, ${\dot{f}}_{coll}$, denotes the changes in particle velocity due to inter-particle collisions. Note that Sikovsky (Reference Sikovsky2014) starts from an equation similar to (3.3), but in that work the fluid velocity seen by the particle (which contributes to the acceleration term) is modelled by a stochastic forcing term. In this paper, no stochastic modelling assumptions are made about the fluid flow seen by the particle. Instead, the present analysis seeks to understand the consequences of the exact conservation equations for the particle phase.

3.2 Particle mass conservation

The particle number density (concentration) is obtained at any $y$ location by integrating $f$ over all possible particle velocities,

(3.4)$$\begin{eqnarray}\displaystyle C(y;t)=C_{0}(t)\int _{-\infty }^{\infty }f(y,v_{y};t)\,\text{d}v_{y}, & & \displaystyle\end{eqnarray}$$

where $C_{0}(t)=\int C(y;t)\,\text{d}y$ is the bulk particle concentration. Integrating (3.3) over all velocities and multiplying by $C_{0}$, we obtain a mass conservation equation for the particle phase,

(3.5)$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(\langle v_{y}|y\rangle C)}{\unicode[STIX]{x2202}y}=0. & & \displaystyle\end{eqnarray}$$

The collisional term exactly vanishes because each collision individually conserves mass. Assuming statistical stationarity, $\unicode[STIX]{x2202}C/\unicode[STIX]{x2202}t=0$. It follows, then, that $\langle v_{y}|y\rangle =0$. This simply states that the net wall-normal flux of particles must vanish at all $y$ locations for the PDF to remain stationary in time.

3.3 Particle wall-normal momentum conservation

The particle wall-normal momentum, which is the same as the particle mass flux introduced in (3.5), is obtained as a first-order moment of $f$,

(3.6)$$\begin{eqnarray}\displaystyle \langle v_{y}|y\rangle C(y;t)=C_{0}\int _{-\infty }^{\infty }v_{y}f(y,v_{y};t)\,\text{d}v_{y}. & & \displaystyle\end{eqnarray}$$

A conservation equation for the wall-normal momentum is obtained by multiplying (3.3) by $C_{0}$ and $v_{y}$ and integrating over all particle velocities. The result is,

(3.7)$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}(\langle v_{y}|y\rangle C)}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(\langle v_{y}^{2}|y\rangle C)}{\unicode[STIX]{x2202}y}-\langle a_{y}|y\rangle C=0. & & \displaystyle\end{eqnarray}$$

As with the mass conservation equation, the collisional term makes no contribution to the wall-normal momentum of the particle ensemble because each inter-particle collision individually conserves momentum. The third term on the left, the acceleration term, is obtained using integration by parts. At steady state, $\unicode[STIX]{x2202}(\langle v_{y}|y\rangle C)/\unicode[STIX]{x2202}t=0$, and the particle wall-normal momentum balance reduces to,

(3.8)$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}y}(\langle v_{y}^{2}|y\rangle C)=\langle a_{y}|y\rangle C. & & \displaystyle\end{eqnarray}$$

Because $\langle v_{y}|y\rangle =0$ from mass conservation, $\langle v_{y}^{2}|y\rangle$ represents the particle wall-normal velocity variance as a function of distance from the wall. This would not be the case for a spatially developing flow such as a turbulent boundary layer, in which case $\langle v_{y}^{2}|y\rangle$ could be split into mean-squared and variance components. Note that the derivative in $y$ on the left-hand side is now a total derivative since time dependence has been removed.

3.3.1 Mechanisms affecting non-uniform concentration

Using the product rule and assuming the Stokes drag formula, $a_{y}=(u_{y}-v_{y})/\unicode[STIX]{x1D70F}_{p}$, the particle wall-normal momentum balance at steady state, equation (3.8), may be rewritten as,

(3.9)$$\begin{eqnarray}\displaystyle \langle v_{y}^{2}|y\rangle \frac{\text{d}C}{\text{d}y}=\left(\frac{\langle u_{y}|y\rangle }{\unicode[STIX]{x1D70F}_{p}}-\frac{\text{d}\langle v_{y}^{2}|y\rangle }{\text{d}y}\right)C. & & \displaystyle\end{eqnarray}$$

The two terms on the right-hand side, proportional to the local concentration level, must balance with the left-hand side term proportional to the concentration gradient. Thus, two mechanisms for generating non-uniform particle concentrations may be deduced from (3.9).

Figure 2. From a point-particle DNS snapshot of $Re_{\ast }=150$ (details in § 4): location of $St^{+}=8$, $\unicode[STIX]{x1D6F7}_{V}=3\times 10^{-5}$ particles with centroid $5\leqslant y^{+}\leqslant 10$, white spheres, for particles superposed with contours at $y^{+}=10$ of the (a) streamwise velocity, (b) wall-normal velocity. The size of each particle is exaggerated for visualization purposes. The velocities are in friction units, $\boldsymbol{u}/u_{\ast }$.

The first right-hand side term in (3.9) represents the average drag force on the particles at a given distance from the wall. For Stokes drag, this is proportional to, $\langle u_{y}|y\rangle$, the average wall-normal fluid velocity as sampled by particles at a given wall-normal location. In arriving at (3.9), it is recognized that the $\langle v_{y}|y\rangle$ contribution from the Stokes drag relation vanishes due to mass conservation, leaving only potential contributions from $\langle u_{y}|y\rangle$. For particles randomly distributed in the flow, this drag term vanishes. However, it is known that heavy particles tend to accumulate in low-speed streaks associated with ejection events in near-wall turbulence (Rashidi, Hetstroni & Banerjee Reference Rashidi, Hetstroni and Banerjee1990; Eaton & Fessler Reference Eaton and Fessler1994; Marchioli & Soldati Reference Marchioli and Soldati2002). Indeed, as shown in figure 2, the particles tend to prefer ‘ejection’ regions, which have lower streamwise velocity ($u_{x}^{\prime }<0$) and positive wall-normal velocity ($u_{y}^{\prime }>0$). As a result, the biased sampling term is positive ($\langle u_{y}|y\rangle >0$) near the wall as particles sample ejection events more often than sweeps. This means that the biased sampling provides a net force which pushes particles away from the wall.

The second term on the right-hand side is the turbophoresis pseudo-force, which creates the migration of particles down gradients in wall-normal velocity variance. The no-slip, no-penetration boundary conditions are enforced on the fluid phase at the wall, and the tendency of the particles to relax toward local fluid velocities will give lower $\langle v_{y}^{2}|y\rangle$ in the immediate vicinity of the wall. This leads to $\text{d}\langle v_{y}^{2}|y\rangle /\text{d}y>0$, i.e. increasing particle wall-normal velocity fluctuations with increasing distance from the wall. As a result, turbophoresis causes a significant net migration of particles toward the wall which is only partially offset by the biased sampling force.

It is noteworthy that the analysis of Guha (Reference Guha1997, Reference Guha2008) has omitted the biased sampling term, but for low Stokes numbers, this term significantly affects the concentration profile. The difference between Eulerian-averaged fluid velocities and average fluid velocities seen by (non-randomly spaced) particles is included as a drift velocity in stochastic models (Minier, Chibbaro & Pope Reference Minier, Chibbaro and Pope2014). Near the wall, this drift velocity is driven by interaction with coherent structures, leading to the biased sampling effect of a net force away from the wall. The importance of this term will be further illustrated when considering the $\unicode[STIX]{x1D70F}_{p}\rightarrow 0$ limit below. At higher Stokes numbers, the $\unicode[STIX]{x1D70F}_{p}$ in the denominator shows that the biased sampling force will diminish and likely become negligible compared to turbophoresis. The prediction of the turbophoresis force may not require detailed knowledge of particle interactions with turbulent coherent structures since it only relies on the gradient of particle wall-normal velocity variance. In contrast, turbulent coherent structures play a direct and unmistakable role in establishing the sampling bias term. For this reason, one objective of the numerical demonstrations below is to quantitatively assess at what Stokes number one may safely neglect the effect of the biased sampling term in (3.9) on the concentration profile, and hence potentially neglect the details of near-wall coherent structures. While turbophoresis and biased sampling are two distinct effects in the statistical equations, it should be appreciated that both naturally emerge from a true account of the particle dynamics in the presence of coherent structures (i.e. the fluid velocity seen by the particles), so they are not uncorrelated in a dynamical sense.

3.3.2 Formal solution

It is straightforward to construct the formal solution of (3.9),

(3.10)$$\begin{eqnarray}\displaystyle C(y)={\mathcal{N}}\exp \left(\underbrace{\frac{1}{\unicode[STIX]{x1D70F}_{p}}\int ^{y}\frac{\langle u_{y}|\unicode[STIX]{x1D702}\rangle }{\langle v_{y}^{2}|\unicode[STIX]{x1D702}\rangle }\,\text{d}\unicode[STIX]{x1D702}}_{\mathit{biased\,sampling}}-\underbrace{\int ^{y}\frac{\text{d}\ln \langle v_{y}^{2}|\unicode[STIX]{x1D702}\rangle }{\text{d}\unicode[STIX]{x1D702}}\,\text{d}\unicode[STIX]{x1D702}}_{\mathit{turbophoresis}}\right), & & \displaystyle\end{eqnarray}$$

where ${\mathcal{N}}$ is an integration constant which is unimportant for exploring causes of non-uniform concentration. Two phoresis integrals can be defined based on this form, as underscored in (3.10). The first phoresis integral,

(3.11)$$\begin{eqnarray}\displaystyle I_{bias}=\frac{1}{\unicode[STIX]{x1D70F}_{p}}\int _{0}^{y}\frac{\langle u_{y}|\unicode[STIX]{x1D702}\rangle }{\langle v_{y}^{2}|\unicode[STIX]{x1D702}\rangle }\,\text{d}\unicode[STIX]{x1D702}, & & \displaystyle\end{eqnarray}$$

quantifies the effect of biased sampling on the resulting concentration profile. For positive $\langle u_{y}|y\rangle$, this term alone would cause increasing concentration with increasing distance to the wall. Meanwhile, the second phoresis integral,

(3.12)$$\begin{eqnarray}\displaystyle I_{turb}=-\int _{0}^{y}\frac{\text{d}\ln \langle v_{y}^{2}|\unicode[STIX]{x1D702}\rangle }{\text{d}\unicode[STIX]{x1D702}}\,\text{d}\unicode[STIX]{x1D702}, & & \displaystyle\end{eqnarray}$$

quantifies the influence of turbophoresis on the final concentration profile. In these two definitions, the wall location ($y=0$) is used for the lower bound. These two definitions enable a straightforward, fair comparison of how biased sampling and turbophoresis affect the concentration profile and when one effect may be safely neglected. A relationship similar to (3.10) was considered by Capecelatro, Desjardins & Fox (Reference Capecelatro, Desjardins and Fox2016) in the context of two-fluid equations, but only after the biased sampling term was found to be negligible and removed from the analysis.

While (3.10) provides a formal expression for the concentration profile, it should be appreciated that it involves unclosed statistical terms depending, most notably, on the statistics of fluid velocities seen by particles at various wall-normal distances. At present, this expression is not necessarily advanced as a specific starting point for model development, but is rather used here to demonstrate qualitative features and provide a framework for elucidating the performance of models for the fluid velocities seen by particles.

3.3.3 The $St\rightarrow \infty$ limit

Seeing that the turbophoresis integral, equation (3.12), may be formally integrated, then (3.10) may also be written as

(3.13)$$\begin{eqnarray}\displaystyle C(y)=\frac{{\mathcal{N}}}{\langle v_{y}^{2}|y\rangle }\exp \left(\frac{1}{\unicode[STIX]{x1D70F}_{p}}\int ^{y}\frac{\langle u_{y}|\unicode[STIX]{x1D702}\rangle }{\langle v_{y}^{2}|\unicode[STIX]{x1D702}\rangle }\,\text{d}\unicode[STIX]{x1D702}\right). & & \displaystyle\end{eqnarray}$$

This form emphasizes that, at large Stokes numbers ($\unicode[STIX]{x1D70F}_{p}^{-1}\rightarrow 0$) when the biased sampling becomes negligible, the concentration profile becomes inversely proportional to the wall-normal particle velocity variance. From a modelling perspective, this means that predicting concentration profiles for high Stokes number particles,

(3.14)$$\begin{eqnarray}\displaystyle C(y)\approx \frac{{\mathcal{N}}}{\langle v_{y}^{2}|y\rangle }, & & \displaystyle\end{eqnarray}$$

requires only an accurate representation of particle wall-normal velocity fluctuation magnitudes. The spatio-temporal details of interactions between particles and near-wall coherent structures becomes less important in this limit, simplifying the modelling task. In a WMLES simulation, it is possible to reproduce fairly accurate wall-normal velocity variances for the fluid (Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2018) in the resolved region of the flow. However, it is less clear whether the variance profiles for particles, including Lagrangian sampling and memory effects, below the first grid point can be reproduced in the WMLES modelling approach.

3.3.4 The $St\rightarrow 0$ limit

In the opposite limit $\unicode[STIX]{x1D70F}_{p}\rightarrow 0$, the particle velocity may be written as (Maxey Reference Maxey1987)

(3.15)$$\begin{eqnarray}\displaystyle v_{y}=u_{y}-\unicode[STIX]{x1D70F}_{p}\frac{\text{D}u_{y}}{\text{D}t}+O(\unicode[STIX]{x1D70F}_{p}^{2}). & & \displaystyle\end{eqnarray}$$

This means that the biased sampling force becomes

(3.16)$$\begin{eqnarray}\displaystyle \lim _{\unicode[STIX]{x1D70F}_{p}\rightarrow 0}\frac{\langle u_{y}|y\rangle }{\unicode[STIX]{x1D70F}_{p}}=\left\langle \frac{\text{D}u_{y}}{\text{D}t}\right\rangle =-\frac{1}{\unicode[STIX]{x1D70C}_{f}}\frac{\unicode[STIX]{x2202}\langle p\rangle }{\unicode[STIX]{x2202}y}. & & \displaystyle\end{eqnarray}$$

Note that the viscous force, $\unicode[STIX]{x1D708}_{f}\unicode[STIX]{x1D6FB}^{2}\langle u_{y}\rangle$, vanishes because $\langle u_{y}\rangle =0$ at steady state in the limit $\unicode[STIX]{x1D70F}_{p}\rightarrow 0$. At $\unicode[STIX]{x1D70F}_{p}=0$, the average over the particle ensemble becomes the same as the fluid ensemble average, the concentration profile becomes flat (incompressibility) and the wall-normal RANS equation is obtained from (3.9),

(3.17)$$\begin{eqnarray}\displaystyle 0=-\frac{\unicode[STIX]{x2202}\langle p\rangle }{\unicode[STIX]{x2202}y}+\frac{\text{d}\langle u_{y}^{2}\rangle }{\text{d}y}, & & \displaystyle\end{eqnarray}$$

as in, for example, equation (5.2.2) of Tennekes & Lumley (Reference Tennekes and Lumley1972). This convergence of the $\unicode[STIX]{x1D70F}_{p}\rightarrow 0$ limit of the biased sampling to the pressure gradient was previously pointed out in passing by Bragg & Collins (Reference Bragg and Collins2014) during a discussion of the analogy between preferential concentration in homogeneous turbulence and turbophoresis in wall-bounded turbulence. The approach to (3.17) illustrates that the sampling bias becomes equal and opposite to turbophoresis in small $St$ number limit, which highlights the importance of preferential sampling of near-wall structures at low Stokes numbers. This analysis confirms that accurately predicting the concentration profile of low Stokes number particles will require a much more detailed representation of the near-wall turbulent fluctuations than is present in WMLES. The numerical results in §§ 4 and 5 explore both low and high $St^{+}$ cases in more quantitative detail.

3.3.5 The $y\rightarrow 0$ limit

The above analysis may also be leveraged to demonstrate the existence of a power-law concentration profile in the viscous sublayer. The derivation stems from a low Stokes number expansion, but the reasons a power law can also be observed at higher Stokes numbers are discussed below. Consider (3.15) in the limit $y\rightarrow 0$, where a Taylor expansion of the fluid velocity field prevails (i.e. in the viscous sublayer),

(3.18)$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}u_{x}(x,y,z,t)=A_{x}(x,z,t)y+O(y^{2}),\\ u_{y}(x,y,z,t)=B_{y}(x,z,t)y^{2}+O(y^{3}),\\ u_{z}(x,y,z,t)=A_{z}(x,z,t)y+O(y^{2}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Note that the linear term ${\sim}y$ for the wall-normal fluid velocity is exactly zero due to the divergence-free condition for incompressible flows (Kim, Moin & Moser Reference Kim, Moin and Moser1987; Pope Reference Pope2000). Substituting (3.18) into (3.15) for particles ‘almost’ following fluid trajectories, i.e.

(3.19)$$\begin{eqnarray}\displaystyle v_{y}=u_{y}-\unicode[STIX]{x1D70F}_{p}\left(\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}t}+u_{x}\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}x}+u_{y}\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}y}+u_{z}\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}z}\right)+O(\unicode[STIX]{x1D70F}_{p}^{2}), & & \displaystyle\end{eqnarray}$$

one obtains an expansion for the particle wall-normal velocities,

(3.20)$$\begin{eqnarray}\displaystyle v_{y}=A_{y}y^{2}-\unicode[STIX]{x1D70F}_{p}\left(y^{2}\frac{\unicode[STIX]{x2202}A_{y}}{\unicode[STIX]{x2202}t}+y^{3}A_{x}\frac{\unicode[STIX]{x2202}B_{y}}{\unicode[STIX]{x2202}x}+2y^{3}B_{y}^{2}+y^{3}A_{z}\frac{\unicode[STIX]{x2202}B_{y}}{\unicode[STIX]{x2202}z}\right)+O(\unicode[STIX]{x1D70F}_{p}^{2},y^{3}). & & \displaystyle\end{eqnarray}$$

From this expression, the particle wall-normal velocity variance can be assessed,

(3.21)$$\begin{eqnarray}\displaystyle \langle v_{y}^{2}|y\rangle =\underbrace{\langle A_{y}^{2}\rangle }_{\unicode[STIX]{x1D6FC}}y^{4}+O(\unicode[STIX]{x1D70F}_{p},y^{5}), & & \displaystyle\end{eqnarray}$$

where statistically stationary flow has been assumed, i.e. $\langle \unicode[STIX]{x2202}_{t}A_{y}\rangle =0$ and $\langle A_{y}\unicode[STIX]{x2202}_{t}A_{y}\rangle =0$. The biased sampling term can also be evaluated,

(3.22)$$\begin{eqnarray}\displaystyle \frac{\langle u_{y}|y\rangle }{\unicode[STIX]{x1D70F}_{p}}=\frac{\langle (u_{y}-v_{y})|y\rangle }{\unicode[STIX]{x1D70F}_{p}}=\underbrace{\left\langle A_{x}\frac{\unicode[STIX]{x2202}B_{y}}{\unicode[STIX]{x2202}x}+2B_{y}^{2}+A_{z}\frac{\unicode[STIX]{x2202}B_{y}}{\unicode[STIX]{x2202}z}\right\rangle }_{\unicode[STIX]{x1D6FD}}y^{3}+O(\unicode[STIX]{x1D70F}_{p},y^{4}). & & \displaystyle\end{eqnarray}$$

Together, equations (3.21) and (3.22) yield

(3.23)$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\langle v_{y}^{2}|y\rangle \approx \unicode[STIX]{x1D6FC}y^{4},\\ \langle u_{y}|y\rangle \approx \unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70F}_{p}y^{3}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Although this has been derived specifically in the $\unicode[STIX]{x1D70F}_{p}\rightarrow 0$ limit, it is reasonable to expect that the scalings in (3.23) hold over a decent range of finite Stokes numbers because the particle ensemble averages are biased toward particles with larger residence time in the near-wall region (trapped particles). These particles with longer residence times by definition have more time to adjust to the local near-wall fluid scalings. In fact, the numerical results below show that scalings for the particle velocity variance and biased sampling terms in (3.23) are observed over a generous range of Stokes numbers.

Substituting (3.23) into (3.13),

(3.24)$$\begin{eqnarray}\displaystyle C(y)=\frac{{\mathcal{N}}}{\unicode[STIX]{x1D6FC}y^{4}}\exp \left(\int ^{y}\frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D702}^{3}}{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D702}^{4}}\,\text{d}\unicode[STIX]{x1D702}\right)=\frac{{\mathcal{N}}}{\unicode[STIX]{x1D6FC}}y^{\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC}-4}. & & \displaystyle\end{eqnarray}$$

At $\unicode[STIX]{x1D70F}_{p}=0$, $\unicode[STIX]{x1D6FD}=4\unicode[STIX]{x1D6FC}$ because the incompressibility of fluid particles implies a uniform concentration profile with biased sampling and turbophoresis balanced according to (3.17). For $\langle u_{y}|y\rangle$ to remain bounded as $\unicode[STIX]{x1D70F}_{p}$ becomes large, it must be that $\unicode[STIX]{x1D6FD}$ becomes small as Stokes number increases. (Note however that the expansion (3.20) will not apply at sufficiently high Stokes number.) At finite Stokes numbers, provided the near-wall scalings given by (3.23) hold, $0<\unicode[STIX]{x1D6FD}<4\unicode[STIX]{x1D6FC}$, then the concentration profile in the viscous sublayer has the form

(3.25)$$\begin{eqnarray}\displaystyle C(y)\sim y^{-\unicode[STIX]{x1D6FE}}, & & \displaystyle\end{eqnarray}$$

with $0<\unicode[STIX]{x1D6FE}<4$. Thus, the conservation of momentum in the wall-normal direction gives a simple understanding for near-wall power laws in particle concentration profiles, with power-law exponent bounded by $\unicode[STIX]{x1D6FE}<4$, provided the biased sampling coefficient $\unicode[STIX]{x1D6FD}$ behaves monotonically with $\unicode[STIX]{x1D70F}_{p}$ (this is confirmed in § 4). This bound is in agreement with previous stochastic models and observed DNS trends (Sikovsky Reference Sikovsky2014). As volume fraction increases, inter-particle collisions can interrupt this power-law behaviour by energizing near-wall particles and causing deviations from the $\langle v_{y}^{2}|y\rangle$ scaling in (3.23).

3.4 Higher-order moments

Higher-order moments,

(3.26)$$\begin{eqnarray}\displaystyle \langle v_{y}^{n}|y\rangle C(y;t)=C_{0}\int _{-\infty }^{\infty }v_{y}^{n}f(y,v_{y};t)\,\text{d}v_{y}, & & \displaystyle\end{eqnarray}$$

can be analysed following the same procedure. The right-hand side term in (3.3) describing the effects of particle–particle collisions no longer vanishes exactly in the evolution equation for higher-order moments. For simplicity, then, only the zero volume fraction limit will be considered here. After explicitly neglecting collisional effects, the $n$th moment of (3.3) gives

(3.27)$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}(\langle v_{y}^{n}|y\rangle C)}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(\langle v_{y}^{n+1}|y\rangle C)}{\unicode[STIX]{x2202}y}-n\langle a_{y}v_{y}^{n-1}|y\rangle C=0. & & \displaystyle\end{eqnarray}$$

At steady state for particles experiencing Stokes drag,

(3.28)$$\begin{eqnarray}\displaystyle \langle v_{y}^{n+1}|y\rangle \frac{\text{d}C}{\text{d}y}=\left(\frac{n}{\unicode[STIX]{x1D70F}_{p}}\langle (u_{y}-v_{y})v_{y}^{n-1}|y\rangle -\frac{\text{d}\langle v_{y}^{n+1}|y\rangle }{\text{d}y}\right)C, & & \displaystyle\end{eqnarray}$$

which has the formal solution

(3.29)$$\begin{eqnarray}\displaystyle C(y)=\frac{{\mathcal{N}}^{\prime }}{\langle v_{y}^{n+1}|y\rangle }\exp \left(\frac{n}{\unicode[STIX]{x1D70F}_{p}}\int ^{y}\frac{\langle (u_{y}-v_{y})v_{y}^{n-1}|\unicode[STIX]{x1D702}\rangle }{\langle v_{y}^{n+1}|\unicode[STIX]{x1D702}\rangle }\,\text{d}\unicode[STIX]{x1D702}\right). & & \displaystyle\end{eqnarray}$$

Dividing this expression by (3.13) raised to the $(n+1)/2$ power and rearranging yields,

(3.30)$$\begin{eqnarray}\displaystyle \frac{\langle v_{y}^{m}\rangle }{\langle v_{y}^{2}\rangle ^{m/2}}={\mathcal{N}}^{\prime \prime }C(y)^{m/2-1}\exp \left(\frac{m-1}{\unicode[STIX]{x1D70F}_{p}}\int ^{y}\frac{\langle (u_{y}-v_{y})v_{y}^{m-2}|\unicode[STIX]{x1D702}\rangle }{\langle v_{y}^{m}|\unicode[STIX]{x1D702}\rangle }\,\text{d}\unicode[STIX]{x1D702}-\frac{m}{2\unicode[STIX]{x1D70F}_{p}}\int ^{y}\frac{\langle u_{y}|\unicode[STIX]{x1D702}\rangle }{\langle v_{y}^{2}|\unicode[STIX]{x1D702}\rangle }\,\text{d}\unicode[STIX]{x1D702}\right), & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $m=n+1$ has been substituted.

Therefore, when the concentration profile has a power law given by (3.25), and similar scaling arguments to those in § 3.3.5 hold, this expression gives a power law for the skewness, flatness, and other higher-order hyper-flatness values. In particular,

(3.31)$$\begin{eqnarray}\displaystyle \frac{\langle v_{y}^{m}\rangle }{\langle v_{y}^{2}\rangle ^{m/2}}\sim y^{-\unicode[STIX]{x1D6FE}(m/2-1)+\unicode[STIX]{x1D6FF}}, & & \displaystyle\end{eqnarray}$$

where the $\unicode[STIX]{x1D6FF}$ comes from the exponential term in (3.30) using (3.18) and (3.20) and following similar steps as before. Note that as $\unicode[STIX]{x1D70F}_{p}$ increases, this exponential term becomes small, leading to $\unicode[STIX]{x1D6FF}\rightarrow 0$ as $St$ becomes large. This limit is in agreement with the stochastic model of Sikovsky (Reference Sikovsky2014), namely, that flatness and hyper-flatness profiles have a power-law form near the wall with exponent $\unicode[STIX]{x1D6FE}(m/2-1)$. That model, however, apparently misses the correction term ($\unicode[STIX]{x1D6FF}$ in (3.31)) coming from the exponential term in (3.30), which is significant for smaller Stokes numbers. This is further apparent in the fact that their results for flatness, $m=4$, give better agreement at $St^{+}=25$ than at $St^{+}=5$.

4.1 Numerical methods and simulation details

The continuum equations for the fluid, equation (2.1), are discretized on a staggered Cartesian grid, and second-order central differencing is employed. Trilinear interpolation is used to compute the flow quantities (e.g. velocity) at particle locations for the drag law, equation (2.2). While previous studies have recommended higher-order interpolation (Yeung & Pope Reference Yeung and Pope1988), we find that trilinear interpolation is sufficient for studying the concentration profile in this work, see appendix F. The Schiller–Naumann form, $\boldsymbol{a}_{SN}$, is used for particle drag. In the present study, we consider the limit that $g\unicode[STIX]{x1D70F}_{p}/u_{\ast }$ is small, and thus neglect gravitational forces on the particles. The presence of a strong gravity force could alter the dynamics of the particles in a way that would depend on the orientation of the gravitational field with respect to the channel geometry. For instance, Lavezzo et al. (Reference Lavezzo, Soldati, Gerashchenko, Warhaft and Collins2010) found that acceleration statistics are significantly altered by gravity in the wall-normal direction.

A fractional step method for time advancement for the fluid and particles is done with Huen’s second-order method (RK2). Particle-particle collisions are computed using a hard-sphere collision model with a specified restitution coefficient. For this treatment, an efficient algorithm for detecting binary collisions using their collision cylinders within a given time step is implemented. The collision outcome is computed based on angle of incidence and the restitution coefficient. Unless otherwise given, a restitution coefficient of $1$ is used (kinetic energy preserving collisions). Particle–wall collisions are similarly treated, with a unity restitution coefficient so that they conserve kinetic energy. The walls are considered smooth for both the continuum carrier fluid and for the particle–wall collisions, although it should be noted that wall roughness can play an important role (Sommerfeld Reference Sommerfeld1992; Kussin & Sommerfeld Reference Kussin and Sommerfeld2002; Benson, Tanaka & Eaton Reference Benson, Tanaka and Eaton2005; Vreman Reference Vreman2007; Konan, Simonin & Squires Reference Konan, Simonin and Squires2011; Milici et al. Reference Milici, De Marchis, Sardina and Napoli2014). Future work could extend this present study to consider rough walls or other effects such as adhesion.

The computational domain for the turbulent channel flow is periodic in $x$ and $z$ with domain size $L_{x}=4\unicode[STIX]{x03C0}h$, $L_{y}=2h$, and $L_{z}=2\unicode[STIX]{x03C0}h$. The mean pressure gradient is imposed by a uniform body force to match the specified friction Reynolds number of $Re_{\ast }=u_{\ast }h/\unicode[STIX]{x1D708}$, where $u_{\ast }^{2}=-h(\text{d}p/\text{d}x)/\unicode[STIX]{x1D708}_{f}$ is the friction velocity. The results in this section are shown for $Re_{\ast }=150$, though appendix C shows equivalent results for $Re_{\ast }=300$ and $600$; the conclusions of this section hold for higher Reynolds numbers as well. The relative influence of particle inertia is represented by the friction Stokes number, $St^{+}=\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{\ast }=u_{\ast }^{2}\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D708}$. When particle–particle collisions are considered, the other dimensionless parameter varied is the volume fraction, $\unicode[STIX]{x1D6F7}_{V}=\unicode[STIX]{x03C0}d_{p}^{3}N_{p}/(6L_{x}L_{y}L_{z})$, where $N_{p}$ is the total number of particles in the domain. The diameter of the particles is held constant at $d_{p}^{+}=0.5$ while the density ratio is changed to vary $St^{+}$. For some of the higher volume fractions shown in this section, the mass fraction is significant although two-way coupling effects are ignored. Appendix D explores two-way coupling effects and why neglecting them is justifiable in the present context. Further, appendix E explores the effect of restitution coefficient for particle–particle collisions. A restitution coefficient of $e=1.0$ is used for the results in this section.

The (DNS) grid resolution used in the homogeneous directions is $\unicode[STIX]{x0394}x^{+}\approx 11$, $\unicode[STIX]{x0394}z^{+}\approx 7$. The grid is stretched in the wall-normal direction using a hyperbolic tangent to yield $\unicode[STIX]{x0394}y_{min}^{+}\approx 0.5$ for the first grid point at the wall (wall-parallel velocities at $y^{+}\approx 0.25$) and $\unicode[STIX]{x0394}y_{max}^{+}\approx 7$ in the centre of the channel. The resulting number of grid points was $172\times 86\times 128$ in the streamwise, wall-normal and spanwise directions, respectively. Sensitivity of the particle concentration profile to further refinement was explored and found to be small, so this grid resolution may be considered sufficient for the present purposes while keeping computational costs low. The particles are initialized with a uniform random distribution, and the simulation proceeds until the particles obtain a stationary distribution before statistics are computed. Aside from the cases shown in this paper, we also verified that the present numerical approach produced results consistent with the benchmark results from Marchioli et al. (Reference Marchioli, Soldati, Kuerten, Arcen, Tanière, Goldensoph, Squires, Cargnelutti and Portela2008c). Appendix F briefly explores the impact of grid resolution and interpolation scheme on the results of this section.

4.2 Simulation results without inter-particle collisions

Figure 3 shows the main results for the statistics of particle ensembles without inter-particle collisions at $Re_{\ast }=150$ for a range of $0\leqslant St^{+}\leqslant 512$, using $d_{p}^{+}=0.5$ and varying the particle density $144\leqslant \unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}\leqslant 36\,864$ to change $St^{+}$. The results from these particle ensembles represent the limit of $\unicode[STIX]{x1D6F7}_{V}\rightarrow 0$ since collisions are neglected. As can be seen from figure 3(a), in the absence of inter-particle collisions, the concentration near the wall can reach values hundreds of times larger than the mean concentration level. The concentration profiles are computed on a uniform grid with spacing of $0.5\unicode[STIX]{x1D6FF}_{\ast }$. Note that the $St^{+}=0$ tracer particles have a flat distribution with slight interpolation and time discretization errors leading to an almost imperceptibly reduced concentration in the viscous sublayer.

Figure 3. Particles statistics from turbulent channel flow at $Re_{\ast }=150$ in the zero volume fraction limit (no inter-particle collisions) at various Stokes numbers: (a) concentration profiles; (b) root-mean-square wall-normal particle velocity, where the dashed grey curve indicates $v_{y}\sim y^{2}$ asymptotic behaviour near the wall; (c) sampling bias $\langle u_{y}|y\rangle$ for the particle ensembles with dashed grey line indicating ${\sim}y^{3}$ behaviour near the wall; (d) sampling bias and turbophoresis integrals, see (3.11) and (3.12).

The particle root-mean-square velocity fluctuations, $v_{y,rms}=\sqrt{\langle v_{y}^{2}|y\rangle }$, shown in figure 3(b) reveal that the asymptotic $v_{y}\sim y^{2}$ behaviour near the wall persists even up to relatively high $St^{+}$. Meanwhile, figure 3(c) shows the sampling bias term $\langle u_{y}|y\rangle$, demonstrating that the ${\sim}y^{3}$ behaviour of that term also extends to relatively large values of $St^{+}$. Taken together, these two figures justify the scaling behaviour, equation (3.23), used in § 3.3.5 well outside of the $St^{+}\ll 1$ range. Further investigation into this observation (not shown) elucidates that this happens because particle statistics near the wall are dominated by trapped particles with long enough residence times to adjust to viscous sublayer fluid fluctuations despite nominally large $St^{+}$.

The resulting phoresis integrals, equations (3.11) and (3.12), are shown in figure 3(d). The integrals are both positive in all cases and are computed on a uniform grid with spacing $0.5\unicode[STIX]{x1D6FF}_{\ast }$ using a trapezoidal rule integration. The exponential of their difference is the concentration profile, as verified in figure 4. At $St^{+}=2$, the sampling bias term mostly cancels the turbophoresis term. As the Stokes number is increased, however, the sampling bias term decreases sharply, leading to the more extreme near-wall enhancements of the concentration profile seen in figure 3. At a large enough Stokes number, say $St^{+}\geqslant 128$, the sampling bias integral is negligible compared to the turbophoresis integral, showing that the concentration profile is simply inversely proportional to the particle velocity variance as in (3.14). It should be noted, however, that for these high Stokes numbers, the particle wall-normal velocity variance deviates significantly from that of the fluid at any given distance from the wall. Therefore, it is still non-trivial that a WMLES designed to reproduce fluid velocity variances in the near-wall region would necessarily produce accurate particle velocity variances and concentration profiles, even in the high Stokes number limit.

Figure 4. For $\unicode[STIX]{x1D6F7}_{V}=0$, the comparison of left (continuous lines) and right sides (dashed lines) of (a) (3.10) and (b) (3.14). Colours and symbols indicate $St^{+}$ following the legend in figure 3.

For completeness, equation (3.10) for the concentration profile is directly verified by comparing left and right sides in figure 4(a). In figure 4(b), the same comparison is done using (3.14) instead, i.e. neglecting the sampling bias effect. This further emphasizes that for $St^{+}\geqslant 128$, the details of the interactions between particles and turbulent structures near the wall are not as important. Instead, the concentration profile can be predicted with good accuracy simply by the inverse of the particle wall-normal velocity variance.

4.3 Simulation results with inter-particle collisions

In the near-wall region, the combination of higher local volume fractions due to turbophoresis with sharp flow gradients enhances the role of inter-particle collisions even at relatively modest bulk volume fractions (Kuerten & Vreman Reference Kuerten and Vreman2015). Inter-particle collisions, therefore, become non-negligible at much lower bulk volume fractions in wall-bounded flows compared to what one might expect using traditional heuristics which do not consider the near-wall region (Elghobashi Reference Elghobashi1994). The main effect of the collisions is to reduce the near-wall concentration levels by flattening the wall-normal fluctuation profiles (Li et al. Reference Li, McLaughlin, Kontomaris and Portela2001; Yamamoto et al. Reference Yamamoto, Potthoff, Tanaka, Kajishima and Tsuji2001). Caraman, Borée & Simonin (Reference Caraman, Borée and Simonin2003) emphasized that particle–particle collisions enhance the wall-normal transport of particles by transfer of streamwise fluctuation energy into wall-normal fluctuations. Because of turbophoresis and enhanced relative velocities, most particle–particle collisions happen very close to the wall (Kuerten & Vreman Reference Kuerten and Vreman2016). As volume fraction increases, inter-particle collisions tend to bring the concentration profile toward a uniform profile, attenuating the effects of turbophoresis. Point-particle methods, therefore, must account for inter-particle collisions to produce accurate concentration profiles.

Figure 5. Statistics for $St^{+}=32$ particles in a turbulent channel flow at $Re_{\ast }=150$ at various volume fractions. Descriptions of panels are the same as in figure 3.

Figure 5 shows the same results as figure 3, but for particle ensembles including inter-particle collisions. For brevity, only $St^{+}=32$ is shown ($\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}=2304$), although the main observations made here apply to other Stokes numbers as well. The most striking observation to be made from figure 5(a) is that the enhanced concentration near the wall due to turbophoresis is suppressed by the effect of inter-particle collisions even at seemingly innocuous bulk volume fractions. With volume fraction increased to $\unicode[STIX]{x1D6F7}_{V}=10^{-4}$, the peak concentration profile at the wall is less than five times the bulk concentration. At the highest volume fraction, the mass fraction is certainly appreciable, but it is shown in appendix D that the two-way coupling force does not strongly impact the concentration profile compared to the particle–particle collisions for the volume fractions considered here. Other quantities such as turbulence modulation or particle acceleration statistics may be more sensitive to the two-way coupling force, but we do not consider these at the present. Collisional effects on the concentration profile are seen even for volume fractions as low as $\unicode[STIX]{x1D6F7}_{V}=10^{-6}$.

The main cause of this change is demonstrated in figure 5(b), where the particle wall-normal fluctuations near the wall increase dramatically above the fluid fluctuation levels as volume fraction is increased. The increased levels of fluctuation can be attributed to a more ballistic behaviour of particles as they collide more frequently and redistribute streamwise fluctuations into wall-normal fluctuations. This also impacts the sampling bias, as shown in figure 5(c). However, in terms of the resulting concentration profile, the dominant effect of collisions is the decrease in the turbophoresis integral, figure 5(d), leading to the decreased near wall (relative) concentration. The increased fluctuations near the wall, far in excess of local fluid fluctuation levels, break the scaling behaviours of (3.23) which are responsible for the near-wall power law in concentration. In fact, figure 5 shows that the concentration profiles for higher volume fraction cases no longer display convincing power laws near the wall. It is also apparent from figure 5(d) that the biased sampling effect is also significantly attenuated by particle–particle collisions as volume fraction increases. This may be attributed to the enhanced collision rates limiting the tendency of particles to preferentially concentrate in low-speed streaks (figure 2).

The verification of (3.10) is included in figure 6(a) for the cases with inter-particle collisions. The conservation of momentum by particle collisions ensures that the ${\dot{f}}_{coll}$ term does not contribute, and (3.13) accurately describes the concentration profiles in both cases. As the volume fraction increases, the sampling bias integral becomes more negligible compared to the turbophoresis integral, meaning that the higher volume fraction cases have concentration profiles very close to the inverse of their wall-normal particle velocity variance profiles, equation (3.14).

Figure 6. For $St^{+}=32$, the comparison of left (continuous lines) and right sides (dashed lines) of (a) (3.10) and (b) (3.14). Colours indicate $\unicode[STIX]{x1D6F7}_{V}$ following the legend in figure 5.

5.1 Modelling and simulation details

The friction Reynolds number for the simulations in this section is $Re_{\ast }=600$, which is high enough for WMLES to make sense while keeping the cost of DNS at a reasonable level. The number of grid points in each direction for the DNS is $682\times 342\times 512$, in keeping with the grid spacings given in the previous section. The WMLES grid is just $128\times 20\times 96$, a reduction by a factor of almost $500$. Unlike the DNS grid described in the previous section, the WMLES grid has uniform spacing in the wall-normal direction, so the first grid point is at $y_{1}^{+}=60$ for the wall-normal velocities and at $y_{1/2}^{+}=30$ for the streamwise and spanwise velocities. The details of the discretization are the same, but the WMLES solves the filtered Navier–Stokes equations,

(5.1)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\widetilde{\boldsymbol{u}}+\widetilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\widetilde{\boldsymbol{u}}=-\frac{1}{\unicode[STIX]{x1D70C}_{f}}\unicode[STIX]{x1D735}\widetilde{p}+\unicode[STIX]{x1D708}_{f}\unicode[STIX]{x1D6FB}^{2}\widetilde{\boldsymbol{u}}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{ij}=\widetilde{u_{i}u_{j}}-\widetilde{u}_{i}\widetilde{u}_{j}$ is the subgrid stress tensor computed using the dynamic Smagorinsky model (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991), where the dynamic constant is found using the least-squares approximation of Lilly (Reference Lilly1992). The algebraic equilibrium wall model is used in lieu of the no-slip condition because the wall is not resolved. This wall model computes the wall shear stress locally (at each wall-adjacent grid point) from,

(5.2)$$\begin{eqnarray}\displaystyle u_{\Vert }(y_{1/2})=u_{\ast }\left[\frac{1}{\unicode[STIX]{x1D705}}\ln \left(\frac{y_{1/2}u_{\ast }}{\unicode[STIX]{x1D708}_{f}}\right)+B\right], & & \displaystyle\end{eqnarray}$$

using the wall-parallel velocity magnitude, $u_{\Vert }(y_{1/2})=\sqrt{u_{x}^{2}+u_{z}^{2}}$, at the first grid point, $y_{1/2}$, to find $\unicode[STIX]{x1D70F}_{w}=\unicode[STIX]{x1D70C}_{f}u_{\ast }^{2}$ via Newton iterations. Constants $\unicode[STIX]{x1D705}=0.41$ and $B=5.2$ are used. The coarse-grained WMLES velocity field captures the large-scale motions and mean flow, but is lacking small-scale features throughout the channel, most critically below the first grid point, $y<y_{1/2}$, where much of the physics related to biased sampling occurs.

The discretization schemes for WMLES are the same as the DNS described above. The same equations for advancing particle trajectories are used, equation (2.2), but with the filtered velocity, $\widetilde{\boldsymbol{u}}$, used for the fluid velocity seen by the particle instead of the (unknown) full velocity, $\boldsymbol{u}$. While deconvolution schemes for approximately reconstructing $\boldsymbol{u}$ from $\widetilde{\boldsymbol{u}}$ (Kuerten Reference Kuerten2006) may work well with WRLES, it is unlikely to produce significant improvements in WMLES because of the much more severe lack of resolution in the near-wall region. Instead, no attempt is presently made to model the SGS velocity fluctuations (as seen by the particle) in the near-wall region. The Smagorinsky subgrid stress model that is used does imply subgrid fluctuations, so in that sense, the fluid velocity seen by the particles is not fully consistent with the fluid momentum equation solved (Minier Reference Minier2016). Nevertheless, the purpose of this investigation is to first focus on how the velocity seen by the particle can be interpolated from the first grid point off the wall, seeing that standard linear or higher-order interpolation schemes assume the flow variation is captured on the grid. Particle–particle collisions are computed deterministically using the same hard-sphere collision model so that the impact of volume fraction, $\unicode[STIX]{x1D6F7}_{V}$, can be demonstrated as well. The particle dynamics has no stochastic input, and so remains deterministic as in the case of DNS.

In the general case, it should be appreciated that the lack of subgrid fluctuations in the present WMLES approach can impact the collision rates predicted by the simulation. The collision cylinder approach used at present is based on smooth particle trajectories within a given time step. Additional sub time step variation in the trajectory (e.g. from a stochastic model) could change the collision rate, most likely increasing it. For the case of inertialess particle–wall collisions in a stochastic modelling framework, Dreeben & Pope (Reference Dreeben and Pope1998) proposed a method for dealing with this effect. Henry et al. (Reference Henry, Minier, Mohaupt, Profeta, Pozorski and Tanière2014) proposed a model for collisions between inertial particles in a stochastic framework. In the present work, because of the relatively moderate $Re_{\unicode[STIX]{x1D70F}}$ (time-scale separation) along with the $St^{+}\gg 1$ values used, $\unicode[STIX]{x0394}t_{WMLES}/\unicode[STIX]{x1D70F}_{p}\ll 1$ so that these effects are not crucial at present.

5.2 Interpolation near the wall

While no attempt is made in the WMLES cases to model subgrid coherent structures, the interpolation scheme can have an important impact on the results. The baseline interpolation scheme is simply the trilinear interpolation used in DNS, which is consistent with a no-slip boundary condition. Because the viscous and buffer layers are not resolved by the WMLES, even the mean streamwise velocity seen by the particles below the first grid point is not well represented by linear interpolation. For that reason, a wall-parallel velocity interpolation consistent with the equilibrium stress model, denoted by ${\mathcal{I}}_{\Vert }$, is developed. This interpolation scheme computes the wall-parallel velocity at a particle location as,

(5.3)$$\begin{eqnarray}\displaystyle \widetilde{u}_{\Vert }(y_{p})=F_{\Vert }(y^{+})\widetilde{u}_{\Vert }(y_{1/2}), & & \displaystyle\end{eqnarray}$$

where $F_{\Vert }(y^{+})=\langle u_{x}\rangle (y_{p}^{+})/\langle u_{x}\rangle (y_{1/2}^{+})$ is specified for the present purposes by a universal velocity profile from Liakopoulos (Reference Liakopoulos1984) which reproduces viscous, buffer, and log-law layers of a typical turbulent wall-bounded flow. An equivalent model could solve the RANS equations on a one-dimensional (1-D) grid below the first grid point, as often done for the equilibrium stress model (Cabot & Moin Reference Cabot and Moin2000; Bose & Park Reference Bose and Park2018), and use this velocity to advect particles.

Furthermore, the analysis in § 3 showed that the concentration profile for high Stokes number particles could potentially be predicted without knowledge of the details of near wall flow structures, provided that the profile of wall-normal velocity variance is accurately represented. To test this idea, the linear interpolation of the wall-normal velocity can be replaced by a more detailed interpolation scheme, denoted by ${\mathcal{I}}_{\bot }$, which is designed to reproduce a realistic profile of wall-normal fluid velocity variance. To do this, the vertical velocity at a particle location is interpolated using,

(5.4)$$\begin{eqnarray}\displaystyle \widetilde{u}_{\bot }(y_{p})=F_{\bot }(y^{+})\widetilde{u}_{\bot }(y_{1}), & & \displaystyle\end{eqnarray}$$

where $F_{\bot }(y^{+})=\sqrt{\langle u_{\bot }^{\prime 2}\rangle }(y_{p})/\sqrt{\langle u_{\bot }^{\prime 2}\rangle }(y_{1})$ is given by a wall-normal velocity root-mean-square profile. For the present work, this profile was simply extracted from DNS and fitted with a sixth-order polynomial function with a $0.9\,\%$ relative error with respect to the DNS profile. This strategy (${\mathcal{I}}_{\bot }$) enforces the right profile of the sampled wall-normal fluid velocity variance that the particles experience near the wall. As with the Liakopoulos profile, this approach is chosen for simplicity of implementation, since this work is viewed more as a conceptual test rather than as a specific modelling proposal. One could imagine constructing $F_{\bot }(y^{+})$ by using an existing Reynolds stress transport (Durbin Reference Durbin1993) or $v^{2}-f$ (Durbin Reference Durbin1991) model to solve for $\langle u_{\bot }^{\prime 2}\rangle$ on a 1-D grid similar to what is done for the mean velocity in equilibrium wall models. Pursuing the details and implementation of the model is beyond the scope of the present work, rather we seek to demonstrate how well such a model could be expected to perform.