Inhomogeneous preferential concentration of inertial particles in turbulent channel flow

TL;DR

This paper extends the theory of inertial particle clustering from homogeneous to inhomogeneous turbulent flows, specifically in channel flows, revealing how flow inhomogeneity enhances particle clustering.

Contribution

It develops a new theoretical framework for particle clustering in inhomogeneous turbulence and validates it with direct numerical simulations in turbulent channel flow.

Findings

Inhomogeneity significantly increases particle clustering.

Clustering is strongest near the viscous sublayer and buffer layer.

Inhomogeneous effects can amplify clustering compared to homogeneous turbulence.

Abstract

Turbophoresis leading to preferential concentration of inertial particles in regions of low turbulent diffusivity is a unique feature of inhomogeneous turbulent flows, such as free shear flows or wall-bounded flows. In this work, the theory for clustering of weakly inertial particles in homogeneous turbulence is extended to the inhomogeneous case of a turbulent channel flow. The inhomogeneity contributes to the cluster formation in addition to clustering in homogeneous turbulence. A space-dependent rate for the creation of inhomogeneous particle concentration is derived in terms of local statistics of turbulence. We provide the formula for the pair-correlation function of concentration that factorizes in product of time and space-dependent average concentrations and time-independent factor of clustering that obeys a power-law in the distance between the points. This power-law characterizes inhomogeneous multifractality of the particle distribution. A unique demonstration and quantification of the combined effects of turbophoresis and fractal clustering in a direct numerical simulation of particle motion in a turbulent channel flow is performed according to the presented theory. The strongest contribution to clustering coming from the inhomogeneity of the flow occurs in the transitional region between viscous sublayer and the buffer layer. Further the ratio of homogeneous and inhomogeneous term depends on the wall distance. The inhomogeneous terms may significantly increase the preferential concentration of inertial particles, thus the overall degree of clustering in inhomogeneous turbulence is potentially stronger compared to particles with the same inertia in purely homogeneous turbulence.

Tables

Table 1. Table 1: Simulation, flow and grid parameters - L i subscript 𝐿 𝑖 L_{i} : domain size in all directions i = x , y , z 𝑖 𝑥 𝑦 𝑧 i=x,y,z , N i subscript 𝑁 𝑖 N_{i} : number of grid points in all directions, δ ​ t 𝛿 𝑡 \delta t : time step, ν 𝜈 \nu : viscosity, U c subscript 𝑈 𝑐 U_{c} : centerline velocity, u τ subscript 𝑢 𝜏 u_{\tau} : friction velocity, R ​ e τ 𝑅 subscript 𝑒 𝜏 Re_{\tau} : friction Reynolds number, S ​ t + 𝑆 superscript 𝑡 St^{+} : viscous Stokes number, Δ ​ x + Δ superscript 𝑥 \Delta x^{+} : grid spacing streamwise direction, Δ ​ y 1 + Δ superscript subscript 𝑦 1 \Delta y_{1}^{+} : grid spacing wall-normal direction (first point), Δ ​ y c + Δ superscript subscript 𝑦 𝑐 \Delta y_{c}^{+} : grid spacing wall-normal direction (central point), Δ ​ z + Δ superscript 𝑧 \Delta z^{+} : grid spacing spanwise direction

$L_x\times L_y\times L_z$	$N_x\times N_y\times N_z$	$\delta t$	$\nu$	$U_c$	$u_{\tau }$	$R{e}_{\tau }$	$S{t}^{+}$	$\Delta x^{+}$	$\Delta y_1^{+}$	$\Delta y_c^{+}$	$\Delta z^{+}$
$8h\pi \times 2h\times 3h\pi$	$2048\times 512\times 1536$	$0.0065$	$5\times {10}^{-5}$	$1.13$	$0.0499$	$999.4$	$1$	$12.26$	$0.0165$	$6.16$	$6.13$

Equations

$$⟨\sum {\lambda }_i(y)⟩=-{\tau }^2{\int }_0^{\infty }⟨{\nabla }^{2}p(0){\nabla }^{2}p(t){⟩}_{c}d{t}^{\prime }+{\partial }_y{v}_{eff}(y),v_{eff}(y)=-\tau {\partial }_y⟨u_y^2⟩+\tau {\int }_0^{\infty }⟨u_y(0){\nabla }^{2}p(t){⟩}_{c}dt.$$

$$\frac{d\mathbf{x}}{dt}=\mathbf{v},\frac{d\mathbf{v}}{dt}=-\frac{\mathbf{v}-\mathbf{u}[t,\mathbf{x}(t)]}{\tau }.$$

$$\mathbf{v}=\mathbf{u}-\tau [\partial \mathbf{u}/\partial t+(\mathbf{u}\cdot \nabla )\mathbf{u}],$$

$$⟨v_y⟩=-\tau ⟨\nabla \cdot (\mathbf{u}{u}_y)⟩=-\tau {\partial }_y⟨u_y^2⟩⟨\nabla \cdot \mathbf{v}⟩=\frac{\partial ⟨v_y⟩}{\partial y}=-\tau \frac{{\partial }^2⟨u_y^2⟩}{\partial y^2}.$$

$$\frac{d\ln V}{dt}=w(t,\mathbf{x}(t)),w(t,\mathbf{x})=\nabla \cdot \mathbf{v}(t,\mathbf{x}),$$

$$\sum _{i=1}^3{\lambda }_i(t)=\frac{1}{t}\ln \left(\frac{V(t)}{V(0)}\right)=\frac{1}{t}{\int }_0^{t}w(t^{\prime },\mathbf{x}(t^{\prime }))d{t}^{\prime }.$$

$${\left(\sum {\lambda }_i\right)}_{bulk}=\underset{t\to \infty }{\lim }\sum {\lambda }_i(t)=\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^t⟨w(t^{\prime },\mathbf{x}(t^{\prime }))⟩d{t}^{\prime }=\underset{t\to \infty }{\lim }⟨w(t,\mathbf{x}(t))⟩,$$

$${\left(\sum {\lambda }_i\right)}_{bulk}=-{\int }_0^{\infty }⟨w(0)w(t)⟩dt,$$

$${\left(\sum {\lambda }_i\right)}_{bulk}=-{\tau }^2{\int }_0^{\infty }⟨{\nabla }^{2}p(0){\nabla }^{2}p(t)⟩dt\propto S{t}^2,$$

$${\lambda }_1(t)+{\lambda }_2(t)=\frac{1}{t}\ln \left(\frac{S(t)}{S(0)}\right),$$

$$⟨{\left(\sum {\lambda }_i(t)\right)}^2⟩-{⟨\sum {\lambda }_i(t)⟩}^2=\frac{1}{t^2}{\int }_0^{t}d{t}_{1}d{t}_2⟨w(t_1,\mathbf{x}(t_1))w(t_2,\mathbf{x}(t_2)){⟩}_c\sim \frac{1}{t}{\int }_0^{\infty }⟨w(0)w(t)⟩dt\sim \frac{|\sum {\lambda }_i|}{t},t\gg {\tau }_{\eta },$$

$$\frac{\sqrt{⟨{\left(\sum {\lambda }_i(t)\right)}^2⟩-{⟨\sum {\lambda }_i(t)⟩}^2}}{|\sum {\lambda }_i(t)|}\sim \frac{1}{\sqrt{|\sum {\lambda }_i|t}}\sim \frac{1}{St\sqrt{t/{\tau }_{\eta }}}.$$

$${\left(\sum {\lambda }_i\right)}_{bulk}(y)\sim -{\tau }^2{\int }_0^{\infty }⟨{\nabla }^{2}p(0,y){\nabla }^{2}p(t)⟩dt,$$

$$⟨{\nabla }^{2}p(0,y){\nabla }^{2}p(t)⟩=\underset{S\to \infty }{\lim }\int \frac{dxdz}{S}{\nabla }^{2}p(0,x,y,z)){\nabla }^{2}p[t,\mathbf{q}(t,x,y,z)],$$

$${\partial }_t\mathbf{q}(t|t_0,\mathbf{x})=\mathbf{v}[t,\mathbf{q}(t|t_0,\mathbf{x})],\mathbf{q}(t|t_0,\mathbf{x})=\mathbf{x},\mathbf{q}(t,\mathbf{x})=\mathbf{q}(t|t_0=0,\mathbf{x}).$$

$$⟨{\nabla }^{2}p(0,y){\nabla }^{2}p(t){⟩}_t=\underset{S\to \infty }{\lim }\int \frac{dxdz}{S}⟨{\nabla }^{2}p(0){\nabla }^{2}p(t){⟩}_t$$

$$⟨f(t,y)⟩=\int \frac{dxdz}{S}f(t,\mathbf{q}(t,x,y,z)).$$

$${\partial }_{t_0}\mathbf{q}(t|t_0,\mathbf{x})+[\mathbf{v}(t_0,\mathbf{x})\cdot \nabla ]\mathbf{q}(t|t_0,\mathbf{x})=0.$$

$$⟨f(t,t_0,y)⟩=\underset{S\to \infty }{\lim }\int \frac{dxdz}{S}f[t,\mathbf{q}(t|t_0,x,y,z)].$$

$$\frac{\partial ⟨f⟩}{\partial t_0}=-\underset{S\to \infty }{\lim }\int \frac{dxdz}{S}{v}_i(t_0,x,y,z){\nabla }_{i}f[t,\mathbf{q}(t|t_0,x,y,z)]$$

$$=-\underset{S\to \infty }{\lim }\int \frac{dxdz}{S}{\nabla }_i\left(v_i(t_0,x,y,z)f[t,\mathbf{q}(t|t_0,x,y,z)]\right)$$

$$+⟨w(0)f(t-t_0)⟩,$$

$$⟨w(0)f(t)⟩=\underset{S\to \infty }{\lim }\int \frac{dxdz}{S}w(0,x,y,z))f[t,\mathbf{q}(t,x,y,z)].$$

$$\frac{\partial ⟨f(t,t_0,y)⟩}{\partial t_0}=⟨w(0)f(t-t_0)⟩-{\partial }_y⟨v_y(0)f(t-t_0)⟩,$$

$$⟨v_y(0)f(t)⟩=\underset{S\to \infty }{\lim }\int \frac{dxdz}{S}{v}_y(0,x,y,z))f[t,\mathbf{q}(t,x,y,z)].$$

$$⟨f(t,y)⟩-⟨f(y)⟩={\int }_0^t\left[{\partial }_y⟨v_y(0)f(t^{\prime })⟩-⟨w(0)f(t^{\prime })⟩\right]d{t}^{\prime }.$$

$$⟨v_y(0)f(t^{\prime })⟩=⟨v_y(0)f(t^{\prime }){⟩}_c+⟨v_y(0)⟩⟨f(t^{\prime })⟩,$$

$$⟨w(0)f(t^{\prime })⟩=⟨w(0)f(t^{\prime }){⟩}_c+⟨w(0)⟩⟨f(t^{\prime })⟩,$$

$$⟨f(t,y)⟩=⟨f(y)⟩+{\int }_0^t\left[{\partial }_y⟨v_y(0)f(t^{\prime }){⟩}_c-⟨w(0)f(t^{\prime }){⟩}_c\right]d{t}^{\prime }-\tau ({\partial }_y⟨u_y^2⟩){\partial }_y{\int }_0^t⟨f(t^{\prime })⟩d{t}^{\prime }.$$

$$⟨\sum {\lambda }_i(t,y)⟩=\frac{1}{t}{\int }_0^t⟨w(t^{\prime },y)⟩d{t}^{\prime }=\frac{1}{t}{\int }_0^{t}d{t}^{\prime }\int \frac{dxdz}{S}w(t^{\prime },\mathbf{q}(t^{\prime },x,y,z)).$$

$$⟨w(t,y)⟩=\tau {\partial }_y{\int }_0^t⟨u_y(0){\nabla }^{2}p(t^{\prime }){⟩}_{c}d{t}^{\prime }-{\tau }^2{\int }_0^t⟨{\nabla }^{2}p(0){\nabla }^{2}p(t^{\prime }){⟩}_{c}d{t}^{\prime }-{\tau }^2({\partial }_y⟨u_y^2⟩){\partial }_y{\int }_0^t⟨{\nabla }^{2}p(t^{\prime })⟩d{t}^{\prime }-\tau {\partial }_y^2⟨u_y^2⟩.$$

$$⟨w(t,y)⟩=-{\tau }^2{\int }_0^{\infty }⟨{\nabla }^{2}p(0){\nabla }^{2}p(t){⟩}_{c}d{t}^{\prime }+{\partial }_y{v}_{eff}(y),v_{eff}(y)=⟨v(y)⟩+{\int }_0^{\infty }⟨u_y(0)w(t){⟩}_{c}dt,t\gg {\tau }_{\eta }(y).$$

Full text

I Introduction

Particulate matter suspended in fluids is a common feature of turbulent flows encountered in industrial devices as well as the environment. Typical examples range from combustion processes in diesel engines and gas turbines Sirignano (1999); Lee et al. (2002), ash or aerosols expelled from volcanic eruptions or chemical or nuclear accidents Crowe et al. (1998); Hidy (2012), to liquid rain droplets in clouds Falkovich et al. (2002); Shaw (2003). Knowledge about the distribution of these particles can be essential for the working reliability and efficiency of engines Williams and Crane (1979), climate predictions Seinfeld and Pandis (2012); Reeks (2014) and the health of living organisms Flagan and Seinfeld (2013). Therefore, understanding of particle-turbulence interaction causing heterogeneous particle distributions and identifying potential high-concentration regions is a crucial aspect. Clustering of inertial particles in homogeneous turbulence conditions has been subject to intensive research for a long time and is thus generally well explained by the phenomenon of small-scale or fractal clustering, see for example Bec et al. (2007a); Sundaram and Collins (1997); Bec et al. (2007b); Jung et al. (2008); Balkovsky et al. (2001); Cencini et al. (2006); Calzavarini et al. (2008); Monchaux et al. (2010); Eaton and Fessler (1994).

Many flows in the environment or industrial applications are physically bounded by walls, which creates a heterogeneity in the flow. In these non-uniform flows, inertial particles are subject to turbophoretic forces that create a strongly inhomogeneous distribution of particles along the direction of non-uniformity. Turbophoresis is generally associated with a net flux of particles towards regions of lower turbulence diffusivity, thus leading to increased particle concentration in near-wall regions of wall-bounded flows. Since inhomogeneous turbulent flows containing non-passive particles are a common occurrence in nature and industry, turbophoresis is an ubiquitous phenomenon that deserves special attention. The turbophoretic mechanism was initially described from a theoretical point of view in Caporaloni et al. (1975) and Reeks (1983). The work of Young and Leeming (1997) developed a theoretical description of the particle transfer to the wall caused by turbulent sweep and ejection events. Subsequently, several experimental investigations (e.g. Kaftori et al. (1995a, b); Righetti and Romano (2004)) have confirmed the accumulation of particles in near-wall regions due to turbophoresis. Increasing computational resources allowed for detailed investigations of turbophoresis based on direct numerical simulations of inhomogeneous turbulent flows. The early work of Eaton & Fessler Eaton and Fessler (1994) found heavy particles to be localized preferably in regions close to the wall that feature a low instantaneous velocity. Further investigations Marchioli and Soldati (2002); Picciotto et al. (2005) indicated that coherent wall structures are very likely responsible for high particle concentrations in low velocity regions in wall-proximity. Furthermore, the spatial development of particle concentrations in a turbulent pipe flow was investigated numerically in Picano et al. (2009). These investigations mainly focused on the increasing concentration of particles with relatively large inertia in the near-wall region and the possibility to derive information about particle surface deposition. Preferential concentration arising due to turbophoresis as well as small-scale clustering has been investigated in a channel flow for heavy particles by Sardina et al. (2012). De Lillo et al. (2016) studied the effects of turbophoresis and small-scale clustering in a shear-flow without walls for a wide range of particle inertia. They showed that turbophoresis is stronger for particles with moderate inertia, whereas small-scale clustering dominates at week inertia.

Multifractality was derived originally for spatially uniform turbulence. It is a small-scale phenomenon holding at scales smaller than Kolmogorov scale. Recently, multifractality was generalized to the case of inhomogeneous turbulence Schmidt et al. (2016). The derivation assumed that during the characteristic time of formation of fractal structures the motion of small parcels of particles is confined in a not too large region of space where statistics of flow gradients can be considered uniform. It was demonstrated that the pair-correlation function of concentration of particles factorizes in the product of (possibly time-dependent) average local concentrations and a geometrical factor of fractal increase of probability of two particles to be close. Thus we can separate dependencies in the particle distribution. The geometry of the time-dependent multifractal to which the particles are confined is statistically stationary. It is space-dependent because of the inhomogeneity of turbulence. However, the overall number of particles, that distribute over the multifractal locally, is determined by the average local concentration that can be time-dependent. For instance in the case of turbophoresis there is depletion of local average concentration because of the particle flux to the wall.

The present work provides an extensive theoretical and numerical investigation of the clustering degree of weakly inertial particles in a turbulent channel flow. The specific case of interest here is the turbophoretic behavior of weakly inertial particles in flows where the turbulent diffusivity varies in one direction (i.e. wall-normal) but the turbulence is homogenous in wall-parallel planes. The universal framework of weakly compressible flow Falkovich et al. (2002); Fouxon (2012); Fouxon et al. (2015) is therefore modified to describe the combined clustering effects occurring due to the inhomogeneity of the flow and fluctuational clustering present in homogeneous turbulence. The theory is used as basis for the analysis of inertial particles in a direct numerical simulation (DNS) of a turbulent channel flow.

Our main theoretical result in this paper is the formula for the space-dependent rate $\left\langle\sum\lambda_{i}(y)\right\rangle$ of creation of inhomogeneities of concentration of particles in terms of local statistics of turbulence,

[TABLE]

In this formula $y$ is the distance to the wall, $\tau$ is the Stokes time of the particles, $\bm{u}$ and $p$ are the turbulent flow velocity and pressure, respectively. The angular brackets designate averaging over local statistics of turbulence and brackets with $c$ standing for the dispersion. Besides the average Eulerian velocity of the particles $-\tau\partial_{y}\langle u_{y}^{2}\rangle$, the rest of the terms come from distinction between Lagrangian and Eulerian averages in the formation of particle inhomogeneities. The first term is the local form of the formula known for spatially uniform turbulence that holds in the bulk. Inhomogeneity of turbulence produces an effective velocity term $v_{eff}(y)$ that in contrast with the bulk term is proportional to the first power of $\tau$. Thus clustering of weakly inertial particles can be much stronger because of inhomogeneity effects.

The numerical investigation of this work is based on DNS data obtained from the Johns Hopkins University Turbulence Database (JHTDB). We implement a Lagrangian tracking algorithm for inertial particles based on time-resolved Eulerian DNS results of a turbulent channel flow. This allows to evaluate a wide range of flow quantities along particle trajectories which serves as basis for the quantification of preferential concentration of inertial particles in inhomogeneous turbulence.

This paper commences with a theoretical analysis (Sec. II) and the equations used for the analysis of the numerical data. This is followed by a description of the channel flow simulation details, including the implementation of the particle tracking algorithm and a brief characterization of the flow quantities of interest (Sec. VI). Subsequently, the performed analysis and results based on the direct numerical simulation (Sec. VII) are illustrated. A summary of these results and the corresponding conclusions are presented in Sec. VIII.

II Sum of Lyapunov exponents

In this Section, we introduce the local rate of production of inhomogeneities of concentration of particles. This rate is the sum of Lyapunov exponents $\sum\lambda_{i}$. We observe that the local rate of particle density increase is determined by the local divergence of velocity. However that divergence is fluctuating in turbulence. The sum of Lyapunov exponents describes the rate of accumulated growth obtained by proper averaging over the fluctuations and it differs both from fluctuating and average divergences. We will see in the coming Sections that it is the average divergence plus the coherent contribution of the fluctuations.

We consider particles with weak inertia suspended in incompressible turbulent flow in a channel. The strength of inertia of a small spherical particle is quantified by the particle Stokes time $\tau$. This is given by $\tau=2a^{2}\rho_{p}/(9\nu\rho)$, with $a$ being the particle radius, $\rho$ and $\rho_{p}$ as the fluid and particle density respectively and $\nu$ is the fluid viscosity. We consider the case of dense particles with $\rho_{p}\gg\rho$ so effects such as memory and added mass can be neglected Maxey and Riley (1983). Assuming that the Reynolds number $Re_{p}$ of the flow perturbation caused by the particles is small ($Re_{p}\ll 1$) we can use the linear law of friction finding the equation of motion,

[TABLE]

Here, $\bm{x}(t)$ is the particle coordinate, $\bm{v}(t)$ is the velocity and $\bm{u}(t,\bm{x})$ is the turbulent flow velocity. Since we consider particles with small inertia, the velocity may be approximated as Maxey (1987),

[TABLE]

where $\tau[\partial\bm{u}/\partial t+(\bm{u}\cdot\bm{\nabla})\bm{u}]$ represents the inertial particle drift.

Despite that the fluid velocity is divergence-free, the inertial drift of the particle velocity results in a ’weakly compressible’ particle flow because $\nabla\cdot\bm{v}\neq 0$. For the case where particles are seeded into a statistically stationary, incompressible channel flow, with two homogeneous $(x,z)$ and one inhomogeneous ($y$) flow direction the mean Eulerian divergence of the particle velocity field reduces to,

[TABLE]

The formula for $\langle v_{y}\rangle$, obtained by averaging Eq. (3), is the well-known result for the average Eulerian particle velocity towards the wall. Since $\langle u_{y}^{2}\rangle$ monotonously grows away from the wall reaching a maximum in the log-law region of the channel, then $\langle v_{y}\rangle$ is directed toward the walls constituting the simplest demonstration of turbophoretic motion of particles towards the regions with smaller intensity of turbulence. In contrast, $\langle\nabla\cdot\bm{v}\rangle$ is not sign-definite: near the wall where $\langle u_{y}^{2}\rangle$ has minimum we have positive $\partial_{y}^{2}\langle u_{y}^{2}\rangle$ (right on the wall $\partial_{y}^{2}\langle u_{y}^{2}\rangle=0$ but it is positive nearby) and negative $\langle\nabla\cdot\bm{v}\rangle$, conversely we have $\partial_{y}^{2}\langle u_{y}^{2}\rangle<0$ and $\langle\nabla\cdot\bm{v}\rangle>0$ already in the buffer layer. Thus there is a point $y_{*}$ where $\langle\nabla\cdot\bm{v}\rangle=0$ so $\langle\nabla\cdot\bm{v}\rangle$ is negative for $y<y_{*}$ and positive otherwise.

Based on Eulerian divergence we could conclude that at $y<y_{*}$ there are positive correlations in positions of particles separated by distance smaller than the viscous scale $\eta$. These particles are in the same divergence Frisch (1995) that is negative on average and thus approach each other. In contrast, for $y>y_{*}$ there would seem to be negative correlations of positions. In reality the particles move and correlate in reaction not to Eulerian divergence of the flow but to the divergence of the flow in the frame of reference that moves with the particle. This divergence determines the evolution of infinitesimal volumes $V$ of particles that obey Batchelor (1959),

[TABLE]

where $\bm{x}(t)$ is the trajectory of some particle located inside the considered infinitesimal volume. This equation holds provided the largest linear size of $V$ is much smaller than $\eta$. The solution of this equation for the logarithmic increase rate of $V(t)$ defines the finite-time sum of Lyapunov exponents $\sum\lambda_{i}(t)$ as,

[TABLE]

The RHS of this equation has the form that occurs in the ergodic theorem Sinai (1972). Indeed, in spatially uniform turbulence ergodicity implies that the RHS converges to the deterministic limit in the limit of $t\to\infty$ (we remark that the convergence holds for almost every point with respect to the stationary measure but in this context it holds for almost every point in space that is with possible exception of points with zero total volume). This limit is independent of the initial position of $V(t)$. Thus the limit equals to the average over the initial position defining the sum of Lyapunov exponents of spatially uniform turbulence Devenish et al. (2012); Falkovich and Fouxon (2004),

[TABLE]

where angular brackets stand for averaging over $\bm{x}(0)$. Here the subscript signifies that this formula holds in the bulk of the flow, far from the boundaries, where the spatially uniform turbulence framework applies. It was demonstrated in Falkovich and Fouxon (2004), and will be proved below somewhat differently, that this can be written as,

[TABLE]

which clarifies that $\sum\lambda_{i}<0$ see details in Falkovich and Fouxon (2004). Taking the divergence of the Navier-Stokes equations obeyed by $\bm{u}$ we find that $w=\tau\nabla^{2}p$ where $p$ is the turbulent pressure. Thus we can write

[TABLE]

where we defined the Stokes number $St=\tau\sqrt{\epsilon/\nu}$. Here $\epsilon$ is the average rate of energy dissipation per unit mass.

In practice, the infinite time limit must be understood as near convergence beyond a certain convergence time $t_{c}$ when for most of the trajectories $\bm{x}(t)$, the LHS of Eq. (6) becomes constant at $t\sim t_{c}$. This time is strongly different for $\sum\lambda_{i}$ and other combinations of Lyapunov exponents. These are defined very similarly where $\lambda_{1}$ and $\lambda_{1}+\lambda_{2}$ are logarithmic increase rates of infinitesimal line and surface elements, respectively Devenish et al. (2012). For instance

[TABLE]

where $S(t)$ is the infinitesimal area of fluid particles Devenish et al. (2012); Schmidt et al. (2016). For $\lambda_{1}$ and $\lambda_{2}$ the convergence time $t_{c}$ is a few Kolmogorov times $\tau_{\eta}$. Here $\tau_{\eta}=\sqrt{\nu/\epsilon}$, which is usually defined as typical time-scale of turbulent eddies at the viscous scale Frisch (1995), is also the correlation time of flow gradients (and thus $w$) in the fluid particle frame Devenish et al. (2012); Falkovich and Fouxon (2004). In contrast, for $\sum\lambda_{i}$ the convergence time is much longer as seen considering the dispersion the of $\sum\lambda_{i}(t)$,

[TABLE]

where $c$ stands for cumulant (dispersion) and we used that the correlation time $\tau_{\eta}$ of $w(t,\bm{x}(t))$ is much smaller than $t$. We observe that the normalized rms deviation of $\sum\lambda_{i}(t)$ obeys,

[TABLE]

Thus fluctuations of $\sum\lambda_{i}(t)$ are small only at quite large times $t\gg\tau_{\eta}/St^{2}$. This is because non-zero $\sum\lambda_{i}$ is due to fluctuations, as evident from Eq. (8), so dispersion of $\sum\lambda_{i}(t)$ has the same order in $St$ as the average. Thus the convergence time for the sum of Lyapunov exponents obeys $t_{c}\sim\tau_{\eta}/St^{2}$. This slow convergence of long-time limit $\sum\lambda_{i}(t)$ at small $St$, seems to be unobserved previously.

Clearly the long-time convergence of $t^{-1}\ln V(t)/V(0)$ to a deterministic limit makes it a very useful quantity because it allows a deterministic prediction despite the randomness of turbulence. We determine $\sum\lambda_{i}$ once and then we can predict it for any arbitrary trajectory. The sum of the Lyapunov exponents, made dimensionless by dividing by a factor of order of $\sqrt{\epsilon/\nu}$ determines the strength of the clustering by giving fractal dimensions of the random attractor formed by particles in space, see e. g. Fouxon (2012) and Section V.

Fortunately, we can extend the notion of the sum of Lyapunov exponents to the inhomogeneous case with which we can describe clustering closer to the walls. We consider the situation when $\bm{x}(t_{c})$ does not deviate from $\bm{x}(0)$ by the characteristic scale of inhomogeneity of turbulent statistics which is the distance to the wall $y$. This demands that the $y-$dependent inhomogeneity time $t_{in}(y)=y/\langle u_{y}^{2}\rangle^{1/2}$, where $\langle u_{y}^{2}\rangle^{1/2}$ is the typical transversal velocity, is much larger than $t_{c}$. Here $t_{in}(y)$ is the time during which the trajectory passes a distance comparable with $y$ over which the statistics changes. If $t_{c}\ll t_{in}(y)$ holds then during the convergence time $\bm{x}(t)$ in Eq. (6) stays in the region where the statistics is roughly uniform. Thus we can use the results for spatially uniform turbulence locally. We could then guess that,

[TABLE]

where the RHS depends on the $y$ coordinate. Here the local correlation function can be defined with the help of temporal averaging. In the case of the channel instead of time averaging we can use averaging over the symmetry plane of the statistics,

[TABLE]

where $S$ is the area of the $x-z$ plane and we introduced ’Lagrangian’ trajectories of particles labeled by their initial positions at $t=t_{0}$,

[TABLE]

The temporal and ”planar” averages are identical because the temporal average

$\langle\nabla^{2}p(0,y)\nabla^{2}p(t)\rangle_{t}$ is independent of $x$ and $z$ so that,

[TABLE]

where interchanging the order of averages in the last term proves that time and plane averages coincide (time average of the plane average is the plane average). Below we will use planar averaging designating it by angular brackets and demonstrate that the guess given by Eq. (13) is incomplete.

III Identity for Lagrangian averages

In this Section we derive an identity for Lagrangian average $\langle f(t,y)\rangle$ of an arbitrary stationary random function $f(t,\bm{x})$ in the frame of particles released at the same distance $y$ from the wall,

[TABLE]

Here this quantity is of interest in the case of $f=w(t,\bm{x})$ when it provides the average sum of Lyapunov exponents studied in the next Section. However other cases of this quantity can be of interest in future studies so we keep arbitrariness of $f$. The derivation is a changed line of thought that appeared in Falkovich and Fouxon (2004). We observe that $\bm{q}(t|t_{0},\bm{x})$ as a function of the initial time $t_{0}$ obeys,

[TABLE]

This expresses that changing initial time and position so that the initial position stays on the same trajectory does not change that trajectory: $\bm{q}(t|t_{0}+\epsilon,\bm{q}(t_{0}+\epsilon|t_{0},\bm{r}))=\bm{q}(t|t_{0},\bm{r})$ (the trajectory that passes through $\bm{r}$ at time $t_{0}$ is the same trajectory that passes through $\bm{q}(t_{0}+\epsilon|t_{0},\bm{r})$ at time $t_{0}+\epsilon$). Differentiating over $\epsilon$ and setting $\epsilon=0$ one finds the equation above.

We introduce the two-time version of average of $w$,

[TABLE]

Because of stationarity $\langle f(t,t_{0},y)\rangle$ depends on $t$, $t_{0}$ only through the difference of the time arguments $t-t_{0}$. We consider the time derivative of Eq. (19) over $t_{0}$ using Eq. (18),

[TABLE]

where we used stationarity and defined the correlation function of $w$,

[TABLE]

Finally, observing that derivatives over $x$ and $z$ in Eq. (20) give zero as integrals of complete derivative and taking the $y-$derivative outside the integration we find

[TABLE]

where we defined the correlation function,

[TABLE]

We find integrating Eq. (22) over $t_{0}$ from $t_{0}=0$ up to $t_{0}=t$ and using $\langle f(t,y)\rangle=\langle f(t,t_{0}=0,y)\rangle$ that

[TABLE]

Here $\langle f(y)\rangle=\langle f(t=0,y)\rangle$ is the Eulerian average of $f(\bm{x})$ over the horizontal plane. The RHS describes the difference between Lagrangian and Eulerian averages that holds because of preferential concentration of particles and thus quantifies the strength of the clustering. This integrals in Eq. (24) do not necessarily converge in the long-time limit. We separate the possibly divergent term introducing cumulants (dispersion)

[TABLE]

where the angular brackets stand for average over $x-z$ coordinates. We find using $\partial_{y}\langle v_{y}(0)\rangle=\langle w(0)\rangle$ that

[TABLE]

where we used $\langle v_{y}(0)\rangle$ provided by Eq. (4). The last term does not necessarily converge in $t\to\infty$ limit. For instance if $\langle f(t)\rangle$ has a finite long-time limit then this term grows linearly with time. We consider this identity in the case of our interest.

IV Space-dependent sum of Lyapunov exponents for channel turbulence

In this Section we study the average sum of Lyapunov exponents. That describes the average logarithmic rate of growth of infinitesimal volumes that start at the same distance from the wall $y$. Performing averaging over initial $x-z$ coordinates of the volume we find from Eq. (6),

[TABLE]

The identity given by Eq. (26) gives in the leading order in $St$ for $f=w$,

[TABLE]

where we used $\langle w\rangle$ given by Eq. (4). In different time correlation functions in this formula we can use trajectories of the tracers in the leading order in $St$. The first term in the integrand is of order $St$ while the second is of order $(St)^{2}$. However both terms have to be kept because the former vanishes in the bulk due to spatial uniformity. In contrast, the third term in the integrand which is of order $(St)^{2}$ is smaller than the first term and can be neglected. The remaining integrals have finite $t\to\infty$ limit. We make the plausible assumption that the integrals converge over the local Kolmogorov time-scale $\tau_{\eta}(y)=\sqrt{\nu/\epsilon(y)}$. Then we find that,

[TABLE]

Finally, we find Eq. (1) using Eq. (27). This formula is valid at not too large times because the third term in Eq. (28) grows with time with $\tau\langle\nabla^{2}p\rangle\approx\sum\lambda_{i}=const$ in the integrand. This term becomes non-negligible at times $t$ obeying $\tau t\partial_{y}\langle u_{y}^{2}\rangle\partial_{y}\ln|\sum\lambda_{i}(y)|\gtrsim 1$. Since the scale of variations of the involved quantities is $y$ then this gives equality $\tau t\langle(u_{y}/y)^{2}\rangle\gtrsim 1$. For spatially uniform turbulence, applicable in the bulk, using Kolmogorov theory we find that $\langle(u_{y}/y)^{2}\rangle\sim\epsilon^{2/3}y^{-4/3}$ is smaller than $\langle(\nabla u)^{2}\rangle$ by a factor of the Reynolds number $Re$. This would give $t\gtrsim(Re/St)\tau_{\eta}$. If we assume that this time is much larger than $t_{c}$ which demands $Re\gg St^{-1}$ then Eq. (1) describes $\sum\lambda_{i}(t)$ for one trajectory so we can use $\sum\lambda_{i}(t)\approx\langle\sum\lambda_{i}(t)\rangle$. A similar consideration can be made in the near wall region of small $y$ where $\langle u_{y}^{2}\rangle\propto y^{4}$. Below we assume that the correction term is negligible in regions of interest (as confirmed by results of a DNS in section VII).

The formula given by Eq. (1) has reductions in the bulk and in the turbulent boundary layer. In the bulk the statistics is uniform and we recover Eq. (9) where there is no dispersion sign in the average because $\langle\nabla^{2}p\rangle=0$ in the bulk. In spatially uniform turbulence the sum of Lyapunov exponents is of order $St^{2}$. In contrast, the inhomogeneous terms are proportional to $St$ and dominate regions of strong inhomogeneity,

[TABLE]

We see that the RHS has the form of the divergence of an effective velocity which is the average Eulerian velocity of particles plus the correction. We can interpret $v_{eff}$ as the average velocity of turbophoretic particles to the wall which differs from the average Eulerian velocity $\langle v_{y}\rangle$ because of the difference between Lagrangian and Eulerian averages. The reason $\left(\left\langle\sum\lambda_{i}\right\rangle\right)_{inhom}$ scales linearly with $St$ is that in the bulk non-zero $\sum\lambda_{i}$ appears because of fluctuations and thus is proportional to $St^{2}$ but in the boundary layer there is an average effect proportional to $St$.

The complete equation (28) is an identity that holds for arbitrary $y$ including those close to the wall. In the passages between this identity and Eq. (1) we introduced the assumption that the third term in Eq. (28) can be neglected and the integrals converge over the local Kolmogorov time-scale $\tau_{\eta}(y)$. In the case of inhomogeneous statistics the integrands can have non-trivial time-dependence because the trajectory samples regions with statistics different from that at the initial point. Nevertheless the convergence seems reasonable.

We consider the question of how well $\left\langle\sum\lambda_{i}(y)\right\rangle$ approximates the fluctuating finite-time Lyapunov exponent $\sum\lambda_{i}(t)=t^{-1}\ln V(t)/V(0)$ for volumes whose initial vertical position is $y$ and horizontal position is arbitrary. Proceeding as we did in studying the similar question in the spatially uniform situation, see Eq. (13), we consider the dispersion at $t\gg\tau_{\eta}$,

[TABLE]

where we used Eq. (6) and consider $t\ll t_{in}(y)$ so the statistics of $w(t)$ does not change over the considered time interval. The difference from the spatially uniform case is that the last integral is no longer necessarily $\sum\lambda_{i}$, see Eq. (1). For $y$ which are not too far from the bulk so we have $\tau^{2}\int_{0}^{\infty}\langle\nabla^{2}p(0)\nabla^{2}p(t)\rangle_{c}dt^{\prime}\gtrsim|\partial_{y}v_{eff}(y)|$, see Eq. (1), we can use Eq. (12). Thus for these $y$ we can use for one trajectory $\sum\lambda_{i}(t)\approx\langle\sum\lambda_{i}(t)\rangle$ for $t\gtrsim\tau_{\eta}(y)/St^{2}(y)$ where $St(y)=\tau/\tau_{\eta}(y)$. In contrast, when $y$ is closer to the wall region where $\tau^{2}\int_{0}^{\infty}\langle\nabla^{2}p(0)\nabla^{2}p(t)\rangle_{c}dt^{\prime}\ll|\partial_{y}v_{eff}(y)|$ we have in that region,

[TABLE]

cf. Eq. (12). In this case convergence is faster.

V Preferential concentration and distinction between Kaplan-Yorke and correlation dimensions

The theory described in previous Sections provided the local rate of production of inhomogeneities of inertial particles in channel turbulence. This rate provides the growth of concentration of particles $n(t,\bm{x})$. Solving the continuity equation at $t>0$, we have

[TABLE]

where the RHS is $\sum\lambda_{i}$ up to the sign. This formula is implied by mass conservation in the particle’s frame, $n(t)V(t)=const$. We see from Eq. (33) that concentration grows at large times as $\exp\left[|\langle\sum\lambda_{i}\rangle|t\right]$. We consider the history of creation of fluctuation of concentration at scale $r\ll\eta$, see Falkovich et al. (2002); Fouxon et al. (2015). This starts from compression of the volume of particles whose initial size is the correlation length $\eta$ of $w$. Over the initial volume the concentration is effectively uniform Fouxon (2012). The smallest dimension of the compressed volume decreases with time $t$ as $\eta\exp[-|\lambda_{3}|t]$ where $\lambda_{3}$ is the third Lyapunov exponent Schmidt et al. (2016); Devenish et al. (2012). Since $\lambda_{3}$ is non-zero for tracers then, in the leading order in Stokes number, we can use $\lambda_{3}$ of the fluid particles. However for fluid particles volumes are conserved. Thus the smallest dimension decreases so that its product with the growing area, giving the volume, stays constant. This gives $|\lambda_{3}|=\lambda_{1}+\lambda_{2}$ where $\lambda_{1}+\lambda_{2}$ is the growth exponent of areas defined in Eq. (10). The fluctuation of concentration grows until the time $t_{*}=|\lambda_{3}|^{-1}\ln(\eta/r)$ when the smallest dimension becomes equal to $r$. Beyond this time there is no growth of correlated fluctuations of concentration Fouxon et al. (2015). We find that the factor of increase of concentration is $\exp[|\sum\lambda_{i}|t_{*}]=(\eta/r)^{D_{KY}}$ where,

[TABLE]

We observe that $D_{KY}$ has the structure of the reduced formula for the Kaplan-Yorke fractal codimension in the case of weak compressibility Kaplan and Yorke (1979); Fouxon (2012). We defined

[TABLE]

where $S(t,x,y,z)$ is an arbitrarily oriented infinitesimal area element located near $(x,y,z)$ initially. This average is independent of $t$ in the considered time interval. Further, the orientation of the surface, that can be defined by the normal, is irrelevant for the long-time limit despite the anisotropy of the statistics of turbulence. This is because orientation reaches (anisotropic) steady state quite fast Devenish et al. (2012) (this is not the case for sedimenting inertial particles where relaxation of orientation is long Fouxon et al. (2015)). There is a significant difference of convergence time for $\sum\lambda_{i}$ and $\lambda_{3}$ as remarked previously: we have $\lambda_{3}(t,y)\approx\langle\lambda_{3}(y)\rangle$ at $t\gg\tau_{\eta}(y)$, cf. spatially uniform case Balkovsky et al. (2001); Devenish et al. (2012).

In the case of spatially uniform turbulence $D_{KY}$ describes all the fractal codimensions Fouxon (2012). For instance the pair-correlation function of concentration scales as $(\eta/r)^{\Delta}$ where $\Delta=2D_{KY}$ is the correlation codimension Balkovsky et al. (2001); Fouxon (2012). However in the case of inhomogeneous turbulence this is no longer true. This can be seen in the simplest context by considering the exponential growth of moments of concentration which in spatially uniform case is determined by $\sum\lambda_{i}$ completely. We have,

[TABLE]

We can use the formula $\langle\exp[X]\rangle=\exp[\langle X\rangle+\langle X^{2}\rangle_{c}/2]$ for averaging of Gaussian random variable $X$ doing averaging of the last term. This can be proved using cumulant expansion and smallness of compressibility Fouxon (2012). We find at $t\gg\tau_{\eta}(y)$,

[TABLE]

In the spatially uniform case described in Section II $|\sum\lambda_{i}|=\int_{0}^{\infty}\langle w(0)w(t)\rangle_{c}dt$ so the growth exponents are $k(k+1)|\langle\sum\lambda_{i}\rangle|$. The exponent is zero at $k=-1$ because of the conservation of the number of particles Balkovsky et al. (2001). The exponents are determined by $\sum\lambda_{i}$ completely. In contrast, for inhomogeneous turbulence factors near $k$ and $k^{2}$ in the formula above become independent. The average concentration can increase or decrease locally without contradicting the global conservation of the number of particles as in turbophoresis. The growth exponents of the moments of concentration are no longer determined uniquely by $\sum\lambda_{i}$. The formulas derived above hold for concentration in the particle’s frame. Similar formulas can be written for the growth of moments of concentration at a fixed spatial point Fouxon (2011).

We consider statistics of particle distribution in space after transients. If we seed particles in the channel then, after transients that at scale $r$ have typical time-scale $t_{*}$, they distribute over a multifractal structure in space Balkovsky et al. (2001); Falkovich et al. (2002); Fouxon (2012). The statistics of the distribution can be obtained by averaging over the $x-z$ plane as in the previous Sections. It was demonstrated in Schmidt et al. (2016) that pair-correlation function of concentration of particles $n(t,\bm{x})$ factorizes in product of (possibly time-dependent) average local concentrations and geometrical factor of fractal increase of probability of two particles to be close. The obvious change of the formula for $w$ in Schmidt et al. (2016) gives,

[TABLE]

We see that the correlation codimension is not $2D_{KY}$. It scales proportionally with $St^{2}$ when $D_{KY}$ has both the term that scales linearly and the terms that scales quadratically with $St$. The reason why the $\partial_{y}v_{eff}(y)$ term in $\sum\lambda_{i}$ drops from $\Delta$ is that this term originates in the average velocity that affects equally the average $n$ and its fluctuations disappearing from the ratio $n/\langle n\rangle$. The form of $\langle n(t,\bm{x})\rangle$ in Eq. (39) is determined by initial and boundary conditions on the concentration and is problem-dependent in contrast with the power-law factor. Further the correlations do not depend on the direction of $\bm{r}$ despite anisotropy of the statistics of turbulence. This is restoration of isotropy that originates in independence of Lyapunov exponents on the initial orientations.

Finally we consider the coarse-grained concentration $n_{l}(\bm{x})$,

[TABLE]

which is mass $m_{l}(\bm{x})$ in small volume of radius $l\ll\eta$ divided by the volume. We can find $m_{l}(\bm{x})$ using the consideration of Fouxon (2012); Schmidt et al. (2016) by tracking the ball of the particles back in time to time $t=-t_{*}$ where $t_{*}=-|\lambda_{3}(y)|^{-1}\ln(\eta/l)$. Since there are no fluctuations of concentration at scale $\eta$ then the mass of the ball at that time is volume $(4\pi l^{3}/3)\exp[-\int_{-t_{*}}^{0}w(t^{\prime},\bm{q}(t^{\prime},\bm{x}))dt^{\prime}]$ times the average concentration $\langle n(-t_{*},\bm{q}(-t_{*},\bm{x}))\rangle$ see details in Fouxon (2012). Comparing the resulting formula for $n_{l}(\bm{x})$ with the formula for $\langle n(\bm{x},t)\rangle$ from Schmidt et al. (2016) we find

[TABLE]

This formula was provided (with a typo) in Schmidt et al. (2016) where the result of the averaging was presented,

[TABLE]

which for $k=2$ corresponds with the previously derived formula for the pair-correlation of concentration. Thus if we use the scaling of $n_{l}/\langle n(\bm{x},t)\rangle$ for defining fractal dimensions Fouxon (2012) then none of the dimensions is $D_{KY}$.

VI Numerical Simulation

We use the direct numerical simulation (DNS) of a turbulent channel flow provided by the JHTDB. A large variety of time and space dependent, Eulerian simulation results are stored on a cluster of databases, which is made accessible to the public. The functionality of this database systems and details on the simulations available, as well as confirmations of their validity, are described in Perlman et al. (2007); Li et al. (2008); Yu et al. (2012); Graham et al. (2016) and other references therein. All details on the DNS computation, specifically numerical schemes, discretization methods and further simulation details of the turbulent channel flow are extensively described in Graham et al. (2016).

The turbulent channel flow with a friction Reynolds number $Re_{\tau}\approx 1000$, considered in this work is a wall bounded flow with no-slip conditions at the top and bottom walls ($y/h=\pm 1$, where $h$ corresponds to half of the channel height) and periodic boundary conditions in the longitudinal and transverse directions. In this channel flow DNS, the streamwise direction $x$ and the transverse direction $z$ can be considered homogeneous, whereas the wall-normal direction $y$ serves as inhomogeneous direction for the purpose of generating turbophoretic drift of inertial particles. The domain spans over the three directions as follows: $L_{x}\times L_{y}\times L_{z}=8\pi h\times 2h\times 3\pi h$, where $h=1$ in dimensionless units. Quantities normalized by the friction velocity $u_{\tau}$, the viscous length $\nu/u_{\tau}$ ($\nu=$ viscosity) or the viscous time scale $\nu/u_{\tau}^{2}$ are presented with the superscript $+$. The wall of the channel is located at $y^{+}=0$, the center of the channel is at $y^{+}=1000$. An overview over the main simulation, flow and grid parameters is given in Table 1.

The time step, $\delta t$, at which the Eulerian flow data can be extracted from the database is $0.0065$ and the total available flow time is approximately $26$ (non-dimensional time units), which corresponds to approximately one flow through time. The total duration of the simulation is thus $t_{end}^{+}=2.6\times 10^{4}$.

We have used the Eulerian results of the channel flow DNS provided by the JHTDB to perform Lagrangian tracking of inertial particles in the channel flow. At total number of, $4\times 10^{6}$ point-particles are randomly seeded across the entire channel domain. The JHTDB allows to extract velocity, velocity gradients and the Hessian of pressure at any arbitrary particle position. We use Eq. (3) to determine the inertial particle velocity. We compute the second term on the RHS of Eq. (3) $(\partial_{t}\bm{u}+(\bm{u}\cdot\nabla)\bm{u})$ based on the material derivative $D\bm{u}/Dt$ along a tracer particle trajectory. This is done by applying a simple finite difference scheme using the tracer particle velocity at the tracer particle position of two consecutive time steps. The inertial particles are advected with time step $\delta t$ applying a second order Adams-Bashforth method for temporal integration.

A verification of the channel flow DNS, including a comparison with previous works, has been done in Graham et al. (2016). We show the Reynolds stresses $u^{\prime}_{i}u^{\prime}_{j}$ as well as the turbulent kinetic energy, normalized by the friction velocity $u_{\tau}$ versus the viscous wall distance $y^{+}$ in Fig. 1(a) and plot the wall-distance logarithmically in order to focus the visualization on the near-wall behavior of the flow. The dashed red line, indicating the Reynolds stresses based on the wall-normal velocity $u^{\prime}_{y}$, shows a maximum at around $y^{+}=100$ dropping quite steeply towards the wall. We highlight this term in Fig. 1, since its second derivative is responsible for the turbophoretic migration of the particles as shown previously in Eq. (4).

The ratio between the particle response time $\tau$ and the Kolmogorov time scale $\tau_{\eta}$, defines the particle Stokes number $St$. Due to the strong dependence of the turbulent kinetic energy on the wall-normal direction (red solid line in Fig. 1(a)), also the dissipation of the turbulent kinetic energy and thus the Kolmogorov microscales, depend on $y$. Fig. 1(b) shows $\tau_{\eta}$ versus the inhomogeneous spatial direction $y^{+}$. For the presented theory to be valid we choose a rather weak particle inertia by setting the Stokes number averaged over the whole channel to $\langle St\rangle=0.1$, based on the averaged Kolmogorov time $\langle\tau_{\eta}\rangle=0.2133$. This yields $\tau=0.0213$ and thus a viscous Stokes number of $St^{+}\approx 1$. The viscous Stokes number is defined in terms of the friction velocity $u_{\tau}$ as $St^{+}=\tau u_{\tau}^{2}/\nu$. Due to the variation of $\tau_{\eta}$ along wall-normal direction (Fig. 1) also the local $St$ depends strongly on $y$, reaching the largest values in the vicinity of the wall.

VII Results and Discussion

VII.1 Results

Before going into the analysis of cluster formation, we start with discussing the effect of turbophoresis. In a channel flow inertial particle migration due to turbophoresis changes the initially random distribution of particles by driving particles towards the wall. A snapshot of the distribution of tracer and inertial particles at $t^{+}=2.5\times 10^{4}$ in the $x-y$ plane near the wall (Fig. 2) shows the increased inertial particle concentration near the wall qualitatively. The plane chosen here is located at $z=1.5\pi$ and covers the full extension of the channel in streamwise direction but only the near wall region between $y^{+}=0$ and $y^{+}=15$. Particles are plotted as black points on top of the time-averaged $\partial_{y}^{2}\langle u_{y}^{2}\rangle$ field. It is visible that the tracer particles shown in Fig. 2(a) distribute randomly in space and appear uniformly distributed while inertial particles (Fig. 2(b)) accumulate near the wall.

The temporal evolution of wall-normal particle concentration profiles are depicted in Fig. 3, which compares the particle concentration $n$ at the initial time step $t^{+}=0$ (circles), as well as $t^{+}=10^{4}$ (triangle), $t^{+}=2\times 10^{4}$ (triangle upside down) and $t^{+}=2.5\times 10^{4}$ (squares), normalized by the initial, random particle concentration $n_{0}$. It is visible that with time more particles accumulate in the vicinity of the wall. In the region below $y^{+}=10$ a significant increase of the particle concentration is observable. In the vicinity of the wall, the ratio $n/n_{0}$ rises almost up to factor three in the considered time span, despite the relatively weak inertia of the particles. Above $y^{+}=10$ the particle concentration decreases below the initial concentration where a minimum of $n/n_{0}=0.95$ is reached at $y^{+}=20$. This leads to a specific interest of the degree of clustering in both, the region where the number of particles is high (below $y^{+}=10$) and the region where the number of particles is low (between $y^{+}=10-70$).

The degree of clustering is not only determined by the particle concentration itself but instead rather by a high probability of finding particles with very small inter-particle distance. As explained in section II, particles approach each other and thus form clusters if $\nabla\cdot\bm{v}<0$ along particle trajectories. This effect is quantified by the previously introduced sum of Lyapunov exponents, $\sum\lambda_{i}$. Negative values of $\sum\lambda_{i}$ correspond to a compression of infinitesimal volumes formed by particles, whereas positive values of $\sum\lambda_{i}$ correspond to diverging volumes. For a precise quantification of the degree of clustering arising from the combined action of inhomogeneous and homogeneous clustering effects, it is necessary to compute the finite-time Lyapunov exponents. We do this in the following by the computation of each individual term of the RHS in Eq. (1) or Eq. (28), respectively.

Since all these terms depend on time, it is important to look at the convergence time of the individual terms. Theoretically, convergence within a few Kolmogorov time scales of the first term ($\tau\partial_{y}\int_{0}^{t}\langle u_{y}(0)\nabla^{2}p(t^{\prime})\rangle_{c}dt^{\prime}$) as well as the second term ($-\tau^{2}\int_{0}^{t}\langle\nabla^{2}p(0)\nabla^{2}p(t^{\prime})\rangle_{c}dt^{\prime}$) on the RHS of Eq.(28) is expected. These two integrals use the cumulant terms as described in Eq. (25). As discussed in section IV the third integral, which is $\int_{0}^{t}\langle\nabla^{2}p\rangle$, does not converge. However, at the time where the other terms have converged the full term $\tau^{2}(\partial_{y}\langle u_{y}^{2}\rangle)\partial_{y}\int_{0}^{t}\langle\nabla^{2}p(t^{\prime})\rangle dt^{\prime}$ (including this integral) remains small and can therefore be excluded from the computation of $\sum\lambda_{i}$ as shown below. For the approximations to be valid the inhomogeneity time $t_{in}$ (Fig. 1)(b) has to be larger than the convergence time of these integrals.

The temporal evolution of the three integrals of Eq. (1) in Fig. 4. The color shading is darker for increasing $y^{+}$, i.e. light curves refer to regions near the wall and dark curves in the bulk of the channel, respectively. The integral of the first term on the RHS of Eq. (1) reaches relatively low values in the viscous sublayer and in the bulk, whereas it becomes more significant in the intermediate (log-layer) regions where convergence takes longer. In the bulk as well as in the viscous sublayer the curves converge to relatively low values. The second integral of Eq. (1) is shown in Fig. 4(b), where it is seen that all curves converge fast. The largest values are found in the regions between $y^{+}=6-25$, where turbulence intensity peaks (see Fig. 1(a)). Figure  4(c) shows $\int_{0}^{t}\langle\nabla^{2}p(t^{\prime})\rangle dt^{\prime}$ that as expected does not converge. Generally, the convergence time is smaller than $t_{in}$. However, as one can see from Fig. 4(a) in the buffer layer region the convergence of the first term is rather slow and an upper limit can be estimated at about $100\langle\tau_{\eta}\rangle$, which is the time we choose to evaluate the mean Lyapunov exponents.

Now we evaluate the individual terms on the RHS of Eq. (28) in Fig. 5(a). Circles show the wall-normal profile of the first term which is the sole term causing small-scale clustering of inertial particles in homogeneous turbulence.

This correlation is zero at the wall but its absolute value increases further away from the wall reaching the largest negative values in the range between $y^{+}=5-20$. Beyond $y^{+}=20$ the curve approaches small magnitudes as turbulence becomes more homogeneous but will not vanish.

The second term on the RHS of Eq. (28) ($\tau\nabla_{y}\int_{0}^{T}\langle u_{y}(0)\nabla^{2}p(t)\rangle$ - x symbols) is negative and contributes to clustering in a part of the buffer layer, whereas it is positive (and counteracting clustering) outside that region.

The third term of Eq. (28) ($-\tau\nabla_{y}\langle u_{y}^{2}\rangle$), shown as $+$ symbols, accounting for the turbophoretic effect solely is negative below $y^{+}=15$. This indicates a compression of infinitesimal volumes and contribution towards clustering (this location has been defined as $y*$ in section II). The curve changes sign contributing against clustering above $y^{+}=15$. Above $y^{+}=15$ the term becomes slightly positive before converging towards [math] at $y^{+}=100$.

The fourth term of Eq. (28) shown as filled black points in Fig. 4(a) is rather small at the time where the other terms have reached convergence as predicted and will be neglected in the remaining analysis below.

In Fig. 5(b) the terms of Eq. (29) are displayed. We divide the total $\sum\lambda_{i}$ (diamond symbols) in a homogeneous or bulk component (circles) and an inhomogeneous component (squares). Despite the linear dependence of the inhomogeneous component on $St$ the clustering due to the homogenous contribution generally exceeds the inhomogeneous component. Both components add up and reach the maximal negative $\sum\lambda_{i}\approx-1$ at $y^{+}=6$.

Since we aim to quantify the preferential concentration via $D_{KY}$ and the correlation codimension $\Delta$ according to Eq. (34) and (39) respectively, in dependence on the inhomogeneous flow direction, the third Lyapunov exponent $\lambda_{3}$ has to be computed. The calculation of $\lambda_{3}$ is performed via the finite-time Lyapunov exponents that provide a measure of the cumulative deformation of the particles Meneveau et al. (2016). This requires an estimation of the instantaneous deformation rate of the particle trajectory, which is done via the instantaneous Lyapunov exponents $\lambda^{\prime}_{ii}$. These instantaneous Lyapunov exponents can be computed by the alignment of the eigenvectors of the Cauchy-Green tensor $C_{ij}$ of a particle trajectory and the velocity gradient tensor. More precisely the instantaneous Lyapunov exponents at each particle location can be found, using

[TABLE]

where $P_{j}$ are the eigenvalues of the strain rate tensor and $\theta_{ij}$ is the angle between eigenvector $i$ of the Cauchy-Green tensor and eigenvector $j$ of the strain rate tensor. Averaging those instantaneous Lyapunov exponents along a Lagrangian path enables us to determine $\lambda_{3}$ depending on the wall distance $y$ of the channel. The corresponding result is presented in the inset of Fig. 6(a). The maximum $|\lambda_{3}|$ is reached at $y^{+}=10$, the curve drops quite fast in both directions. As apposed to spatially uniform turbulence the estimate that $\lambda_{3}\approx\tau_{\eta}^{-1}$ is not true along the entire channel. Below $y^{+}=100$ the product of $\tau_{\eta}\lambda_{3}$ drops significantly, whereas it is constant above $y^{+}=100$.

The clustering introduced by $(\langle\sum\lambda_{i}\rangle)_{bulk}$ is quantified using $(D_{KY})_{hom}$. The resulting values are shown as blue circles in Fig. 6(a). With $(D_{KY})_{hom}=0.1$ the maximum level is moderately high but interestingly it stays almost constant across a large region ranging from $y^{+}=3$ to $y^{+}=30$. The values obtained for $\Delta_{hom+inhom}$, combining homogeneous and inhomgeneous clustering, are relatively high in regions below $y^{+}=10$, see Fig. 6(a). At the wall $\Delta$ is [math] but starts increasing rapidly. The curve peaks around $y^{+}=4$ with $\Delta_{hom+inhom}\approx 0.3$. It then decreases slowly to almost [math] at $y^{+}=200$. In contrast to $(D_{KY})_{hom}$ the correlation codimension takes on very large values in a much narrower range from $y^{+}=2-10$.

One key difference between the clustering in homogeneous and inhomogeneous turbulence is that for weakly inertial particles there is a linear dependence on St (for the first term on the RHS of Eq. (28)). The ratio $\left|\frac{\langle\sum\lambda_{i}\rangle_{i}nhom}{\langle\sum\lambda_{i}\rangle_{hom}}\right|$ changes not only significantly throughout the channel height but also for different Stokes numbers. One would expect this term to dominate in inhomogeneous regions of the flow. In the case examined in this study so far, this is not the case. The homogeneous term is generally higher or at least of the same magnitude as the inhomogeneous term. Therefore, we want to extend this study to other $St$. In Fig. 6(a) we show the ratio between the homogeneous and the inhomogeneous contributions to clustering additionally for the case where $\langle St\rangle=1$ or $St^{+}=1$ (empty circles) and $St^{+}=0.1$ (filled squares) and $St^{+}=10$ (filled triangles). We find that the Stokes number has a strong impact on which term dominates the clustering. However, predicting the behavior is not trivial due to the variations of the different terms along the $y-$ direction. As mentioned before, for $St^{+}=1$ the ratio of $\left|\frac{\langle\sum\lambda_{i}\rangle_{i}nhom}{\langle\sum\lambda_{i}\rangle_{hom}}\right|$ stays below $1$ throughout the entire channel. Interestingly, for particles with small inertia the inhomogeneous term will dominate, as can be seen for the case of $St^{+}=0.1$ in Fig. 6(b). For particles with larger inertia the behavior becomes more complex. The inhomogeneous terms will dominate or be of the same order as the homogeneous terms right at the wall but become less important throughout the rest of the channel. The homogeneous terms dominate even more strongly than in the case of $St^{+}=1$ further away from the wall.

VII.2 Discussion

The results presented in the previous section indicate a strong dependence of particle concentration and clustering on the wall-normal direction. Turbophoresis in a turbulent channel flow drives particles towards the wall. After $t^{+}=2.5\times 10^{4}$ the initial concentration at the wall is exceeded by almost factor 3 within the considered time span. The turbophoretic particle migration is driven by $-\tau\nabla_{z}^{2}\langle u_{z}^{2}\rangle$ shown in Fig. 5(a). The change of sign in this term at $y^{+}\approx 15$ determines the wall-distance below which particles start accumulating and above which the particle concentration decreases.

We analyze the space-dependent rate of creation of inhomogeneous particle concentration to quantify inertial particle clustering in inhomogeneous turbulence. All three terms investigated according to Eq. (1) and Eq. (28) determining $\langle\sum\lambda_{i}\rangle$, depend differently on the wall distance. The term $\tau^{2}(\partial_{y}\langle u_{y}^{2}\rangle)\partial_{y}\int_{0}^{\infty}\langle\nabla^{2}p(t^{\prime})\rangle dt^{\prime}$ was small enough to be neglected for times within the convergence time of the other terms. Particle clustering for the case of $St^{+}=1$ is dominated by the homogeneous fractal clustering, represented by $-\tau^{2}\int_{0}^{t}\langle\nabla^{2}p(0)\nabla^{2}p(t^{\prime})\rangle_{c}dt^{\prime}$ even close to the wall where inhomogeneity is strongest. Its dependence on $y$ (Fig. 5) is similar to the one of the turbulent kinetic energy, shown in Fig. 1(a). In regions of strong turbulence, $\tau\nabla^{2}p$ takes on large values and causes stronger clustering. Despite the linear dependence on $\tau$ the inhomogeneous contribution $(\langle\sum\lambda_{i}\rangle)_{inhom}$ to the overall clustering degree is smaller than $(\langle\sum\lambda_{i}\rangle)_{bulk}$ in the case of $St^{+}1$. Though $(\langle\sum\lambda_{i}\rangle)_{inhom}$ behaves similarly to $(\langle\sum\lambda_{i}\rangle)_{bulk}$ its absolute value peaks already at $y^{+}=4-5$ and is below the peak of $(\langle\sum\lambda_{i}\rangle)_{bulk}$. As can be seen in Fig. 6 the significance of the inhomogeneous terms decreases with larger wall-distance. We use the correlation codimension of the multifractal formed by inertial particles to quantify the strength of particle clustering. The values found for the the correlation codimension ($\Delta$) for the case under investigation here are relatively large for particles with such small inertia. In the region between $y^{+}=2$ to $y^{+}=10$, $\Delta$ reaches up to $0.3$. In the region between $y^{+}=10-70$ — where the particle concentration drops below the initial concentration — the correlation codimension decreases rapidly. Whereas, the local homogeneous turbulence contribution to clustering $(D_{KY})_{hom}$ stays almost constant at $(D_{KY})_{hom}=0.1$ from $y^{+}=3-30$. Studies with particles of similar inertia find lower values of $\Delta$ (between $10^{-2}$ and $10^{-1}$) in homogeneous isotropic turbulence, e.g. Saw et al. (2012); Chun et al. (2005). This shows that turbulent inhomogeneity enhances the clustering degree particularly in the viscous sublayer and the onset of the buffer layer for $St^{+}=1$. The difference between $\Delta$ and $D_{KY}$ shows also that the inhomogeneous terms enhance clustering but also affect the regions where clustering occurs by the peak of $\Delta$ towards the wall. This can be explained by the maximal absolute values of $(\langle\sum\lambda_{i}\rangle)_{inhom}$ and $(\langle\sum\lambda_{i}\rangle)_{bulk}$ in Fig. 5(b). Therefore, a complete investigation, unifying the clustering effects of both mechanism is essential for a precise quantification of the preferential concentration in inhomogeneous turbulent flows.

Due to the linear dependence of $(\langle\sum\lambda_{i}\rangle)_{inhom}$ on $\tau$ and the quadratic dependence of $(\langle\sum\lambda_{i}\rangle)_{bulk}$ on $\tau$, the situation changes for different particle inertia. In Fig. 6(a) we saw that for particles with even smaller inertia $St^{+}=0.1$ the inhomogeneous terms dominate clustering throughout the entire channel. Instead for particles with large inertia ($St^{+}=10$) the inhomogeneous terms dominate near the wall but then become outweight by the homogeneous clustering term. The complementary effects of inhomogeneous and homogeneous clustering cause a different degree of clustering depending on inertia and the wall-normal direction.

VIII Conclusion

In this study the combined effects of turbophoresis and small-scale fractal clustering of weakly inertial particles have been analyzed theoretically as well as numerically. A novel theoretically approach allows to describe the preferential concentration resulting from homogeneous and inhomogeneous particle clustering in inhomogeneous turbulence. We determine a space-dependent rate ($\langle\sum\lambda_{i}(y)\rangle$) that create inhomogeneities of concentration of particles. We find that $\sum\lambda_{i}$ depends linearly on the particle Stokes number, as opposed to homogeneous turbulence where $\sum\lambda_{i}\propto St^{2}$. The theoretical predictions for the creation of particle inhomogeneities have been investigated by direct numerical simulations of a turbulent channel flow using the JHTDB. The results reveal a strong turbophoretic migration of particles, despite the relatively small inertia. We find that the clustering degree depends strongly on the wall-normal direction. The strongest effects of preferential concentration can be observed in the transition of the viscous sublayer to the buffer layer at $y^{+}=2-10$, where local homogeneous terms contribute to clustering and clustering due to the inhomogeneity of the flow is strong. The correlation codimension rises up to $0.3$ which is remarkably high for particles with relatively weak inertia of $St^{+}=1$ or $\langle St\rangle=0.1$. The values found for $\Delta$ in the near-wall region are larger compared to what particles with the same inertia would show in homogeneous turbulence. The contributions to clustering from homogeneous and inhomogeneous terms varies strongly with the wall-distance but as we show also with particle inertia. For particles with very small inertia the inhomogeneous terms outweigh the homogeneous contributions significantly. However, for particles with large inertia homogeneous turbulence will mainly determine clustering except in the vicinity of the wall, where the contribution of inhomogeneous turbulence might be stronger. The findings of this work allow for a precise quantification of inhomogeneous preferential concentration of weakly inertial particles in non-uniform turbulent flows. The case of a turbulent channel flow investigated here, serves as general example for all inhomogeneous turbulent flows. Thus, the presented results can be easily transferred to investigate preferential concentration of weakly inertial particles in other frequently occurring turbulent flows, e.g. pipe or free shear flows.

IX Acknowledgments

L.S. would like to thank Stephen S. Hamilton for his support regarding the work with the database provided by the Johns Hopkins University. Financial support from the Swiss National Science Foundation (SNSF) under Grant No. 144645 is gratefully acknowledged.