Particle transport in a turbulent pipe flow: direct numerical simulations, phenomenological modelling and physical mechanisms

Abstract

In particle-laden turbulent wall flows, transport of particles towards solid walls is phenomenologically thought to be governed by the wall-normal turbulence intensity supporting the underlying particle–eddy interactions that are usually modelled by a combination of turbophoresis and turbulent diffusion. We estimate the turbophoretic and turbulent diffusive coefficients as a function of wall-normal coordinate directly from a generated direct numerical simulation (DNS) database of low volume fraction point particles in a turbulent pipe flow. These coefficients are then used in an advection–diffusion equation to estimate the particle concentration as a function of wall-normal distance and time, with favourable comparison against DNS for smaller Stokes number ($St^+$) particles suggesting a limitation of the common gradient diffusion hypothesis for larger $St^+$ particles. Using DNS we explore the non-trivial effects of $St^+$, pipe wall condition (particle absorbing or elastic) as well as the influence of drag and lift force on the velocity and particle statistics giving rise to different particle concentrations. We then appraise various Eulerian-based models of turbophoretic and turbulent diffusive coefficients and, finally, use physical insights from Lagrangian correlation times, conditional quadrant analysis and flow topology to shed further light on the particle transport as a function of various parameters and the limits of gradient diffusion hypothesis.

1. Introduction

The transport of particles in turbulent pipe flows is of particular relevance in a wide range of applications in technology and science. These include drying of mined resources, as well as mineral ore particles (such as drying of brown coal for energy production; e.g. Okuma et al. 1989) and transport of inhaled airborne medication (e.g. inhaled salbutamol for asthma treatment; Melchor et al. 1993). Particles in a turbulent flow are advected along curved streamlines. Streamline curvature is an essential feature of turbulence, with large curvatures representing small-scale motions (e.g. Schaefer 2012). Inertia associated with a particle limits its behaviour to follow curved streamlines of the surrounding fluid, and this feature can markedly differ as their sluggishness increases from infinitesimal to finite-sized particles (e.g. Haller & Sapsis 2008). Presently, from the multitude of particle–turbulence physics that exists, we restrict our attention to gas–solid flows, and where particles are small compared with all dynamical scales of turbulence. Gravitational settling is not considered but many of the concepts presented are flexible and could be adapted to scenarios where gravitational settling is a factor, such as horizontal pipes. The particles considered in this paper are sub-Kolmogorov in diameter, defined by the ratio of their diameter $d_p$ to the Kolmogorov scale, $d_p / \eta < 1$ which characterises sensitivity to carrier phase gradient scales across their diameter. In addition, the particle Reynolds numbers $Re_p = d_p |\boldsymbol{u}-\boldsymbol{u}_p|/\nu$ (where $\boldsymbol{u}_p$, $\boldsymbol{u}$ and $\nu$ are the particle velocity and velocity and kinematic viscosity of the fluid, respectively) in this study are typically low. Vortex shedding only occurs when particle Reynolds number $Re_p > 400$ (Elghobashi 1994) and with sufficient particle loading is then capable of augmenting turbulence.

In wall-bounded flows, the particle relaxation time $\tau _p = \rho _p d_p^2 / (18 \rho _f \nu )$ (with $\rho _p$ and $\rho _f$ the particle and fluid densities, respectively) is non-dimensionalised with the viscous time scale $\tau _f = \nu /U_\tau ^2 \ (U_\tau = \sqrt {\tau _w/\rho _f}$ is the friction velocity and $\tau _{w}$ is the wall shear stress), and this parameter is known as the viscous Stokes number $St^+ = \tau _p/\tau _f$. The two main effects associated with finite-sized particles in wall flows well known in the literature are particle clustering (Fessler, Kulick & Eaton 1994), also termed preferential concentration, and turbophoresis (e.g. Sardina et al. 2012). For a recent review of particles in wall flows (as well as in homogeneous turbulence) see Brandt & Coletti (2022), and as noted in that particle-laden turbulence review, there is a paucity of experimental data featuring particles in wall-bounded turbulent flows. Further, when point-particle direct numerical simulation (PP-DNS) of pipes (Picano, Sardina & Casciola 2009) is compared with channels (Sardina et al. 2012), an augmentation of near-wall particle concentration for pipes is observed. This greater wall accumulation for pipes was directly contrasted by Zahtila et al. (2020), suggesting that the two geometric distinguishing features of a pipe, an axisymmetric restriction on centreline and curvature in the solid surface at the wall, are not negligible. A pipe was chosen for this study to maintain intimate connection with the preceding deposition experiments and applied mathematical modelling. In the present study, our focuses are, turbophoresis (the migration of particles towards solid walls due to variation in turbulence intensity) specifically in pipe flows with varying $St^+$ particles, influence of particle–wall collision boundary condition type and the role of forces experienced by the inertial particles.

The phenomenon of turbophoresis is particularly relevant for particles depositing on a surface. The chief illustration of this phenomena has been experimental observations in turbulent vertical pipes in the seminal work of Liu & Agarwal (1974). Olive oil droplets from the output of an aerosol generator were connected to a steady air flow through a pipe and they found that the quantity of particulate deposited onto the inside pipe surface increased by orders of magnitude with increasing $St^+$. This dramatic increase in deposition arising from increased particle relaxation time (i.e. $St^+$) gave clear indication of an inertial effect due to interaction between finite size particles and turbulence. Peak deposition occurred at $St^+ \approx 30$, thereafter the deposition quantity remained more-or-less constant or slightly declined indicating a saturation limit that prevents further growth in the deposition rate.

Particle deposition is typically characterised by the deposition velocity, $V_{dep} \equiv J_{wall}/c_b$, which is the particle transfer rate on the wall, where the particle flux $J_{wall} \equiv N_{p,t}/({A}_{pipe}\Delta t)$; $N_{p,t}$ is the number of particles transferred to ${A}_{pipe}$ the surface area of the pipe in the time interval $\Delta t$, and the bulk concentration $c_b$ is the ratio of the total number of particles in the pipe at any given instant ($N_p$) and the pipe volume ($V_{pipe}$), i.e. $c_b \equiv N_p/V_{pipe}$ (e.g. Zhang & Ahmadi 2000), which can change with time depending on the wall condition (absorbing or not). Here the deposition velocity is non-dimensionalised by the friction velocity, $V^{+}_{dep} = V_{dep}/U_{\tau }$. We use the superscript $+$ to denote quantities scaled with the viscous variables ($U_\tau$ and $\nu$). This definition is well suited to particulates that adhere to surfaces upon collision. In this paper we consider both fully absorbing and elastically rebounding particles, because real pipes are somewhere in between. For the elastic case, we define particles to have ‘transferred to the wall’ if they reside within a single viscous unit of the wall. Note that $V_{dep}$ is a function of time, but with increasing time $V_{dep}$ reaches a nominally steady value (as described later).

Figure 1 presents from the literature experimental data and analytical estimates of the deposition velocity $V_{dep}^{+}$ in the nominal ‘steady-state’ scenario along with values computed in this study from DNS simulation sets named ‘Set 1’ through to ‘Set 4’. (Note that Sets $2$ and $4$ are ‘absorbing’ walls, consistent with the experimental conditions, whereas the other two sets are ‘elastic’ or ‘reflecting’ walls. The details are discussed later.) The dramatic increase in particle deposition at around $St^+ \simeq O(1)$ arises because low $St^+$ particles behave more like fluid elements which are less likely to deposit owing to their inability to reach the wall (as a consequence of the impermeability condition on the wall). On the other hand, due to their inertia, larger $St^+$ particles when wall-directed can evade the no-slip condition effects enforced on a fluid continuum and can reach the wall, hence resulting in increased deposition. In figure 1, there is agreement in the overall trend between each of these datasets. Figure 1 also shows three distinct regimes of deposition: (i) for low $St^+$ turbulent diffusivity dominates where $V_{dep}^{+}$ is nominally constant or there is slight decrease (due to Brownian diffusion that we ignore in this paper); (ii) turbulent diffusion by eddy impaction, in which deposition dramatically increases due to an increasing turbophoretic velocity, and (iii) particle inertia-moderated regime, where $V_{dep}^{+}$ is again nominally constant or there is a slight decline (due to an Eulerian acceleration term, cf. § 4). Importantly, however, we point out that as time passes, spatial concentration gradients form in the pipe, and a single value of $V_{dep}^{+}$ to characterise particle concentration may be misleading.

Figure 1. Deposition velocity variation with Stokes number, $St^+$, current DNS data collected for $t^+ = 720$ viscous time units ($t^+ = t/(\nu /U_\tau ^2)$) after initial uniform concentration reproduces corresponding experimental and analytical studies. Current data symbols corresponding to constituent simulation sets detailed in table 1; experimental work (Liu & Agarwal 1974) represented by open circles and analytical reference represented by a solid blue line (Guha 1997). Values averaged over the entire simulation interval are included (with small right shift for clarity).

Starting from a uniform concentration of particles in a realistic pipe (of finite length), it is not clear, a priori, whether a steady state will be reached or whether the steady-state value of $V_{dep}^{+}$ would be still applicable if the particle concentration within the pipe falls to very low values after continuous deposition. The unsteady evolving concentration gradients remain important in understanding the particle migration. In addition, note that figure 1 is only concerned about particle concentration at the wall, whereas one could be interested in the particle concentration as a function of distance from the wall. Furthermore, delineating the effects of wall conditions (absorbing or elastic) and $St^+$ requires additional investigation. The present paper attempts to clarify our understanding of unsteady particle concentration variations in pipes as well as the pertinent phenomenological modelling and fluid dynamics analysis. Before presenting the specific aims of this paper, we briefly describe some relevant aspects of particle migration.

1.1. The phenomenological and fluid dynamical aspects of particle migration

From a phenomenological perspective, particles without inertia migrate primarily owing to ‘turbulent diffusion’, but inertial particles have an additional mode of transport, i.e. turbophoresis, which is owing to variation in turbulence along the wall-normal direction. These two modes can simplistically be illustrated by a kinetic theory argument following Reeks (1983). Here the inertial particle movement is not due to random molecular motion, rather owing to turbulent eddies. Consider a one-dimensional case where the particle concentration $c$ is a function of distance $r$, i.e. $c(r)$. Particles with relaxation time $\tau _p$ and velocity $v$ will take a distance $a \sim v \tau _p$ (the ‘particle stop distance’) in order to reach a particular location $r$. Thus, the net particle flux $J$ at $r$ (in the positive direction of $r$) is the difference between fluxes coming towards $r$ from distances $(r-a)$ and from $(r+a)$. Note that the flux from $r+a$ to $r$ is in the negative $r$ direction. As flux is the product of concentration and velocity,

(1.1)\begin{equation} \left.\begin{aligned} J & \sim c(r-a)v(r-a) - c(r+a)v(r+a) , \\ & \approx \left( c - a\frac{\partial c}{\partial r} \right) \left( v - a\frac{\partial v}{\partial r} \right) - \left( c + a\frac{\partial c}{\partial r} \right) \left( v + a\frac{\partial v}{\partial r} \right) , \\ & = \left({-}2a \frac{\partial v}{\partial r}\right) c - ( 2av) \frac{\partial c}{\partial r},\ \text{and with } a \sim v \tau_p, \\ J & \sim \left( -\tau_p \frac{\partial v^2}{\partial r}\right) c - ( 2v^2 \tau_p) \frac{\partial c}{\partial r}. \end{aligned}\right\} \end{equation}

Here, $v$ is usually replaced by the particle root-mean-square (r.m.s.) velocity in the $r$ direction. In (1.1) the first and second terms are the contributions due to turbophoresis and turbulent diffusion, respectively. The coefficient of $c$ in the first term is related to the turbophoretic velocity $D_1$ and the negative coefficient of ${\partial c}/{\partial r}$ is akin to the turbulent diffusivity $D_2$. Therefore, phenomenological considerations suggest

(1.2)\begin{equation} {J}(r) = D_1(r) c(r) - D_2(r) \frac{\partial c(r)}{\partial r}.\end{equation}

The negative sign in front of $D_2$ is simply related to the fact that conventionally when only turbulent (or Fickian) diffusion is present and $c$ increases with $r$, the flux is in the negative $r$ direction. It is important to note that gradient diffusion models such as (1.2) are strictly valid only where the concentration is essentially constant over the distance $a$, or when the turbulent motions posses scales that are small compared with the distance over which the concentration varies (e.g. Corrsin 1974).

Turbulent diffusion of inertial particles given by $D_2$ have been considered by many authors (e.g. Tchen 1947; Friedlander 1957; Csanady 1963) along the lines of diffusion of passive scalars. The term $D_1 \sim -\tau _p {\partial v^2}/{\partial r}$, however, was independently derived by Caporaloni et al. (1975) and Reeks (1983), who named this migration phenomenon ‘turbophoresis’.

The turbophoretic term is essential to explain the deposition dependency on particle size that has been observed experimentally and yields better agreement with measurements than diffusion-only models (e.g. Davies 1966), which underpredict deposition by orders of magnitude. There has been considerable effort spent in modelling $D_1$ and $D_2$, which are fundamentally related to Lagrangian particle and fluid motions, in terms of Eulerian fluid statistics that are more easily available (and we comment on some aspects of it in § 4). In this regard, a simple unified theory of deposition was constructed by Guha (1997) starting from Reynolds averaging of the Eulerian particle continuity and momentum equations (which are discussed further in § 4). Notwithstanding the modelling efforts, we still do not know the distribution of the ‘actual’ $D_1$ and $D_2$ in a turbulent wall flow. We use our DNS data to evaluate these distributions (in § 4).

Among others, fluid dynamics aspects of particle migration in a boundary layer were studied through the means of DNS by Marchioli & Soldati (2002) in a turbulent channel. Particle transfer was found to be dominated by the same coherent structures that control momentum transfer to the wall. For example, particles of relaxation timescales comparable to the viscous time scale ($St^+ = 3$) with wall-directed velocity were highly correlated to advection by turbulent sweeps. The effectiveness of these transfer mechanisms declined with increase in particle relaxation time due to the larger timescales associated with heavy inertial particles. These particles are able to decorrelate their trajectory from fluid events and to rebound elastically from a solid wall, which introduces wall effects (e.g. Chan et al. 2021). In addition, the role of gravity, that we do not consider in this paper, is emerging as an additional complication when the flow is horizontal, as additional attention is required for the role of preferential sampling of fluid turbulent events, turbophoretic velocity and the Stokes settling velocity (e.g. Bragg, Richter & Wang 2021).

1.2. Present work

In this paper, one of our aims is to determine the turbophoretic and diffusive flux directly from the PP-DNS of a particle-laden turbulent pipe flow. The time-evolving contributions of $D_1$ and $D_2$ are particularly unclear in the case of rebounding particle–wall collisions in conjunction with sharp concentration gradients that establish close to the wall. We vary the coefficient of restitution to delineate the bounds of deposition as influenced by particle–wall collisions being absorbing or elastic to investigate effects of wall absorption featuring prominently in the theory of deposition (e.g. Johansen 1991; Young & Leeming 1997). This is of crucial importance as the coefficient of restitution is not widely known for many scenarios. We also investigate the effect of the Saffman–Mei lift force $\boldsymbol {F}_l$ (along with the compulsory dominant drag force $\boldsymbol {F}_d$). Particles across four decades of $St^{+}$ are considered and this constitutes particles spanning all distinct deposition regimes (e.g. Guha (1997), and cf. figure 1).

In § 2 we provide the details of the DNS set-up for a one-way coupled particle-laden turbulent pipe flow. Per the usual meaning, one-way coupling implies that the particles are affected by the fluid but the fluid flow is not affected by the particles, implying we consider low volume fraction of particles. The simulations are carried out for a single Reynolds number but four different $St^+$, absorbing/elastic particle–wall collision types and combinations of $\boldsymbol {F}_d$ and $\boldsymbol {F}_l$, resulting in a total of $16$ different simulation datasets (and an additional 2 datasets to address the effects of density ratio $\rho _p/\rho _f$, if any). Subsequently, we present basic particle statistics (mean and r.m.s. velocities) comparing them to fluid statistics, as well as time-varying particle concentration (and deposition) statistics in § 3. The procedure to determine turbophoretic and diffusive coefficients from DNS as well as accompanying $D_1$ and $D_2$ wall-normal profiles are presented in § 4. Here we also present existing turbulence-based models for $D_1$ and $D_2$ and compare them with our DNS results. Furthermore, we use these $D_1$ and $D_2$ to solve the one-dimensional phenomenological advection-diffusion equation for particle concentration $c$, and compare the time-varying $c$ with DNS. In § 5 we explore physical fluid dynamics aspects behind the particle migration using flow visualisation, Lagrangian particle auto-correlation, Lagrangian time scales and conditional particle statistics based on coherent motions, identified by fluid ejections/sweeps as well as by the invariants of the fluid velocity-gradient tensor. Finally, we conclude in § 6.

2. Computational set-up

The Euler–Lagrangian (e.g. Balachandar & Eaton 2010) framework is used to carry out the particle-laden flow simulations in a pipe, ignoring the effect of gravity (figure 2). The simulations are carried out using the open-source code, OpenFOAM (Weller et al. 1998), which has been validated as capable of conducting particle-laden DNS-type calculations (e.g. Chan et al. 2021), provided there is sufficient time and grid resolution to capture the full range of dynamical scales. The carrier (fluid) phase is fully developed and statistically steady prior to the injection of particles into the domain. The particles are stationary (zero velocity) when released and evenly distributed throughout, which yields a uniform concentration.

Figure 2. Sketch of the numerical smooth-wall pipe flow set-up with dispersed particles. Eulerian simulation of the fluid with Lagrangian particle tracking.

The Navier–Stokes equations for incompressible flows are solved for the fluid phase:

(2.1)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} ={-}\frac{1}{\rho_f}\boldsymbol{\nabla} p + \nu\boldsymbol{\nabla}^2 \boldsymbol{u} + F_{x}\boldsymbol{i}, \end{gather}$$

where $\boldsymbol {u}$ is the fluid velocity vector, $p$ is the pressure fluctuation and ${F_{x}}(t)$ is a spatially uniform forcing term to ensure a constant mass flux in the streamwise direction. Equations (2.1) and (2.2) are solved with the conventional finite-volume method approach and for the pressure–velocity coupling consisting of the merged Pressure Implicit Splitting Operator (PISO) and Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithms, PISO–SIMPLE (PIMPLE). Further, OpenFOAM is a collocated grid solver and to remove spurious pressure oscillations, careful discretisation is required of the pressure gradient and Laplacian of pressure that appears in the pressure equation. Avoiding oscillations on non-staggered grids was first suggested by Rhie & Chow (1983), and constitutes a remedy for the error term that arises from an inconsistent stencil for the gradient and Laplacian operators. This requirement to avoid spurious oscillations is satisfied in OpenFOAM by applying Gauss’ theorem to the Laplacian term that arises in the pressure equation and by interpolating variables stored at cell centers to cell faces. A detailed account of Rhie–Chow correction in OpenFOAM is available from Kärrholm (2006). An unstructured hybrid Cartesian ‘O-grid’ discretisation of hexahedral elements was used in this study, which follows the meshing details outlined in Chan et al. (2015) and Cheng et al. (2020). The cell dimensions at the wall are, $\Delta y^{+} \approx 0.22$, $R\Delta \theta ^{+} \approx 4.71$ and $\Delta x^{+} = 5.89$. A geometric stretching is applied in the wall-normal direction with a cell-to-cell expansion ratio of less than 1.05; elements are equidistant in the streamwise direction and in the radial–azimuthal plane follow an ‘O-grid’. As is documented in an earlier work (Chan et al. 2021), these cell dimensions can be compared with the computed Kolmogorov scale $\eta$, taken at the wall. We report Kolmogorov-scaled grid spacing of $\Delta y / \eta \approx 0.14$ (at the wall), $R\Delta \theta / \eta \approx 2.94$ and $\Delta x / \eta = 3.68$. This is an acceptable resolution for DNS, first because the main contribution to dissipation arises from wall-normal gradients and second because average dissipation occurs on only $O(\eta )$ rather than $\eta$ itself (Moin & Mahesh 1998). Further, the convective terms are obtained by linear interpolation to surface quantities, and diffusive terms include non-orthogonal correction for the Laplacian operator; in addition, there is no employment of flux limiters or upwinding in this study. The OpenFOAM Eulerian fluid phase is verified against simulations conducted with the spectrally accurate DNS code Nek5000 (Fischer, Lottes & Kerkemeier 2008), and good agreement is found (see Appendix A).

The simulation was conducted in a periodic smooth wall pipe with a domain length of $L_x = 4 {\rm \pi}R$ where $R$ is the radius of the pipe. This domain length was found to be sufficiently long to obtain converged second-order fluid statistics in a smooth wall pipe (e.g. Chin et al. 2010). All simulations are carried out at an approximately fixed friction Reynolds number $Re_\tau = {R U_\tau }/{\nu } = 180$. Fluid quantities are interpolated to particle coordinates by means of inverse distance weighting (e.g. Shepard 1968) from the Eulerian grid.

We assume particles as rigid spheres, with high density ratio ($\rho _p / \rho _f= 1000$). The general equation of motion for rigid spherical particles is given by Maxey & Riley (1983), and under our circumstances of high particle-to-fluid density ratio, many of the forces on a small rigid sphere become negligible, such as the fluid acceleration, added mass and Basset terms. The last of these will ostensibly decay and is omitted to keep the physical processes as simple as possible and to allow comparison with other researchers. Gravity is omitted in order to preserve the generality of flow direction statistics; however, researchers have found that for heavy particles, the inclusion of gravity has a non-monotonic effect on velocity profiles and deposition coefficients (e.g. Marchioli, Picciotto & Soldati 2007). Further, an important parameter in particle-laden flows is the volume fraction of particles. For each particle species being studied, there is a corresponding sufficiently high particle volume fraction beyond which the presence of particles is classified as modifying the turbulence, and these studies are termed two-way coupled. Our one-way coupling simulations are only applicable for dilute loading ratios of particles, which would correspond to a particle volume fraction of $\varPhi _v$ less than $\approx 10^{-6}$ (e.g. Elgobashi 2006).

In this study, the point-particle dynamics are primarily driven by drag force considerations and of secondary importance are models for the shear and lift forces acting on the particle. The drag correlation of Putnam (1961) is used to correct Stokes drag for finite Reynolds number effects. This requires particle-flow slip velocity information obtained by interpolation of the fluid velocity from the Eulerian mesh to the instantaneous particle location, yielding the drag force $\boldsymbol {F}_d$. The lift force $\boldsymbol {F}_l$ is modelled by using the Saffman–Mei lift force (Mei 1992), which empirically corrects for finite Reynolds number effects the original expression for the lift on a sphere in a shear flow (Saffman 1965).

Among the catalogue of expressions for shear-induced lift force, Saffman–Mei was chosen as its ability to capture near-wall statistics is corroborated by state-of-the-art particle-resolved DNS (Costa, Brandt & Picano 2020b). As such, equations for particle motion are

(2.3)$$\begin{gather} \frac{{\rm d} \boldsymbol{x}_p}{{\rm d} t} = \boldsymbol{u}_p, \end{gather}$$

(2.4)$$\begin{gather}m_p \frac{{\rm d} \boldsymbol{u}_p}{{\rm d} t} = \underbrace{\frac{m_p}{\tau_p}\left(1+\frac{1}{6} Re_p^{2/3}\right)(\boldsymbol{u}-\boldsymbol{u}_p)}_{\boldsymbol{F}_d} +\underbrace{ m_p \frac{\rho_f}{\rho_p}C_L ((\boldsymbol{u} - \boldsymbol{u_{p}}) \times \boldsymbol{\omega})}_{\boldsymbol{F}_l}, \end{gather}$$

where

(2.5a,b)\begin{equation} C_L = \frac{3}{2{\rm \pi} \sqrt{Re_s}} C_{LS} , \quad \text{and} \quad C_{LS} = \begin{cases} 6.46f & Re_p < 40 ,\\ 6.46 \cdot 0.0524 \sqrt{\beta Re_p} & Re_p \geq 40. \end{cases} \end{equation}

Here $\boldsymbol {x}_p$ and $\boldsymbol {u}_p$ are particle position and velocity vectors, $m_p$ is the spherical particle mass and $\boldsymbol {\omega } = \nabla \times \boldsymbol {u}$ is the vorticity of the fluid at the location of the point particles, $Re_p = d_p |\boldsymbol{u}-\boldsymbol{u}_p|/\nu$ and $Re_s = {d_p^{2} |\boldsymbol {\omega }|}/{\nu }$, $\beta = Re_s/(2 Re_p)$, $f = (1 - \alpha ) \exp (-0.1 Re_p) + \alpha$ and $\alpha = 0.3314 \sqrt {\beta }$. The barycentric tracking algorithm tracks particles on the mesh and information is interpolated to particle coordinates. The ordinary differential equations (ODEs) associated with particle motion are integrated in time through an implicit Euler scheme. These details are the same as our earlier contribution in Chan et al. (2021).

Adhesiveness of solid walls determines whether upon impact a particle will absorb or rebound. Particle–wall collision requires attention as analytical (e.g. Guha 1997) and experimental treatments (e.g. Liu & Agarwal 1974) of deposition have thus far considered only the absorbing case and the validity of extending these treatments to rebounding walls remains unknown. The particle–wall collision model is studied by considering the extremes in the normal and tangential coefficients of restitution, $e_{n}$ and $e_{t}$, respectively. Simulations are conducted for fully absorbing ($e_n = 0$, $e_t = 1$) and perfectly elastic ($e_n = 1$, $e_t = 0$) walls. The particle tracking algorithm in OpenFOAM follows a particle until a cell face crossing occurs (Macpherson, Nordin & Weller 2009), and when this cell face represents a solid wall, an interaction model is required. The particle velocity is updated upon crossing a wall face by

(2.6)$$\begin{gather} \boldsymbol{u}_{p,n}^{u} ={-}e_n \boldsymbol{u}_{p,n}^{o}, \end{gather}$$

(2.7)$$\begin{gather}\boldsymbol{u}_{p,t}^{u} = (1- e_t)\boldsymbol{u}_{p,t}^{o}, \end{gather}$$

where $\boldsymbol {u}_{p,n}$ and $\boldsymbol {u}_{p,t}$ are the particle normal and tangential velocity, respectively, and superscripts ‘$u$’ and ‘$o$’ denote ‘updated’ and ‘old’ quantities, respectively. Particle–particle collisions are neglected.

We present four core sets of DNSs (tables 1 and 2). In Set 1, particle dynamics are solely informed by $\boldsymbol {F}_d$, the drag force and particle–wall collisions are fully elastic ($e_n = 1$, $e_t = 0$), consistent with recent turbulent channel multiphase DNS datasets (e.g. Sardina et al. 2012; Johnson, Bassenne & Moin 2020). In Set 2, again particle dynamics are solely informed by $\boldsymbol {F}_d$, the drag force, but the particle–wall collisions are treated as fully absorbing ($e_n = 0$, $e_t = 1$), which is more consistent with the experiments on particle deposition (e.g. Liu & Agarwal 1974) and the analytical work presented in § 1.1. Sets 3 and 4, are similar to Sets 1 and 2, respectively, except that the Saffman–Mei lift force $\boldsymbol {F}_l$ is included in the particle dynamics, which was found to yield better agreement with fully resolved particle simulations (Costa et al. 2020b). The purpose of these four sets is to test the robustness of particle–wall interaction and force treatments on particle dynamics. Each of these sets are run with $St^+ = 0.1,1,10$ and $100$. Table 2 provides a summary of datasets generated along with the symbols used to denote each of the $16$ datasets. In addition, simulation Set 5 run at $\rho _p / \rho _f= 10{\,}000$ (compared with the other sets at $\rho _p / \rho _f= 1000$) was collected to support the key findings reported, and can be directly contrasted to simulation Set 3. For this set, only $St^{+} = 10$ and $100$ particles are studied as the Saffman–Mei lift force is not found to be prominent for the smaller diameters.

Table 1. Summary of different particle types used in simulation Sets 1-4, with constant $\rho _p / \rho _f= 1000$, and Set 5 featuring $\rho _p / \rho _f= 10{\,}000$. All simulation sets contain particle number $N_p = 2 \times 10^6$ for each species.

Table 2. Particle–wall interaction types and forces studied: $\boldsymbol {F}_d$ and $\boldsymbol {F}_l$ are the drag and Saffman–Mei lift forces, cases run at $Re_\tau = 180$, corresponding to centreline $Re_{cl} = 7043$ and bulk $Re_D = 5357$, based on the centreline and bulk velocity, respectively, and pipe diameter. Markers decrease lightness with increasing $St^{+}$ number. Simulation Set 5 is partitioned separately as it is briefly addresses density ratio effects that underpin the findings in this paper.

3. Particle velocity statistics and concentration

3.1. Statistical moments of particle phase

Instantaneous particle velocity is decomposed into mean and fluctuating components, e.g. $u_{px} = U_{px} + u_{px}^{\prime }$, and the mean and r.m.s. of the fluctuating particle velocities are shown in figure 3 for Set 1 (see table 2). Here $y = R - r$, where $R$ is the pipe radius. The data are temporally averaged over $\Delta t^+ = 3600$ (which is equal to 1190 integral time units, i.e. ratio of pipe radius to bulk velocity) from the beginning of particle injection. As expected, the $St^+ = 0.1$ particles behaving like tracers do not show a shift from the fluid in their mean velocity nor r.m.s. profiles. Note that fluid statistics are shown by solid grey lines. As $St^+$ increases, deviation between the streamwise velocity of the particle ($U_{px}^{+}$) and the fluid ($U_{x}^{+}$) is minor, except for $St^+ = 100$ particles, where the $U_{px}^{+}$ profile is significantly flatter. This is more clearly observed in the inset of figure 3(a), where the difference between the particle and fluid velocity $\Delta U^{+} = U_{px}^{+} - U^{+}$ is plotted. This behaviour has also been observed by Mortimer, Njobuenwu & Fairweather (2019) previously in channel flows. Although the precise details of this mechanism remains unknown, $St^+ = 100$ particles lag the fluid in the bulk flow but retain their inertia as they approach the near-wall region and, consequently, have higher-than-fluid streamwise velocity for $y^{+} \lessapprox 10$.

Figure 3. Turbulent statistics for simulations in Set 1 (see table 2 for symbols). Mean particle velocities in the (a) streamwise and (b) radial directions; r.m.s. of particle velocity fluctuations in the (c) streamwise, (d) radial directions and (e) the particle Reynolds shear stresses. Note that fluid statistics are shown by solid grey lines. Also plotted in (d) with the same symbols as in table 2 but grey-coloured is analytical prediction of the r.m.s. of radial particle velocity fluctuation, obtained from the DNS determined eddy viscosity and r.m.s. profile for the fluid via (4.7) and (4.8a,b) from § 4.2.

The particle mean radial velocity ($U_{pr}^{+}$) profiles show a peak at approximately $y^+ = 36$, indicating that particle drift towards walls is a feature of wall turbulence (recall that radial velocity is positive when directed towards the wall and $U_r^{+}$ for the fluid is zero). When particle–wall collisions are modelled as elastic, there are competing effects contributing to the $U_{pr}^{+}$; and some of them are, particles bouncing against the wall, particle sluggishness and preferential sampling of fluid events by the particles. The $St^+ = 100$ particles are qualitatively observed to collide with the wall and flee the near-wall region, re-entraining into the pipe turbulent core, whereas $St^+ = 10$ particles that begin near the wall and are wall-fleeing often show high curvatures in their trajectory of motion in the buffer region and return to the viscous sublayer.

The results of an international benchmark effort in turbulent channels presented in Marchioli et al. (2008) reveal considerable scatter between different numerical schemes in the prediction of the r.m.s. velocity fluctuations of particles. Results in their benchmark effort however consistently found, compared with the fluid, an increase in the peak for streamwise particle velocity fluctuations for particle populations $St^+ = 5$ and $25$, which is similar to figure 3(c) for pipes. These velocity fluctuations near the wall are significantly enhanced for the largest particles as they are able to retain their velocity from the buffer region while entering the viscous sublayer. Furthermore, consistent with Mortimer et al. (2019), figure 3(c) shows that at the wall particle velocity fluctuations do not go to zero for large $St^+$ values. However, we recall a limitation in our PP-DNS study: particle radius effects are not incorporated into the particle–wall collision model. The particle-tracking algorithm follows a particle until cell face crossing occurs and in the case of a boundary, applies the particle–wall boundary condition, (2.6) and (2.7). Having recognised this limitation, the largest particles in core simulation sets of this study are $d_p^{+} = 1.341$, and so $r_p^{+} = 0.67$ would be the collision radius but the high particle r.m.s. profiles persist until $y^{+} = 10$. In addition, in Appendix B, figure 18 reports that the near-wall profiles are characterised by the $St^{+}$ number as this finding persists with higher-density $\rho _p / \rho _f= 10\ 000$ particles.

The $u_{px,rms}^{\prime +}$ profiles in the outer region approximately collapse for all particle populations. In figure 3(d) the $u_{pr,rms}^{\prime +}$ profiles in the radial direction monotonically decrease in the outer region with increasing $St^+$ number, and all except the largest $St^+$ particles collapse onto the fluid profile very close to the wall. Turning attention to the particle Reynolds shear stresses $\overline {u_{px}^{\prime }u_{pr}^{\prime }}^{+}$ in figure 3(e), we observe a slight gain in the peak of the profile for $St^{+} = 0.1,1$ and $10$ particles and then collapse onto the fluid Reynolds shear stress away from the wall. The situation is different for $St^{+} = 100$ particles, where there is a slight gain close to the wall but thereafter, dampening in the profile compared with fluid. This is highly likely to be driven by the substantial decrease in the radial velocity fluctuations observed in figure 3(d).

The profiles presented in figure 3 for Set 1 are supplemented by comparison with the corresponding profiles for the remaining simulation sets in Appendix B, figures 18 and 19. Inclusion of the Saffman–Mei lift force for the elastic simulations yields similar statistics ($U_{px}^{+},u_{px,rms}^{\prime +},u_{pr,rms}^{\prime +}$) for particles with $St^+ = 0.1,1 \ \textrm {and} \ 10$ but $St^+ = 100$ particles show a striking difference in their near-wall quantities. This is most likely because the Saffman–Mei lift force for the latter plays a crucial role by assisting in wall-fleeing, and providing one mechanism for lowering the residence time of the particles in the viscous sublayer. Note that the ‘wall-fleeing’ effect of larger $St^+$ particles is due to the lift force (see (2.4)) where the wall-normal shear and lower particle velocity compared with fluid constitute a lift force away from the wall. Comparison between the elastic and absorbing wall simulation sets (cf. figures 18 and 19) shows increased $U_{pr}^{+}$ for absorbing walls because of the consequent lower wall-fleeing. Again absorbing walls highlight the distinctive $St^+ = 100$ particles, and classifies these particles as ballistic rather than merely deviating from flow behaviour. The $St^+ = 100$ particles also show the most significant deviation for $U_{pr}^{+}$ profiles between simulation sets, which suggests the importance of the wall and Saffman–Mei lift force for heavier or larger particles in future modelling efforts. The peculiar characteristics of $St^+=100$ particles will become apparent later when we consider modelling aspects.

Furthermore, we observe the lift force gains importance with increasing $St^{+}$ in the near wall-region, and approaches in the mean sense $20\,\%$ of the drag force magnitude close to the wall for the largest particles studied. In Appendix B, we further quantify the relative importance of the lift force by comparing its magnitude to that of the drag force as a function of distance from the wall. The role of the lift force is not immediately obvious, as typically in wall flows it aids particles to migrate away from walls but Vreman (2007) found inclusion of the lift force causes the near-wall particle concentration to increase by a factor $1.5$ or $2$, dependent on the particle size and overall particle loading, which is also transiently observed in this paper prior to strong concentration gradients being established (see discussion towards the end of Appendix B).

3.2. Time-evolving concentration and spatial distribution profile

The time-varying number fraction of inertial particles ($N_p / N_t$, where $N_p$ is the instantaneous number of particles residing in a region and $N_t$ is the total number of particles in the simulation) that reside within each of the turbulent regions, viscous sublayer, buffer, logarithmic and core, is presented in figure 4. This region-by-region analysis has been previously employed by Picano et al. (2009) where particles were modelled by drag force $\boldsymbol {F}_d$ alone and released on the centreline; their results for ${St^+ = 10}$ are also plotted in figure 4 by solid (grey) lines. In addition, reference is made to cross-sectional slices of instantaneous particle distributions in figure 5 plotted over pseudo-colours of fluid streamwise vorticity. Profiles of particle wall-normal number fraction and spatial preferential concentrations are jointly presented in this section as they are intimately connected (Sardina et al. 2012); e.g. biased sampling of fluid events in turn leads to ejection mechanisms near the wall that control particle transfer. Figure 5 highlights the preferential concentration in particle-laden flows (e.g. Squires & Eaton 1991) as well as the diffusional deposition that preferentially occurs in streamwise oriented streaks (e.g. Narayanan et al.'s (2003) observations in channel flows).

Figure 4. Time-evolving averaged local particle number ($N_p/N_t$) in regions: $(\textit {a})$ viscous sublayer, $(\textit {b})$ buffer, $(\textit {c})$ logarithmic and $(\textit {d})$ core. Solid grey line from Picano et al. (2009) for $St^+=10$, which injected particles on the centreline. See table 2 for symbols. Green lines are plotted to represent expected fraction in the case of uniform dispersion.

Figure 5. Streamwise slices of thickness $\approx R/12$ for instantaneous particle positions for (a) $t^+ = 36$, (b) $t^+ = 900$, (c) $t^+ = 1800$, (d) $t^+ = 2700$ and (e) $t^+ = 3600$ coloured by streamwise vorticity $\omega _x$. Results from simulation Set 1.

The key feature that we see in figures 4 and 5 for $St^+ = 0.1$ particles is very low levels of accumulation in the viscous sublayer with time. Even after $\approx 300$ integral time units ($\equiv U_b t / R$), wall-normal region-by-region particle number fraction shows only minor deviation due to accumulation in the viscous sublayer. These particles are able to trace flow structures of all time scales (see figure 5), thus remaining evenly distributed throughout the domain, both in terms of wall-normal distance and their local organisation with respect to flow structures.

Particles with $St^+ = 1$ accumulate at a moderate rate in the viscous sublayer suggesting that deposition cannot be neglected for this population. Accumulation of $St^+ = 1$ particles did not appreciably vary between simulation Sets 1–4 indicating that deposition modelling for this species does not require inclusion of $\boldsymbol {F}_l$, the lift force and particle–wall interactions did not have any significant effect. These particles are unable to trace high-wavenumber fluid motions and accumulate outside regions of high vorticity magnitude which can be observed in figure 5; particle–turbulence interactions highlighted by preferential concentration can be observed in all regions of the cross-section.

The situation is more complex for larger particles. For simulations in Set 1, $St^+ = 10$ and $St^+ = 100$ particles both wall-accumulate at a rate that is similar but orders of magnitude greater than the $St^+ = 0.1$ and $St^+ = 1$ species. To delineate the effects of $\boldsymbol {F}_l$ and wall-boundary conditions separately, let us first consider the absorbing wall case (so that particle ejection from the viscous sublayer effects are absent). For $St^+=10$ with absorbing walls, figure 4(a) shows an increased $N_p / N_t$ with $\boldsymbol {F}_l$ included (back diamonds) than without $\boldsymbol {F}_l$ (green diamonds). Note that $\boldsymbol {F}_l$ direction is towards (or away from) the wall when particle velocity is higher (or lower) than the fluid velocity. As the inertial $St^+=10$ particles approach the wall they retain their higher velocity and, hence, experience a wall-ward lift force increasing the $N_p / N_t$. On the other hand, $N_p / N_t$ distribution of $St^+=100$ does not seem to be significantly affected by $\boldsymbol {F}_l$. For elastic walls, in general, $N_p / N_t$ is lower because of particle rebound. Interestingly, for $St^+=10$ particles with elastic walls, although initially $\boldsymbol {F}_l$ helps to increase $N_p / N_t$ (similar to the absorbing wall case), as the time increases particles show longer residence times in the wall region. Next, they may be ejected by fluid motions or rebound against the wall or their now lower-than-fluid velocity directs $\boldsymbol {F}_l$ away from the wall. Therefore, over long time intervals particles with $\boldsymbol {F}_l$ included have a reduced $N_p / N_t$. This effect is accentuated for $St^+=100$ particles with elastic walls (compare red and blue pentagram symbols in figure 4(a).

Increased (decreased) $N_p / N_t$ within the near-wall region implies a decreased (increased) number fraction away from the wall. This is reflected for the buffer region shown in figure 4(b) especially in the case of $St^{+} = 100$ particles where the number fraction even increases relative to the initial distribution. This is consistent with the inertia-moderated regime described by Guha (1997) (see also figure 1), where deposition velocity decreases slightly with increasing $St^+$. Simulations in Set 3 reveal a distinctive phenomena; rather than wall-accumulate, particles converge to a steady state. Convergence is faster in $St^+ = 100$ particles and an order of magnitude decrease in viscous sublayer concentration compared with results from Set 1 suggests the importance of $\boldsymbol {F}_l$. Spatially, $St^+ = 10$ particles show major deviation from randomness, accumulating outside regions of intense vorticity (see again figure 5). However, $St^+ = 100$ particles with their significant inertia appear ballistic in nature, as they filter fluid events and remain randomly distributed. Note, it should be clear that $St^+ = 1$ and $St^+ = 100$ particles both remain randomly distributed but for totally different reasons.

A key finding in this section is the identification of relative importance of the lift force $\boldsymbol {F}_l$ and particle–wall interaction type for non-dimensional $St^+$ number. The effect of $\boldsymbol {F}_l$ is significant for $St^+ = 10$ particles and dramatic for $St^+ = 100$ particles when particle–wall interaction is elastic. The effect of $\boldsymbol {F}_l$ is less important when particle–wall interaction is absorbing type, and is unimportant for $St^+ = 0.1$ and $St^+ = 1$ particles. As such, figure 4 is useful in decision-making on modelling of particle dynamics and particle–wall interactions. In addition, in Appendix B, we investigate the relevance of density ratio $\rho _p/\rho _f$ and whether $St^{+}$ number alone is sufficient to characterise particle wall-accumulation by comparing velocity and particle concentration for Sets 3 and 5. There, it is reported that although the velocity statistics are similar, the density ratio may influence the ratio of Saffman–Mei lift force to drag force acting on the particle (because increased $\rho _p / \rho _f$ for some $St^{+}$ implies reduced $d_p$ and, hence, reduced $F_l$), and this in turns modifies the particle wall-normal concentration.

3.3. Deposition velocity

Figures 6(a), 6(b), 6(c) and 6(d) corresponding to $St^+=0.1,1,10$ and $100$ show the time series of $V_{dep}^+$ revealing varying levels of agreement between Sets $1$–$4$ for different Stokes numbers. Simulations with absorbing walls (green and black symbols) quickly attain a steady deposition rate. The slight exception is the $St^+ = 10$ particles in Set 2, where $V_{dep}^+$ gradually increases with time before becoming approximately constant. The presence of elastic walls results in either a gradual increase in the deposition velocity when only the drag force $\boldsymbol {F}_d$ is included, or with the inclusion of the Saffman–Mei lift force $\boldsymbol {F}_l$, a fast equilibrium concentration is achieved by $t^+ \approx 600$ and, thereafter, fluctuations occur in particles departing/re-accumulating in the near-wall region.

Figure 6. Deposition velocity time series for Sets 1–4: (a) $St^+ = 0.1$; (b) $St^+ = 1$; (c) $St^+ = 10$; and (d) $St^+ = 100$. Note that results from simulation Set 3, in which particles each time-interval may resuspend into the flow after wall-accumulating yielding ‘negative’ deposition velocity are marked by absolute valued filled circles. For symbols, see table 2.

These evolving $V_{dep}^+$ profiles in Sets 1 and 3 that feature elastically rebounding walls can be attributed to evolving wall-normal concentration gradients. As particle number fraction increases in the viscous sublayer (reported in § 3.2), there is a growing role of ejection mechanisms by vortical motions (e.g. Marchioli, Picciotto & Soldati 2003). In addition, the importance of the lift force $\boldsymbol {F}_l$ becomes clear as this additional mechanism assists particles to flee the viscous sublayer. In particular, for rebounding walls, time-evolving deposition velocity changes by several orders of magnitude for $St^+ = 10$ and $100$ particles when lift force $\boldsymbol {F}_l$ is included. The findings in this section indicate that the analytical modelling that prominently features absorbing walls and steady deposition rates (Guha 1997; Young & Leeming 1997) is on good footing, but the situation involving rebounding walls is more complex, where near-wall fleeing mechanisms are of vital importance.

The large variation in $V_{dep}^+$ for changing $St^+$ can be better illustrated by taking the average of $V_{dep}^+$ across time. These average values of $V_{dep}^+$ as a function of $St^+$ is shown in figure 1 for different Sets. The symbols denote an averaging time from $t^+=0$ to $720$, which corresponds to a time based on bulk flow, $t U_b ( = R \,(t U_\tau ^2/\nu ) (U_b R/\nu )(\nu ^2/(U_\tau ^2 R^2))=R t^+ Re_b/Re_{\tau }^2 \approx 0.165 R t^+) \approx 119 R$ that is based on our case where $Re_b=5357$ and ${Re_\tau = 180}$. These distances are not unusual for typical experiments; for example, the length of deposition pipe in Liu & Agarwal (1974) is $\approx 160R$. The slight right-shifted symbols in figure 1 correspond to long average up to $t^+=3600$. It is clear that the DNS results for the absorbing wall with one or two particle forces (Sets $2$ and $4$) are in closest agreement with the experimental data of Liu & Agarwal (1974). As discussed in the previous paragraph, reflecting walls show a slightly different behaviour especially for the large $St^+$ cases where particle bouncing or ejection from the wall is common, thus leading to a reduced $V_{dep}^+$.

Recall that deposition velocity is simply the normalised particle flux $J$ at the wall. Knowledge of $J$ as function of $r$, on the other hand, would not only enable us to estimate concentration $c$ at the wall (i.e. $V_{dep}^+$) but also the full distribution of $c$ with $r$ as well as its variation in time. In the next section we address several aspects of this modelling that will provide further insight into particle migration and later also allow us to estimate $c(r;t)$.

4. Transport equation modelling

Once the general form for flux of particles $J$ (number of particles crossing per unit area per unit time) is given by (1.2), conservation of particle number leads to the general continuity equation:

(4.1)\begin{equation} \frac{\partial c}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{J}} = 0.\end{equation}

Although ${\boldsymbol {J}}$ is expressed in a generalised vector form here, in our case, the flux is only considered in $r$ direction, i.e. ${\boldsymbol {J}}$ has only one component $J$. (The vector form (4.1) is convenient for numerical finite-volume methodology used later in the paper.) Using our expression for flux (1.2) and expanding (4.1) in cylindrical coordinates,

(4.2)\begin{equation} \frac{\partial c}{\partial t} = \overbrace{\underbrace{-D_{1} \frac{1}{r}\frac{\partial}{\partial r} (rc)}_\text{Term C1} \underbrace{-c\frac{\partial D_1}{\partial r}}_\text{Term C2} \\ }^{\text{Turbophoretic contribution}} \overbrace{\underbrace{+ D_2\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial c}{\partial r}\right)}_\text{Term G1} \underbrace{+ \frac{\partial D_2 }{\partial r}\frac{\partial c}{\partial r}}_\text{Term G2} }^{\text{Diffusive contribution}}.\end{equation}

Therefore, (4.2) establishes a partial differential equation that contains two spatially varying parameters, $D_1$ and $D_2$, encapsulating the tendency for the inhomogeneous turbulence to transport particles towards the wall and the effect on transport due to wall-normal average concentration gradients, respectively. Interpretation of these Eulerian terms are as follows: $C1 \textrm { and } G1$ represent particle migration associated with fixed coefficients whereas $C2 \textrm { and } G2$ represent the effect of first-order variations in $D_1$ and $D_2$. As mentioned before, one of our aims is to estimate the efficacy of the model for particle migration constructed in this paper, i.e. (4.2). The first step towards this is a direct determination of the coefficients $D_1$ and $D_2$ from the current DNS. Later, we compare the DNS determined $D_1$ and $D_2$ with the existing turbulence-based models.

4.1. Determining turbophoretic and diffusive coefficients with the DNS database

To determine the turbophoretic and diffusive coefficients $D_1$ and $D_2$, respectively, the particle fluxes are measured across imagined cylindrical surfaces in the domain. The radii of these cylinders vary from a minimum of $r^{+} = 4$ to a maximum of $r^{+} = 178$ with an increment of $\Delta r^{+} = 6$. The expression for particle flux is given by

(4.3)\begin{equation} J = \frac{n_p}{A\Delta t}, \end{equation}

where $n_p$ is the number of particles that crosses the surface, $A$ is the area of the cylindrical surface and $\Delta t$ is the time interval over which subsequent snapshots of particle interval are recorded. Note that $J$ is positive in the positive $r$ direction, i.e. when particles migrate towards the wall. Here the concentration $c(r)$ is evaluated by $c(r) = N_{p,l} / ({\rm \pi} L_z (r_{o}^2 - r_{i}^2))$ where $N_{p,l}$ is the local number of particles in a thin annulus with length, $L_z$ (the pipe length), with outer and inner radius $r_o$ and $r_i$, the thickness of which approach each other in the limit. The concentration is evaluated at twice as many points as the flux is measured, for a typical calculation, $r_o^+ - r_i^+ \approx 3$, and the concentration gradient is obtained through finite differencing. To obtain $D_1$ and $D_2$, (1.2) is arranged, such that

(4.4)\begin{equation} \frac{J}{c(r)} = D_1(r) - D_2(r) \frac{1}{c(r)} \frac{\partial c(r)}{\partial r}.\end{equation}

Determination of $D_1$ and $D_2$ follows a statistical approach through the simultaneous observations of particle flux $J(r)$, concentration $c$ and concentration slope $\partial c/\partial r$. To evaluate $J$, the radial positions of individual particles are recorded at sequential snapshots in time. Particles whose linearly interpolated trajectory cross a surface are determined as either wall-directed or wall-fleeing, and a net particle number flux is determined. Now that we know $J$, $c$ and $\partial c / \partial r$, we linearly regress (4.4) through ensembles of the known terms to determine the two remaining unknowns $D_1$ and $D_2$. One thousand equally spaced instances are collected spanning $3600$ viscous time units; measured instances are graphed, and a line of regression drawn.

Figure 7 left and right panels respectively show $D_1$ and $D_2$ (with symbols as in table 2). In general, value of $D_1$ increases with increasing $St^+$, which is consistent with the increased role of turbophoresis in particle deposition with increasing $St^{+}$. The values of $D_2$ on the other hand show only mild dependency on $St^{+}$, again consistent with our notion of turbulent diffusive action. Before further discussion, we report some analytical estimates for $D_1$ and $D_2$ from the properties of the underlying turbulent fluid. Subsequently, we compare the analytical estimates with DNS in figure 7.

Figure 7. Empirically determined coefficients $D_1$ (left panel) and $D_2$ (right panel), corresponding respectively to turbophoretic drift and gradient-driven diffusivity. (a,c,e,g) Solid grey lines, $\widetilde {D_1}$ from the analytical study of Guha (1997), i.e. (4.9); blue solid and dashed lines, (4.6) (with ${u}^{\prime 2}_{pr, rms}$ from DNS) and (4.6) using (4.7) to estimate ${u}^{\prime 2}_{pr, rms}$ (with $u_{r, rms}^{\prime 2}$ from DNS). (b,d, f,h) Solid grey lines, $\tilde {D}_2/10$ from Guha's (1997) estimates, where $\widetilde {D_2} = \nu _t$; dash-dotted grey lines, $\tilde {D}_2/10$ from (4.10) with the $\nu _t$ expression from Guha (1997); and dashed lines, $\tilde {D}_2/10$ from (4.10) with $\nu _t$ from present DNS. All grey lines are divided by a numerical value of $10$ to be shown on the same vertical axis limit. Blue lines indicate $\tilde {D}_2$ from (4.11), whereas dark/light blue are where $f(St)$ is taken for elastic/absorbing walls. Inset included for $St^{+} = 0.1$ and $1$ to aide judgement of agreement between the models and present data.

4.2. Analytical estimation of turbophoretic and diffusive coefficients, and comparisons with DNS results

Theoretical efforts to determine the coefficients in (4.4) (e.g. Caporaloni et al. 1975; Guha 1997; Zhao & Wu 2006) have broadly suggested that the expression for wall-directed flux can be written as

(4.5)\begin{equation} J = \widetilde{D_1} c - \widetilde{D_2}\frac{\partial c}{\partial r},\end{equation}

where $\widetilde {D_2}$ is the particle eddy diffusivity, and $\widetilde {D_1}$ is the turbophoretic velocity. Note that (1.2) and (4.5) are similar, however, in (1.2) the coefficients ($D_1,D_2$) are evaluated from a DNS whereas in (4.5), which follows the work of others, the coefficients ($\widetilde {D_1},\widetilde {D_2}$) are connected to the underlying turbulent velocity field. Additional terms that can be included in (4.5) (but do not form a part of this study) are, for example, the effect of Brownian diffusion, gravitational settling or electrostatic force.

A model for turbophoresis, $\widetilde {D_1}$ was first given by Caporaloni et al. (1975) as (also see (1.1))

(4.6)\begin{equation} \widetilde{D_1} ={-} \tau_p \frac{d {u}^{\prime 2}_{pr, rms}}{d{r}}. \end{equation}

As the particle properties are difficult to model, the interest is to relate them to the underlying fluid properties. The first step is to estimate the parameter, $\mathfrak {R}$, which is defined as the ratio of mean square of the particle velocity fluctuations to the mean square of fluid velocity fluctuations:

(4.7)\begin{equation} \mathfrak{R} = \frac{u^{\prime 2}_{pr, rms}}{u_{r, rms}^{\prime 2}}.\end{equation}

If we know $\mathfrak {R}$, and given that estimates of $u_{r, rms}^{\prime 2}$ are easier to obtain, we can find $u_{pr, rms}^{\prime 2}$ and hence $\widetilde {D_1}$. An experimental correlation for the $\mathfrak {R}$ parameter given by Binder & Hanratty (1991) (which is inspired by the original works of Friedlander 1957) is

(4.8a,b)\begin{equation} \mathfrak{R} = \frac{1}{1 + 0.7(\tau_p / \tau_L)} \quad \text{and} \quad \tau_L = \frac{\nu_t}{u_{r, rms}^{\prime 2}}, \end{equation}

where $\tau _L$ is the Lagrangian time scale for the fluid (see Johansen 1991) and $\nu _t$ is the eddy viscosity (defined in the usual manner as $\nu _t \equiv -\overline {u_r^{\prime }u_x^{\prime }}/(\textrm {d} U_x/\textrm {d} r)$). It can be observed that in the limit of small $\tau _p$ that represents the case where particles faithfully follow the fluid eddies, $\mathfrak {R}$ is indeed equal to unity. This completes one empirical formulation for $\widetilde {D_1}$.

Guha (1997), on the other hand, starts with Eulerian equations for particle continuity and momentum equation, and then uses the Reynolds decomposition for mean and turbulence; after neglecting a number of terms mostly related with triple correlations, they arrive at an equation for the mean wall-normal particle velocity, identified with $\widetilde {D_1}$. The resulting ODE is

(4.9)\begin{equation} \widetilde{D_1} \frac{{\rm d} \widetilde{D_1}}{{\rm d} r} + \frac{\widetilde{D_1}}{\tau_p} ={-}\frac{{\rm d} u_{pr,rms}^{\prime 2}}{{\rm d} r}.\end{equation}

Note the similarity of this with (4.6), except that the particle acceleration term is absent in (4.6). For $u_{pr,rms}^{\prime 2}$, Guha (1997) uses (4.7) and (4.8a,b) along with the analytically fitted experimental correlations for $\nu _t$ and $u_{r,rms}^{\prime 2}$ from Davies (1966) and Kallio & Reeks (1989), respectively.

For comparison of $\widetilde {D_1}$ with DNS, we first write the ODE (4.9) in finite difference form and we then iteratively solve. The necessary boundary conditions for the solution are no flux through the pipe wall and at the centreline. The solution to (4.9) is plotted in figure 7 (left panel) using solid grey lines. The more basic expression arising from (4.6) is plotted in ‘blue’ solid lines with ${u}^{\prime 2}_{pr, rms}$ from DNS, whereas ‘blue’ dashed lines correspond to (4.6) but now using (4.7) to estimate ${u}^{\prime 2}_{pr, rms}$ (with $u_{r, rms}^{\prime 2}$ from DNS). Significant differences are only apparent for $St^+ = 100$ particles.

The empirical findings obtained from DNS and our flux-splitting method (using (4.4)) strongly corroborate the analytically obtained turbophoretic velocity for smaller particles of $St^+ = 0.1$ and $St^+ = 1$. This comparison validates the analytical treatments via (4.9) and (4.6) that have been used for many decades in modelling particle-laden flows. The empirical and analytical profiles for $St^+ = 10$ are in modest agreement but there is a significant gap between the profiles for $St^+ = 100$ particles. These results suggest that the expression for $D_1$ from (4.6) that comes basically from (1.1) is essentially in agreement with the DNS, except of course for $St^+ = 100$ particles. Possible causes for the anomalous behaviour of $St^+ = 100$ particles are discussed in § 5, but briefly, they are related to the fact that Lagrangian auto-correlation based integral time scale for $St^+ = 100$ particles in the radial direction is much longer compared with the corresponding time scale for smaller $St^+$ particles.

Furthermore, it is intriguing to compare the $D_1$ profiles from the DNS with absorbing and elastically rebounding boundary conditions. Physically, this corresponds to a pipe with and without an adhesive surface. Particles of $St^+ = 0.1,1$ and $10$ show little difference between their turbophoretic velocity with these boundary conditions, suggesting that the wall collision effects are not important. However, for the case of $St^+ = 100$ particles (cf. figure 7 bottom left panel), there is a fairly dramatic reduction in the magnitude of $D_1$ for elastically rebounding particles. We now turn to the discussion of $D_2$.

It appears that unlike $\widetilde {D_1}$, there is relatively less understanding with regard to modelling $\widetilde {D_2}$. For example, Guha (1997) simply takes the particle eddy diffusivity $\widetilde {D_2}$ to be equal to the turbulent eddy viscosity $\nu _t$. Zhao & Wu (2006), however, incorporate the effect of finite particle relaxation time and use (a relation attributed by them to Hinze 1975)

(4.10)\begin{equation} \frac{\widetilde{D_2}}{\nu_t} = \left(1 + \frac{\tau_p}{\tau_L} \right)^{{-}1}. \end{equation}

The DNS and the analytical estimates (in grey lines) for $D_2$ are compared in figure 7 (right panel). Here the analytical estimates of $D_2$ are divided by a factor of $10$ so that they can be shown in figures with same vertical axes. Solid grey lines are from Guha's (1997) estimates $\widetilde {D_2} = \nu _t$, whereas dash-dotted and dashed grey lines come from (4.10) with the $\nu _t$ expression from Guha (1997) and with $\nu _t$ from present DNS, respectively. It is clear that $\widetilde {D_2} = \nu _t$ is not appropriate, and even though (4.10) furnishes a slightly better comparison with the DNS (for larger $St^+$) in terms of magnitude, they are still off, and more importantly the peak location of $D_2$ in the analytical estimates is not at the correct wall-normal location.

To delve a little deeper into this issue of $D_2$, note that the form (4.10) can be obtained by starting from (1.1), where $D_2 \sim v^2 \tau _p \sim u_{pr,rms}^{\prime 2} \tau _p$. Using (4.7) and (4.8a,b) leads to $D_2 \sim \mathfrak {R} u_{r,rms}^{\prime 2} \tau _p$, and if one takes (in absence of any realistic data) the eddy viscosity $\nu _t \sim u_{r,rms}^{\prime 2} \tau _p$, it will result in $D_2 \sim \nu _t \mathfrak {R}$, which is the functional form of (4.10). Note that the least realistic assumption here seems to be $u_{r,rms}^{\prime 2} \tau _p \sim \nu _t$, where one relates the turbulent flow property $\nu _t$ directly to the particle time-scale $\tau _p$. We, therefore, revert to the basic expression: $D_2 \sim u_{pr,rms}^{\prime 2} \tau _p$, or $D_2^+ (= D_2/\nu ) \sim u_{pr,rms}^{\prime 2 +}\, St^+$ and, more generally,

(4.11)\begin{equation} D_2^+= u_{pr,rms}^{\prime2 +} f(St^+),\end{equation}

where $f(St^+)$ is a function of $St^+$. It is evident from (4.11) that, for a fixed $St^+$ particle, the form of wall-normal distribution of $D_2^+$ is completely specified by $u_{pr,rms}^{\prime 2 +}$ except for a constant $f(St^+)$, and this seems appropriate when we compare $D_2^{+}$ data in figure 7 (right panel) and $u_{pr,rms}^{\prime 2 +}$ profiles in figures 3(d) and figures 18 and 19. Furthermore, by examining the DNS data in figure 7 (right panel) it is obvious that both the absorbing-wall sets as well as the elastic ones behave similarly, and the variation is only between absorbing and elastic cases. Therefore, using our DNS database we find that $f(St^+)$ follows approximately a power law with $St^+$. Specifically, for absorbing walls, $f(St^+) \approx 0.62 {St^+}^{2/3}$, whereas for elastic walls $f(St^+) \approx 0.7 {St^+}^{3/4}$. The light/dark blue solid lines in figure 7 (right panel) correspond to (4.11) with these $f(St^+)$ for absorbing/elastic walls. At present these power-law variations are merely empirical. The simplistic (4.11) is by no means perfect, but it does capture the peak position in $D_2^+$ and the wall-normal variation of lower $St^+$ cases reasonably well. The $St^+=100$ case is problematic, as has been the case in modelling $D_1^+$ also, and it is likely due to similar reasons.

Even if we have perfect models for $D_1$ and $D_2$, it is not clear whether the time-dependent concentration $c(t)$ can be predicted with the phenomenological model (4.2). If $c(t)$ from (4.2) compares well with the DNS, then we can use the reduced (4.2) (compared with the full Navier–Stokes, continuity and particle equations) to shed light on the key features of particle migration.

4.3. Time-dependent solution of the phenomenological concentration equation (4.2) and comparisons with the DNS

For the numerical solution procedure of (4.2) the physical pipe domain is represented as a 1D function of the distance from the wall, $y^+ = R^+ - r^+$; then ($D_1, D_2$) coefficients determined in § 4.1 are substituted into (4.2) and time marched. Details of the solution procedure are described in Appendix C. The simulations for all $St^{+}$ are initialised with a uniform concentration, i.e. $c/c_b = 1$ at $t = 0$. (Recall that $c_b$ is the bulk concentration at any instant in time.) The time evolution of $c(r)$ are presented in figures 8(a) to 8(d) for $St^+=0.1$, $1$, $10$ and $100$, respectively. Here we restrict our focus to simulation Set 1 (drag force only, elastic walls), but additional results are presented in Appendix C.

Figure 8. Predictive capacity for the simplified model against the observable from the DNS, corresponding to simulation Set 1: (a) $St^{+}=0.1$, (b) $St^{+}=1$, (c) $St^{+}=10$ and (d) $St^{+}=100$. Solid blue line indicates model prediction and light blue shading corresponds to the region bounding the 95$\%$ confidence interval for determined $D_1$ and $D_2$. The five different time instances correspond to $t^{+} = 0\unicode{x2013} 3600$ with intervals of ${\Delta t^{+} = 720}$.

In general, the major qualitative features observed from the DNS (in symbols) are well captured in this model (solid blue lines), particularly the negligible and low wall accumulation at lower $St^+$ values and dramatic increase in wall accumulation for particles with significant inertia. The shaded (light-blue) regions around the model profiles (in solid blue lines) are uncertainty estimates based on the regressed values of $D_1$ and $D_2$, and further details of the error estimation are presented in Appendix C1. Favourable agreement between our time-marched solution and the time-varying DNS concentration profiles allows us to probe further into the contributing factors to particle migration at the phenomenological level.

As shown in (4.2), $C_1$ and $C_2$ represent two turbophoretic terms whereas $G_1$ and $G_2$ the diffusive terms. Terms $C_2$ and $G_2$ show the effects of radial variations in $D_1$ and $D_2$, respectively. These terms are shown in figures 9(a) and 9(b) for $St^+=1$ and $10$, respectively. The terms are evaluated adjacent to the wall until $y^+=3$, and they are interpreted in terms of accumulation–depletion depending on the flux term signs. Further, the split contributions to the flux divergence are non-dimensionalised by the pipe radius $R$, deposition velocity $V_{dep}$ determined in § 3.3 and bulk concentration $c_b$.

Figure 9. Decomposition of contribution to normalised divergence of fluxes in numerical solution adjacent to the wall for (a) $St^+ = 1$ and (b) $St^+ = 10$. Black terms correspond to convective contributions and blue terms to diffusive contributions, and solid lines correspond to $C_1,G_1$ with dashed lines corresponding to $C_2,G_2$.

We observe that the relative magnitudes of the convective and diffusive contributions for $St^+ = 1$ and $St^+ = 10$ are similar but the magnitude of each term is larger for the more inertial particles suggesting a $St^{+}$ dependence. Higher magnitudes are expected as $St^+ = 10$ particles set up sharper concentration gradients close to the wall and are more likely to enter the viscous sublayer. This analysis suggests that close to the wall, the signs of the individual convection–diffusion mechanisms (and so their role) are similar for $St^+ = 1$ and $St^+ = 10$ particles, but their magnitude grows with $St^{+}$. This near-wall mechanism acts in conjunction with the high-turbophoretic force (via $D_1$) at $y^+ = 15$ which is much larger for $St^+ = 10$ than $St^+ = 1$ particles. Interestingly, the main balance is between the turbophoretic term $C_2$ (and a smaller effect from $C_1$ and $G_1$) that is driving accumulation and the diffusion term $G_2$ depleting the particles concentration. As such, it is the gradients in coefficients $D_1$ and $D_2$ that dominate the net particle migration. This implies that a constant turbophoretic velocity or diffusive coefficient will not capture the particle migration physics appropriately. Away from the wall we observe (but not shown here) that the turbophoretic components $C_1$ and $C_2$ are both negative, which implies that they are responsible for depletion of concentration in that region. Furthermore, the sum of all the fluxes at these locations are also negative, which implies a net reduction in concentration, and this is expected because most particles migrate from these regions towards the wall. This is, of course, in contrast to the near-wall region (see figure 9) where the total flux convergence is positive, suggesting an increase in particle concentration.

5. Physics informing transport equation modelling

Although the preceding analysis constitutes purely Eulerian modelling and facilitates comparison with applied mathematical modelling, particles are perhaps more naturally analysed in a Lagrangian frame of reference. The properties of individual particles can be studied through a Lagrangian description, and subsequently related back to Eulerian modelling. In § 4 we discussed various phenomenological modelling aspects, whereas in § 3 we considered different fluid and particle velocity statistics. In the following we explore three analysis methods that further clarify the dependence of particle dynamics on the underlying coherent turbulent motion allowing us to better understand the concentration statistics and modelling.

5.1. Lagrangian velocity correlation

Lagrangian statistics in homogeneous isotropic turbulence have been studied extensively (e.g. Yeung & Pope 1989), where particles are statistically identical, regardless of their initial position, but the case of wall turbulence is more complicated. Particles in wall flows have statistical dependency on their wall-normal distance. Lagrangian trajectories from the Set 1 DNS in this study are visualised in figure 10(a) (i.e. top panel), where it is apparent that the motion of inertial particles diverges from tracers. Particles are tagged and tracked once they reside within a small interval of the target heights of $y^+ = 5$, 30 and 118. Those tagged from the core, buffer and viscous regions (shown in figure 10(a) in the left, middle and right panels, respectively) show markedly different pathlines between $St^+ = 0.1$ and $St^+ = 100$ (identified by their symbols). In the near-wall region, there is a dramatic reduction in the pathline curvature of $St^+ = 100$ particles.

Figure 10. (a) Trajectories, coloured by streamwise distance traversed, of particles tracked after entering the turbulent core (left panel), buffer (middle panel) and viscous regions (right panel). Particles of different $St^+$ are identified by their respective symbols as indicated in table 2. (b) Wall-normal trajectories of inertial particles, starting from $y^{+} = 5$. Initial position are marked by black dots and final positions by red or blue dots. Trajectories plotted for a duration of $t^{+} = 360$, at which point, blue and red dots indicate particles that are respectively nearer to or further away from the wall than $y^{+} = 30$.

We make the additional observation that particle motion in the buffer region exhibits high curvature, which is due to the high fluid streamline curvature in this region (e.g. Perven, Philip & Klewicki 2021). Particles with low $St^+$ enter this region and with even probability may journey to the pipe center or return to the wall. However, particles with higher $St^+$ are unable to follow these curvy motions and are flung wall-wards resulting in a reduced population that journeys to the pipe centre.

To further illustrate visually the particle movement, in figure 10(b) (i.e. bottom panel) we present wall-normal trajectories of particles (for $St^+=0.1, 1, 10$ and $100$ from left to right panels, respectively). Particles with $St^+=100$ tend to show motions compared with other $St^+$ with less curvature and diminished total pathline distance. These visual features of particle paths can be quantified by using Lagrangian correlations. In general, Lagrangian correlations can be defined by

(5.1)\begin{equation} \rho_{ij}(\tau;y_0) = \frac{\langle v_{pi}^{\prime}(t_0;y_0) v_{pj}^{\prime}(t_0 + \tau;y_0)\rangle}{\langle v_{pi}^{\prime 2}(t_0;y_0) \rangle^{1/2} \langle v_{pj}^{\prime 2}(t_0 + \tau;y_0)\rangle^{1/2}}, \end{equation}

where $i,j$ represent the three coordinate directions, $t_0$ and $y_0$ are the initial time and wall-normal location of the particles, and the fluctuating particle velocity $v_{pi}^{\prime }(t_0 + \tau,y_0) = v_{pi}(t_0 + \tau,y_0) - \langle v_{pi}(t_0 + \tau,y_0)\rangle$ is the particle velocity fluctuation relative to the Lagrangian average over the total averaging time.

In figure 11 the panels in three columns (from left to right) show Lagrangian correlation functions where the particle populations commence their motion at the wall normal heights of $y^+_0=3$, $30$ and $118$, whereas the three rows (top to bottom) correspond to autocorrelation of particle velocity in the streamwise, radial and azimuthal directions $\rho _{xx}$, $\rho _{rr}$ and $\rho _{\theta \theta }$. Of particular interest for understanding the one-dimensional model discussed in § 4.3 is the radial component $\rho _{rr}$. Here, the profiles of the particle populations $St^+ = 0.1$ and $St^+ = 1$ show little difference, which is expected. Increasing in size, the $St^+ = 10$ particle population remains correlated for appreciably longer; but more notable is the order of magnitude increase in correlation period for the $St^+ = 100$ particles, as well as their flatter profiles in the neighbourhood of $\tau ^+=0$. This would suggest that an Eulerian gradient transport model for particle transport is less suitable for $St^{+} = 100$ particles, and this can be further quantified by evaluating several timescales associated with the Lagrangian correlations.

Figure 11. Lagrangian auto-correlation of streamwise ($\rho _{xx}$), wall-normal ($\rho _{rr}$) and azimuthal ($\rho _{\theta \theta }$) particle fluctuating velocity. For symbols see table 2.

In order to quantify the time interval over which particle velocity is correlated with itself, the Lagrangian integral scale associated to each velocity component is estimated according to

(5.2)\begin{equation} T_{ii} = \int^{\tau_c}_0 \rho_{ii} (\tau,y_0)\,{\rm d} \tau, \end{equation}

where $\tau _c$ is the time lag at which the auto-correlation first crosses 0.1 (e.g. Stelzenmuller et al. 2017). This definition is chosen because the usual definition of the integration over an infinite time interval cannot be applied to all profiles, as some profiles cross zero and others do not. Results of $T_{ii}$ for the correlations presented in figure 11 are shown in figure 12. Note that $St^+=100$ data are markedly different to other $St^+$ cases. Of particular interest are values of timescale in the radial direction $T_{rr}$ (middle panel), wherein $St^+=100$ has $T_{rr}^+ = T_{rr}U_\tau (U_\tau /\nu ) \approx 100$. Note that the pipe flow is for $Re_\tau =R (U_\tau /\nu ) = 180$, which implies that $T_{rr}U_\tau \approx 0.55R$, i.e. $St^+=100$ particle motions are correlated over about half $R$. This is perhaps not surprising for these particle which skip most of the small-scale fluid motions. These large correlation times also point towards the difficultly (we encountered in § 4) in modelling $St^+=100$ particles using gradient hypothesis (e.g. Corrsin 1974).

Figure 12. Lagrangian integral scales ($T_{ii}$) corresponding to data in figure 11: (a) $T_{xx}^+$, (b) $T_{rr}^+$ and (c) $T_{\theta \theta }^+$. For symbols see table 2.

Notwithstanding the fact that timescale discussion sheds light on the modelling efforts, the physical movement of the particles is governed by the turbulent motions specific to wall flows. These coherent turbulent motions are discussed in the following section focusing on different wall regions, which can also help to clarify some observations, such as in the buffer region the empirically determined $D_1$ reaching a peak value, the dramatic difference for $St^+=100$ data and the effect of elastic/absorbing wall conditions.

5.2. Coherent motions and flow topology

We spatially correlate the particle wall-ward velocity in the buffer region with the turbulent coherent ‘sweep’ and ‘ejection’ motions (e.g. Willmarth & Lu 1972; Brodkey, Wallace & Eckelmann 1974; Robinson 1991). These turbulent coherent motions are made up of strongly coherent ($Q2$ and $Q4$) events and additionally interaction-type ($Q1$ and $Q3$) events which remain correlated over significantly less time. Classification into quadrants for turbulent events constitutes signage of turbulent fluctuations $u^{\prime }$ and $v^{\prime }$ where $Q2$ ejections are low-speed streaks ($u^{\prime } < 0$) moving away from the wall ($v^{\prime } > 0$), and $Q4$ sweeps are high-speed streaks moving towards the wall. Note that in pipe flows the wall-normal velocity $v = -u_{r}$. Particle transfer has previously been found to be controlled by strongly coherent sweeps and ejections which entrain particles (e.g. Marchioli & Soldati 2002). Visualisation of coherent motions (as line contours) with entrained particles (as dots) from simulation Set 1 are plotted in figures 13 and 14 for $Q2$ and $Q4$ events, respectively. These plots answer the question: if a particle is wall-directed (or wall-fleeing), how likely is it that they are entrained within a strongly coherent turbulent event? Particle points within the $Q2$/$Q4$ events indicate that the particles are correlated to, and hence likely transported by, these events. Strongly coherent $Q2$ and $Q4$ motions are able to transfer $St^+ = 1$ particles towards the wall, but also capable are $Q1$ and $Q3$ events. Inertial $St^+ = 10$ particles particles, however, require strongly coherent motions to be transferred towards the wall as the momentum of these particles is sufficient to filter $Q3$ events. Coherent motions are still able to transfer $St^+ = 100$ particles wall-ward but this mechanism is less effective.

Figure 13. Instantaneous overlay at $y^+ = 30$ between coherent $Q2$ ejections ($u_x^{\prime } < 0,u_r^{\prime } < 0$) and particles: (a) $St^+ = 0.1$, (b) $St^+ = 1$ (c) $St^+ = 10$ and (d) $St^+ = 100$. Contours levels for $Q2$ range from $(u_x^{\prime } < 0) \times (u_r^{\prime } < 0) = 0.001$–$5.0$. Particles plotted with instantaneously wall-fleeing velocity ($u_{pr} < 0$) and residing at $y^+ = [29,31]$.

Figure 14. Instantaneous overlay at $y^+ = 30$ between coherent $Q4$ sweeps ($u_x^{\prime } > 0,u_r^{\prime } > 0$) and particles: (a) $St^+ = 0.1$, (b) $St^+ = 1$, (c) $St^+ = 10$ and (d) $St^+ = 100$. Contours levels range from $(u_x^{\prime } > 0) \times (u_r^{\prime } > 0) = 0.001$–$5.0$. Particles plotted with instantaneously wall-ward velocity ($u_{pr} > 0$) and residing at $y^+ = [29,31]$.

To quantify the connection between particles and $Q2$/$Q4$ events illustrated in figures 13 and 14, we calculate probability of particle radial velocity direction ($u_{pr}$) conditioned on $Q2$ (where fluid velocity $u_x^{\prime } < 0,u_r^{\prime } < 0$) as well as on $Q4$ ($u_x^{\prime } > 0,u_r^{\prime } > 0$). (Recall that wall-normal velocity $v=-u_r$.) These conditional probabilities (as percentages) for varying $St^+$ are presented in figures 15(a) and 15(b) for $Q2$ and $Q4$ event, respectively. The statistics are collected over the entire time of the simulation within the wall-normal range of $y^+ = [25,35]$. The large conditional probability in figure 15(a) for all particles except $St^+=100$ suggest that indeed the coherent ejection motion is responsible for upward particle velocity (and wall-fleeing effect of particles). Furthermore, figure 15(a) also shows that modelling of the wall as either elastic or absorbing has little effect on these statistics when we consider the $St^+ = 0.1,1$ and $10$ particles, but a dramatic effect on the sluggish $St^+ = 100$ particles. This suggests a memory effect for $St^+ = 100$ pertaining to their approach towards $y^+ = [25,35]$ from either the wall, or from the pipe center. For the $St^+ = 100$ Stokes-drag-only cases, $Pr(u_{pr} < 0 | Q2) = 0.71$ for elastic walls and decreases to $Pr(u_{pr} < 0 | Q2) = 0.53$ for absorbing walls, both of which are smaller than for the probabilities for smaller $St^+$ suggesting these particles retain their attained inertia and can coast across strong vortical regions. To re-emphasise the point of conducting four distinct simulation sets, with toggled inclusion of the Saffman–Mei lift force $\boldsymbol {F}_l$, and wall-boundary type, we note that in the buffer region, for $St^+ = 0.1,1 \ \textrm {and} \ 10$ particles, there is negligible ($<1\,\%$) dependence of the conditional statistics on $\boldsymbol {F}_l$. Our findings are consistent with the general finding in wall flows where the Saffman–Mei lift force $\boldsymbol {F}_l$, is most prominent close to walls and decays approaching the centreline (e.g. Wang et al. 1997). For the $St^+ = 100$ particles, however, our simulations show that in the case of absorbing walls, inclusion of the Saffman–Mei lift force yields a moderate ($9\,\%$) reduction in the $Pr(u_{pr} < 0 | Q2)$ statistic between Set 2 and Set 4. In the case of elastic walls, we again see a reduction of $9\,\%$ in $Pr(u_{pr} > 0 | Q4)$ between Set 1 and Set 3. For sweep events, figure 15(b) shows characteristics similar to ejection events in figure 15(a). It seems, however, that sweep events (compared with ejection) events are more efficient in bringing even larger $St^+$ particles to the wall.

Figure 15. Probability of $(\textit {a})$ particle fleeing given it is sampling a Q2 event ($Pr(u_{pr} < 0 | Q2)$) or $(\textit {b})$ directed towards ($Pr(u_{pr} > 0 | Q4)$) pipe walls given sampling of a Q4 event, conditioned on the strongly coherent event types: $Q2 (u_x^{\prime } < 0,u_r^{\prime } < 0)$ and $Q4 (u_x^{\prime } > 0,u_r^{\prime } > 0)$. Statistics collected at $y^+ = [25,35]$ from simulation Set 1. For symbols see table 2.

Although the preceding coherent motions analysis sufficiently explains particle transfer mechanisms, study of local organisation of particles is needed to give insight into biased sampling of fluid events that lead to wall-ward transfer. We employ the flow classification of Chong, Perry & Cantwell (1990) that maps flow topology with respect to invariants of the velocity gradient tensor, $\boldsymbol {A} = \boldsymbol {\nabla } \boldsymbol {u}$. The invariants are $P = A_{ii}=0$ (owing to continuity), $Q = -(1/2)A_{ij}A_{ji}$ and $R = -(1/3)A_{ij}A_{jk}A_{ki}$. For data presented on a $Q$–$R$ plot, large positive and negative $Q$ indicate vorticity and strain dominated regions, respectively; in addition, for $Q<0$ a positive $R$ suggests bi-axial strain or sheet-like structures (e.g. Davidson 2015). In figure 16, we plot four decades of the joint probability density function (PDFs) for fluid events (in contour lines) separated by wall-normal region in the pipe (in three columns of figure 16) and then overlay the events sampled by inertial particles from Set 1 (shown by coloured scatter plots). In the $Q$–$R$ plane, there is little deviation between the conditional PDFs of particles able to trace flow events, $St^+ = 0.1$ and the fluid. There is a slight ‘down-shift’ in the sampling of $Q$ events for particles of comparable time scale to the small fluid scales, $St^+ = 1$, when they reside in the buffer and logarithmic regions (cf. figures 16e and 16f). This is suggestive of particles avoiding vortical regions in preference for strain-dominated regions. Another feature that we observe is a dramatic reduction for $St^+ = 10$ and $100$ events to sample more extreme events in $Q$–$R$ space. See also particle-laden channel DNS by Rouson & Eaton (2001). This is more prominent in the viscous sublayer than in the buffer/logarithmic regions. The consequences, first, are that once particles reach the viscous sublayer, they organise outside of the streamwise vortices, and prefer regions of low streamwise velocity rather than the high-speed streaks and, second, particles that do rebound from the wall must pass through the vortex cores present in the buffer region if they are to continue towards the centreline of the pipe. As they approach the vortex cores, however, they are likely to be flung back towards the wall region.

Figure 16. Joint PDF of $Q$–$R$, conditionally sampled by the fluid (contour lines) and particles (coloured scatter points). Solid grey line corresponds to the discriminant, $D = \frac {27}{4}R^2 + Q^3 = 0$. Left-to-right columns, viscous sublayer, buffer region and logarithmic layer. Top-to-bottom rows, $St^+ = 0.1,1,10$ and $100$.

6. Summary and conclusions

Four sets of particle-laden turbulent flow DNS have been generated in a smooth-wall pipe with particles spanning four decades of $St^{+}$ number ($0.1$, $1$, $10$ and $100$). The datasets first address common point-particle modelling of forces and particle–wall-type interaction. It was observed that inclusion of Saffman–Mei lift force $\boldsymbol {F}_l$ and presence of absorbing or elastically rebounding walls significantly influenced particle dynamics for the largest $St^{+}$ number studied. The variance between datasets contextualises primary findings of this paper that establish a simple representation of turbophoresis physics. This representation was formed by employing a phenomenological flux model informed by DNS-driven statistical determination of turbophoresis and diffusive contributions to particle migration. These findings support the modelled turbophoretic velocity arising originally from Caporaloni et al. (1975) and Reeks (1983) and further quantified by the Reynolds averaging of the particle mass and momentum equations of Guha (1997). Specifically there is broad agreement in the modelled and DNS-determined peak turbophoretic velocity; in addition, DNS-driven profiles of turbophoresis for $St^{+} = 100$ do not continue to dramatically increase beyond the $St^{+} = 10$ profiles, supporting Guha (1997) which included a particle acceleration term within modelling to yield an inertia-moderated deposition regime. The diffusive terms, however, were found to be an order of magnitude lower than that suggested by all available models, and based on the original equations of Reeks (1983), we suggest another model.

Subsequently, turbophoretic and diffusive fluxes are employed in a 1D advection– diffusion equation for particle concentration, which when time-marched demonstrates reasonable agreement with DNS data. Expanding the terms in this simple representation of flux facilitates further delineating turbophoretic and diffusive contributions responsible for accumulation of particle concentration at walls. The dominant contribution was found to arise from gradients in turbophoretic profiles dispelling conjecture that turbophoresis may in smooth wall pipes be adequately modelled by a wall effect or constant turbophoretic term alone.

Targeted investigation of relevant particle-flow interaction physics that underpins turbophoresis gave insight into $St^{+}$ number dependency for trajectories and correlations, Reynolds stress turbulent events and flow topology. The specific linkage of these flow physics to turbophoresis are respectively: (i) suitability of Eulerian models that assume instantaneously local fluid events determine particle statistics; (ii) distinguishes whether mechanisms responsible for migration are agnostic to wall type and inclusion of Saffman–Mei lift force $\boldsymbol {F}_l$; and (iii) sheds light on preferential concentration throughout the domain and in near-wall ejection. Characterising particles by comparison with tracer-like $St^{+} = 0.1$ particles revealed persistent Lagrangian velocity correlations for $St^{+} = 100$, which suggest possible breakdown of gradient diffusion hypothesis in modelling particles in turbulent wall flows for larger $St^+$. Concerning particle accumulation at walls, chiefly important is radial auto-correlation, larger particles were distinct as their Lagrangian integral scale was an order of magnitude greater than tracer scales. Unsurprisingly, for $St^{+} = 100$ particles a dramatic reduction was observed in the efficacy of 3D mechanisms in wall turbulence and Reynolds stress turbulent events in migrating particles towards the wall. This stresses an important complexity for turbophoresis modelling, for moderate $St^{+}$ there is a diminished efficacy of Reynolds stress momentum fluxes in transporting particles wall-wards but an increased turbophoretic velocity. Turbophoresis is, at least in part, informed by an aggregate of motions rather than sampling of fluid events at any single wall-normal location.

Lastly, focusing on quantification of particle–eddy interaction through classification of biased flow topology sampling revealed particles with increased inertia organise close to the wall preventing efficient re-entrainment into the core. Two themes emerge, although turbophoretic particle migration reaches a peak value in the buffer region, particle deposition fate is inherently influenced by sequential events sampled along Lagrangian trajectories and, second, future particle modelling efforts may require delineation of particles that are of $St^{+} = O(100)$ as these large ballistic particles sufficiently differ in their fluid–particle interaction from the remaining $St^{+}$ particles in this study.

For the present findings which are determined in a pipe flow, similarity of turbulent fluctuations profiles with other geometries, e.g. channels, should guide relevancy of the results presented here. However, it remains as a separate study to gauge the difference in mechanisms pertaining to channels and pipes, particularly in horizontal flows where gravity will drive particles wall-ward. Of particular interest would be modification to ejection mechanisms that depend on over-sampling of boundary layer low-speed strata, and the ability of a particle to traverse streaks.