Hostname: page-component-7b7567c8c-7pwfn Total loading time: 0 Render date: 2026-11-14T21:57:23.158Z Has data issue: false hasContentIssue false
  • English
  • Français

Clustering of rapidly settling, low-inertia particle pairs in isotropic turbulence. Part 2. Comparison of theory and DNS

Published online by Cambridge University Press:  22 May 2019

Sarma L. RaniOpen the ORCID record for Sarma L. Rani[Opens in a new window],
Rohit DhariwalOpen the ORCID record for Rohit Dhariwal[Opens in a new window] and
Donald L. KochOpen the ORCID record for Donald L. Koch[Opens in a new window]
Sarma L. Rani*
Affiliation:
Department of Mechanical and Aerospace Engineering, University of Alabama in Huntsville, Huntsville, AL 35899, USA
Rohit Dhariwal
Affiliation:
Department of Civil and Environmental Engineering, Duke University, Durham, NC 27708, USA
Donald L. Koch
Affiliation:
Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, NY 14853, USA
*
†Email address for correspondence: sarma.rani@uah.edu
Rights & Permissions [Opens in a new window]

Abstract

Part 1 (Rani et al. J. Fluid Mech., vol. 871, 2019, pp. 450–476) of this study presented a stochastic theory for the clustering of monodisperse, rapidly settling, low-Stokes-number particle pairs in homogeneous isotropic turbulence. The theory involved the development of closure approximations for the drift and diffusion fluxes in the probability density function (p.d.f.) equation for the pair relative positions $\boldsymbol{r}$. In this part 2 paper, the theory is quantitatively analysed by comparing its predictions of particle clustering with data from direct numerical simulations (DNS) of isotropic turbulence containing particles settling under gravity. The simulations were performed at a Taylor micro-scale Reynolds number $Re_{\unicode[STIX]{x1D706}}=77.76$ for three Froude numbers $Fr=\infty ,0.052,0.006$, where $Fr$ is the ratio of the Kolmogorov scale of acceleration and the magnitude of gravitational acceleration. Thus, $Fr=\infty$ corresponds to zero gravity, and $Fr=0.006$ to the highest magnitude of gravity among the three DNS cases. For each $Fr$, particles of Stokes numbers in the range $0.01\leqslant St_{\unicode[STIX]{x1D702}}\leqslant 0.2$ were tracked in the DNS, and particle clustering quantified both as a function of separation and the spherical polar angle. We compared the DNS and theory values for the exponent $\unicode[STIX]{x1D6FD}$ characterizing the power-law dependence of clustering on separation. The $\unicode[STIX]{x1D6FD}$ from the $Fr=0.006$ DNS case are in reasonable agreement with the theoretical predictions obtained using the second drift closure (referred to as DF2). To quantify the anisotropy in clustering, we calculated the leading–order coefficient in the spherical harmonics expansion of the p.d.f. of pair relative positions. The coefficients predicted by the theory (DF2) again show reasonable agreement with those calculated from the DNS clustering data for $Fr=0.006$. However, we note that in spite of the high magnitude of gravity, the clustering is only marginally anisotropic both in DNS and theory. The theory predicts that the spherical harmonic coefficient scales with $\unicode[STIX]{x1D6FD}(=\unicode[STIX]{x1D6FD}_{2}St_{\unicode[STIX]{x1D702}}^{2})$, where $\unicode[STIX]{x1D6FD}_{2}$ is the ratio of the drift and diffusion flux coefficients. Since the drift flux, and thereby $\unicode[STIX]{x1D6FD}_{2}$, is seen to decrease with gravity for $St_{\unicode[STIX]{x1D702}}<1$, the anisotropy is also correspondingly diminished.

JFM classification

Turbulent Flows: Isotropic turbulence

Information

Type
JFM Papers
Information
Journal of Fluid Mechanics , Volume 871 , 25 July 2019, pp. 477 - 488
Copyright
© 2019 Cambridge University Press 

1 Introduction

This paper presents part 2 of the current work on a stochastic theory for the clustering of monodisperse, low-inertia particle pairs that are settling rapidly in homogeneous isotropic turbulence. The theory is developed in the limits of $Fr\ll St_{\unicode[STIX]{x1D702}}\ll 1$ , where $Fr$ is the Froude number and the Stokes number $St_{\unicode[STIX]{x1D702}}$ is the ratio of the particle response time to the Kolmogorov time scale. The part 1 paper (Rani, Gupta & Koch Reference Rani, Gupta and Koch2019) presented the derivation of: (1) closure approximations for the drift and diffusion fluxes in the probability density function (p.d.f.) equation for pair relative positions, and (2) the analytical solution to the p.d.f. $\langle P\rangle (r,\unicode[STIX]{x1D703})$ , where $r$ is the separation, and $\unicode[STIX]{x1D703}$ is the spherical polar angle. Part 2 focuses on the quantitative analysis of the theory by comparing its predictions of particle clustering with the results from DNS of isotropic turbulence containing particles settling under gravity.

Gravitational settling modifies particle dynamics in important ways. One of the main effects of gravity is that it introduces anisotropy into the particle sampling of the underlying turbulence, and thereby in the spatial clustering of particles. Gravity also alters the correlation times of fluid velocity gradients along particle trajectories. The altered time scales, in turn, modulate the path-history effects that play a key role in determining particle clustering (part 1 presents a more detailed discussion of the path-history effects and their role in clustering). The theory developed in this study incorporates the effects on clustering of settling-induced anisotropy, as well as of settling-modulated flow time scales.

Recently, Ireland, Bragg & Collins (Reference Ireland, Bragg and Collins2016) performed a detailed DNS investigation of the effects of gravity on the dynamics of single particles, as well as particle pairs. In their study, the Froude number $Fr=0.052$ , which is representative of fluid accelerations in cumulus clouds. They considered a wide range of Taylor micro-scale Reynolds numbers ( $88\leqslant Re_{\unicode[STIX]{x1D706}}\leqslant 597$ ) and Stokes numbers ( $0\leqslant St_{\unicode[STIX]{x1D702}}\leqslant 56.2$ ). For $St_{\unicode[STIX]{x1D702}}<1$ , they showed that the principal effect of gravity on particle clustering was to decrease the inward (radial) drift, thereby reducing the radial distribution function (r.d.f.). They also found that gravity mitigates the preferential concentration mechanism by reducing the interaction time between particles and the underlying turbulence. Specifically, gravity reduces the Lagrangian time scales of strain rate and rotation rate along particle trajectories. As shown in Chun et al. (Reference Chun, Koch, Rani, Ahluwalia and Collins2005), the drift flux is proportional to the time integral of the two-time correlation of $[\unicode[STIX]{x1D61A}^{2}(t)-\unicode[STIX]{x1D619}^{2}(t)]$ along the trajectory of the primary particle, where $\unicode[STIX]{x1D61A}^{2}$ and $\unicode[STIX]{x1D619}^{2}$ are the second invariants of the strain-rate and rotation-rate tensors, respectively. Ireland et al. (Reference Ireland, Bragg and Collins2016) also quantified the anisotropy in particle clustering due to gravity through the use of spherical harmonic functions to represent the r.d.f. dependence on the polar angle $\unicode[STIX]{x1D703}$ .

Bec, Homann & Ray (Reference Bec, Homann and Ray2014) performed a DNS study of the effects of settling on inertial particle clustering in isotropic turbulence. They considered both low- and high-Stokes-number particles ( $St_{\unicode[STIX]{x1D702}}\lesssim 1$ and $St_{\unicode[STIX]{x1D702}}>1$ , respectively) that are settling under low- and high-gravity conditions ( $Fr>1$ and $Fr<1$ , respectively). Bec et al. (Reference Bec, Homann and Ray2014) observed that when $Fr\ll 1$ , the clustering of $St_{\unicode[STIX]{x1D702}}\lesssim 1$ particles decreased and became anisotropic in that particles formed streaks along the vertical (or gravity) direction. The opposite effect is seen for $St_{\unicode[STIX]{x1D702}}>1$ particles, whose clustering increased when compared to the zero-gravity case. The opposing effects of gravity on the clustering of low- and high-Stokes-number particles was also seen by Ireland et al. (Reference Ireland, Bragg and Collins2016). Other DNS studies of settling particles (Ayala et al. Reference Ayala, Rosa, Wang and Grabowski2008; Onishi, Takahashi & Komori Reference Onishi, Takahashi and Komori2009; Woittiez, Jonker & Portela Reference Woittiez, Jonker and Portela2009; Parishani et al. Reference Parishani, Ayala, Rosa, Wang and Grabowski2015) focused primarily on the collision rates of droplets in isotropic turbulence.

In this part 2 paper, we present a comparison of the theory-predicted particle clustering with the corresponding DNS data. For Stokes numbers $St_{\unicode[STIX]{x1D702}}<1$ , we compare the exponent $\unicode[STIX]{x1D6FD}$ characterizing the power-law dependence of the p.d.f. $\langle P\rangle (r,\unicode[STIX]{x1D703})$ on separation $r$ . We also compare the degree of anisotropy of clustering obtained from the theory with that from DNS. In the case of DNS, particle clustering is quantified at $Re_{\unicode[STIX]{x1D706}}=77.76$ for three Froude numbers $Fr=\infty ,0.052,0.006$ (in the order of increasing gravity), and six Stokes numbers in the range $0.01\leqslant St_{\unicode[STIX]{x1D702}}\leqslant 0.2$ . While $Fr=0.052$ was chosen from the Ireland et al. (Reference Ireland, Bragg and Collins2016) study, the choice of $Fr=0.006$ had a twofold motivation. First, $Fr=0.006$ is smaller than the Stokes numbers considered, and thereby is in accordance with the parametric regime in which the theory is applicable. Second, too small an $Fr$ would have meant that the domain size along the settling direction would have to be significantly increased to account for the shorter particle-settling times through the domain. A larger domain size that satisfies the spatial resolution and Courant number requirements of DNS would lead to substantially higher computational costs.

For the highest gravity case, $Fr=0.006$ , the particle settling time through the periodic box of length $2\unicode[STIX]{x03C0}$ is close to the integral time scale (we will elaborate on this subsequently). As a result, numerical errors may occur since particles that have exited the domain and been reintroduced into it may again encounter the same correlated eddies that they have already seen previously on their way down the box (Woittiez et al. Reference Woittiez, Jonker and Portela2009; Ireland et al. Reference Ireland, Bragg and Collins2016). To eliminate this numerical artifact, we performed DNS with a bigger domain size of $4\unicode[STIX]{x03C0}$ along the vertical direction for the $Fr=0.006$ case.

The organization of the paper is as follows. Section 2 presents the computational details of the DNS runs, as well as the particle evolution algorithm. Quantification of the two-time correlations of dissipation rate and enstrophy is discussed in § 3. The development of a model energy spectrum that closely matches the DNS energy spectrum is presented in § 4. The model spectrum then serves as an input to calculating the drift and diffusion fluxes. In § 5, we present the comparison of theory predictions of particle clustering with the DNS data. Section 6 summarizes the key findings of this part 2 paper.

2 Computational method

2.1 Fluid phase

Direct numerical simulations of forced isotropic turbulence were performed using a pseudo-spectral method involving the discrete Fourier expansions of flow variables. The simulation domain, consisting of a cube of length $2\unicode[STIX]{x03C0}$ , is discretized into $N^{3}$ grid points, with periodic boundary conditions along all three Cartesian directions.

The governing equations for the flow are the Navier–Stokes equations in rotational form and the continuity equation (Brucker et al. Reference Brucker, Isaza, Vaithianathan and Collins2007; Ireland et al. Reference Ireland, Vaithianathan, Sukheswalla, Ray and Collins2013)

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D74E}\times \boldsymbol{u}=-\unicode[STIX]{x1D735}(p/\unicode[STIX]{x1D70C}_{f}+\boldsymbol{u}^{2}/2)+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D735}\times \boldsymbol{u}$ is the vorticity, $\unicode[STIX]{x1D70C}_{f}$ is the fluid density, and $p$ is the pressure.

Transforming (2.1) and (2.2) into Fourier space and eliminating pressure using the spectral form of continuity yields

(2.3) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}\right)\widehat{\boldsymbol{u}}=-\left(\unicode[STIX]{x1D644}-\frac{\unicode[STIX]{x1D73F}\unicode[STIX]{x1D73F}}{\unicode[STIX]{x1D705}^{2}}\right)\boldsymbol{\cdot }\widehat{\unicode[STIX]{x1D74E}\times \boldsymbol{u}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D73F}$ is the wavenumber vector, and $\unicode[STIX]{x1D705}^{2}=\unicode[STIX]{x1D73F}\boldsymbol{\cdot }\unicode[STIX]{x1D73F}$ . Direct evaluation of the convolution $\widehat{\unicode[STIX]{x1D74E}\times \boldsymbol{u}}$ is extremely computationally intensive. Hence, a pseudo-spectral approach is adopted wherein $\unicode[STIX]{x1D74E}\times \boldsymbol{u}$ is first computed in physical space, and then transformed into the spectral space.

Since the time-derivative and viscous stress terms on the left-hand side of (2.3) are linear in $\widehat{\boldsymbol{u}}$ , one may evolve these terms in time exactly by multiplying (2.3) with the integrating factor, $\exp (\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}t)$ . This yields the following equation (in index notation):

(2.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}[\exp (\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}t)\widehat{u}_{i}]=\text{RHS}_{i}\times \exp (\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}t),\end{eqnarray}$$

where $\text{RHS}_{i}=(-\unicode[STIX]{x1D6FF}_{im}+\unicode[STIX]{x1D73F}_{i}\unicode[STIX]{x1D73F}_{m}/\unicode[STIX]{x1D73F}^{2})\unicode[STIX]{x1D716}_{mjk}{\mathcal{F}}\{\unicode[STIX]{x1D714}_{j}u_{k}\}$ represents the right-hand side of (2.3), and $\unicode[STIX]{x1D716}_{mjk}{\mathcal{F}}\{\unicode[STIX]{x1D714}_{j}u_{k}\}$ represents the convolution $\widehat{\unicode[STIX]{x1D74E}\times \boldsymbol{u}}$ , and $\unicode[STIX]{x1D716}_{mjk}$ is the Levi-Civita tensor.

Equation (2.4) is then discretized in time using the second-order Runge–Kutta (RK2) method giving

(2.5) $$\begin{eqnarray}\widehat{u_{i}}^{n+1}=\widehat{u_{i}}^{n}\exp (-\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}t)+\{\text{RHS}_{i}^{n}\exp (-\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}t)+\text{RHS}_{i}^{n+1}\},\end{eqnarray}$$

where $n$ is the the previous time-step level and $h$ is the time-step size. To prevent convective instabilities, time-step size $h$ is chosen such that the Courant–Friedrichs–Lewy (CFL) number ${\leqslant}0.5$ . The pseudospectral algorithm introduces aliasing errors which are removed by zeroing the fluid velocities in spectral space for wavenumbers satisfying $\unicode[STIX]{x1D705}\geqslant \unicode[STIX]{x1D705}_{max}$ , where $\unicode[STIX]{x1D705}$ is the wavenumber magnitude, $\unicode[STIX]{x1D705}_{max}=\sqrt{2}N/3$ , and $N$ is the number of grid points along each dimension.

To achieve statistically stationary turbulence, we employ the deterministic forcing method developed by Witkowska, Juvé & Brasseur (Reference Witkowska, Juvé and Brasseur1997), wherein the turbulent kinetic energy dissipated during a time step is added back to the flow at the low wavenumbers. It may be noted that in this method, there is no explicit forcing term $\boldsymbol{f}$ added to the Navier–Stokes equations. Instead, one scales the velocity components in the wavenumber band $[\unicode[STIX]{x1D705}_{min},\unicode[STIX]{x1D705}_{max}]$ by a factor such that the energy dissipated during a given time step is resupplied, as follows:

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\boldsymbol{u}}(\unicode[STIX]{x1D73F},t+\unicode[STIX]{x0394}t)=\widehat{\boldsymbol{u}}(\unicode[STIX]{x1D73F},t+\unicode[STIX]{x0394}t)\sqrt{1+\frac{\unicode[STIX]{x0394}E_{diss}(\unicode[STIX]{x0394}t)}{\displaystyle \int _{\unicode[STIX]{x1D705}_{min}}^{\unicode[STIX]{x1D705}_{min}}E(\unicode[STIX]{x1D705},t+\unicode[STIX]{x0394}t)\,\text{d}\unicode[STIX]{x1D705}}}\quad \forall \unicode[STIX]{x1D705}\in [\unicode[STIX]{x1D705}_{min},\unicode[STIX]{x1D705}_{max}], & \displaystyle\end{eqnarray}$$
(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}E_{diss}(\unicode[STIX]{x0394}t)=\int _{\unicode[STIX]{x1D705}_{min}}^{\unicode[STIX]{x1D705}_{min}}E(\unicode[STIX]{x1D705},t)\,\text{d}\unicode[STIX]{x1D705}-\int _{\unicode[STIX]{x1D705}_{min}}^{\unicode[STIX]{x1D705}_{min}}E(\unicode[STIX]{x1D705},t+\unicode[STIX]{x0394}t)\,\text{d}\unicode[STIX]{x1D705}, & \displaystyle\end{eqnarray}$$

where $\widehat{\boldsymbol{u}}(\unicode[STIX]{x1D73F},t)$ is the spectral velocity, $\unicode[STIX]{x0394}E_{diss}$ is the total energy dissipated during $\unicode[STIX]{x0394}t$ , and $E(\unicode[STIX]{x1D705},t+\unicode[STIX]{x0394}t)$ is the spectral turbulent kinetic energy in a wavenumber shell with magnitude $\unicode[STIX]{x1D705}$ at time $t+\unicode[STIX]{x0394}t$ . In the current study, the velocity components in the range $\unicode[STIX]{x1D705}\in (0,\sqrt{2}]$ are forced using (2.6).

2.2 Particle phase

The governing equations of motion for a heavy spherical particle, whose diameter is much smaller than the Kolmogorov length scale, may be written as

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{x}_{p}}{\text{d}t}=\boldsymbol{v}_{p}, & \displaystyle\end{eqnarray}$$
(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{v}_{p}}{\text{d}t}=\frac{\boldsymbol{u}(\boldsymbol{x}_{p},t)-\boldsymbol{v}_{p}}{\unicode[STIX]{x1D70F}_{v}}+\boldsymbol{g}, & \displaystyle\end{eqnarray}$$

where we assumed Stokes drag to be the principal force on the particle, $\boldsymbol{x}_{p}$ and $\boldsymbol{v}_{p}$ are the particle position and velocity, respectively, and $\unicode[STIX]{x1D70F}_{v}$ is the particle viscous relaxation time. In (2.9), $\boldsymbol{u}(\boldsymbol{x}_{p},t)$ is the fluid velocity at the particle’s location. We neglect two-way coupling effects, as well as particle collisions. In order to solve (2.8) and (2.9) numerically, $\boldsymbol{u}(\boldsymbol{x}_{p},t)$ needs to be evaluated. This is achieved by interpolating, to the particle position, fluid velocities at a stencil of grid points surrounding the particle. We use the eighth-order Lagrange interpolation method that is based on a stencil of $8\times 8\times 8$ fluid velocities.

Temporal update of particle motion is achieved through a modified second-order Runge–Kutta (RK2) method in which the standard RK2 weights are replaced by exponential integrators as follows Ireland et al. (Reference Ireland, Vaithianathan, Sukheswalla, Ray and Collins2013):

(2.10) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{p}(t_{0}+h)=\text{e}^{-h/\unicode[STIX]{x1D70F}_{v}}\boldsymbol{v}_{p}(t_{0})+w_{1}\boldsymbol{u}[\boldsymbol{x}_{p}(t_{0})]+w_{2}\boldsymbol{u}[\boldsymbol{x}_{p}(t_{0})+\boldsymbol{v}_{p}(t_{0})h]+(1-\text{e}^{-h/\unicode[STIX]{x1D70F}_{v}})\unicode[STIX]{x1D70F}_{v}\boldsymbol{g},\qquad & & \displaystyle\end{eqnarray}$$

where $h$ is the time step, and the exponential integrators $w_{1}$ and $w_{2}$ are given by

(2.11a,b ) $$\begin{eqnarray}\displaystyle w_{1}\equiv \left(\frac{h}{\unicode[STIX]{x1D70F}_{v}}\right)\left[\unicode[STIX]{x1D719}_{1}\left(\frac{-h}{\unicode[STIX]{x1D70F}_{v}}\right)-\unicode[STIX]{x1D719}_{2}\left(\frac{-h}{\unicode[STIX]{x1D70F}_{v}}\right)\right],\quad w_{2}\equiv \left(\frac{h}{\unicode[STIX]{x1D70F}_{v}}\right)\unicode[STIX]{x1D719}_{1}\left(\frac{-h}{\unicode[STIX]{x1D70F}_{v}}\right), & & \displaystyle\end{eqnarray}$$
(2.12a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{1}(z)\equiv \frac{\text{e}^{z}-1}{z},\quad \unicode[STIX]{x1D719}_{2}(z)\equiv \frac{\text{e}^{z}-z-1}{z^{2}}. & & \displaystyle\end{eqnarray}$$

For small values of $Fr$ , the periodic box length $L$ is an important consideration since it can artificially influence the motion of settling inertial particles. Specifically, the use of periodic boundary conditions is problematic if the time it takes the settling particles to traverse the length $L$ is ${\leqslant}O(T_{E})$ , where $T_{E}$ is the large eddy turnover time. Several studies (Woittiez et al. Reference Woittiez, Jonker and Portela2009; Ireland et al. Reference Ireland, Bragg and Collins2016; Dhariwal & Bragg Reference Dhariwal and Bragg2018) considered this issue in detail and found that box sizes larger than $L=2\unicode[STIX]{x03C0}$ may be needed, particularly when $St_{\unicode[STIX]{x1D702}}\gtrsim 1$ . In this paper, the smallest Froude number is $Fr=0.006$ . Considering $St_{\unicode[STIX]{x1D702}}=0.1$ , $Fr=0.006$ and $T_{E}=1.568$ , it can be shown that the settling times of particles through a domain length $L=2\unicode[STIX]{x03C0}$ is $T_{settle}^{L=2\unicode[STIX]{x03C0}}\approx 1.75$ . Since $T_{settle}^{L=2\unicode[STIX]{x03C0}}$ is only marginally greater than $T_{E}$ , we considered a box dimension of $4\unicode[STIX]{x03C0}$ along the direction of gravity for the $Fr=0.006$ case.

3 Correlation times in second drift closure: $T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}$ , $T_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ , $T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}$ , $T_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D716}}$

The drift flux $q_{i}^{d}(\boldsymbol{r},t)$ in the transport equation for the p.d.f. $\langle P\rangle$ is given by

(3.1) $$\begin{eqnarray}\displaystyle q_{i}^{d}(\boldsymbol{r},t)=-\langle P\rangle (\boldsymbol{r};t)\frac{St_{\unicode[STIX]{x1D702}}^{2}}{\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D702}}^{2}}r_{k}\int _{-\infty }^{t}d_{ik}\,\text{d}t^{\prime }, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D702}}$ is the inverse of the Kolmogorov time scale. In the Part 1 paper, we derived two closure forms, referred to as DF1 and DF2, for the integral on the right-hand side of (3.1). DF1 is based on the assumption that the fluid velocity gradient along particle trajectories has a Gaussian distribution. In DF2, we regard the strain-rate and rotation-rate tensors scaled by the turbulent dissipation rate and enstrophy, respectively, as normally distributed.

In the second closure of drift flux (DF2), this integral is given by

(3.2) $$\begin{eqnarray}\displaystyle \int _{-\infty }^{t}d_{ik}\,\text{d}t^{\prime } & = & \displaystyle \frac{1}{4\unicode[STIX]{x1D708}^{2}}\left\{\frac{1}{3}\unicode[STIX]{x1D6FF}_{ik}[\langle \unicode[STIX]{x1D716}^{2}\rangle T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}+\langle \unicode[STIX]{x1D716}\unicode[STIX]{x1D701}\rangle T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}-\langle \unicode[STIX]{x1D701}\unicode[STIX]{x1D716}\rangle T_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D716}}-\langle \unicode[STIX]{x1D701}^{2}\rangle T_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}]\right.\nonumber\\ \displaystyle & & \displaystyle +\,2\langle \unicode[STIX]{x1D716}^{2}\rangle \int _{-\infty }^{t}\exp \left(-\frac{t-t^{\prime }}{T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}}\right)\langle \unicode[STIX]{x1D70E}_{ij}(t)\unicode[STIX]{x1D70E}_{lm}(t^{\prime })\rangle \langle \unicode[STIX]{x1D70E}_{jk}(t)\unicode[STIX]{x1D70E}_{lm}(t^{\prime })\rangle \,\text{d}t^{\prime }\nonumber\\ \displaystyle & & \displaystyle -\left.2\langle \unicode[STIX]{x1D701}^{2}\rangle \int _{-\infty }^{t}\exp \left(-\frac{t-t^{\prime }}{T_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}}\right)\langle \unicode[STIX]{x1D70C}_{ij}(t)\unicode[STIX]{x1D70C}_{lm}(t^{\prime })\rangle \langle \unicode[STIX]{x1D70C}_{jk}(t)\unicode[STIX]{x1D70C}_{lm}(t^{\prime })\rangle \,\text{d}t^{\prime }\right\}.\end{eqnarray}$$

Equation (3.2) contains the auto- and cross-correlation times of the dissipation rate $\unicode[STIX]{x1D716}$ and enstrophy $\unicode[STIX]{x1D701}$ along particle trajectories – $T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}$ , $T_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ , $T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}$ , and $T_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D716}}$ . Since the particles are settling rapidly ( $Sv_{\unicode[STIX]{x1D702}}=g\unicode[STIX]{x1D70F}_{v}/u_{\unicode[STIX]{x1D702}}\gg 1$ ), we may regard the surrounding turbulence as essentially frozen during a particle response time $\unicode[STIX]{x1D70F}_{v}$ . Therefore, based on Taylor’s hypothesis, the Lagrangian time scales may be expressed as the respective spatial correlation lengths divided by the particle settling velocity. For instance, $T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}=L_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}/g\unicode[STIX]{x1D70F}_{v}$ and $T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}=L_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}/g\unicode[STIX]{x1D70F}_{v}$ . Due to isotropy, $L_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}=L_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D716}}$ , so that the terms $\langle \unicode[STIX]{x1D716}\unicode[STIX]{x1D701}\rangle T_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}}$ and $(-\langle \unicode[STIX]{x1D701}\unicode[STIX]{x1D716}\rangle T_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D716}})$ cancel on the right-hand side of (3.2). Therefore, the remaining unknown length scales are $L_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}$ and $L_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ . We now discuss the procedure for computing the length scale $L_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}$ through DNS (an analogous process is used to compute $L_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ ).

The length scale $L_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}$ is defined as

(3.3) $$\begin{eqnarray}\displaystyle L_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}=\frac{\displaystyle \int R_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}(r)\,\text{d}r}{\langle \unicode[STIX]{x1D716}^{2}\rangle }, & & \displaystyle\end{eqnarray}$$

where $R_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}(r)=\langle \unicode[STIX]{x1D716}(\boldsymbol{x};t)\unicode[STIX]{x1D716}(\boldsymbol{x}+\boldsymbol{r};t)\rangle$ is the spatial correlation of dissipation rate. The correlation $R_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}$ is evaluated using Fourier transforms as

(3.4) $$\begin{eqnarray}\displaystyle R_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}(r)=\int \text{d}\unicode[STIX]{x1D73F}\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D73F})\text{e}^{\text{i}\unicode[STIX]{x1D73F}\boldsymbol{\cdot }\boldsymbol{r}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}$ is the Fourier coefficient of $R_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}$ , and $\unicode[STIX]{x1D73F}$ and $\boldsymbol{r}$ are the wavenumber and relative position vectors, respectively. Expressing the integral $\int \!\text{d}\unicode[STIX]{x1D73F}$ in spherical coordinates, and performing the integrations in polar and azimuthal angles, we have

(3.5) $$\begin{eqnarray}\displaystyle R_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}(r)=\int _{1}^{\unicode[STIX]{x1D705}_{max}}\text{d}\unicode[STIX]{x1D705}D_{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D705})\frac{\sin (\unicode[STIX]{x1D705}r)}{\unicode[STIX]{x1D705}r}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{max}=\sqrt{2}N/3$ is the maximum resolved $\unicode[STIX]{x1D705}$ in the DNS ( $\unicode[STIX]{x1D705}=|\unicode[STIX]{x1D73F}|$ ), and

(3.6) $$\begin{eqnarray}\displaystyle D_{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D705})=\left\langle \mathop{\sum }_{|\unicode[STIX]{x1D73F}|=\unicode[STIX]{x1D705}}\widehat{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D73F},t)\widehat{\unicode[STIX]{x1D716}}^{\ast }(\unicode[STIX]{x1D73F},t)\right\rangle . & & \displaystyle\end{eqnarray}$$

In the above equation, $\widehat{(\cdot )}$ denotes the Fourier coefficient, and the superscript $\ast$ the complex conjugate. We evaluate $D_{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D705})$ using DNS, while the integral in (3.5) is calculated through numerical quadrature. In (3.6), $\langle \cdots \rangle$ denotes averaging over an ensemble of temporal snapshots of the statistically stationary turbulent velocity field. The length scales thus determined from the current DNS are given in table 1.

Table 1.

Correlation length scales of $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D701}$ for $Re_{\unicode[STIX]{x1D706}}=77.76$ .

4 DNS inputs to theory

Table 2.

Simulation parameters for the DNS study. All dimensional parameters are in arbitrary units. $Re_{\unicode[STIX]{x1D706}}\equiv u_{rms}\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}$ is the Taylor micro-scale Reynolds number, $u_{rms}\equiv \sqrt{(2k/3)}$ is the fluid root-mean-square fluctuating velocity, $k$ is the turbulent kinetic energy, $\unicode[STIX]{x1D708}$ is the fluid kinematic viscosity, $\unicode[STIX]{x1D716}\equiv 2\unicode[STIX]{x1D708}\int _{0}^{\unicode[STIX]{x1D705}_{max}}\unicode[STIX]{x1D705}^{2}E(\unicode[STIX]{x1D705})\,\text{d}\unicode[STIX]{x1D705}$ is the dissipation rate of turbulent kinetic energy, $L\equiv 3\unicode[STIX]{x03C0}/(2k)\int _{0}^{\unicode[STIX]{x1D705}_{max}}E(\unicode[STIX]{x1D705})/\unicode[STIX]{x1D705}\,\text{d}\unicode[STIX]{x1D705}$ is the integral length scale, $\unicode[STIX]{x1D706}\equiv u_{rms}\sqrt{(15\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716})}$ is the Taylor micro-scale, $\unicode[STIX]{x1D702}\equiv \unicode[STIX]{x1D708}^{3/4}/\unicode[STIX]{x1D716}^{1/4}$ is the Kolmogorov length scale, $T_{E}\equiv L/u_{rms}$ is the large-eddy turnover time, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\equiv \sqrt{(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716})}$ is the Kolmogorov time scale, $\unicode[STIX]{x1D705}_{max}$ is the maximum resolved wavenumber, $\unicode[STIX]{x0394}t$ is the time step, and $N_{p}$ is the number of particles per Stokes number.

To be able to consistently compare theory and DNS, it is important that the theory use the same turbulence parameters as those in statistically stationary DNS. Hence, inputs to the theory such as the Kolmogorov and integral length scales, dissipation rate and $Re_{\unicode[STIX]{x1D706}}$ are all identical to those in table 2, which lists the turbulence quantities in our DNS runs. In particular, one also has to ensure that the energy spectrum $E(\unicode[STIX]{x1D705})$ needed in the theory (to compute the drift and diffusion flux coefficients) closely matches the DNS energy spectrum. This was achieved by suitably selecting the parametric inputs to the model spectrum provided in Pope (Reference Pope2000), as follows:

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle E(\unicode[STIX]{x1D705})=C\unicode[STIX]{x1D716}^{2/3}\unicode[STIX]{x1D705}^{-5/3}f_{L}(\unicode[STIX]{x1D705}L)f_{\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D705}\unicode[STIX]{x1D702}), & \displaystyle\end{eqnarray}$$
(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle f_{L}(\unicode[STIX]{x1D705}L)=\left(\frac{\unicode[STIX]{x1D705}L}{[(\unicode[STIX]{x1D705}L)^{2}+c_{L}]^{1/2}}\right)^{5/3+p_{0}}, & \displaystyle\end{eqnarray}$$
(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle f_{\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D705}\unicode[STIX]{x1D702})=\text{exp}\{-\unicode[STIX]{x1D6FD}([(\unicode[STIX]{x1D705}\unicode[STIX]{x1D702})^{4}+c_{\unicode[STIX]{x1D702}}^{4}]^{1/4}-c_{\unicode[STIX]{x1D702}})\}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}=5.2$ and $p_{0}=2$  (Pope Reference Pope2000). The parameters $c_{L}$ and $c_{\unicode[STIX]{x1D702}}$ are determined from the following constraints:

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{3}{2}u_{rms}^{2}=\int _{1}^{\unicode[STIX]{x1D705}_{max}}E(\unicode[STIX]{x1D705})\,\text{d}\unicode[STIX]{x1D705}, & \displaystyle\end{eqnarray}$$
(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}=2\unicode[STIX]{x1D708}\int _{1}^{\unicode[STIX]{x1D705}_{max}}\unicode[STIX]{x1D705}^{2}E(\unicode[STIX]{x1D705})\,\text{d}\unicode[STIX]{x1D705}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D716}$ is the dissipation rate, and the wavenumber limits $[1,\unicode[STIX]{x1D705}_{max}]$ are the same as in DNS. These wavenumber limits are also used in calculating the drift and diffusion flux coefficients. The parameters $c_{L}$ and $c_{\unicode[STIX]{x1D702}}$ are numerically evaluated using the DNS values of $u_{rms}$ , $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D708}$ from table 2. The resulting values of $c_{L}$ and $c_{\unicode[STIX]{x1D702}}$ are shown in table 3. In figure 1, the model spectrum calculated from (4.1)–(4.3) is compared with the DNS energy spectrum for $Re_{\unicode[STIX]{x1D706}}=77.76$ . Good agreement is seen between the model and DNS spectra.

Figure 1.

Comparison of the DNS and model energy spectra at $Re_{\unicode[STIX]{x1D706}}=77.76$ .

Table 3.

Parameters for the model energy spectrum at $Re_{\unicode[STIX]{x1D706}}=77.76$ . After determining $c_{L}$ and $c_{\unicode[STIX]{x1D702}}$ , the parameter $C$ was adjusted to match the DNS energy spectrum. Pope (Reference Pope2000) suggested $C=1.5$ .

5 Results

We have shown in the part 1 paper that the p.d.f. of pair separation is given by

(5.1) $$\begin{eqnarray}\displaystyle \langle P\rangle (r,\unicode[STIX]{x1D707})=r^{\unicode[STIX]{x1D6FD}_{2}St_{\unicode[STIX]{x1D702}}^{2}}\frac{1}{4\unicode[STIX]{x03C0}}\left[1+St_{\unicode[STIX]{x1D702}}^{2}\unicode[STIX]{x1D6FD}_{2}\left(\frac{1}{2}\ln (1-\unicode[STIX]{x1D707}^{2})-(\ln 2-1)\right)\right], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D707}=\cos \unicode[STIX]{x1D703}$ , with $\unicode[STIX]{x1D703}\in (-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0})$ being the spherical polar angle. The coefficient $\unicode[STIX]{x1D6FD}_{2}$ in the power-law exponent ( $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{2}St_{\unicode[STIX]{x1D702}}^{2}$ ) is given by

(5.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{2}=\frac{\unicode[STIX]{x1D706}_{2}+2\unicode[STIX]{x1D706}_{1}}{\unicode[STIX]{x1D6FC}_{2}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{1}$ and $\unicode[STIX]{x1D706}_{2}$ are the drift flux coefficients, and $\unicode[STIX]{x1D6FC}_{2}$ is the diffusion flux coefficient. These coefficients are determined through a numerical quadrature process, discussed in Part 1. It may be recalled that there are two forms of $\unicode[STIX]{x1D706}_{1}$ and $\unicode[STIX]{x1D706}_{2}$ , corresponding to the two drift closures DF1 and DF2. Thus, we will compare two theoretical values of $\unicode[STIX]{x1D6FD}_{2}$ with the corresponding DNS value.

DNS were performed at the Reynolds number $Re_{\unicode[STIX]{x1D706}}=77.76$ for three Froude numbers $Fr=\infty ,0.052,0.006$ , where $Fr$ is the ratio of the gravitational and Kolmogorov accelerations. Thus, $Fr=\infty$ corresponds to zero gravity, and $Fr=0.006$ to the highest magnitude of gravity considered in the current DNS. For each $Fr$ , particles of six Stokes numbers $St_{\unicode[STIX]{x1D702}}=0.01,0.02,0.04,0.08,0.15,0.2$ were tracked. The spatial clustering of particles is quantified both through the radial distribution function (r.d.f.) that is independent of the polar angle ( $\unicode[STIX]{x1D703}$ ), and the angular distribution function (a.d.f) that includes the dependence on $\unicode[STIX]{x1D703}$ (Ireland et al. Reference Ireland, Bragg and Collins2016). The r.d.f. $g(r)$ scales with separation as $g(r)\sim r^{\unicode[STIX]{x1D6FD}}$ for $r\lesssim \unicode[STIX]{x1D702}$ (Reade & Collins Reference Reade and Collins2000). The DNS value of the exponent $\unicode[STIX]{x1D6FD}$ is then determined through a least squares curve fit of the r.d.f. for $0.6\unicode[STIX]{x1D702}\leqslant r\leqslant 2.5\unicode[STIX]{x1D702}$ .

In figure 2, we compare the exponent $\unicode[STIX]{x1D6FD}$ obtained from DNS and theory. The theoretical $\unicode[STIX]{x1D6FD}$ for DF1 and DF2 (both valid in the $Fr\ll 1$ limit), and the DNS $\unicode[STIX]{x1D6FD}$ for the three $Fr$ are presented in figure 2. A key feature of DF2 is that it accounts for the two-time autocorrelations and cross-correlations of dissipation rate and enstrophy (or, equivalently of the second invariants of strain-rate and rotation-rate tensors). Chun et al. (Reference Chun, Koch, Rani, Ahluwalia and Collins2005) showed that these correlations quantify, as well as illustrate, the mechanisms driving the clustering of non-settling particles. It is evident from figure 2 that the DF1 $\unicode[STIX]{x1D6FD}$ are significantly lower than those of DNS for all three Froude numbers. However, the DF2 $\unicode[STIX]{x1D6FD}$ are in reasonable agreement with the DNS $\unicode[STIX]{x1D6FD}$ for $Fr=0.006$ . The improved performance of DF2, compared to DF1, is because DF2 accounts for the correlations of dissipation rate and enstrophy, but DF1 does not.

Figure 2.

Comparison of the power-law exponent $\unicode[STIX]{x1D6FD}$ obtained from theory and DNS. Results obtained using both DF1 and DF2 are shown (referred to as Theory 1 and Theory 2, respectively). There is an uncertainty of ${\sim}$ 8 % in the DNS values of $\unicode[STIX]{x1D6FD}$ .

Next, we consider the anisotropy in particle clustering due to settling. Ireland et al. (Reference Ireland, Bragg and Collins2016) quantified the anisotropy through the use of the angular distribution function (a.d.f.), $g(\boldsymbol{r})$ , which is defined as (Ireland et al. Reference Ireland, Bragg and Collins2016)

(5.3) $$\begin{eqnarray}g(\boldsymbol{r})=\frac{N(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})/V(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})}{N/V},\end{eqnarray}$$

where $N(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ is the number of particle pairs in a solid-angle region of volume $V(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ , while $N$ is the total number of pairs in the domain of volume $V$ . The differential volume $V(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ is defined as (Ireland et al. Reference Ireland, Bragg and Collins2016)

(5.4) $$\begin{eqnarray}V(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\sin (\unicode[STIX]{x1D703})\unicode[STIX]{x0394}\unicode[STIX]{x1D703}\unicode[STIX]{x0394}\unicode[STIX]{x1D719}[(r+\unicode[STIX]{x0394}r)^{3}-(r+\unicode[STIX]{x0394}r)^{3}]/3.\end{eqnarray}$$

They then expressed $g(\boldsymbol{r})$ in terms of the Legendre spherical harmonic functions, as

(5.5) $$\begin{eqnarray}\frac{g(\boldsymbol{r})}{g(r)}=\mathop{\sum }_{l=1}^{\infty }\frac{{\mathcal{C}}_{2l}^{0}(r)}{{\mathcal{C}}_{0}^{0}(r)}Y_{2l}^{0}(\cos \unicode[STIX]{x1D703}),\end{eqnarray}$$

where

(5.6) $$\begin{eqnarray}g(r)={\mathcal{C}}_{0}^{0}(r)=\int _{0}^{\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}g(\boldsymbol{r}).\end{eqnarray}$$

Applying the orthogonality of Legendre polynomials to (5.5), we get

(5.7) $$\begin{eqnarray}\displaystyle \frac{{\mathcal{C}}_{2}^{0}(r)}{{\mathcal{C}}_{0}^{0}(r)}=\frac{5}{2}\frac{\displaystyle \int _{0}^{\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}g(\boldsymbol{r})Y_{2}^{0}(\cos \unicode[STIX]{x1D703})}{g(r)}. & & \displaystyle\end{eqnarray}$$

The theoretical value of the coefficient ratio is

(5.8) $$\begin{eqnarray}\displaystyle \left[\frac{{\mathcal{C}}_{2}^{0}(r)}{{\mathcal{C}}_{0}^{0}(r)}\right]_{theory} & = & \displaystyle \frac{5}{2}\frac{\displaystyle \int _{0}^{\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}\langle P\rangle (r,\unicode[STIX]{x1D703})Y_{2}^{0}(\cos \unicode[STIX]{x1D703})}{\displaystyle \int _{0}^{\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}\langle P\rangle (r,\unicode[STIX]{x1D703})}\nonumber\\ \displaystyle & = & \displaystyle \frac{5\unicode[STIX]{x1D6FD}}{12}=\frac{5\unicode[STIX]{x1D6FD}_{2}St_{\unicode[STIX]{x1D702}}^{2}}{12}.\end{eqnarray}$$

In table 4, we compare the values of ${\mathcal{C}}_{2}^{0}(r)/{\mathcal{C}}_{0}^{0}(r)$ obtained using theory and DNS, the latter for $Fr=0.006$ . We see that the degree of anisotropy predicted by the theory is in reasonable agreement with that computed using DNS, particularly for $St_{\unicode[STIX]{x1D702}}\leqslant 0.10$ . However, the theory generally overpredicts the coefficient ratio as compared to DNS.

Table 4.

Comparison of ${\mathcal{C}}_{2}^{0}(r)/{\mathcal{C}}_{0}^{0}(r)$ obtained using theory (DF2 only) and current DNS for $Fr=0.006$ . There is an uncertainty of ${\sim}$ 10 % in the DNS values of the coefficient ratio.

6 Conclusions

Part 2 of this study focuses on the quantitative analysis of the theory through a direct comparison of theory predictions with the data obtained in our DNS runs. While the theory is derived in the $Fr\ll 1$ regime, DNS were performed for three Froude numbers $Fr=\infty ,0.052,0.006$ at $Re_{\unicode[STIX]{x1D706}}=77.76$ . In the DNS runs, the $Fr=0.006$ case represented the run with the highest magnitude of gravitational acceleration. For this case, a domain size of $4\unicode[STIX]{x03C0}$ was used in the direction of gravity to mitigate the numerical errors arising from particles spuriously sampling the same eddies more than once as they settle through the domain. A model energy spectrum was derived that closely matched the DNS energy spectrum. The correlation length scales of dissipation rate and enstropy, $L_{\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}}$ and $L_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}}$ , needed for DF2 were also obtained using DNS. The model energy spectrum, and the two correlation lengths were then used in the calculation of the drift and diffusion coefficients needed for the power law exponent, as well as the coefficient ratio (for quantifying anisotropy). In terms of the dependence of clustering on separation, we see that the DF1 $\unicode[STIX]{x1D6FD}$ are significantly lower than the DNS $\unicode[STIX]{x1D6FD}$ for the three Froude numbers. But, the DF2 $\unicode[STIX]{x1D6FD}$ are in reasonable quantitative agreement with the DNS values for $Fr=0.006$ . These results demonstrate that the two-time correlations of dissipation rate and enstropy constitute an important mechanism driving the drift flux responsible for clustering. We also quantified the anisotropy in particle clustering due to the settling. We see that the values of the coefficient ratio ${\mathcal{C}}_{2}^{0}(r)/{\mathcal{C}}_{0}^{0}(r)$ obtained using DF2 are in reasonable agreement with the DNS values at $Fr=0.006$ , particularly for $St_{\unicode[STIX]{x1D702}}\leqslant 0.1$ . Despite the high magnitude of gravity in both theory and DNS, we observe only marginal anisotropy in clustering, which can be explained based on the theory. We see from the theory (5.8) that the spherical harmonic coefficient scales with the exponent $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{2}St_{\unicode[STIX]{x1D702}}^{2}$ , where $\unicode[STIX]{x1D6FD}_{2}$ is the ratio of the drift and diffusion flux coefficients. We know that $\unicode[STIX]{x1D6FD}_{2}$ decreases monotonically with gravity, as the drift flux driving clustering is also diminished as gravity is increased. Therefore, at higher gravity, the spherical harmonic coefficients quantifying anisotropy are also correspondingly smaller. The current two-part study presents the development and analysis of an analytical theory for the clustering of low-inertia particle pairs that are settling rapidly in isotropic turbulence. We see that the theory accurately captures the qualitative trends in particle clustering. It also shows reasonable quantitative agreement with DNS data for low Stokes and Froude numbers.

Acknowledgement

S.L.R. gratefully acknowledges NSF support through the grant CBET-1436100.

References

Ayala, O., Rosa, B., Wang, L.-P. & Grabowski, W. W. 2008 Effects of turbulence on the geometric collision rate of sedimenting droplets. Part 1. Results from direct numerical simulation. New J. Phys. 10 (7), 075015.Google Scholar
Bec, J., Homann, H. & Ray, S. S. 2014 Gravity-driven enhancement of heavy particle clustering in turbulent flow. Phys. Rev. Lett. 112 (18), 184501.10.1103/PhysRevLett.112.184501Google Scholar
Brucker, K. A., Isaza, J. C., Vaithianathan, T. & Collins, L. R. 2007 Efficient algorithm for simulating homogeneous turbulent shear flow without remeshing. J. Comput. Phys. 225, 20–32.10.1016/j.jcp.2006.10.018Google Scholar
Chun, J., Koch, D. L., Rani, S. L., Ahluwalia, A. & Collins, L. R. 2005 Clustering of aerosol particles in isotropic turbulence. J. Fluid Mech. 536, 219–251.10.1017/S0022112005004568Google Scholar
Dhariwal, R. & Bragg, A. D. 2018 Small-scale dynamics of settling, bidisperse particles in turbulence. J. Fluid Mech. 839, 594–620.10.1017/jfm.2018.24Google Scholar
Ireland, P. J., Bragg, A. D. & Collins, L. R. 2016 The effect of Reynolds number on inertial particle dynamics in isotropic turbulence. Part 2. Simulations with gravitational effects. J. Fluid Mech. 796, 659–711.10.1017/jfm.2016.227Google Scholar
Ireland, P. J., Vaithianathan, T., Sukheswalla, P. S., Ray, B. & Collins, L. R. 2013 Highly parallel particle-laden flow solver for turbulence research. Comput. Fluids 76, 170–177.10.1016/j.compfluid.2013.01.020Google Scholar
Onishi, R., Takahashi, K. & Komori, S. 2009 Influence of gravity on collisions of monodispersed droplets in homogeneous isotropic turbulence. Phys. Fluids 21 (12), 125108.10.1063/1.3276906Google Scholar
Parishani, H., Ayala, O., Rosa, B., Wang, L.-P. & Grabowski, W. W. 2015 Effects of gravity on the acceleration and pair statistics of inertial particles in homogeneous isotropic turbulence. Phys. Fluids 27 (3), 033304.10.1063/1.4915121Google Scholar
Pope, S. B. 2000 Turbulent Flows. Cambridge University Press.10.1017/CBO9780511840531Google Scholar
Rani, S. L., Gupta, V. K. & Koch, D. L. 2019 Clustering of rapidly settling, low-inertia particle pairs in isotropic turbulence. Part 1. Drift and diffusion flux closures. J. Fluid Mech. 871, 450–476.10.1017/jfm.2019.204Google Scholar
Reade, W. C. & Collins, L. R. 2000 Effect of preferential concentration on turbulent collision rates. Phys. Fluids 12 (10), 2530–2540.10.1063/1.1288515Google Scholar
Witkowska, A., Juvé, D. & Brasseur, J. G. 1997 Numerical study of noise from isotropic turbulence. J. Comput. Acoust. 5 (03), 317–336.10.1142/S0218396X97000186Google Scholar
Woittiez, E. J. P., Jonker, H. J. J. & Portela, L. M. 2009 On the combined effects of turbulence and gravity on droplet collisions in clouds: a numerical study. J. Atmos. Sci. 66 (7), 1926–1943.10.1175/2005JAS2669.1Google Scholar
Figure 0

Table 1.Correlation length scales of $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D701}$ for $Re_{\unicode[STIX]{x1D706}}=77.76$.

Figure 1

Table 2.Simulation parameters for the DNS study. All dimensional parameters are in arbitrary units. $Re_{\unicode[STIX]{x1D706}}\equiv u_{rms}\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}$ is the Taylor micro-scale Reynolds number, $u_{rms}\equiv \sqrt{(2k/3)}$ is the fluid root-mean-square fluctuating velocity, $k$ is the turbulent kinetic energy, $\unicode[STIX]{x1D708}$ is the fluid kinematic viscosity, $\unicode[STIX]{x1D716}\equiv 2\unicode[STIX]{x1D708}\int _{0}^{\unicode[STIX]{x1D705}_{max}}\unicode[STIX]{x1D705}^{2}E(\unicode[STIX]{x1D705})\,\text{d}\unicode[STIX]{x1D705}$ is the dissipation rate of turbulent kinetic energy, $L\equiv 3\unicode[STIX]{x03C0}/(2k)\int _{0}^{\unicode[STIX]{x1D705}_{max}}E(\unicode[STIX]{x1D705})/\unicode[STIX]{x1D705}\,\text{d}\unicode[STIX]{x1D705}$ is the integral length scale, $\unicode[STIX]{x1D706}\equiv u_{rms}\sqrt{(15\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716})}$ is the Taylor micro-scale, $\unicode[STIX]{x1D702}\equiv \unicode[STIX]{x1D708}^{3/4}/\unicode[STIX]{x1D716}^{1/4}$ is the Kolmogorov length scale, $T_{E}\equiv L/u_{rms}$ is the large-eddy turnover time, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\equiv \sqrt{(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716})}$ is the Kolmogorov time scale, $\unicode[STIX]{x1D705}_{max}$ is the maximum resolved wavenumber, $\unicode[STIX]{x0394}t$ is the time step, and $N_{p}$ is the number of particles per Stokes number.

Figure 2

Figure 1.Comparison of the DNS and model energy spectra at $Re_{\unicode[STIX]{x1D706}}=77.76$.

Figure 3

Table 3.Parameters for the model energy spectrum at $Re_{\unicode[STIX]{x1D706}}=77.76$. After determining $c_{L}$ and $c_{\unicode[STIX]{x1D702}}$, the parameter $C$ was adjusted to match the DNS energy spectrum. Pope (2000) suggested $C=1.5$.

Figure 4

Figure 2.Comparison of the power-law exponent $\unicode[STIX]{x1D6FD}$ obtained from theory and DNS. Results obtained using both DF1 and DF2 are shown (referred to as Theory 1 and Theory 2, respectively). There is an uncertainty of ${\sim}$8 % in the DNS values of $\unicode[STIX]{x1D6FD}$.

Figure 5

Table 4.Comparison of ${\mathcal{C}}_{2}^{0}(r)/{\mathcal{C}}_{0}^{0}(r)$ obtained using theory (DF2 only) and current DNS for $Fr=0.006$. There is an uncertainty of ${\sim}$10 % in the DNS values of the coefficient ratio.