1. Introduction
Wall-bounded turbulent flows laden with inertial particles abound in environmental and industrial contexts, such as the transport of particulate matter in the atmosphere, sediment transport in rivers, the separation of fine particles within industrial cyclones, and fluidized bed reactors. The chaotic and multiscale nature of the fluid turbulence coupled with the particle dynamics results in fascinating phenomena, which, however, are still challenging to understand and model (Balachandar & Eaton Reference Balachandar and Eaton2010; Brandt & Coletti Reference Brandt and Coletti2022).
In the past decades, many studies have been devoted to the dynamics of turbulent wall transport of small inertial particles, driven by their prevalent nature and rich physics. Traditionally, these systems have been classified by how the dispersed phase influences the overall flow behaviour (Elghobashi Reference Elghobashi1994). Precisely, one-way coupling (1WC) corresponds to very small particle loadings, so small that their influence on the dynamics of the suspending fluid turbulence is negligible; two-way coupling (2WC) denotes flow regimes where mass loading is high enough such that the particles modify flow observably, but particle–particle interactions are negligible; finally, regimes where both particle–particle and particle–fluid interactions influence the flow dynamics are grouped into the four-way coupling category.
In first-principles, fully resolved direct numerical simulations (DNS) of a particle-laden flow, the flow around each particle needs to be resolved explicitly (Balachandar & Eaton Reference Balachandar and Eaton2010; Maxey Reference Maxey2017) (so-called particle-resolved DNS, PR-DNS). While this approach is free from modelling assumptions for the dispersed phase dynamics, it is computationally expensive due to the explicit imposition of no-slip and no-penetration boundary conditions at the surface of many particles moving in a turbulent medium. This is incredibly challenging when particles are tiny and there is a scale separation between the particle size and the smallest (Kolmogorov) turbulence scale due to extreme resolution requirements. Fortunately, in this case, one may be able to resort to the so-called point-particle approximation (PP-DNS), where interphase coupling is considered to be localized to a point.
Indeed, the point-particle approximation has been the method of choice for simulating the particle dynamics in turbulent wall flows. In these cases, it is assumed that the local properties of an undisturbed flow at the particle position drive the dispersed-phase dynamics (Gatignol Reference Gatignol1983; Maxey & Riley Reference Maxey and Riley1983). In the case of highly inertial particles (i.e. large particle-to-fluid density ratios), the so-called Maxey–Riley–Gatignol equations simplify to a drag force term (Arcen, Tanière & Oesterlé Reference Arcen, Tanière and Oesterlé2006), which nevertheless yields highly non-trivial particle dynamics, even in isotropic turbulence (Toschi & Bodenschatz Reference Toschi and Bodenschatz2009). In wall-bounded turbulent flows, the inhomogeneous turbulence results in even richer dynamics, with the particle distribution driven by the interplay between small-scale clustering, turbophoresis, and the interaction between the particles and near-wall turbulence structures (Reeks Reference Reeks1983; Soldati & Marchioli Reference Soldati and Marchioli2009; Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012), resulting in very inhomogeneous particle concentrations peaking at the wall, with strong preferential sampling of low-speed regions, as reproduced in a plethora of numerical studies, such as Fessler, Kulick & Eaton (Reference Fessler, Kulick and Eaton1994), Uijttewaal & Oliemans (Reference Uijttewaal and Oliemans1996), Marchioli et al. (Reference Marchioli, Giusti, Salvetti and Soldati2003, Reference Marchioli, Soldati, Kuerten, Arcen, Taniere, Goldensoph, Squires, Cargnelutti and Portela2008), Kuerten (Reference Kuerten2006), Soldati & Marchioli (Reference Soldati and Marchioli2009), Bernardini, Pirozzoli & Orlandi (Reference Bernardini, Pirozzoli and Orlandi2013) and Jie et al. (Reference Jie, Cui, Xu and Zhao2022).
When 2WC effects are important, point-particle DNS must describe the back-reaction of the dispersed phase in the flow. This is a major challenge, as the point-particle dynamics is driven by a local undisturbed fluid velocity, while the local flow field is being disturbed by the particles (see e.g. Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2015). The classical approach, known as the particle-in-cell method, was developed by Crowe, Sharma & Stock (Reference Crowe, Sharma and Stock1977) and is widely used – even in the present study – but it requires a sufficiently high number of particles per grid cell, and cannot reproduce simple benchmarks of a sedimenting sphere in a quiescent medium. Indeed, approaches for a consistent treatment are being actively investigated; see Gualtieri et al. (Reference Gualtieri, Picano, Sardina and Casciola2015), Horwitz & Mani (Reference Horwitz and Mani2016, Reference Horwitz and Mani2020), Ireland & Desjardins (Reference Ireland and Desjardins2017) and Horwitz et al. (Reference Horwitz, Iaccarino, Eaton and Mani2022). Investigations of particle-laden turbulent flows in the two-way-coupling regime are found in e.g. Vreman et al. (Reference Vreman, Geurts, Deen, Kuipers and Kuerten2009), Zhao, Andersson & Gillissen (Reference Zhao, Andersson and Gillissen2010), Kuerten (Reference Kuerten2016), Richter & Sullivan (Reference Richter and Sullivan2014), Capecelatro, Desjardins & Fox (Reference Capecelatro, Desjardins and Fox2018), Wang & Richter (Reference Wang and Richter2019) and Battista et al. (Reference Battista, Gualtieri, Mollicone, Salvadore and Casciola2023).
While employing the point-particle approximation for small inertial particles is physically sound, validating the fidelity of the approximation in 1WC and 2WC conditions remains a challenge. Experimental data are available (Eaton & Fessler Reference Eaton and Fessler1994; Kaftori, Hetsroni & Banerjee Reference Kaftori, Hetsroni and Banerjee1995). Still, there are few parameter-matched numerical studies due to limitations in terms of Reynolds number and well-controlled experiments, and only in recent years have efforts in this direction started to appear (Wang et al. Reference Wang, Fong, Coletti, Capecelatro and Richter2019). Fortunately, PR-DNS of small inertial particles in turbulence has become possible, thanks to the continuous growth in available computer power and development of efficient numerical methods, with the first direct comparisons between point-particle models and particle-resolved simulations starting to appear for forced homogeneous isotropic turbulence and decaying HIT with moving particles (Schneiders, Meinke & Schröder Reference Schneiders, Meinke and Schröder2017a,Reference Schneiders, Meinke and Schröderb; Fröhlich et al. Reference Fröhlich, Schneiders, Meinke and Schröder2018; Mehrabadi et al. Reference Mehrabadi, Horwitz, Subramaniam and Mani2018), and turbulent channel flow (Horne & Mahesh Reference Horne and Mahesh2019; Costa, Brandt & Picano Reference Costa, Brandt and Picano2020a).
These results from PR-DNS of particles in wall-bounded turbulence confirm that a sole drag force may not suffice to accurately reproduce the particle dynamics even for relatively small particles with a large density ratio under 1WC conditions. Where the shear rate is high near the wall, lift forces are also important. This has already been suggested in early works using PP-DNS, since the work by McLaughlin (Reference McLaughlin1989) (see also Botto et al. Reference Botto, Narayanan, Fulgosi and Lakehal2005). In this regard, while Arcen et al. (Reference Arcen, Tanière and Oesterlé2006) reported that the lift force has a negligible impact on the dispersed phase statistics, Marchioli, Picciotto & Soldati (Reference Marchioli, Picciotto and Soldati2007) and Shin et al. (Reference Shin, Portela, Schaerer and Mangiavacchi2022) found that the inclusion of the lift force in PP-DNS can lead to weakened near-wall particle accumulation in upward and horizontal channels, respectively. McLaughlin (Reference McLaughlin1989) showed that the inclusion of the lift force resulted in a higher deposition rate; Wang et al. (Reference Wang, Squires, Chen and McLaughlin1997) found that neglecting the lift force resulted in a slight reduction in the deposition rate.
While these findings are not necessarily contradictory, as there are some variations in the governing parameters in the different studies and lift force models, many questions remain elusive. (1) Which form of lift force model is appropriate for reproducing with high fidelity the dynamics of small inertial particles in a turbulent wall flow? (2) Under which conditions do lift forces matter in the particle dynamics? (3) What are the consequences of choosing an inaccurate lift force model in dispersed phase dynamics and near-wall turbulence modulation?
In the recent direct comparison between PR-DNS and PP-DNS for small particles in turbulent channel flow by Costa et al. (Reference Costa, Brandt and Picano2020a) (see also Costa, Brandt & Picano Reference Costa, Brandt and Picano2020b), it was shown that lift models are vital for reproducing the near-wall particle dynamics, at least for inertial particles with a size of the order of one viscous wall unit. Surprisingly, a modified Saffman lift model predicted the particle dynamics perfectly, with the Saffman force scaled by a normalized shear rate, while conventional lift models yield poorer predictions (Costa et al. Reference Costa, Brandt and Picano2020b). In addition to reducing the near-wall concentration peak, it was also seen that lift force causes a large increase in correlated streamwise–wall-normal particle velocity fluctuations. This quantity has dramatic drag-increasing effects for sufficiently high mass loading, as it modulates the streamwise momentum budget; see Costa, Brandt & Picano (Reference Costa, Brandt and Picano2021). This direct link between lift force and drag increase makes their accurate modelling crucial.
In the present work, we address the three questions presented above by performing DNS of turbulent channel flow laden with small inertial particles, using the point-particle approximation, with three different lift models, from the classical Saffman lift model to the one that perfectly predicts previous PR-DNS data. We consider different inner-scaled particle sizes 
$D^+\sim 1\unicode{x2013}0.1$ in 1WC and 2WC conditions for a Stokes number that is known to feature strong wall accumulation and preferential concentration in low-speed regions. Our analysis shows that currently available lift models are bound to fail near the wall for small inertial particles due to a vanishing particle-to-fluid slip velocity. We then use 2WC simulations to illustrate how different lift force models can result in qualitatively different turbulence modulation (i.e. turbulence attenuation versus turbulence enhancement). Near-wall accumulation is still significantly reduced for the smallest particle size (
$D^+=0.1$). Still, the lift-induced increase of correlated velocity fluctuations near the wall is negligible, thus lift force has few consequences in turbulence modulation.
This paper is organized as follows. Next, in § 2, the numerical method, lift force models and computational set-ups are described. Then in § 3, the effects of lift force models on the particle dynamics and near-wall accumulation are investigated using 1WC DNS, based on which we try to explore the reason for the failure of conventional lift force models. Following this, we demonstrate qualitatively the effect of lift force models on turbulence modulation and momentum transfer with 2WC DNS. Finally, conclusions are drawn in § 4.
2. Methodology
2.1. Governing equations and numerical method
The fluid phase is governed by the incompressible Navier–Stokes equations,
\begin{equation} \left.\begin{gathered} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \\ \rho(\partial_t \boldsymbol{u} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u})={-}\boldsymbol{\nabla} p +\mu\,\nabla^2\boldsymbol{u} + \boldsymbol{f}, \end{gathered}\right\} \end{equation}
where 
$\boldsymbol {u}$ denotes the fluid velocity vector, 
$p$ is the fluid pressure, 
$\mu$ is the dynamic viscosity, 
$\rho$ is the fluid density, and 
$\boldsymbol {f}$ is the particle feedback force to the fluid phase in the case of 2WC point-particle simulations, here computed by using a standard particle-in-cell approach, which spreads the particle force to the nearest Eulerian grid points with a linear kernel (Boivin, Simonin & Squires Reference Boivin, Simonin and Squires1998; Lee & Lee Reference Lee and Lee2019; Zhang, Cui & Zheng Reference Zhang, Cui and Zheng2023).
We consider 
$x$, 
$y$ and 
$z$ as the streamwise, wall-normal and spanwise directions, respectively. No-slip and no-penetration boundary conditions are specified at the domain walls, while periodicity is applied along the streamwise and spanwise directions. These equations are discretized in space using a pseudo-spectral approach along 
$x$ and 
$z$ and second-order finite differences along 
$y$. Wray's low-storage third-order Runge–Kutta scheme is employed for time marching (see Wray Reference Wray1990).
Spherical particles with density 
$\rho _p \gg \rho$ and particle diameter 
$D$ in the absence of gravity are tracked with a standard Lagrangian point-particle method, with their dynamics governed by
\begin{equation} m_p \dot{\boldsymbol{U}}_p = \boldsymbol{F}_d+\boldsymbol{F}_l,\quad \dot{\boldsymbol{X}}_p = \boldsymbol{U}_p, \end{equation}with the Schiller–Naumann drag force
\begin{equation} \boldsymbol{F}_d ={-}3{\rm \pi}\mu D \boldsymbol{U}_s(1+0.15\,{\textit{Re}}_p^{0.687}), \end{equation}
where 
$\boldsymbol {U}_p$, 
$\boldsymbol {X}_p$, 
$\boldsymbol {F}_d$ and 
$\boldsymbol {F}_l$ are the particle velocity, position, drag and lift forces, 
$m_p$ is the particle mass, 
$\boldsymbol {U}_s=\boldsymbol {U}_p-\boldsymbol {u}|_{\boldsymbol {x}=\boldsymbol {X}_p}$ is the local slip velocity evaluated at the particle location, and 
${\textit {Re}}_p$ 
$=|\boldsymbol {U}_s|\,D / {\nu }$ is the particle Reynolds number. Since the particle-to-fluid density ratio is high, other dynamic effects such as added mass, fluid acceleration, and the Basset history force are negligible. While in most practical scenarios of particle-laden wall transport, gravity is important at high density ratios, we neglect it in the present work. We do this to isolate the interplay between more intricate wall accumulation mechanisms (e.g. turbophoresis) and lift forces, in a flow that could otherwise feature significant settling effects. Finally, for the sake of simplicity, a perfectly elastic hard-sphere rebound is employed for particle–wall interactions.
2.1.1. Lift force models
We consider standard shear-induced lift force models where lift force 
$\boldsymbol {F}_l$ acts perpendicular to the local shear sampled by the particle and points along 
$\boldsymbol {\omega }\times \boldsymbol {U}_s$, where 
$\boldsymbol {\omega } = (\boldsymbol {\nabla } \times \boldsymbol {u} )|_{\boldsymbol {x}=\boldsymbol {X}_p}$ is the undisturbed flow vorticity evaluated at the particle position. For convenience, let us express the lift force in terms of the dimensionless lift coefficient 
$C_L$, defined as
\begin{equation} C_L = \frac{\boldsymbol{F}_l}{\tfrac{1}{8}{\rm \pi}\rho D^2\,|\boldsymbol{U}_s|^2}\boldsymbol{\cdot} \frac{\boldsymbol{\omega}\times\boldsymbol{U}_s}{|\boldsymbol{\omega}\times\boldsymbol{U}_s|}. \end{equation}
Assuming that the Oseen length 
$\ell _u=\nu /{| \boldsymbol {U}_s |}$ is much larger than the Saffman length 
$\ell _\omega =\sqrt {\nu /{| \boldsymbol {\omega } |}}$, i.e. 
$\varepsilon \gg 1$ with 
$\varepsilon =\sqrt {| \boldsymbol {\omega }|\,\nu } /| \boldsymbol {U}_s |$, an explicit lift solution can be derived:
\begin{equation} C_L = \frac{18}{{\rm \pi}^2}\,\varepsilon\,J(\varepsilon), \end{equation}with
\begin{equation} J(\varepsilon\gg1) =J^\infty=2.255. \end{equation}
Hereafter, the lift expression (2.5) together with the 
$J$ function by (2.6) will be referred to as the Saffman model, which is expected to be valid in the double limits 
${\textit {Re}}_p\to 0$ and 
${\textit {Re}}_\omega \to 0$ provided that 
$\varepsilon \gg 1$, where 
${\textit {Re}}_\omega =| \boldsymbol {\omega } |\,D^2/\nu$ is the shear Reynolds number.
Still in the double limits 
${\textit {Re}}_p\to 0$ and 
${\textit {Re}}_\omega \to 0$, the 
$J$ function in (2.5) at finite 
$\varepsilon$ turns out to be a volume integral in Fourier space. Its value cannot be put in closed form, but was estimated numerically in Asmolov (Reference Asmolov1989), McLaughlin (Reference McLaughlin1991) and more recently, Shi & Rzehak (Reference Shi and Rzehak2019). Based on these numerical data, various empirical correlations of 
$J(\varepsilon )$ were proposed (see Shi & Rzehak (Reference Shi and Rzehak2019) for a comprehensive review), with the most commonly used being the one proposed by Mei (Reference Mei1992), i.e.
\begin{equation} J(\varepsilon) =0.3\,J^{\infty} \left(1+\tanh \left[\tfrac{5}{2}(\log _{10} \varepsilon+0.191)\right]\right) \left(\tfrac{2}{3}+\tanh (6 \varepsilon-1.92)\right). \end{equation}
Hereafter, the lift expression (2.5) together with the 
$J$ function by (2.7) will be referred to as the Mei model, which is expected to be valid irrespective of 
$\varepsilon$, provided that 
${\textit {Re}}_p$ and 
${\textit {Re}}_\omega$ are small.
The flow is assumed to be unbounded for the two lift models introduced above. In particle-laden channel flows, near-wall accumulation of particles is often observed, as discussed above. Typically, the peak in the particle concentration appears within the viscous sublayer, i.e. for 
$y^+ \leq 5$, wherein the wall effect on the lift is crucial if 
$D^+ \geq 1$ (Balachandar & Eaton Reference Balachandar and Eaton2010; Shi et al. Reference Shi, Rzehak, Lucas and Magnaudet2021). In this context, the explicit lift solution can be derived if the wall lies in the ‘inner region’ of the flow disturbance in the low-
${\textit {Re}}_p$ limit, specifically, if the separation between the particle and the wall, 
$\ell$, is much smaller than the inertial lengths (
$\ell \ll \min {(\ell _u, \ell _\omega )}$). The corresponding lift solutions (see Shi & Rzehak (Reference Shi and Rzehak2020) for a comprehensive review) take the general form
\begin{equation} C_L = A + B\,Sr + C\,Sr^2, \end{equation}
where 
$Sr= | \boldsymbol {\omega } |\,D/| \boldsymbol {U}_s |$ is the normalized shear rate, and 
$A$, 
$B$ and 
$C$ are pre-factors that are independent of 
$Sr$. Hereafter, the terms proportional to 
$Sr$ and 
$Sr^2$ will be referred to as the linear and quadratic inner contributions 
$C_{L, {in}}^{linear}$ and 
$C_{L,{in}}^{quad}$, respectively.
Situations where the wall is in the ‘outer region’ of the flow disturbance (
$\ell \gg \max {(\ell _u, \ell _\omega )}$), while still in the limit 
${\textit {Re}}_p\to 0$, were considered in Asmolov (Reference Asmolov1990), McLaughlin (Reference McLaughlin1993) and Takemura, Magnaudet & Dimitrakopoulos (Reference Takemura, Magnaudet and Dimitrakopoulos2009). Assuming that 
$Sr={\textit {Re}}_\omega ^{1/2}\,\varepsilon \leq O(1)$ (i.e. the same condition for the shear rate as considered in Saffman Reference Saffman1965), the outer-region lift approaches the Saffman solution (2.5), namely 
$C_{L,{out}}^{linear}\propto \varepsilon$. Together with the scaling in the inner region that 
$C_{L,{in}}^{linear}\propto Sr$, it appears that the inner-region lift contribution transitions into the outer region by a pre-factor 
$\varepsilon /Sr$. If one assumes that the quadratic contribution follows the same transition, then it may be speculated that 
$C_{L,{in}}^{quad}$ scales as 
$\varepsilon \,Sr$ in the outer region, i.e. 
$C_{L,{out}}^{quad} \propto \varepsilon \,Sr$. For the particle-laden channel flow,
\begin{equation} Sr\approx\frac{(u_\tau^2/\nu)D}{u_s}=\left(\frac{u_s}{u_\tau}\right)^{{-}1}D^+, \end{equation}
where 
$u_s$ is the streamwise mean slip velocity, namely 
$u_s=\langle \boldsymbol {U}_s\boldsymbol {\cdot }\boldsymbol {e}_x\rangle$, with 
$\boldsymbol {e}_x$ the unit vector along the streamwise direction, and 
$u_\tau$ is the conventional wall friction velocity. As will be demonstrated in figure 5 in § 3.1, 
$u_s$ might change its sign in the buffer layer, leading to extremely large values of 
$Sr$ in the near-wall region.
The analysis above implies that there is also a quadratic lift contribution 
$C_{L}^{quad} \propto \varepsilon \,Sr$, which might dominate the lift generation in the inner wall region. This is consistent with the form of the lift force proposed in Costa et al. (Reference Costa, Brandt and Picano2020a,Reference Costa, Brandt and Picanob), where
\begin{equation} C_L = \frac{18}{{\rm \pi}^2}\,\varepsilon\,Sr\,J^\infty, \end{equation}
with 
$J^\infty =2.255$ according to Saffman (Reference Saffman1965). Despite its simplified form, this correlation aligns well with their PR-DNS results. Hereafter, the lift model by (2.10) will be referred to as the CBP model.
Finally, it should be noted that in practice, to avoid singularities in the numerical calculation, we implemented a slight variant of the lift coefficient as described in (2.4), where the quotient on the right-hand side to define the unit vector is modified to 
$\boldsymbol {\omega }\times \boldsymbol {U}_s/|\boldsymbol {\omega }|\,|\boldsymbol {U}_s|$. In practice, the results are not very sensitive to this choice, as 
$\boldsymbol {\omega }\times \boldsymbol {U}_s$ is nearly aligned with the wall-normal direction.
2.2. Validation
The present code has demonstrated successful applications in wall-bounded particle-laden flows, including open channel flow (Wang & Richter Reference Wang and Richter2019; Gao, Samtaney & Richter Reference Gao, Samtaney and Richter2023) and planar Couette flow (Richter & Sullivan Reference Richter and Sullivan2014; Richter Reference Richter2015) loaded with inertial particles, but the lift force is neglected. Hence, for completeness, we present the validation of the 1WC PP-DNS code, including the lift force models. Specifically, we reproduce in figure 1 the PP-DNS of particle-laden turbulent channel flow in the 1WC regime using these three lift models, as reported in Costa et al. (Reference Costa, Brandt and Picano2020a,Reference Costa, Brandt and Picanob) in the 1WC regime. This corresponds to case CL1 in table 1. The agreement is excellent, which validates our implementation.
Figure 1. (a) Normalized particle volume fraction 
$\phi /\varPhi$. Outer-scaled second-order moments of particle velocity: (b) streamwise velocity root mean square, (c) wall-normal velocity root mean square, (d) Reynolds shear stress profile. The points denote the corresponding lift-included 1WC PP-DNS from Costa et al. (Reference Costa, Brandt and Picano2020b). Here, 
$U_b$ is the flow bulk velocity.
Table 1. Physical parameters used in the present point-particle DNS campaign: 
$N_p$ denotes the total number of particles, and 
$D^+$ is the inner-scaled particle diameter. The particle density ratio was adjusted to keep the same value of the Stokes number 
$St^+=\tau _{p}/(\nu /u_\tau ^2) = 50$ and bulk solid mass fraction 
$\varPsi _m = 0.337$ for all cases.
2.3. Computational parameters
We perform point-particle DNS of channel flow at 
${\textit {Re}}_{\tau }=180$ in a computational domain 
$(L_x,L_y,L_z) = (6h,2h,3h)$, where 
$h$ is the half-channel height. The domain is discretized on 
$(N_x,N_y,N_z) = (160,320,160)$ grid points; the grid is slightly stretched to refine the near-wall resolution, corresponding to grid spacing 
$(\varDelta _x^+,\varDelta _y^+,\varDelta _z^+) = (6.75,0.9,3.375)$, with the conventional ‘
$+$’ superscript denoting viscous wall scaling. Note that the grid spacing near the wall is comparable to that of the largest particle sizes. All simulations are carried out at a constant time-marching step fixed at 
$\Delta t^+ = 0.1$ (normalized by viscous unit 
$\nu /u^2_{\tau }$), which corresponds to 
$\textrm {CFL} \approx 0.4$. The DNS results confirm this time step to be sufficiently small that particles and fluid elements could not pass through a grid cell per time step (Zheng, Feng & Wang Reference Zheng, Feng and Wang2021). The total simulation time is 
$T \approx 450h/u_\tau$ (approximately 
$80\,000$ viscous time scales 
$\nu /u_\tau ^2$), which is long enough to guarantee converged statistics (Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012).
We consider three different set-ups for varying particle diameter, while 
$St^+ = \tau _{p}/(\nu /u_\tau ^2)=50$ is fixed by varying the particle density. These parameters are described in table 1. This target particle Stokes number was chosen since it is known to feature highly inhomogeneous particle distributions in wall turbulence; see Sardina et al. (Reference Sardina, Schlatter, Brandt, Picano and Casciola2012). Moreover, it should be noted that a value 
$D^+ =3$ corresponds to about one Kolmogorov length scale in the channel bulk, an order of magnitude that is often investigated in the literature using PP-DNS (see e.g. Bernardini Reference Bernardini2014; Motoori, Wong & Goto Reference Motoori, Wong and Goto2022; Zhang et al. Reference Zhang, Nathan, Tian and Chin2022). We performed 1WC and 2WC simulations at a fixed mass fraction for different particle sizes and lift force models. Naturally, in the 1WC simulations, the mass fraction is not a governing parameter, and the high number of particles ensures statistical convergence of the results. In the 2WC PP-DNS, the bulk mass fraction is fixed to 
$\varPsi _m \approx 30\,\%$, to ensure significant turbulence modulation and highlight the effects of lift force model choice in the flow statistics.
3. Results and discussion
3.1. Particle dynamics
We start by describing the dynamics of the dispersed phase under the effect of lift force models using a simple 1WC approach. This allows us to measure observables sampled by the particle that dictate the validity region of the lift models (e.g. 
$Sr$, 
$\varepsilon$ and 
${\textit {Re}}_p$). As we will see, current lift models are bound to fail in the buffer layer.
Figure 2 shows the profiles of the normalized mean particle concentration for all particle sizes. When the lift force is not taken into account, the preferential accumulation of particles in the wall region is significant, which has been observed widely (Picano, Sardina & Casciola Reference Picano, Sardina and Casciola2009; Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012; Lee & Lee Reference Lee and Lee2015; Gao et al. Reference Gao, Samtaney and Richter2023; Gualtieri et al. Reference Gualtieri, Battista, Salvadore and Casciola2023), and termed as turbophoresis (Reeks Reference Reeks1983; Johnson, Bassenne & Moin Reference Johnson, Bassenne and Moin2020). However, the particles show less wall accumulation when subjected to lift forces, and this effect is more pronounced the larger the particle size. The Saffman model yields a much smaller peak, while the other two models feature smaller deviations, with the highest peak corresponding to the Mei model. When the particles are very small (
$D^+=0.1$), the differences in the lift force models become almost indistinguishable, with all the models robustly showing the same reduction of wall accumulation. Yet the difference in accumulation in the absence of lift is clear.
Figure 2. Local particle volume fraction 
$\phi$, normalized by the bulk value 
$\varPhi$, as a function of the wall-normal distance (1WC): (a) CL1, (b) CM1, (c) CS1. The insets show the normalized particle volume fraction versus the inner-scaled wall distance.
To qualitatively describe the lift force effect on the spatial localization of wall accumulation, figure 3 shows instantaneous snapshots of streamwise velocity fluctuations in a wall-parallel plane within the buffer layer (
$y^+=12$) and in a plane of constant streamwise location (
$x^+=540$), along with the corresponding particle positions for the case with largest particles, CL1. As expected, the particles show strong spatial localization and inhomogeneous distribution, with larger local density occurring in elongated clusters (Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012). Effects of the lift force model on the gross characteristics of the particle distribution are readily discernible. One notable manifestation is the difference in particle numbers in the wall-parallel (
$x$–
$z$) plane. Consistent with figure 2, the lift-free model yields the highest number of particles, while the Saffman model yields the lowest. Conversely, the Mei and CBP models fall somewhere between these two. This indicates the weakening accumulation of particles near the wall by the lift force, consistent with observations by Marchioli et al. (Reference Marchioli, Picciotto and Soldati2007) and Shin et al. (Reference Shin, Portela, Schaerer and Mangiavacchi2022). The other manifestation concerns small-scale particle clustering. At first glance, lift-free, Mei and CBP simulations show a more significant tendency to over-sample low-speed regions than the Saffman one. This tendency is checked in figure 4, showing the probability density functions (p.d.f.s) of streamwise fluid velocity fluctuations sampled by the particles within the 
$y^+ \le 12$ region. For larger particle sizes, particles under the Saffman lift force sample a less focused range of streamwise velocity fluctuations than other models for larger particle sizes. In particular, for 
$D^+=3$ and 
$1$, the CBP and Mei cases show approximately the same tendency to sample low-speed regions, differently than the Saffman model and comparable to the lift-free case for 
$D^+=3$. Preferential accumulation gets more focused with decreasing particle size, with the lift-free case producing the strongest preferential accumulation for 
$D^+=1$ and 
$0.1$.
Figure 3. Inner-scaled instantaneous streamwise velocity fluctuation and the corresponding instantaneous snapshots of particle locations in the 
$x$–
$z$ plane (
$y^+=12$, left-hand column) and 
$y$–
$z$ plane (
$x^+=540$, right-hand column) (case CL1, 1WC). The wall-parallel plane shows all particles with inner-scaled wall normal position 
$Y_p^+ \le 12$. (a) Lift-free, (b) Saffman model, (c) Mei model, (d) CBP model.
Figure 4. Probability density functions (p.d.f.s) of inner-scaled streamwise fluid velocity fluctuations conditioned at the particle location within the 
$y^+ \le 12$ region (1WC): (a) CL1, (b) CM1, (c) CS1.
Figure 5 shows the inner-scaled profiles of streamwise mean particle velocity 
$\langle u_p \rangle$ and local slip velocity 
$\langle u_s \rangle$. Compared with the lift-free particles with 
$D^+ \ge 1$, lift forces increase 
$\langle u_p \rangle$ in the viscous sublayer and buffer layer, while the differences between cases become negligible beyond 
$y^+ \sim 30$. Once more, particles under the Saffman lift force show the highest variation with respect to the lift-free case, yielding the highest 
$\langle u_p \rangle$; particles under the Mei lift force yield the lowest, and the ones under the CBP lift force fall in between. Consistently, for the smallest particle sizes, the overall modification of 
$\langle u_p \rangle$ by lift forces is quite small. This tendency of lift forces becoming less important with decreasing particle size is expected – as we will discuss later in more detail (§ 3.3), the relative importance of lift to drag forces close to the wall approximately scales with the inner-scaled particle diameter 
$D^+$.
Figure 5. Inner-scaled streamwise mean particle velocity and slip velocity profiles (1WC): (a,d) CL1, (b,e) CM1, (c,f) CS1. The horizontal dashed lines denote 
$\langle u_s \rangle = 0$.
At the wall, while the fluid velocity must vanish, particles can feature a mean apparent slip where they flow faster than the fluid due to their inertia (Zhao, Marchioli & Andersson Reference Zhao, Marchioli and Andersson2012). This is shown clearly in the particle mean slip (figures 5d–f). At larger wall distances, the tendency of particles to over-sample slower-than-average fluid velocity (Kiger & Pan Reference Kiger and Pan2002; Baker & Coletti Reference Baker and Coletti2021) is also clearly reflected in the mean particle slip velocity. A similar trend can be observed in the corresponding streamwise mean slip velocity profiles; however, the difference in 
$\langle u_s \rangle$ is larger than 
$\langle u_p \rangle$, and this is attributed to the negative apparent slip velocity (
$\langle u |_{\boldsymbol {x}=\boldsymbol {X}_p} \rangle - \langle u \rangle < 0$), which reflects preferential sampling of slower-than-average fluid. Interestingly, lift forces increase the apparent slip velocity near the wall. This can be understood from the sign of the lift force in this region, which points towards the wall due to the positive slip velocity. Particles sampling higher-momentum regions will be driven towards the wall by the lift force, and their inertia results in a higher mean slip.
The sign of 
$\langle u_s \rangle$ changes in approximately the same location of the buffer layer, irrespective of the particle size and lift force model. This mechanism seems relatively robust and has been observed in numerous other studies with inertial particles at different Stokes numbers and different flow Reynolds numbers (see e.g. Mortensen et al. Reference Mortensen, Andersson, Gillissen and Boersma2008; Zhao et al. Reference Zhao, Marchioli and Andersson2012; Wang & Richter Reference Wang and Richter2020; Gao et al. Reference Gao, Samtaney and Richter2023). It will be shown that this change in sign makes lift force modelling in wall turbulence extremely challenging (see discussions in § 3.3).
Figure 6 shows the 
$x$–
$y$ component inner-scaled velocity covariances for the particle and fluid phases, respectively. The spanwise velocity variances are hardly affected by lift forces, and thus not shown. As expected, away from the wall where the mean shear is low, the particle statistics closely follow those of the lift-free cases. Near the wall, instead, the lift force enhances the streamwise and wall-normal particle velocity variances (
$\langle u'^2_p \rangle$ and 
$\langle v'^2_p \rangle$) and Reynolds stress (
$\langle u'_p v'_p \rangle$), even for the smallest particle size (
$D^+=0.1$). Enhancement of wall-normal velocity fluctuations near the wall has indeed been reported in previous numerical and experimental studies (Fong, Amili & Coletti Reference Fong, Amili and Coletti2019; Costa et al. Reference Costa, Brandt and Picano2020a). Indeed, as also shown in Costa et al. (Reference Costa, Brandt and Picano2020a), this enhancement of velocity fluctuations near the wall cannot be reproduced when lift forces are not considered in the particle dynamics. It seems that lift force is very effective at correlating streamwise and wall-normal velocity fluctuations near the wall – under high shear rate, small variations in streamwise velocity naturally translate to changes in vertical acceleration through the lift forces, which in turn induce wall-normal velocity fluctuations. Streamwise velocity fluctuations will also be amplified naturally, since streamwise and wall-normal velocity fluctuations are correlated through the mean shear. As we will see in § 3.4, this amplification of near-wall velocity fluctuations by lift forces has major consequences on turbulent modulation and drag changes in the flow at high particle mass loading.
Figure 6. Second-order moments of mean particle (solid lines) and fluid (dashed lines) velocity (1WC). (a–c) Inner-scaled streamwise and (d–f) wall-normal velocity variances, and (g–i) Reynolds shear stress profiles, where (a,d,g), (b,e,h) and (c,f,i) denote cases CL1, CM1 and CS1, respectively.
3.2. Mechanism for near-wall accumulation
To understand the role of the lift force in altering the near-wall particle distributions, it is essential to examine three primary mechanisms responsible for wall-normal inertial particle transport. The first is the lift-induced migration, which is the direct consequence of the inertial lift force. The second is turbophoresis (Reeks Reference Reeks1983), causing particle migration from regions of higher to lower turbulence intensity. In wall-bounded flows, turbophoresis results in strong particle accumulation in the viscous sublayer, where turbulent fluctuations vanish. Note, however, that in most cases, the peak in the particle concentration appears not at the wall but at a distance of 
$O(D)$ from the wall (Marchioli & Soldati Reference Marchioli and Soldati2002; Costa et al. Reference Costa, Brandt and Picano2020a), owing to the wall–particle collisions as well as hydrodynamic wall–particle interactions (Goldman, Cox & Brenner Reference Goldman, Cox and Brenner1967; Vasseur & Cox Reference Vasseur and Cox1977; Zeng, Balachandar & Fischer Reference Zeng, Balachandar and Fischer2005). Finally, there is a biased sampling effect, due to the tendency of particles to sample ejection regions in the buffer layer featuring negative streamwise velocity and positive (i.e. repelling) wall-normal velocity fluctuations (Marchioli & Soldati Reference Marchioli and Soldati2002; Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012). As a result, biased sampling usually provides a net slow particle drift away from the wall. The interaction among the three mechanisms is sketched in figure 7 and will be elaborated below.
Figure 7. Illustration of the interplay between different mechanisms responsible for the wall-normal particle transport.
Lift-induced migration is directed along 
$\boldsymbol {\omega }\times \boldsymbol {U}_s$, which in the case of a turbulent wall flow is predominantly in the wall-normal direction. According to figures 5(d–f), particles lead the liquid flow for 
$y^+\lesssim 10$, and lag behind the flow further away from the wall. Hence the lift-induced migration is towards the wall in the viscous sublayer, and acts cooperatively with the turbophoresis effect. The way these two mechanisms interact in the viscous sublayer may seem surprising at first glance, as figure 2 clearly indicates a suppression of the near-wall particle accumulation by lift force. The reasoning behind this suppression may be understood by noting the following two issues. First, the lift force is towards the channel centre within the buffer layer and beyond 
$y^+\approx 10$. The resulting outward migration, together with the strong turbulence intensity in the buffer layer, tends to entrain particles in this region outwards, leading to a net decrease in the concentration in the region 
$y^+\lesssim 10$ as revealed by figure 2. Second, the lift force influences the wall-normal particle transport by turbophoresis. This phenomenon is linked to the particle wall-normal velocity variance 
$\langle v'^2_p \rangle$ (Reeks Reference Reeks1983; Johnson et al. Reference Johnson, Bassenne and Moin2020). Indeed, in the limit of vanishing particle Reynolds number, the particle wall-normal momentum balance at steady state yields a turbophoresis pseudo-force proportional to 
$\mathrm {d}\langle v'^2_p \rangle /\mathrm {d}y$, which creates the migration of particles down gradients in particle wall-normal velocity variance (Sikovsky Reference Sikovsky2014; Johnson et al. Reference Johnson, Bassenne and Moin2020). Recall figures 6(d–f), where the variation of 
$\langle v'^2_p \rangle$ with the wall distance is depicted. Clearly, the lift force tends to increase the ‘apparent inertia’ of the particles, enhancing the deviation between the particle velocity variance and that of the fluid within the viscous sublayer. This leads to an attenuation in the corresponding turbophoresis pseudo-force. In particular, for the two cases where 
$D^+\geq 1$ (figures 6d,e), the prediction with the Saffman lift model yields 
$\mathrm {d}\langle v'^2_p \rangle /\mathrm {d}y<0$ in the viscous sublayer, indicating that the turbophoresis pseudo-force is outwards, driving the particles away from the wall. Consequently, the lift force competes with the turbophoresis pseudo-force in these two cases, leading to a highly flattened near-wall peak in the particle fraction, as seen in figures 2(a,b). In contrast, they cooperate in predictions using the Mei and the CBP models, as figure 2(a) shows. There, the lift force compensates the attenuation in turbophoresis, making the negative slopes in 
$\phi (\kern0.7pt y^+)$ from results using Mei and the CBP models closely follow the lift-free case (highlighted in the figure insets). Note, however, that the magnitude of the Mei lift (CBP lift) is proportional to 
${\textit {Re}}_\omega ^{1/2}$ (
${\textit {Re}}_\omega$), which approaches 
$D^+$ (
${D^+}^2$) in the viscous sublayer. Consequently, for the two cases where 
$D^+\leq 1$, the lift-induced migration is weak and incapable of compensating the still-pronounced attenuation in turbophoresis revealed by figure 6(e,f). This is why in these two cases the slopes in 
$\phi (\kern0.7pt y^+)$ from results employing the Mei and CBP models are smaller in magnitude than in the lift-free case. Still, the difference in accumulation with respect to the lift-free case for smaller sizes is clear, due to near-wall particles residing for long times near the wall while being subjected to small but persistent lift forces.
Given the discussion above, it is worth examining in more detail the link between the chosen lift model and the corresponding modulation in the turbophoresis pseudo-force, particularly to what extent this modulation relates to the lift force magnitude. Combining the turbophoresis pseudo-force (figures 6d–f) and lift force as shown in figure 8, a monotonic increase in the suppression of the turbophoresis pseudo-force with increasing lift is observed at 
$D^+ \leq 1$, while no such trend is evident at 
$D^+=3$. In the latter case, the suppression of the turbophoresis pseudo-force is most pronounced in the prediction using the Saffman model, even though the magnitude of the corresponding lift force is only half that of the prediction using the CBP model for 
$y^+\lesssim 5$. Hence the differences in peak concentration are directly linked not to the lift force magnitudes, but rather to the modulation in the turbophoresis pseudo-force induced by these models. In the present work, the turbophoresis pseudo-force is always suppressed by the presence of lift, leading to a decrease in the peak concentration. Hence the precise relationship between the extent of this suppression and the lift force magnitude remains unclear, making it challenging to directly connect the chosen lift models with the changes in the predicted concentration profiles. Finally, it would be worth investigating these dynamics at lower Stokes numbers, as both turbophoresis and preferential sweeping will be less pronounced, and the nature of the interaction between lift and turbulence may change drastically. Given the increased issues with lift force modelling at lower Stokes number, discussed in the next subsection, addressing this question should be the object of a future, dedicated study.
Figure 8. The inner-scaled lift-induced acceleration (i.e. normalized by 
$u^3_{\tau }/\nu$) profiles along the wall-normal direction (1WC): (a) CL1, (b) CM1, (c) CS1. The horizontal dashed lines denote 
$\langle F_{l,y} \rangle = 0$.
To gain insights into biased sampling of sweeps/ejection events under lift forces, we carried out a quadrant analysis for the fluid flow experienced by the particles. In the 
$u^{\prime }\unicode{x2013} v^{\prime }$ plane, with positive 
$v^{\prime }$ directed away from the wall, the flow experienced by the particle is categorized into four types of events: first quadrant events (Q1), characterized by outward motion of high-speed fluid, with 
$u^{\prime } > 0$ and 
$v^{\prime } > 0$; second quadrant events (Q2), characterized by outward motion of low-speed fluid, with 
$u^{\prime } < 0$ and 
$v^{\prime } > 0$, which are usually called ejections; third quadrant events (Q3), characterized by inward motion of low-speed fluid, with 
$u^{\prime } < 0$ and 
$v^{\prime } < 0$; and finally, fourth quadrant events (Q4), which represent motions of high-speed fluid towards the wall, with 
$u^{\prime } > 0$ and 
$v^{\prime } < 0$, and are usually called sweeps. Figure 9 shows the joint p.d.f.s in the 
$u^{\prime }\unicode{x2013} v^{\prime }$ plane of the fluid fluctuation seen by ascending (
$v_p >0$) and descending (
$v_p <0$) particles within the 
$y^+ \le 12$ region. The first observation is that a larger proportion of particles tends to move outwards (ascending), irrespective of the presence of lift force. The proportion of ascending particles without lift force decreases as the particle size decreases, whereas particles subject to lift force exhibit an inverse trend. Second, the lift-free ascending particles are highly likely to sample ejections (Q2), while the lift-free descending ones are prone to sample sweeps (Q4) for particles with 
$D^+ \geq 1$. This agrees with experimental observations in channel flow laden with finite-size particles (Kiger & Pan Reference Kiger and Pan2002; Yu, Vinkovic & Buffat Reference Yu, Vinkovic and Buffat2016; Baker & Coletti Reference Baker and Coletti2021). However, the particles under the effect of lift forces tend not to over-sample sweeps, irrespective of the wall-normal particle velocity sign, which indicates that the wall-normal velocity of the particles is less correlated with the fluid one. This also confirms that inertial particles under the effect of lift forces are not prone to slowly drive towards the wall in low-speed regions, because the competing wall-repelling effect of lift forces seems to dominate their dynamics. In close inspection of the high-intensity region of the joint p.d.f.s on the 
$u'$-axis, we note that particles under the Saffman lift force tend to sample a wider range of negative 
$u'$ regions, compared with particles under the Mei and CBP lift forces. This aligns with the p.d.f.s of 
$u'$ as shown in figure 4. Third, irrespective of the presence of the lift force or wall-normal particle velocity sign, the smallest particles (
$D^+=0.1$) over-sample ejection (Q2) and inward motions of low-speed fluid regions (Q3). Sweep events, bringing high-momentum fluid towards the wall, are experienced only weakly by particles. Overall, the percentage of particles sampling Q2 events only slightly surpasses that for Q3, suggesting that biased sampling may play a minor role in driving high-inertia particles (
$St^+=50$) to depart from the inner wall region. This observation aligns with numerical results from Marchioli & Soldati (Reference Marchioli and Soldati2002), where biased sampling is found to be more efficient for transferring particles with smaller inertia.
Figure 9. Joint p.d.f.s of the streamwise and wall-normal fluid fluctuating velocities seen by ascending (
$v_p >0$) and descending (
$v_p <0$) particles within the 
$y^+ \le 12$ region (1WC): (a) CL1, (b) CM1, (c) CS1. The numbers in the images denote the proportion of ascending particles.
3.3. Emerging challenges in lift force modelling
The aforementioned results allow us to address several issues concerning the effects of lift force in turbulent wall flows, which may have been overlooked previously in PP-DNS. Ideally, such discussions would be based on results from PR-DNS. However, in this instance, we will extrapolate from PP-DNS results, in conjunction with the CBP model, as an analogue to PR-DNS outcomes. This approach is justified by findings from Costa et al. (Reference Costa, Brandt and Picano2020a,Reference Costa, Brandt and Picanob), which demonstrate that PP-DNS are capable of satisfactorily reproducing the corresponding PR-DNS results for the most challenging case considered in this study, with 
$D^+=3$.
Let us first address the circumstances under which the effects of lift force could be considered negligible for large inertial particles, and thereby disregarded in PP-DNS. Given that biased sampling does not play a key role in the wall-normal transport of inertial particles, a criterion based on the ratio of lift to the turbophoresis pseudo-force seems sufficient to describe the relative significance of the lift force. In this context, the dimensionless particle size 
$D^+$ seems a promising candidate, for the following two reasons. First, 
${D^+}^2$ approaches the shear Reynolds number 
${\textit {Re}}_\omega$ in the viscous sublayer; the latter directly measures the magnitude of the shear lift irrespective of the choice of lift models (Saffman Reference Saffman1965; McLaughlin Reference McLaughlin1991; Mei Reference Mei1992; Shi & Rzehak Reference Shi and Rzehak2019; Costa et al. Reference Costa, Brandt and Picano2020b). Second, the turbophoresis effects seem to decay with increasing 
$D^+$. This is clearly revealed by comparing the slopes of 
$\langle v'^2_p \rangle (\kern0.7pt y^+)$ in figures 6(d–f), which correspond to results with increasing 
$D^+$. Gluing these two aspects together, it seems that the force ratio of lift to turbophoresis approximately scales as 
${D^+}^k$ with 
$k>0$. Although it is not feasible to determine the exact value of 
$k$, the results summarized in § 3.1 made it clear that by no means should the lift force be neglected for 
$D^+\geq 1$. On the other hand, it is reasonable for the lift force to be neglected for 
$D^+\leq 0.1$, except for the under-prediction of near-wall accumulation discussed above (recall figure 2).
We should now address a major issue in commonly applied lift force models for PP-DNS of turbulent wall flows. As outlined in § 2.1, the shear lift models originating from the pioneering work of Saffman (Reference Saffman1965) are applicable rigorously only in the double limits 
${\textit {Re}}_p\to 0$ and 
${\textit {Re}}_\omega \to 0$. In the viscous sublayer, 
${\textit {Re}}_p\approx D^+(\langle u_s\rangle /u_\tau )\approx D^+$ for 
$D^+\geq 1$ according to figure 5, whereas 
${\textit {Re}}_\omega \approx {D^+}^2$. Consequently, neither of the two limits is satisfied for 
$D^+\geq 1$. In addition to the double limits above, there is a constraint concerning the dimensionless shear rate 
$Sr$. This is made clear by noting that 
$Sr={\textit {Re}}_\omega ^{1/2}\,\varepsilon$ and hence 
$Sr\leq O(1)$ even in the strong shear limit 
$\varepsilon \geq 1$ considered in Saffman (Reference Saffman1965). This is why in the wall-bounded situation, the ‘outer-region’ lift solutions (Asmolov Reference Asmolov1990; McLaughlin Reference McLaughlin1993; Takemura et al. Reference Takemura, Magnaudet and Dimitrakopoulos2009) following the methodology of Saffman (Reference Saffman1965) match the corresponding ‘inner-region’ solutions (Cherukat & McLaughlin Reference Cherukat and McLaughlin1994; Magnaudet, Takagi & Legendre Reference Magnaudet, Takagi and Legendre2003; Shi & Rzehak Reference Shi and Rzehak2020) only for 
$Sr=O(1)$. For particle-laden channel flows with small inertial particles (Mortensen et al. Reference Mortensen, Andersson, Gillissen and Boersma2008; Zhao et al. Reference Zhao, Marchioli and Andersson2012; Wang & Richter Reference Wang and Richter2020; Gao et al. Reference Gao, Samtaney and Richter2023), this constraint is more seriously violated than that for the double limits above, as the slip velocity reverses at 
$y^+\approx O(10)$ (see e.g. figures 5(d–f) in the present work), making 
$Sr\to \infty$ for particles within the buffer layer. This issue is highlighted in figure 10, where we show the p.d.f.s of the inverse of 
$Sr$ (horizontal axis) as functions of the wall-normal distance 
$y^+$ (vertical axis). Apparently, a large portion of particles experience a strong shear where 
$Sr\geq 1$, in particular for the two cases where 
$D^+\geq 1$. Finally, it is worth noting that the probability of a particle experiencing very high values of 
$Sr$ near the wall will increase with decreasing values of 
$St^+$ as, for the same local shear rate, particles will experience an ever smaller slip velocity. This exacerbates the violation of the 
$Sr=O(1)$ constraint, meaning that current lift force models are even less suitable for turbulent channel flow laden with particles at lower Stokes numbers.
Figure 10. The p.d.f.s of 
$u_s/(\omega _z D)$ along the wall-normal direction (1WC DNS with the CBP model): (a) CL1, (b) CM1, (c) CS1. The vertical dashed lines denote 
$u_s/(\omega _z D) = \pm 1$.
Additionally, the discussion above focuses on the lift force arising from the shear in the ambient flow. In wall-bounded flows, as considered in this work, the lift force may deviate significantly from its unbounded counterpart if the particle is close to the wall. In brief, the presence of the wall leads to a repulsive transverse force in the absence of the ambient shear (Vasseur & Cox Reference Vasseur and Cox1977; Zeng et al. Reference Zeng, Balachandar and Fischer2005), and in the presence of an ambient shear, a suppression of the shear-induced lift (McLaughlin Reference McLaughlin1993; Magnaudet et al. Reference Magnaudet, Takagi and Legendre2003; Takemura et al. Reference Takemura, Magnaudet and Dimitrakopoulos2009). Although analytical solutions for the combined lift force are achievable in some asymptotic limits (see e.g. Wang et al. (Reference Wang, Squires, Chen and McLaughlin1997) and, more recently, Shi & Rzehak (Reference Shi and Rzehak2020) for an overview of these solutions), there is so far no satisfactory way to glue these solutions together to achieve a general treatment of the wall effect (Ekanayake et al. Reference Ekanayake, Berry, Stickland, Dunstan, Muir, Dower and Harvie2020; Ekanayake, Berry & Harvie Reference Ekanayake, Berry and Harvie2021). Here, instead of trying to make any further progress in achieving such a general solution, it is more feasible to discuss under which condition this intricate wall effect may be neglected. Within the 
$O(D)$-depth layer close to the wall, wall–particle collisions govern the wall-normal position of the peak particle concentration. Hence the hydrodynamic wall–particle interaction may play a significant role only beyond 
$y^+\approx D^+$. On the other hand, it is known that (McLaughlin Reference McLaughlin1993) the amount of reduction of the shear-induced lift due to the wall decreases as 
$(\ell /\ell _\omega )^{-5/3}$, and reduces to only 
$20\,\%$ of its magnitude in the limit 
$\ell /\ell _\omega \to 0$ at 
$\ell /\ell _\omega \approx 3$. Hence the wall effect may be considered negligible if the position of the peak particle concentration, which is of 
$O(D)$ from the wall, is larger than approximately 
$3\ell _\omega$, which yields approximately 
$D^+\geq 3$. This finding, together with the previous one where (shear-induced) lift force is non-negligible for 
$D^+\geq 1$, indicates that such a wall effect can be disregarded except in the overlap regime where 
$1\leq D^+ \leq 3$.
In addition to the ambient shear, particle rotation is known to produce a lift force as well (Rubinow & Keller Reference Rubinow and Keller1961). This lift contribution is often neglected in prior PP-DNS and not considered in the present work. Nevertheless, it might be worth discussing its relevant importance in the following. Assuming that the particle inertia is small enough for the particle to always stay in the torque-free state, particle rotation leads to a lift contribution 
$C_L^\varOmega =0.5\,Sr$ in the double limits considered by Saffman (Reference Saffman1965). The ratio between this spin-induced lift contribution and that induced by ambient shear is proportional to 
${\textit {Re}}_\omega ^{1/2}$, indicating that 
$C_L^\varOmega$ is merely a component of the second-order lift contribution from the ambient shear. This is confirmed by the recent work of Candelier et al. (Reference Candelier, Mehaddi, Mehlig and Magnaudet2023), where second-order lift contributions (with respect to 
${\textit {Re}}_\omega ^{1/2}$) were obtained using matched asymptotic expansions. In addition to the spin-induced contribution 
$C_L^\varOmega =0.5Sr$, they also obtained a second-order contribution from the ambient shear, 
$C_L^{'(2)}=-0.505\,Sr$, which counterbalances the spin-induced contribution. In other words, the effect of particle rotation on the lift is negligible in the torque-free state. This conclusion is valid irrespective of the presence of the wall (Cherukat & McLaughlin Reference Cherukat and McLaughlin1994; Magnaudet et al. Reference Magnaudet, Takagi and Legendre2003; Shi & Rzehak Reference Shi and Rzehak2020). However, large-inertia particles with a large density ratio as considered in this work might not always stay torque-free. Specifically, for particles with density ratio 
$\rho _p/\rho \geq 100$, the relaxation time in response to torque is comparable with that to slip (see e.g. Bagchi & Balachandar Reference Bagchi and Balachandar2002). This finding, together with the pronounced slip in the viscous sublayer revealed by figures 5(a,b), indicates that particles with 
$D^+\geq 1$ may experience a strong angular acceleration and possibly a non-negligible lift contribution from forced rotation. Investigating the effects of rotation by accounting for the particles’ conservation of angular momentum would therefore be interesting but is outside the scope of the present work due to the higher complexity associated with drag and lift closures for angular momentum.
Finally, it is worth noting that while the CBP lift force model used in Costa et al. (Reference Costa, Brandt and Picano2020a,Reference Costa, Brandt and Picanob) appears to provide reasonable predictions in some cases, it is not valid universally. Specifically, it does not align with any of the inner-region solutions (Shi & Rzehak Reference Shi and Rzehak2020), which are applicable for particles extremely close to the wall (e.g. within the viscous sublayer), nor does it approach the unbounded solution (Saffman Reference Saffman1965; McLaughlin Reference McLaughlin1991) for particles located within the outer layer of the channel flow. The agreement between the PP-DNS utilizing this lift force model and PR-DNS suggests that the model offers a satisfactory blending of results in the two limits for the specific cases explored in Costa et al. (Reference Costa, Brandt and Picano2020a) (see table 1). A preliminary attempt to establish a more general lift force model can be found in the work of Ekanayake et al. (Reference Ekanayake, Berry, Stickland, Dunstan, Muir, Dower and Harvie2020, Reference Ekanayake, Berry and Harvie2021), where a lift model satisfying the above conditions was proposed in the limit of small shear Reynolds numbers, i.e. 
${\textit {Re}}_\omega \ll 1$. However, the applicability of this model is limited to 
$D^+\ll 1$ because in the viscous layer, 
${\textit {Re}}_\omega \approx (u_\tau ^2/\nu ) D^2/\nu =(D^+)^2$. While this provides some insights, its practicality may be limited since the lift force no longer plays a pivotal role for 
$D^+<1$, as demonstrated here. Hence further work extending Ekanayake et al. (Reference Ekanayake, Berry, Stickland, Dunstan, Muir, Dower and Harvie2020, Reference Ekanayake, Berry and Harvie2021) to finite 
${\textit {Re}}_\omega$ (hence 
$D^+$) would be significant.
3.4. Lift-induced turbulence modulation
Finally in this section, we use 2WC point-particle DNS to illustrate that lift force models can have tremendous consequences in basic integral quantities such as the overall drag. To achieve this, we use the same parameters as in table 1, which fixes the total mass fraction and particle Stokes number while varying the particle diameter. We should reiterate that high-fidelity 2WC point-particle DNS are highly challenging and are currently being investigated actively. Still, the simple particle-in-cell 2WC DNS employed here suffice to illustrate the impact of lift force model choice on important turbulence metrics. In what follows, we illustrate the effects of lift force on the overall drag, and analyse the streamwise momentum fluxes that contribute to it.
Figure 11 shows the friction coefficient for the 2WC DNS, defined as 
$C_f=2 u^2_{\tau }/U^2_b$. While all small-inertia particles have a drag-increasing effect, the difference between lift force models is highly amplified as the particle size increases. At 
$D^+=3$, in particular, particles under the Saffman lift force show a drag reduction compared to the other particle sizes, almost reaching an overall drag decrease. Instead, the other lift models show a monotonic drag increase.
Figure 11. The wall friction coefficient 
$C_f$ versus the particle diameter 
$D^+$ (2WC). The horizontal dashed line denotes the mean 
$C_f$ of the particle-free case.
To better understand and quantify this difference, we follow the analysis in Costa et al. (Reference Costa, Brandt and Picano2021) and investigate the stress budget in the 2WC limit of vanishing volume fraction, while keeping the mass fraction finite:
\begin{equation} \tau = \rho u_\tau^2\left(1-\frac{y}{h}\right) \approx \underbrace{\left(\mu\,\frac{\mathrm{d}\langle u\rangle}{\mathrm{d}y} - \rho\langle u^\prime v^\prime\rangle\right)}_{\tau_\nu + \tau_{T,f}} - \underbrace{\langle \psi\rangle \langle u_p^\prime v_p^\prime\rangle}_{\tau_{T,p}}, \end{equation}
where 
$\langle \psi \rangle = \rho _p\langle \phi \rangle$ is the mean local mass fraction. The terms on the right-hand-side of this equation are the viscous stress, fluid Reynolds stress and particle Reynolds stress (hereafter 
$\tau _\nu$, 
$\tau _{T,f}$ and 
$\tau _{T,p}$). This equation boils down to the stress budget in the single-phase limit (
$\tau _\nu$ and 
$\tau _{T,f}$), corrected by an additional turbulent particle momentum flux 
$\tau _{T,p}$ that involves the local mass fraction and the particle counterpart of the Reynolds shear stress. This term relates to the transfer of streamwise momentum by correlated particle velocity fluctuations, and shows the importance of high particle mass fraction near the wall when combined with correlated streamwise–wall-normal particle velocity fluctuations. We have seen in § 3.1 that while the former decreases in the presence of lift force, the latter may be increased by it.
The profiles of stress terms and the corresponding contributions of each momentum transfer mechanism to the overall friction coefficient, normalized by the overall friction coefficient of the particle-free case, are shown in figure 12. In the left-hand plots, the normalized shear stress profiles are shown, while on the right-hand side we show the so-called FIK (Fukagata–Iwamoto–Kasagi) identity employed to a particle suspension (Fukagata, Iwamoto & Kasagi Reference Fukagata, Iwamoto and Kasagi2002; Yu et al. Reference Yu, Xia, Guo and Lin2021). In the case of a turbulent plane channel flow, this identity decomposes the mean friction coefficient (
$C_f$) into laminar and turbulent contributions. Here, the latter can be further decomposed into a contribution due to the fluid phase and one due to the particles. Each contribution is given by the weighted integral 
$C_{f,i} = (1/h^2)\int _0^h 6 (h-y) \tau _{i}\,\mathrm {d} y/(\rho U_b^2)$, with 
$\tau _i$ being one of the stress mechanisms described above. Note that the residual of this stress budget is virtually zero, which validates the approximation in (3.1).
Figure 12. (a,c,e) Budget of streamwise momentum and (b,d,f) contribution of each stress contribution 
$\tau _i$ to the mean wall friction coefficient 
$C_f=2\tau _w/(\rho U_b^2)$, i.e. 
$(1/h^2)\int _0^h 6 (h-y) \tau _i \,\mathrm {d}y/(\rho U_b^2)$ normalized by that of the particle-free case: (a,b) CL1, (c,d) CM1, (e,f) CS1. Different line colours in (a,c,e) denote different cases, which match the legend in (c).
As anticipated from figure 11, all 2WC DNS show increased drag compared to the unladen flow. Yet, remarkably, different lift models show qualitatively different turbulence modulation for the largest particle size (
$D^+=3$). Particles under the Saffman lift force have a major net turbulence attenuation effect, which results in an almost drag-reducing flow with respect to the unladen case. Conversely, particles under the Mei and CBP lift forces have a net drag increase, caused mostly by an increase in the 2WC stress term 
$\psi \langle u'_pv_p'\rangle$. This results from the amplification of particle Reynolds shear stresses 
$\langle u'_pv'_p\rangle$ discussed in § 3.1, which compensate for the decrease in local mass fraction 
$\psi$. Naturally, as the particle size decreases, the differences between lift models become less pronounced but still noticeable at all particle sizes, with particles under the Saffman lift force consistently showing the highest fluid turbulence attenuation. For the smallest particle size 
$D^+=0.1$, the two turbulent stress (fluid and particle) terms still show noticeable but in practice negligible differences (recall figure 11).
Hence it appears that the most important dynamical effects of lift forces in 2WC conditions can be seen as finite-size effects that directly modify single-point moments of particle velocity for 
$D^+\gtrsim 1$ (recall figure 6). However, the lift force still affects near-wall turbulence modulation when particles are relatively small, as they strongly affect the mean particle concentration (recall figure 2), directly linked to turbulence modulation via the 2WC stress term.
4. Conclusions
We have investigated the importance of lift forces in the dynamics of particle-laden turbulent channel flow. Three relevant shear-induced lift force models are analysed using PP-DNS for varying particle size and fixed Stokes number. Our results confirm that shear-induced lift forces strongly influence particle dynamics. We have analysed the main mechanisms for near-wall particle accumulation under lift forces: lift-induced migration, turbophoresis and biased sampling. Very close to the wall, where the particle-to-fluid slip velocity is positive and particle fluctuations are vanishing, lift forces and turbophoresis cooperate in driving particles towards the wall. Particle–wall interactions and biased sampling balance this effect. Further away from the wall, lift forces and turbophoresis compete, with the former having a wall-repelling effect. Concerning biased sampling, we observe that for larger particle sizes (
$D^+ \geq 1$), the lift-free ascending particles prefer to sample ejections, while the descending ones tend to sample sweeps. Conversely, particles under lift forces tend to sample ejection regions while experiencing a marginal occurrence of sweep events, regardless of the employed lift model or wall-normal particle velocity sign. This is again a consequence of the competing wall-repelling effect of lift force. For the smallest particles (
$D^+=0.1$), all cases exhibit over-sampling in ejection and inward motions of low-speed fluid regions (i.e. Q3 events). Overall, the sampling percentage of the ejections events only slightly surpasses that of Q3, suggesting that biased sampling plays a minor direct role in near-wall particle transport. However, the biased sampling of low-speed regions exposes particles for long times to persistent lift forces, which results in significantly reduced particle accumulation, even for small particle sizes.
These different dynamics have major consequences for turbulence modulation due to a modification of the near-wall values of (1) the local particle concentration, and (2) the correlated streamwise–wall-normal particle velocity fluctuations. While the latter is modified only for sufficiently large particle size and therefore can be seen as a finite-size effect, the former is still non-negligible for relatively small particle sizes, and can still result in small but visible changes in overall drag. Indeed, correlated particle velocity fluctuations show little sensitivity to lift forces for sufficiently small values of inner-scaled particle diameter (
$D^+\lesssim 0.1$). Conversely, for small particle sizes, the modification of near-wall particle accumulation due to lift showed little sensitivity to the flavour of shear-induced lift model. Hence to account accurately for lift-induced turbulence modulation with 
$D^+=O(1)$ particle sizes, reliable shear-induced lift models are needed that correctly reproduce the particle turbulent momentum flux. This could be seen as a necessary correction for finite-size effects, in the same spirit as the Faxén correction in the Maxey–Riley–Gatignol equation.
While such correction is important for reliable particle-modelled simulations of turbulent wall flows, our analysis has shown that current lift force models are bound to lose their validity in the near-wall region for small inertial particles. The reason is that the well-known tendency of particles to over-sample low-speed regions in the buffer layer, along with their tendency to flow faster than the fluid in the viscous sublayer, results in a region near the wall of vanishing particle-to-fluid slip velocity. Naturally, this issue will be even more pronounced at lower Stokes numbers, where the range of very low slip velocities in regions of high shear is wider. Unfortunately, current shear-induced lift force models are available only for non-dimensional shear rates 
$Sr$ of order unity, which robustly tends to infinity within the buffer layer. Hence the condition 
$Sr\leq 1$ required by these models is seriously violated in wall turbulence with small inertial particles.
Acknowledgements
The Cray XC40 Shaheen II at KAUST was used for all simulations reported. P.C. thanks F. Picano from Padova University for insightful discussions.
Funding
This research was partially supported by the KAUST Office of Sponsored Research (OSR) under award no. OSR-2019-CCF-3666 and under baseline research funds of M.P.
Declaration of interests
The authors report no conflict of interest.
Open access
