# Asymptotic closure model for inertial particle transport in turbulent
boundary layers

Y. Zhang  Department of Civil and Environmental Engineering, Duke University,
Durham, North Carolina 27708, USA  G. Wang  Life Science Technologies
Department, Interuniversity Microelectronics Centre(IMEC), Leuven 3000,
Belgium  A.D. Bragg<EMAIL_ADDRESS>Department of Civil and
Environmental Engineering, Duke University, Durham, North Carolina 27708, USA

###### Abstract

Transport equations for heavy inertial particles in turbulent boundary layers
may be derived from an underlying phase-space probability density function
(PDF) equation. These equations, however, are unclosed, and the standard
closure approach is to use a quasi-Normal approximation (QNA) in which the
fourth moments are approximated as behaving as if the velocities were Normally
distributed. Except for particles with weak inertia, the QNA leads to large
quantitative errors, and is not consistent with the known asymptotic
predictions of Sikovsky (Flow Turbulence Combust, vol. 92, 2014, pp. 41-64)
for the moments of the PDF in the viscous sublayer. We derive a new closure
approximation based on an asymptotic solution to the transport equations in
regions where the effect of particle inertia is significant. The new closure
is consistent with the asymptotic predictions of Sikovsky, but applies even
outside the viscous sublayer. Comparisons with direct numerical simulations
(DNS) show that the new closure gives similar results to the QNA (with the QNA
results in slightly better agreement with the DNS) when the viscous Stokes
number is $St<10$, but for $St>10$ the new model is in far better agreement
with the DNS than the QNA. While the predictions from the new model leave room
for improvement, the results suggest that this new closure strategy is a very
effective alternative to the traditional QNA approach, and the closure could
be refined in future work.

##  I Introduction

The classical model for predicting the concentrations of settling inertial
particles in wall-bounded turbulence is that of Rouse [1]. This model is
purely phenomenological, and assumes that the effect of particle inertia is
negligible except for the finite gravitational settling velocity it
introduces. While extensions of this model to particles with small but finite
Stokes number have been considered using regular perturbation theory [2],
developing a model for moderate to large Stokes numbers remains a significant
challenge. In order to understand in detail the additional physical mechanisms
introduced by finite particle inertia, and therefore the additional terms that
an extended Rouse model must capture, in [3], the settling and concentration
profiles of settling inertial particles in wall-bounded turbulence was
explored using theory and direct numerical simulations (DNS). The particles
were assumed to be small and heavy, with their motion governed by a Stokes
drag force and gravity. In contrast to the phenomenological approach of [1], a
rigorous, systematic approach for deriving transport equations for inertial
particles in turbulent flows is to derive them as moments of an underlying
phase-space PDF equation [4, 5, 6, 7, 8]. Therefore, in [3] the exact (but
unclosed) transport equations governing the particle concentration and
settling velocities were constructed as moments of a PDF equation for the
particle position and velocity. An analysis of the equations led to new
insights into the physical mechanisms governing these processes, and how they
differ from those in homogeneous turbulence that was explored in [9, 10]. Data
from DNS was then used to evaluate the various terms appearing in these
equations, providing insights into the role played by the inertial mechanisms
that are absent from the classical model of Rouse [1]. These additional terms
were shown in [3] to be of leading order importance in determining the
particle settling and concentrations, unless the particle Stokes number is
very small, and these terms must therefore be accounted for if a Rouse-type
model is to be extended to the case of finite Stokes number particles.

While the analysis of [3] led to new physical insights, in order to develop a
predictive theoretical model the hierarchy of moment equations derived from
the PDF equation most be closed. The usual closure assumption is to make a
quasi-Normal approximation (QNA) [11] (which can also be derived using a
Chapman-Enskog approach [12]), according to which the hierarchy of moment
equations is closed by assuming that the fourth moment behaves as if the
particle velocities were Normally distributed. Aside from being ad-hoc (in
general), this closure approach is known to lead to significant errors in the
predictions from the closed moment equations, errors that are both
quantitative and qualitative in nature (similar issues also arise when the QNA
is used in models of inertial particle-pair transport in isotropic turbulence
[13, 14]). These issues will be discussed in detail in §II.2.

In this paper we explore an alternative closure approximation that captures
the asymptotic behavior of the relation between the fourth and second moments
of the particle velocity in regimes of the flow where locally the effect of
particle inertia is significant. This closure captures the strong non-
Gaussianity of the particle velocities in the near-wall region of the flow,
and is consistent with the asymptotic behavior of the moments in the viscous
sublayer that was described in [15, 16]. The predictions of the new closed
model are compared with DNS data and we find that while the new model is not
always in full quantitative agreement with the DNS, it does provide far
superior predictions compared to the QNA when the particle inertia is moderate
to strong. Moreover, even when there are quantitative discrepancies, the
predictions are qualitatively consistent with the DNS data, unlike the QNA
model whose solutions also feature a spurious bifurcation near the wall as the
Stokes number is increased beyond a threshold value[17]. Therefore, while
there is still room for improvement, the results suggest that the new closure
approach is promising, and could be further refined in future work.

##  II Transport equations for inertial particles in a turbulent boundary
layer

###  II.1 Hierarchy of moment equations

In this work, we consider the transport of small, heavy inertial particles
subject to the equation of vertical motion (in what follows, all variables are
in wall units, with the usual “$+$” superscript omitted for notational
simplicity)

$\displaystyle\ddot{z}^{p}(t)\equiv\dot{w}^{p}(t)=\frac{1}{St}\Big{(}u^{p}(t)-w^{p}(t)\Big{)},$
(1)

where $St\equiv\tau_{p}/\tau_{*}$ is the particle Stokes number, $\tau_{p}$ is
the particle response time, $\tau_{*}\equiv\nu u_{*}^{-2}$ is the fluid time
scale based on the friction velocity $u_{*}$, $z^{p}(t),w^{p}(t)$ are the
vertical particle position and velocity, and $u^{p}(t)$ is the vertical fluid
velocity at the particle position. The particle volume and mass loadings are
assumed sufficiently small to ignore particle-particle collisions and two-way
coupling. While this is a highly simplified system, the great difficulties in
developing statistical transport equations even for this simple system mean
that incorporating additional complexities are best left for future work.
Moreover, while as discussed in the introduction, our ultimate interest is in
a model for settling inertial particles, we are here focusing on developing an
improved closure method for the simpler case of non-settling particles. An
extension of the new model proposed in this paper to the case of settling
particles will be the subject of a future study.

The joint probability density function (PDF) for $z^{p}(t),w^{p}(t)$ in a 2D
phase-space with time-independent coordinates $z,w$ is defined as

$\displaystyle\mathcal{P}(z,w,t)\equiv\Big{\langle}\delta(z^{p}(t)-z)\delta(w^{p}(t)-w)\Big{\rangle},$
(2)

where $\delta(\cdot)$ is the Dirac delta distribution, and the evolution
equation is [16, 3]

$\displaystyle\partial_{t}\mathcal{P}=-\nabla_{z}\Big{(}\mathcal{P}w\Big{)}-\nabla_{w}\Big{(}\mathcal{P}\langle\dot{w}^{p}(t)\rangle_{z,w}\Big{)},$
(3)

where the operator $\langle\cdot\rangle_{z,w}$ denotes an ensemble average
conditioned on $z^{p}(t)=z$, $w^{p}(t)=w$, and $\nabla_{z}$ and $\nabla_{w}$
denote differentiation with respect to $z$ and $w$, respectively.

The $N^{th}$ order moment $\mathcal{M}_{N}(z,t)$ of the PDF is defined as

$\displaystyle\mathcal{M}_{N}(z,t)\equiv\int_{\mathbb{R}}w^{N}\mathcal{P}\,dw=\varrho\mathcal{W}_{N},$
(4)

where $\varrho\equiv\mathcal{M}_{0}$ is the concentration PDF and
$\mathcal{W}_{N}\equiv\langle[{w}^{p}(t)]^{N}\rangle_{z}$ are the moments of
the particle vertical velocity.

The evolution equation for $\mathcal{M}_{N}(z,t)$ can be obtained from (3) and
is given by [16]

$\displaystyle\partial_{t}\mathcal{M}_{N}=-\nabla_{z}\mathcal{M}_{N+1}+N\varrho\Big{\langle}\dot{w}^{p}(t)[w^{p}(t)]^{N-1}\Big{\rangle}_{z}.$
(5)

Both terms on the rhs of this equation are unclosed. Closed expressions for
the second term can be derived once a closure for
$\langle\dot{w}^{p}(t)\rangle_{z,w}$ is specified. Given (1) we have

$\displaystyle\langle\dot{w}^{p}(t)\rangle_{z,w}=\frac{1}{St}\Big{(}\langle
u^{p}(t)\rangle_{z,w}-w\Big{)},$ (6)

and various approximations have been introduced for closing the conditional
average $\langle u^{p}(t)\rangle_{z,w}$ (see [8] for a detailed discussion).
The closure approximation that can be proven to be formally consistent with
the fully-mixed condition in the limit $St\to 0$ is that derived using the
Furutsu-Novikov formula assuming that the fluid velocity field has Gaussian
statistics [6, 8], leading to the closure

$\displaystyle\frac{1}{St}\mathcal{P}\langle
u^{p}(t)\rangle_{z,w}\approx\mathcal{P}\kappa-\nabla_{z}\Big{(}\mathcal{P}\lambda\Big{)}-\nabla_{w}\Big{(}\mathcal{P}\mu\Big{)},$
(7)

where $\kappa(z,w,t)$ is a drift coefficient, and $\lambda(z,w,t)$,
$\mu(z,w,t)$ are dispersion coefficients for diffusion in $z,w$-space,
respectively. The details of these coefficients are not important for the
present discussion, and so will be given later, except to note that under
standard approximations, their dependence on $w$ is neglected, and
$\kappa=\nabla_{z}\lambda$ is assumed [18] (although this is strictly only
valid in the steady-state when $St\to 0$ [8]). These approximations will be
assumed throughout this paper, under which (7) simplifies to

$\displaystyle\frac{1}{St}\mathcal{P}\langle
u^{p}(t)\rangle_{z,w}\approx-\lambda\nabla_{z}\mathcal{P}-\mu\nabla_{w}\mathcal{P}.$
(8)

While the focus of this study is on moment equations derived from a kinetic
PDF equation, higher dimensional PDF equations that also include $u^{p}(t)$ in
the phase-space have been considered extensively. For these, $u^{p}(t)$ is
usually modeled via a generalized Langevin model (GLM) [19], which is an
approach that was first developed in the context of single-phase turbulent
flows [20]. The kinetic and GLM PDF approaches have their own merits, as
discussed in [7, 21]. Since an ultimate goal of our work is to develop a model
that extends that of [1] to the case of finite inertia particles, the kinetic
approach is preferred here, since the use of a GLM based PDF equation
ultimately requires one to construct the solutions to the moment equations via
a Monte-Carlo method, rather than simply as the solution to a set of coupled
PDEs (although [22] explores direct numerical solutions of the PDEs defined
via the GLM PDF model).

###  II.2 Quasi-Normal approximation

The standard approach for closing the first term on the rhs of (5) is to use a
quasi-Normal approximation (QNA) [11] (also derived using a Chapman-Enskog
approach [12]). This approach may be summarized as follows. By specifying the
particle acceleration $\dot{w}^{p}(t)$ that appears in the rhs of (5) using
(1), then the steady-state form of the $N=3$ equation can be re-arranged to
give

$\displaystyle\mathcal{M}_{3}=-(St/3)\nabla_{z}\mathcal{M}_{4}+\varrho\Big{\langle}{u}^{p}(t)[w^{p}(t)]^{2}\Big{\rangle}_{z}.$
(9)

Assuming that $\mathcal{W}_{4}$ behaves as if $w^{p}(t)$ were Normally
distributed leads to $\mathcal{W}_{4}\approx 3\mathcal{W}_{2}^{2}$ and hence
$\mathcal{M}_{4}\approx 3\mathcal{M}_{2}^{2}/\varrho$. Inserting this into
(9), and using (8) then leads to

$\displaystyle\mathcal{M}_{3}\approx-
St\nabla_{z}(\mathcal{M}_{2}^{2}/\varrho)-St\lambda\nabla_{z}\mathcal{M}_{2}.$
(10)

This equation can then be substituted into the transport equation for
$\mathcal{M}_{2}$, and after some manipulation, this leads to a second-order
ODE for $\mathcal{W}_{2}$ [23]

$\displaystyle 0$
$\displaystyle=(\mathcal{W}_{2}+\lambda)\nabla_{z}^{2}\mathcal{W}_{2}+\nabla_{z}\lambda\nabla_{z}\mathcal{W}_{2}-2\mathcal{W}_{2}/St^{2}+2\mu/St,$
(11)

whose solution can be be used to obtain $\varrho$ from the steady-state $N=1$
equation, namely

$\displaystyle 0$
$\displaystyle=-St(\lambda+\mathcal{W}_{2})\nabla_{z}\varrho-
St\varrho\nabla_{z}\mathcal{W}_{2}.$ (12)

Since the model assumes that the fluid velocity field is Gaussian when closing
$\langle u^{p}(t)\rangle_{z,w}$, then the QNA is self-consistent in the regime
$St\ll 1$. However, for $St\geq O(1)$, the statistics of $w^{p}(t)$ are
expected to be strongly non-Gaussian even if the fluid velocity field is
Gaussian [15], owing to the non-local nature of the inertial particle
dynamics. There are also known to be two particular errors introduced by the
QNA, which we now discuss.

First, the QNA leads to behavior for $\mathcal{W}_{N}$ that is inconsistent
with the asymptotic behavior for $z\to 0$ when $St\geq O(1)$ [15]. In
particular, for $St\to 0$, the scaling of the vertical fluid velocity field
for $z\to 0$ implies $\mathcal{W}_{N}\propto z^{2N}$, and the QNA result
$\mathcal{W}_{4}\approx 3\mathcal{W}_{2}^{2}$ is consistent with this.
However, for $St\geq O(1)$, $\mathcal{W}_{N}\propto z^{\gamma}$ [15, 16],
where $\gamma(St)$ is the power-law exponent describing $\varrho$ in the limit
$z\to 0$, namely $\varrho\sim z^{-\gamma}$. The QNA is not consistent with
this because it predicts $\mathcal{W}_{4}\propto z^{2\gamma}$ rather than the
correct behavior $\mathcal{W}_{4}\propto z^{\gamma}$.

The second issue is that the QNA equation for $\mathcal{W}_{2}$ given by (11)
predicts a bifurcation in the solution as $St$ exceeds a threshold value [17],
which through (12) also leads to a bifurcation in the solution for $\varrho$.
This predicted bifurcation is not supported by DNS data and is argued to be
unphysical [17], and will be illustrated in §III.

In view of these serious issues with the QNA for $St\geq O(1)$, an alternative
closure approximation for $\mathcal{M}_{4}$ is desirable that is both
consistent with the known asymptotic behavior of the particle velocities in
the limit $z\to 0$, and also avoids the unphysical bifurcations predicted by
the QNA model.

###  II.3 Asymptotic closure approximation

An alternative closure approximation is motivated by the observation in [16]
that the normalized solutions to the steady-state transport equations for
$\mathcal{W}_{N}$ can be written as

$\displaystyle\begin{split}\mathcal{W}_{N}\Big{/}\mathcal{W}_{2}^{N/2}=\mathcal{C}_{N}\varrho^{N/2-1}\exp&\Bigg{(}\frac{N-1}{St}\int^{z}\mathcal{W}_{N}^{-1}(q)\Big{\langle}\Big{(}{u}^{p}(t)-w^{p}(t)\Big{)}[w^{p}(t)]^{N-2}\Big{\rangle}_{q}\,dq\\\
&\quad-\frac{N}{2St}\int^{z}\mathcal{W}_{2}^{-1}\Big{\langle}u^{p}(t)\Big{\rangle}_{q}\,dq\Bigg{)},\end{split}$
(13)

where $\mathcal{C}_{N}$ are constants with respect to $z$, but will in general
depend on $St$. In view of this result, for $St\gg 1$ the quantity
$\mathcal{W}_{N}/\mathcal{W}_{2}^{N/2}$ behaves asymptotically as

$\displaystyle\mathcal{W}_{N}\Big{/}\mathcal{W}_{2}^{N/2}\sim\mathcal{C}_{N}\varrho^{N/2-1}\Big{(}1+O(1/St)\Big{)}.$
(14)

This asymptotic result is valid for arbitrary $z$, however, it is expected
that very large values of $St$ would be needed in practice in order to observe
this behavior across the entire boundary layer (in particular, it would
require that the Stokes number based on the largest timescale in the flow is
$\gg 1$). The results in [15] also imply that (14) is valid in the viscous
sublayer even for $St=O(1)$ since for $z\to 0$ the result in (14) reduces to
the asymptotic results for the regime $St\geq O(1)$ predicted by [15].

The result in (14) yields the asymptotic closure approximation (ACA)

$\displaystyle\mathcal{M}_{4}\sim\mathcal{C}_{4}\mathcal{M}_{2}^{2},$ (15)

whose most important difference compared to the QNA result
$\mathcal{M}_{4}\approx 3\mathcal{M}_{2}^{2}/\varrho$ is the absence of the
factor $1/\varrho$. It is precisely because the QNA contains the factor
$1/\varrho$ that it leads to a behavior for $\mathcal{W}_{4}$ that is
inconsistent with the asymptotic behavior predicted by [15] in the limit $z\to
0$. In the near-wall region where $\varrho$ can be very large and exhibits a
power-law dependence on $z$ [15, 16], the QNA and ACA for $\mathcal{M}_{4}$
will be radically different, both in terms of their qualitative and
quantitative behavior.

In order to use (15) to close the moment equations, the constant
$\mathcal{C}_{4}$ must be specified. While this will in general depend upon
$St$, the simplest choice is to use $\mathcal{C}_{4}=3/\varrho(z_{b})$, where
$z_{b}$ is the upper boundary of the solution domain and $\varrho(z_{b})$ is
the boundary condition imposed when solving (12). In a flow with friction
Reynolds number $Re_{\tau}\to\infty$ and $St$ large but finite, then provided
that $z_{b}$ is large enough to correspond to a height at which the effects of
the particle inertia are negligible,
$\mathcal{M}_{4}\sim\mathcal{C}_{4}\mathcal{M}_{2}^{2}$ approaches the QNA
result $\mathcal{M}_{4}\approx 3\mathcal{M}_{2}^{2}/\varrho$ as $z\to z_{b}$.
This is a self-consistent choice given that the closure for $\langle
u^{p}(t)\rangle_{z,w}$ assumes that the wall-normal fluid velocities are
Normally distributed, and therefore the PDF of $w^{p}(t)$ should be only
weakly perturbed from a Normal distribution in regions where the effect of
particle inertia is weak.

An important point is that although (15) will not be accurate when $St\ll 1$,
this does not in practice matter. The reason for this is two-fold. First,
since

$\displaystyle\mathcal{M}_{3}=-(St/3)\nabla_{z}\mathcal{M}_{4}+\varrho\Big{\langle}{u}^{p}(t)[w^{p}(t)]^{2}\Big{\rangle}_{z},$
(16)

then in the regime $St\ll 1$ the contribution from the term involving
$\mathcal{M}_{4}$ (whose closure based on (15) is not accurate for $St\ll 1$)
will be very small (noting that $\mathcal{C}_{4}\mathcal{M}_{2}^{2}$ and its
gradient are finite in the limit $St\to 0$). Therefore, errors in the closure
for $\mathcal{M}_{4}$ will only lead to small errors in the overall model
predictions for $\mathcal{M}_{2}$. Second, with the aforementioned choice
$\mathcal{C}_{4}=3/\varrho(z_{b})$, then (15) asymptotes to the QNA for small
$St$ (for which $\varrho$ is almost uniform), and this is know to yield
reasonable predictions for $St\ll 1$.

Using (15) to specify $\mathcal{M}_{4}$ in (16), and substituting the
resulting equation for $\mathcal{M}_{3}$ into the equation for
$\mathcal{M}_{2}$ leads to the second-order ODE for $\mathcal{M}_{2}$

$\displaystyle
0=\mathcal{A}\nabla_{z}^{2}\mathcal{M}_{2}+\nabla_{z}\mathcal{A}\nabla_{z}\mathcal{M}_{2}-2\mathcal{M}_{2}/St^{2}+2\mu\varrho/St,$
(17)

where

$\displaystyle\mathcal{A}\equiv(2\mathcal{C}_{4}/3)\mathcal{M}_{2}+\lambda.$
(18)

Since (17) explicitly contains $\varrho$, then (17) must be solved
simultaneously with the equation governing $\varrho$, namely (12). However, we
have found that the numerical stability of solutions to the coupled equations
for $\varrho$ and $\mathcal{M}_{2}$ is improved if instead a second-order ODE
is solved for $\varrho$. This may be obtained by substituting (12) (which
comes from the equation for $\mathcal{M}_{1}$) into the equation for
$\mathcal{M}_{0}$, yielding

$\displaystyle 0$
$\displaystyle=-\lambda\nabla_{z}^{2}\varrho-\nabla_{z}\lambda\nabla_{z}\varrho-\nabla_{z}^{2}\mathcal{M}_{2}.$
(19)

###  II.4 Boundary conditions & numerical solution

For the QNA model, two boundary conditions must be specified for
$\mathcal{W}_{2}$. A standard choice is to use
$\nabla_{z}\mathcal{W}_{2}|_{z_{a}}=0$ and either
$\nabla_{z}\mathcal{W}_{2}|_{z_{b}}=0$ [17] or
$\mathcal{W}_{2}|_{z_{b}}=St\mu(z_{b})$, where $z_{a},z_{b}$ are the lower and
upper boundary points. The Neumman condition
$\nabla_{z}\mathcal{W}_{2}|_{z_{b}}=0$ is suitable if $z_{b}$ lies in the
quasi-homogeneous region of the wall-bounded flow, or at the centerline of,
e.g. a channel flow. The Dirichlet condition
$\mathcal{W}_{2}|_{z_{b}}=St\mu(z_{b})$ is less restrictive since it is
appropriate provided that the local equilibrium solution to (11) is accurate,
without requiring anything about the gradients of $\mathcal{W}_{2}$. In the
QNA solutions shown later, this Dirichlet boundary condition will be used. The
point $z_{a}$ can be specified as $z_{a}=d_{p}/2$, where $d_{p}$ is the
particle diameter.

Since $\varrho$ is decoupled from $\mathcal{W}_{2}$ in the QNA model, the
solution for $\varrho$ can be obtained after obtaining $\mathcal{W}_{2}$ by
solving (12), for which a Dirichlet boundary condition $\varrho(z_{b})$ can be
used. Given that $\varrho$ is a PDF for $z^{p}(t)$, its integral over the full
flow should be equal to one. Due to the linearity of (12), $\varrho(z_{b})$
can be chosen arbitrarily, and the solution can be subsequently re-normalized
to satisfy this integral condition. However, if the model is only being solved
over a portion of the flow (e.g. the boundary layer), then the absolute values
of $\varrho$ cannot be determined, but only the concentration profile relative
to some reference value. In this case, the choice of $\varrho(z_{b})$ is
arbitrary and may be simply set to one.

In the new ACA model, the equations for $\varrho$ and $\mathcal{M}_{2}$ are
coupled, and the boundary conditions should be chosen to be consistent with
equation (12) which requires

$\displaystyle\nabla_{z}\mathcal{M}_{2}|_{z_{a}}$
$\displaystyle=-\lambda\nabla_{z}\varrho|_{z_{a}}.$ (20)

One choice would be to use
$\nabla_{z}\mathcal{M}_{2}|_{z_{a}}=-\lambda\nabla_{z}\varrho|_{z_{a}}=0$, and
this is the appropriate choice for $St\gg 1$ because
$\lim_{St\to\infty}\lambda=0$. For moderate values of $St$, an alternative is
to specify $\nabla_{z}\varrho|_{z_{a}}$ based on the local equilibrium
solution to $\varrho$. This is obtained by using the local equilibrium
solution $\mathcal{W}_{2}=St\mu$ in (12) yielding $\varrho^{eq}$. With this,
the Neumman boundary condition for $\mathcal{M}_{2}$ is obtained

$\displaystyle\nabla_{z}\mathcal{M}_{2}|_{z_{a}}$
$\displaystyle=-\lambda\nabla_{z}\varrho^{eq}|_{z_{a}}.$ (21)

This is similar to the approach described in [17] to specify
$\nabla_{z}\mathcal{W}_{2}|_{z_{a}}$ as an alternative boundary condition for
the QNA model. However, we found that (21) can lead to numerical instability
of the solution of the ACA model, and therefore we will use
$\nabla_{z}\mathcal{M}_{2}|_{z_{a}}=0$ for all $St$ values considered. Note
that this is consistent with the use of $\nabla_{z}\mathcal{W}_{2}|_{z_{a}}=0$
when solving the QNA model.

For the upper boundary, the local equilibrium solution
$\mathcal{M}_{2}|_{z_{b}}=St\mu(z_{b})\varrho(z_{b})$ may be used. Given the
linearity of the equation for $\varrho$, we may use $\varrho(z_{b})=1$, and
the solution can be subsequently normalized to yield
$\int_{z_{a}}^{z_{b}}\varrho dz=1$ in the case where the interval
$[z_{a},z_{b}]$ spans the height of the whole flow. Since we are using
$\nabla_{z}\mathcal{M}_{2}|_{z_{a}}=0$, then consistent with (20) we use
$\nabla_{z}\varrho|_{z_{a}}=0$ to specify the second boundary condition for
(19).

The QNA and ACA models involve second order, nonlinear ODEs. To solve them,
linearization with Newton-Raphson iteration was used. The local equilibrium
solution $\mathcal{W}_{2}(z)=St\mu|_{z_{b}}$ is used as the initial guess
(using $\mathcal{W}_{2}(z)=St\mu(z)$ leads to numerical issues for larger $St$
values), and the solutions converged rapidly, usually within 3 or 4
iterations.

##  III Comparison between models and DNS

In this section we compare the predictions from the QNA and ACA models for
$\varrho$ and $\mathcal{W}_{2}$ with DNS data of particle transport in an open
channel flow. The DNS data is from the same dataset as that in [3], except
that here there is no gravitational settling, and elastic particle-wall
collisions are used which leads to a steady state with
$\mathcal{M}_{1}(z)=0\forall z$. For the transport equations, the dispersion
coefficients $\lambda$ and $\mu$ must be specified, and for these we use the
standard local approximations [24, 18]

$\displaystyle\lambda(z)$ $\displaystyle\approx\frac{\tau_{L}\langle
uu\rangle}{St(1+St/\tau_{L})},$ (22) $\displaystyle\mu(z)$
$\displaystyle\approx\frac{\lambda}{\tau_{L}},$ (23)

where $u$ is the vertical fluid velocity at a fixed position (in contrast to
$u^{p}(t)$ which is the vertical fluid velocity along a particle trajectory).
In the results that follow, the DNS data for the fluid wall-normal Reynolds
stress $\langle uu\rangle$ is used, while the model discussed in [17] for the
fluid Lagrangian timescale seen by the particle $\tau_{L}$ was used.

The model equations were solved on a domain $z\in[d_{p},z_{b}]$ with
$z_{b}=200$. Regarding this choice of $z_{b}$, in the DNS the open channel
surface is located at $z=312.5$, and the solutions to the model are
insensitive to the choice of $z_{b}$ for the $St$ values considered if it is
chosen in the range $z_{b}\in(150,250)$. For $z_{b}$ significantly outside of
this range, the model predictions are compromised because the Dirichlet
boundary conditions $\mathcal{W}_{2}|_{z_{b}}=St\mu(z_{b})$ and
$\mathcal{M}_{2}|_{z_{b}}=St\mu(z_{b})\varrho(z_{b})$ are no longer
appropriate, given that they are based on a local equilibrium solution to the
equations for $\mathcal{W}_{2}$ and $\mathcal{M}_{2}$.

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{rhozmin_plot.pdf} \put(120.0,-3.0){$St$}
\put(0.0,80.0){\rotatebox{90.0}{$\varrho(z)$}} \put(51.0,151.0){ACA,
$z=z_{a}$} \put(51.0,140.0){ACA, $z=1$} \put(51.0,129.0){QNA, $z=z_{a}$}
\put(51.0,118.0){QNA, $z=1$} \end{overpic} (a)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt 199.16928pt,scale={0.4},clip,tics=20]{Wzmin_plot.pdf}
\put(120.0,-3.0){$St$} \put(-2.0,77.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(51.0,151.0){ACA, $z=z_{a}$} \put(51.0,140.0){ACA, $z=1$}
\put(51.0,129.0){QNA, $z=z_{a}$} \put(51.0,118.0){QNA, $z=1$} \end{overpic}
(b)

Figure 1: Comparison of the predictions from the ACA and QNA models for (a)
$\varrho(z)$ and (b) $\mathcal{W}_{2}(z)$ as a function of $St$ and for
$z=z_{a}$ and $z=1$.

We begin by comparing the QNA and ACA model predictions with each other, in
order to highlight their key differences. In §II.2 it was discussed that the
QNA model leads to a bifurcation in the solution as $St$ exceeds a threshold
value, as also discussed in [17]. In figure 1 we compare the model predictions
for $\varrho$ and $\mathcal{W}_{2}$ as a function of $St$ and for $z=z_{a}$
and $z=1$. The results in figure 1(a) show that the QNA model predicts that
$\varrho(z_{a})$ gradually increases with increasing $St$ until $St\approx
35$, at which point $\varrho(z_{a})$ suddenly reduces. The results for $z=1$
show similar behavior, except the drop in $\varrho$ is more gradual. Figure
1(b) shows the associated behavior of $\mathcal{W}_{2}(z_{a})$, for which the
QNA model predicts that $\mathcal{W}_{2}(z_{a})$ slowly decreases with
increasing $St$ until $St\approx 35$, and then $\mathcal{W}_{2}(z_{a})$
rapidly increases, before slowly increasing with increasing $St$. This is the
bifurcation behavior discussed in [17], which appears to be spurious, and is
not predicted by the asymptotic analysis of [15]. Once this bifurcation
occurs, the solutions for $\mathcal{W}_{2}$ in the viscous sublayer
dramatically switch from exhibiting a power-law dependence on $z$ to becoming
independent of $z$, which in turn causes a corresponding switch in the
behavior of $\varrho$, due to their coupling according to (12). This can be
observed in figure 1 by noting that for $St>35$, the solutions for $\varrho$
and $\mathcal{W}_{2}$ from the QNA model are the same for $z=z_{a}$ and $z=1$.

By contrast, the predictions from the ACA model do not show such abrupt
changes in the behavior of either $\varrho$ or $\mathcal{W}_{2}$, nor does
this model predict that these become independent of $z$. The results do
indicate, however, that the ACA model predicts that $\varrho(z)$ as a function
of $St$ peaks too early, noting that DNS data suggests the near-wall
concentration is strongest somewhere around $St\approx 30$ [25, 16]. This is
not surprising, however, given that the ACA is effectively derived for $St\gg
1$ as the leading order approximation of an asymptotic series.

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{Rho_plot_2p8.pdf} \put(130.0,0.0){$z$}
\put(0.0,70.0){\rotatebox{90.0}{$\varrho(z)/\varrho(z_{b})$}}
\put(163.0,150.0){DNS} \put(163.0,138.0){ACA} \put(163.0,128.0){QNA}
\put(150.0,70.0){$St=2.8$} \end{overpic} (a)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt 199.16928pt,scale={0.4},clip,tics=20]{Rho_plot_4p6.pdf}
\put(130.0,0.0){$z$}
\put(0.0,70.0){\rotatebox{90.0}{$\varrho(z)/\varrho(z_{b})$}}
\put(163.0,150.0){DNS} \put(163.0,138.0){ACA} \put(163.0,128.0){QNA}
\put(150.0,70.0){$St=4.6$} \end{overpic} (b)

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{Rho_plot_9p3.pdf} \put(130.0,0.0){$z$}
\put(0.0,70.0){\rotatebox{90.0}{$\varrho(z)/\varrho(z_{b})$}}
\put(163.0,150.0){DNS} \put(163.0,138.0){ACA} \put(163.0,128.0){QNA}
\put(150.0,70.0){$St=9.3$} \end{overpic} (c)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt 199.16928pt,scale={0.4},clip,tics=20]{Rho_plot_46p5.pdf}
\put(130.0,0.0){$z$}
\put(0.0,70.0){\rotatebox{90.0}{$\varrho(z)/\varrho(z_{b})$}}
\put(163.0,150.0){DNS} \put(163.0,138.0){ACA} \put(163.0,128.0){QNA}
\put(150.0,70.0){$St=46.5$} \end{overpic} (d)

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{Rho_plot_128.pdf} \put(130.0,0.0){$z$}
\put(0.0,70.0){\rotatebox{90.0}{$\varrho(z)/\varrho(z_{b})$}}
\put(163.0,150.0){DNS} \put(163.0,138.0){ACA} \put(163.0,128.0){QNA}
\put(150.0,70.0){$St=128$} \end{overpic} (e)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt 199.16928pt,scale={0.4},clip,tics=20]{Rho_plot_512.pdf}
\put(130.0,0.0){$z$}
\put(0.0,70.0){\rotatebox{90.0}{$\varrho(z)/\varrho(z_{b})$}}
\put(163.0,150.0){DNS} \put(163.0,138.0){ACA} \put(163.0,128.0){QNA}
\put(150.0,70.0){$St=512$} \end{overpic} (f)

Figure 2: Comparison of DNS data for $\varrho$ with the predictions from the
QNA and ACA for different $St$.

We now compare the model predictions with DNS for $\varrho$ and
$\mathcal{W}_{2}$ and for Stokes numbers $St=2.8,4.6,9.3,46.5,128,512$,
spanning particles with relatively weak to strong inertia. Figure 2 shows the
results for $\varrho$, and for $St=2.8,4.6,9.3$, the QNA and ACA models give
similar predictions that are in quite good agreement with the DNS data, with
the ACA performing slightly better for $St=2.8,4.6$, and the QNA performing
slightly better for $St=9.3$. Both models slightly underpredict $\varrho$ in
the range $7\lesssim z\lesssim 70$. While there are various possible
explanations for this, one is that the underpredictions are due to errors
introduced by the local approximation for $\lambda$, which can lead to errors
for particle transport in turbulent boundary layers [18].

Another possibility is that the underpredictions are due to errors in the
closure approximation $(1/St)\varrho\langle
u^{p}(t)\rangle_{z}\approx-\lambda\nabla_{z}\varrho$ that appears in the
transport equation governing $\varrho$. In general, the exact expression for
$(1/St)\varrho\langle u^{p}(t)\rangle_{z}$ would involve contributions from
infinitely many other terms involving diffusion coefficients and derivatives
of $\varrho$ of all orders [3]. That the closure $(1/St)\varrho\langle
u^{p}(t)\rangle_{z}\approx-\lambda\nabla_{z}\varrho$ only involves a gradient
term is a consequence of the assumption that the fluid velocity fluctuations
have Gaussian statistics [18]. In [15] it was demonstrated that a gradient
diffusion closure is asymptotically exact in the viscous sublayer. It might
also be expected to be reasonable in and beyond the log-region of a boundary
layer where deviations of the fluid velocity statistics from being Gaussian
are not expected to be strong. However, in the buffer region where the
turbulent production term peaks and where there are intense gradients, the
higher-order contributions to $(1/St)\varrho\langle u^{p}(t)\rangle_{z}$
arising from non-Gaussian fluid velocity fluctuations could be important. This
then could explain some of the descrepancies between the model and DNS results
for $\varrho$ observed in the region $7\lesssim z\lesssim 70$.

For $St=46.5,128,512$ the QNA has gone past the bifurcation $St$ value
discussed earlier, and its predictions for $\varrho$ are in serious error. Not
only does it drastically underpredict the values of $\varrho$ in the viscous
sublayer, but it also erroneously predicts that $\varrho$ is independent of
$z$ in this region. For $St=46.5$, the ACA model underpredicts the DNS data
for $\varrho$.

However, its predictions are closer to the DNS than the QNA, and most
importantly, the ACA predicts that for this $St$, $\varrho$ exhibits a power-
law dependence on $z$, in agreement with the DNS and asymptotic analysis of
[15], but which the QNA fails to reproduce. For $St=128$, the ACA predictions
for $\varrho$ are in very good agreement with the DNS down to $z\approx 4$,
below which the ACA underpredicts the DNS data. For $St=512$, the ACA is in
excellent agreement with the DNS across the range of $z$ considered. The
improvement of the ACA predictions as $St$ is increased is of course
consistent with the asymptotic nature of its closure approximation.

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{W_plot_2p8.pdf} \put(130.0,0.0){$z$}
\put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}} \put(163.0,62.0){DNS}
\put(163.0,51.0){ACA} \put(163.0,40.0){QNA} \put(163.0,29.0){Fluid}
\put(150.0,90.0){$St=2.8$} \end{overpic} (a)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt 199.16928pt,scale={0.4},clip,tics=20]{W_plot_4p6.pdf}
\put(130.0,0.0){$z$} \put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(163.0,62.0){DNS} \put(163.0,51.0){ACA} \put(163.0,40.0){QNA}
\put(163.0,29.0){Fluid} \put(150.0,90.0){$St=4.6$} \end{overpic} (b)

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{W_plot_9p3.pdf} \put(130.0,0.0){$z$}
\put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}} \put(163.0,62.0){DNS}
\put(163.0,51.0){ACA} \put(163.0,40.0){QNA} \put(163.0,29.0){Fluid}
\put(150.0,90.0){$St=9.3$} \end{overpic} (c)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt 199.16928pt,scale={0.4},clip,tics=20]{W_plot_46p5.pdf}
\put(130.0,0.0){$z$} \put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(163.0,62.0){DNS} \put(163.0,51.0){ACA} \put(163.0,40.0){QNA}
\put(163.0,29.0){Fluid} \put(150.0,90.0){$St=46.5$} \end{overpic} (d)

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{W_plot_128.pdf} \put(130.0,0.0){$z$}
\put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}} \put(163.0,62.0){DNS}
\put(163.0,51.0){ACA} \put(163.0,40.0){QNA} \put(163.0,29.0){Fluid}
\put(150.0,90.0){$St=128$} \end{overpic} (e)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt 199.16928pt,scale={0.4},clip,tics=20]{W_plot_512.pdf}
\put(130.0,0.0){$z$} \put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(163.0,62.0){DNS} \put(163.0,51.0){ACA} \put(163.0,40.0){QNA}
\put(163.0,29.0){Fluid} \put(150.0,90.0){$St=512$} \end{overpic} (f)

Figure 3: Comparison of DNS data for $\mathcal{W}_{2}(z)$ with the predictions
from the QNA and ACA

Taken together, these results show that the ACA model provides good
predictions for $\varrho$ for small and large $St$, but leads to some
underpredictions for intermediate $St$. For both intermediate and large $St$,
however, the new ACA model provides a significant improvement compared with
the traditional QNA, being not only in much better quantitative agreement with
the DNS, but also correctly capturing the power-law asymptotic behavior of
$\varrho$ in the viscous sublayer, which the QNA does not correctly predict
for intermediate and large $St$.

Figure 3 compares the QNA and ACA model predictions with the DNS data for
$\mathcal{W}_{2}(z)$. The DNS data for the fluid vertical Reynolds stress is
also shown for comparison, in order to highlight the extent to which the
models capture the effect of the particle inertia on the velocities. For
$St=2.8,4.6$ the QNA and ACA models are both in very good agreement with the
DNS, with the ACA model predictions in almost exact agreement with the DNS for
$St=4.6$, while the QNA model slightly underpredicts $\mathcal{W}_{2}(z)$ for
$z\lesssim 3$ when $St=4.6$. For $St=9.3$ the QNA and ACA models are both in
very good agreement with the DNS down to $z\approx 3$, but below this the QNA
model underpredicts the DNS, while the ACA model overpredicts the DNS. For
$St=46.5$ the QNA and ACA models are in good agreement with the DNS down to
around $z=10$, with both models capturing the strong effects of particle
inertia in this regime. However, this $St$ value exceeds the bifurcation $St$
value for the QNA model, and related to this is that the QNA predicts that
$\mathcal{W}_{2}(z)$ becomes constant with values that far exceed those of the
DNS at small $z$. The ACA model also significantly overpredicts the DNS at
small $z$, but the values are much closer to the DNS than those of the QNA
model. Most importantly, while the QNA model predicts that
$\mathcal{W}_{2}(z)$ becomes constant for small $z$, the ACA model preserves
the power-law like behavior observed in the DNS. This demonstrates then that
despite the quantiative shortcomings of the ACA model at $St=46.5$, it
significantly improves upon the QNA model in terms of preserving the right
kind of qualitative behavior. For $St=128$, the QNA model is accurate down to
around $z=30$, but below this its predictions are in enormous error compared
with the DNS, both quantitatively and qualitatively. By contrast, the ACA
model is accurate down to around $z=5$. Below this it overpredicts the DNS
data, but again preserves a power-like type behavior in this region, which is
in much better qualitative agreement with the DNS than the QNA model. Finally,
for $St=512$, the QNA model is only accurate down to around $z=100$, and
significantly over predicts the DNS data below this, while again predicting
that $\mathcal{W}_{2}(z)$ becomes constant at small $z$, in stark contrast to
the DNS.

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{Linear_W_plot_2p8.pdf}
\put(130.0,0.0){$z$} \put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(130.0,62.0){DNS} \put(130.0,51.0){ACA} \put(130.0,40.0){QNA}
\put(130.0,29.0){Fluid} \put(160.0,140.0){$St=2.8$} \end{overpic} (a)
\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{Linear_W_plot_4p6.pdf}
\put(130.0,0.0){$z$} \put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(160.0,140.0){$St=4.6$} \end{overpic} (b)

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{Linear_W_plot_9p3.pdf}
\put(130.0,0.0){$z$} \put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(160.0,140.0){$St=9.3$} \end{overpic} (c)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{Linear_W_plot_46p5.pdf}
\put(130.0,0.0){$z$} \put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(160.0,140.0){$St=46.5$} \end{overpic} (d)

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.4},clip,tics=20]{Linear_W_plot_128.pdf}
\put(130.0,0.0){$z$} \put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(160.0,140.0){$St=128$} \end{overpic} (e)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt 199.16928pt,scale={0.4},clip,tics=20]{Linear_W_plot_512.pdf}
\put(130.0,0.0){$z$} \put(0.0,75.0){\rotatebox{90.0}{$\mathcal{W}_{2}(z)$}}
\put(160.0,140.0){$St=512$} \end{overpic} (f)

Figure 4: Comparison of DNS data for $\mathcal{W}_{2}(z)$ with the predictions
from the QNA and ACA

The ACA model slightly overpredicts the DNS for $z>20$, but is in excellent
qualitative and quantiative agreement with the DNS below this. Again, the
improvement of the ACA predictions compared with the DNS as $St$ is increased
is consistent with the asymptotic nature of its closure approximation.

In figure 3 we again compare the QNA and ACA model predictions with the DNS
data for $\mathcal{W}_{2}(z)$, but this time in a linear scale in order to
highlight the behavior at larger values of $z$. At larger $z$, the QNA and ACA
predictions are almost identical, and this is because in these regions
$\varrho$ does not differ strongly from one, and if $\varrho$ were identically
equal to one, then the QNA and ACA closures would be identical. The results
show that at greater distances from the wall, e.g. $z\gtrsim 50$, the QNA and
ACA models are in general in good agreement with the DNS, with some
underpredictions for smaller $St$ that become smaller as $St$ is increased.
Comparing the DNS data for $\mathcal{W}_{2}(z)$ with the fluid Reynolds stress
shows that for the range of $St$ considered there is a strong effect of $St$
on $\mathcal{W}_{2}(z)$, and the models do a very good job of capturing this
effect of the particle inertia.

##  IV Conclusions

We have developed a new closure approximation for the moment equations
describing inertial particle transport in turbulent boundary layers that are
derived from an underlying phase-space PDF equation. Traditionally, a quasi-
Normal approximation (QNA) has been used to close the equations, but while
this yields good results when the particle Stokes number $St$ is sufficiently
small, it leads to significant errors for larger $St$, errors that are both
quantitative and qualitative in nature. We derive a new closure approximation
based on an asymptotic solution to the transport equations in regions where
the effect of particle inertia is significant. This new closure approximation
(referred to as the asymptotoc closure approximation, ACA) differs strongly
from the QNA closure in regions where the particle concentration $\varrho$
deviates strongly from being uniform, but asymptotes to the QNA approximation
when the concentration is uniform.

Comparisons of the model predictions for $\varrho$ and the variance of the
vertical particle velocity $\mathcal{W}_{2}$ with DNS data show that while the
QNA and ACA model make similar predictions at smaller $St$ that are in good
agreement with the DNS, their predictions differ dramatically at larger $St$.
The ACA model predictions are in good agreement with the DNS over a much wider
range of the boundary layer. At smaller distances from the wall, even when the
ACA model predictions are not in quantitative agreement with the DNS, they
correctly preserve the power-law like behavior of $\varrho$ and
$\mathcal{W}_{2}$, unlike the QNA model that erroneously predicts that these
functions become independent of $z$. For very large $St$, the ACA model is in
excellent quantitative agreement with the DNS data. The new ACA model
therefore dramatically improves on the traditional QNA model.

In order to address the remaining quantitative deficiencies of the ACA model,
two possibilities should be explore in future work. First, the coefficient
$\mathcal{C}_{4}$ that appears in equation (15) was obtained by enforcing that
in regions of the boundary layer where the particle inertia is weak that the
ACA closure asymptotes to the QNA closure. This yields a value for
$\mathcal{C}_{4}$ that is independent of $St$, whereas in reality it probably
should depend on $St$. Improving the specification of $\mathcal{C}_{4}$ to
include an appropriate $St$ dependence could improve the accuracy of the ACA
model. Second, the ACA closure is formally obtained as the leading order term
in an asymptotic series for the regime $St\gg 1$. It may be possible to
improve upon this by either incorporating the next term in the expansion, or
perhaps by using a renormalization approach to perform a partial summation of
some of the terms in the series. This will be explored in future work,
together with extensions of the model to include the effect of gravitational
settling.

##  V Acknowledgements

Support from the Army Research Office (award $\\#$ W911NF-22-2-0222) is
gratefully acknowledged. The computing resources and services used in this
work were provided by the VSC (Flemish Supercomputer Center), funded by the
Research Foundation - Flanders (FWO) and the Flemish Government.

## References

  * Rouse [1937] H. Rouse, Modern conceptions of the mechanics of turbulence, Trans. Am. Soc. Civ. Eng. 102, 463 (1937).
  * Richter and Chamecki [2018] D. H. Richter and M. Chamecki, Inertial effects on the vertical transport of suspended particles in a turbulent boundary layer, Boundary Layer Meteorol. 167, 235 (2018).
  * Bragg  _et al._ [2021] A. D. Bragg, D. H. Richter, and G. Wang, Mechanisms governing the settling velocities and spatial distributions of inertial particles in wall-bounded turbulence, Phys. Rev. Fluids 6, 064302 (2021).
  * Reeks [1991] M. W. Reeks, On a kinetic equation for the transport of particles in turbulent flows, Phys. Fluids 3, 446 (1991).
  * Zaichik [1997] L. Zaichik, Modelling of the motion of particles in non-uniform turbulent flow using the equation for the probability density function, J. Appl. Math. Mech. 61, 127 (1997).
  * Swailes and Darbyshire [1997] D. Swailes and K. Darbyshire, A generalized Fokker-Planck equation for particle transport in random media, Physica A 242, 38 (1997).
  * Reeks [2005] M. W. Reeks, On probability density function equations for particle dispersion in a uniform shear flow, J. Fluid. Mech. 522, 263–302 (2005).
  * Bragg  _et al._ [2012a] A. Bragg, D. C. Swailes, and R. Skartlien, Drift-free kinetic equations for turbulent dispersion, Phys. Rev. E 86, 056306 (2012a).
  * Maxey [1987] M. R. Maxey, The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields, J. Fluid Mech. 174, 441 (1987).
  * Tom and Bragg [2019] J. Tom and A. D. Bragg, Multiscale preferential sweeping of particles settling in turbulence, J. Fluid Mech. 871, 244 (2019).
  * Zaichik and Alipchenkov [2010] L. I. Zaichik and V. M. Alipchenkov, Modelling of transport and dispersion of arbitrary-density particles in turbulent flows, International Journal of Heat and Fluid Flow 31, 850 (2010), sixth International Symposium on Turbulence, Heat and Mass Transfer, Rome, Italy, 14-18 September 2009\.
  * Swailes  _et al._ [1998] D. Swailes, Y. Sergeev, and A. Parker, Chapman–enskog closure approximation in the kinetic theory of dilute turbulent gas-particulate suspensions, Physica A: Statistical Mechanics and its Applications 254, 517 (1998).
  * Bragg and Collins [2014a] A. Bragg and L. Collins, New insights from comparing statistical theories for inertial particles in turbulence: I. spatial distribution of particles., New J. Phys. 16, 055013 (2014a).
  * Bragg and Collins [2014b] A. Bragg and L. Collins, New insights from comparing statistical theories for inertial particles in turbulence: II. relative velocities of particles., New J. Phys. 16, 055014 (2014b).
  * Sikovsky [2014] D. P. Sikovsky, Singularity of inertial particle concentration in the viscous sublayer of wall-bounded turbulent flows, Flow, Turbulence and Combustion 92 (2014).
  * Johnson  _et al._ [2020] P. L. Johnson, M. Bassenne, and P. Moin, Turbophoresis of small inertial particles: theoretical considerations and application to wall-modelled large-eddy simulations, J. Fluid Mech. 883, A27 (2020).
  * Sikovsky [2019] D. P. Sikovsky, Particle reynolds stress model for wall turbulence with inertial particle clustering, Journal of Physics: Conference Series 1382, 012099 (2019).
  * Bragg  _et al._ [2012b] A. Bragg, D. C. Swailes, and R. Skartlien, Particle transport in a turbulent boundary layer: Non-local closures for particle dispersion tensors accounting for particle-wall interactions, Phys. Fluids 24, 103304 (2012b).
  * Minier and Peirano [2001] J.-P. Minier and E. Peirano, The pdf approach to turbulent polydispersed two-phase flows, Physics Reports 352, 1 (2001).
  * Pope [2000] S. B. Pope,  _Turbulent Flows_ (Cambridge University Press, New York, 2000).
  * Reeks  _et al._ [2018] M. Reeks, D. C. Swailes, and A. D. Bragg, Is the kinetic equation for turbulent gas-particle flows ill posed?, Phys. Rev. E 97, 023104 (2018).
  * van Dijk and Swailes [2012] P. van Dijk and D. Swailes, Hermite-dg methods for pdf equations modelling particle transport and deposition in turbulent boundary layers, J. Comput. Phys. 231, 4904 (2012).
  * Skartlien [2007] R. Skartlien, Kinetic modeling of particles in stratified flow – evaluation of dispersion tensors in inhomogeneous turbulence, Int. J. Multiph. Flow. 33, 1006 (2007).
  * Zaichik [1999] L. I. Zaichik, A statistical model of particle transport and heat transfer in turbulent shear flows, Phys. Fluids 11, 1521 (1999).
  * Marchioli  _et al._ [2008] C. Marchioli, A. Soldati, J. Kuerten, B. Arcen, A. Tanière, G. Goldensoph, K. Squires, M. Cargnelutti, and L. Portela, Statistics of particle dispersion in direct numerical simulations of wall-bounded turbulence: Results of an international collaborative benchmark test, Int. J. Multiph. Flow 34, 879 (2008).