Turbophoresis in forced inhomogeneous turbulence

We show, by direct numerical simulations, that heavy inertial particles (characterized by Stokes number St) in inhomogeneously forced statistically stationary isothermal turbulent flows cluster at the minima of mean-square turbulent velocity. Two turbulent transport processes, turbophoresis and turbulent diffusion together determine the spatial distribution of the particles. If the turbulent diffusivity is assumed to scale with turbulent root-mean-square velocity, as is the case for homogeneous turbulence, the turbophoretic coefficient can be calculated. Indeed, for the above assumption, the non-dimensional product of the turbophoretic coefficient and the rms velocity is shown to increase with St for small St, reach a maxima for St ≈10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${\rm St} \approx 10$\end{document} and decrease as ∼ St -0.33\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \sim {\rm St}^{-0.33}$\end{document} for large St.


Introduction
The subject of this paper is turbulent flows of Newtonian fluids with small particles. We limit ourselves to spherical particles whose radii are significantly smaller than the viscous scale of the fluid. We further assume that the material density of each particle is significantly larger than that of the fluid and the number density of the particles is quite low. In this parameter regime the equations of motion of each individual particle are well described bẏ where X and V are respectively the position and velocity of the particle, τ is the characteristic relaxation time of the particle and u is the flow velocity, whose magnitude is denoted by u. As the number density of the particles is assumed to be small, the effects of the particles on the flow and particle-particle interactions are both ignored. The flow velocity u is obtained by solving the Navier-Stokes equation with proper boundary conditions. Particles that satisfy (1) are called Heavy Inertial Particles (HIP). This approximation is often used to model various natural and industrial flows. Two important examples are: (a) droplets in warm clouds (see, e.g., [1,2] for a review) and (b) dust grains in astrophysical accretion disks (see, e.g., [3], for a review). The study of heavy inertial particles in homogeneous and isotropic turbulent (HIT) flows is a topic of numerous research publications. In addition to the fluid Reynolds number, such a problem is described by a second dimensionless number, the Stokes number, defined by St ≡ τ /τ f , where τ f is a characteristic time scale of the flow. Recent analytical, numerical, and experimental works have shown that HIPs self-organize in fractal clusters in HIT, (see, e.g., [4][5][6][7][8] and references therein) and at small distances the probability distribution function (PDF) of their relative velocities shows power-law behavior [9,10].
The focus of this paper is the large-scale clustering of heavy inertial particles in inhomogeneous, but unstratified, turbulence. Numerical and experimental studies in wall-bounded flows, e.g., channel and pipe flows, have shown that heavy inertial particles distribute themselves preferentially near the wall -a phenomenon that has been called turbophoresis [11] by analogy with thermophoresis [12][13][14]. Thermophoresis -motion of Brownian particles governed by gradient in temperature [15]-is well understood within the framework of local thermodynamic equilibrium in statistical mechanics [12]. Recent works [13,14] have demonstrated that the flux of molecules in thermophoresis (J ) can be well described by (2) where D is the diffusion coefficient, N the concentration of the molecules, D T the thermophoretic coefficient, and the symbol · th denotes averaging over a state in local thermal equilibrium at temperature T .
To understand turbophoresis, let us consider a statistically stationary inhomogeneous turbulent flow. For simplicity, we consider the turbulence to be inhomogeneous in only one coordinate direction, which we denote by z. This is often the case in most practical applications, e.g., in turbulent channel or pipe flows. We start with the mean flux (J ) of inertial particles in space, where the symbol (•) denotes Reynolds averaging, and N denotes the number density of the particles. We use the standard Reynolds decomposition, where each quantity of interest is decomposed into mean and fluctuating parts. After Reynolds averaging we obtain The second term on the right-hand side (RHS) of this equation is usually modelled through closure by nv = −κ∇N , where κ is the turbulent diffusivity. To calculate the first term in the RHS of (5) we use (1) to write (following ref. [16]) which is an approximation valid for small τ . On Reynolds averaging, we obtain In general, we are interested in the component of the flux J along the direction of inhomogeneity, which in wall bounded flows, e.g., pipe or channel flows, is the wall-normal direction. Projecting (7) along the wall-normal direction we obtain because u z = 0. This demonstrates that although the mean flow velocity in the wall normal direction is zero, the mean particle velocity is not. Putting together, the flux of inertial particles along the direction of inhomogeneity is given by The second term on the RHS is the usual Fickian (turbulent) diffusion, the first one is called the turbophoretic term. There are several ways to obtain (9), here we have followed refs. [17][18][19]. Clearly, (9) is an approximation that holds for small St (small τ ). How can this expression be generalized for arbitrary values of the Stokes number?
In this paper we show by direct numerical simulations (DNS) that the large-scale clustering of heavy inertial particles in inhomogeneous turbulence can be described by a simple generalization of (9): Here κ t is the turbophoretic coefficient, κ is turbulent diffusivity, and K ≡ (u 2 /2). The expression of molecular flux in thermophoresis, (2), can be considered as an inspiration for (10). So far, neither numerical nor experimental studies have measured κ t and κ simultaneously and it is not clear how this can be done. Here, we perform careful numerical experiments to measure their ratio: which we call the turbophoretic ratio. We find that our data is well described by Tu = σ t /(2 √ K). We nondimensionalize σ t by multiplying it with u rms ; the non-dimensionalized σ t is a non-monotonic function of the Stokes number -it increases linearly with St for small St, reaches a maximum value around St ≈ 10 and decreases with St for large St 1 . This is the first numerical measurement of the turbophoretic ratio and its dependence on the Stokes number. Table 1. Parameters for our DNS runs: the runs are grouped by the function F (z) which determines the spatial variation of the mean-square turbulent velocity K, in particular; A) F (z) = exp[−(z/H) 2 ] with H = 1, B) F (z) = sin 2 (z), and C) F (z) = sin 2 (2z). The numbers of Eulerian grid points are N 3 , the number of HIPs is given by Np.

Run
N

Model
So far, most numerical [20][21][22][23] and experimental [24] observations of large-scale clustering and deposition [25,26] of HIP in inhomogeneous flows are in wall-bounded flows (see, e.g., [19] for a review). In such flows it has not been possible to disentangle the individual effects of inhomogeneities of mean-square turbulent velocity, large-scale shear, and physical boundaries. For a deeper understanding of turbophoresis, we choose to do DNS of the Navier-Stokes equation under isothermal conditions, with an inhomogeneous external force. Here ∂ k denotes partial derivatives with respect to coordinate k, D t ≡ ∂ t + u k ∂ k is the advective derivative, u (whose k-th component is u k ), p, and ρ are respectively the velocity, pressure, and density of the flow, μ is the dynamic viscosity, and S kj ≡ ∂ k u j + ∂ j u k − δ jk (2/3)∂ k u k is the viscous stress. The simulations are performed in a three-dimensional periodic box with sides L x = L y = L z = 2π. In addition we use the ideal gas equation of state with a constant speed of sound c s = 1.
The flow attains a statistically stationary state where the average energy dissipation by viscous forces is balanced by the average energy injection by the external force f which is random, Gaussian, white-in-time, with zero-mean, concentrated on a shell of wavenumber with radius k f in Fourier space. The amplitude of the external force is chosen such that the Mach number is always less than 0.1, i.e., the flow is weakly compressible. Crucially, the flow is made inhomogeneous (but statistically isotropic) by the function F (z) which is a function of one coordinate direction, z 2 . We choose the function F (z) to be a slowly varying function of z, i.e., the characteristic wavenumber of F (z), q F k f . In particular, we examine three cases: In the context of the present problem, length scales larger than 1/k f are the large scales.
We use the pencil-code [27], which uses a sixth-order finite-difference scheme for space derivatives and a third-order Williamson-Runge-Kutta [28] scheme for time derivatives. The same setup, with a homogeneous external force, has been used in studies of scaling and intermittency in fluid and magnetohydrodynamic turbulence [29][30][31]. We introduce the particles into the DNS after the flow has reached a statistically stationary state. Then we simultaneously solve the equations of the flow, (12), and the HIPs (1). To solve for the HIPs in the flow we have to interpolate the flow velocity to typically off-grid positions of the HIPs. We use a tri-linear method for interpolation.
We define the Reynolds number by Re ≡ u rms /(νk f ), where u rms is the root-mean-square velocity of the flow averaged over the whole domain and the kinematic viscosity ν = μ/ρ. As we are interested in large-scale clustering, we use the large eddy turnover time, τ f ≡ 1/(u rms k f ), which is the characteristic time scale of the flow at scale 1/k f , to define the Stokes number. In our simulations we use three different values of L x k f /2π = 5, 15, and 30. Here we implement Reynolds averaging by averaging over the x and y coordinate directions. On top of the Reynolds averaging, we also perform time averaging over several thousand large eddy turnover times. The various relevant parameters of our simulations are shown in table 1.

Results
Once our simulations reach a statistically stationary state we plot, in fig. 1, the mean-square turbulent velocity K(z) ≡ (1/2)u 2 and the mean number density of the particles, N , for the three different F (z) and for different values of the Stokes number. Here the Eulerian number density of the particles N (x) ≡ δ 3 (x − X), where δ 3 (·) is the three-dimensional Dirac delta function. Clearly, the particles cluster away from the maxima of K, i.e., turbophoresis is observed.
We now want to quantify turbophoresis using (10). At a stationary state, the mean flux, J , must be zero 3 . Furthermore, as the turbulence is inhomogeneous in only one direction, z, (10) reduces to the following one-dimensional problem: Clearly, it is impossible to measure the turbulent transport coefficients κ t and κ individually. Instead, only their ratio, which we call the turbophoretic ratio, Tu, can be measured. In addition, to integrate (13) and compare with our data we need to make conjectures about the spatial dependence of Tu. The naive one would be to assume that Tu does not depend on space; which would imply N = N 0 exp(−TuK).
In the top panel of fig. 2 we show that this conjecture does not fit our data well even in the best of cases. Let us note now, that although we do not know exactly what κ would be in inhomogeneous turbulence, the turbulent diffusion has been a well studied problem in homogeneous and isotropic turbulence, in which case K is not a function of space and one typically obtains κ ∼ √ K. In the present problem K is a slowly varying function of space, hence we conjecture that κ(z) ∼ K(z) is also a slowly varying function of space, i.e., Tu ≡ σ t /(2 √ K) with σ t a constant in space. With this substitution, (13) can be integrated to obtain which fits our numerical data quite well as we show for two different cases in the middle and the bottom panel of fig. 2. If κ ∼ √ K were true even in inhomogeneous turbulence we would have obtained σ t ∼ κ t , i.e., the coefficient σ t is proportional to the turbophoretic coefficient κ t . To summarize, two results lie at the heart of our paper: A) We quantify large-scale clustering in inhomogeneous turbulent flows through two turbulent transport coefficients, turbulent diffusion and turbophoresis, whose ratio we call the turbophoretic ratio, Tu. The turbophoretic ratio, Tu is a function of space such that the function σ t ≡ (2 √ K)Tu is a constant. B) We find that the σ t (non-dimensionalized by multiplying with u rms ) is a non-monotonic function of the Stokes number with maximum clustering appearing at St ≈ 10.

Discussions
Large-scale clustering of inertial particles in inhomogeneous turbulent flows is commonly observed in geophysical and engineering problems. But in trying to understand such clustering we encounter a difficult problem due to an inhomogeneous distribution of mean-square turbulent velocity, presence of large-scale shear, anisotropy, gravity and boundary layers. To make progress, it is essential that we study each effect in isolation. To this end, in this paper, we have proposed a model in which large-scale clustering is due solely to the inhomogeneous distribution of meansquare turbulent velocity. We find that the minimum of the number density of the particles appears at the maximum of the mean-square turbulent velocity. Furthermore, we find that this clustering is determined by two processes: turbophoresis and turbulent diffusion. We numerically measure the ratio of these two transport coefficients, which we call the turbophoretic ratio, Tu, which is given by Tu = σ t /(2 √ K) with σ t as a universal function of St, increasing with St for small St reaching a peak near St ≈ 10 and then decreasing with St.
Traditionally such problems are studied by first constructing exact, but unclosed, relations between various moments of the probability distribution functions of flow velocity and density and velocity of the particle phase and then imposing arbitrary closures to obtain expressions for transport coefficients, see, e.g., refs. [32][33][34]. In this paper we have taken a different approach. We have assumed that the large scale transport is determined by two undertermined turbulent transport coefficients and then numerically calculated the ratio of the transport coefficients as a function of the Stokes number.
A special mention must be made of ref. [35] where the same clustering phenomenon is studied in turbulent Kolmogorov flows. Although the authors do not interpret their results as a balance between turbophoretic and turbulent diffusive fluxes as we do, they do observe that clustering increases for small St but this trend reverses smoothly at higher values of St.
In the introduction we have provided an approximate analytical theory for turbophoresis which gives a linear dependence of the turbophoretic coefficient on St for small St. To the best of our knowledge, this result was first obtained in refs. [36,37], using the Lagrangian History Direct Interaction Approximation, (see also [38] for a review). Our data does reproduce this trend for very small St although for those values of St the clustering is quite weak and the errorbars on σ t are large.
An alternative way to analytically understand turbophoresis would be to replace the fluid velocity, u, in eq. (1) by a Gaussian white noise whose strength is the function F (z) of one coordinate, z. Treating eq. (1) as a Langevin problem it is then possible to write the corresponding Fokker-Planck equation which is an equation for the probability distribution function p(z, v). The particle number density can be obtained by N (z) = p(z, v)dv. For a general F (z) a stationary solution of this Fokker-Planck equation is not known. If the strength of the noise F (z) = c is taken to be a constant, then it is possible to obtain a flux-free stationary solution of the Fokker-Planck equation to obtain p(z, v) ∼ exp(−τ v 2 /c 2 ). If F (z) is assumed to be a slowly varying function of space then using the method of multiple scales it is possible to calculate a flux that has both a diffusive and a turbophoretic part (see the supplemental material of ref. [39]). But in such a model the turbophoretic ratio is no longer a function of St but a constant. To capture the non-monotonic dependence of Tu on St, as we observe, it may be necessary to use a colored noise instead of white-in-time noise. In ref. [40] it has been shown that for F (z) = z 2 it is possible to find a non-trivial solution of the Fokker-Planck equation, which leads to a localization-delocalization transition. In other words no large-scale clustering should be observed above a critical value of the Stokes number. This is not found in our numerical experiments.
The large-scale clustering observed in turbulent pipe flows [41], which is similar to what has been observed in turbulent channel flows in ref. [23], can be analyzed in a similar manner to extract the turbophoretic ratio whose dependence on the Stokes number is similar to what has been observed by us but for the data very close to the wall and near the center of the pipe. We believe that to understand such departure the effects of the boundary layer and anisotropy have to be taken into account.
In our simulations we have assumed that the motion of the particles is well described by the HIP approximation, (1). For small St the validity of this approximation has been critically examined, by comparing DNS and experimental data for homogeneous and isotropic turbulence, in ref. [43]. To the best of our knowledge no such study exists for large St and for inhomogeneous turbulence. It may turn out that the large St behavior of σ t shown in fig. 3 may not be easily realisable in experiments. Nevertheless, the non-monotonic behavior of σ t as a function of St with a peak at around St ≈ 10 should be experimentally realisable. Note further that we have assumed that the particles are passive, which implies that our results apply to experiments where the mass loading parameter is small. Clustering of particles is also found in simulations where the feedback from the particles to the flow is included, see, e.g., refs. [33,34]. In those cases, gravity may also play an important role, i.e., horizontal and vertical pipes may show different clustering phenomena, see e.g., refs. [44][45][46][47]. A description of such clustering through turbophoresis and turbulent diffusion can also be done, but this is outside the scope of this paper.
In this work, we have considered an isothermal flow which is relevant to most engineering applications. In geophysical problems, it is common to have fully developed turbulent flows with an inhomogeneous distribution of mean temperature. In such cases, the HIPs are found to accumulate in the vicinity of the minima of mean temperature. This large-scale clustering phenomenon, which can be described by turbulent thermal diffusion, was predicted in an analytical study [18,[48][49][50][51][52] and detected in laboratory experiments [53][54][55], DNS [56] and atmospheric observations [57].