Scale Dependence of Multiplier Distributions for Particle Concentration, Enstrophy and Dissipation in the Inertial Range of Homogeneous Turbulence

TL;DR

This study uses DNS data to analyze how cascade multipliers for particle concentration, enstrophy, and dissipation depend on scale and local Stokes number, providing insights into turbulent clustering at high Reynolds numbers.

Contribution

It introduces a method to correct small particle number effects and demonstrates scale-dependent multiplier PDFs that collapse when scaled with local Stokes number, advancing turbulence cascade modeling.

Findings

Multiplier PDFs are scale dependent but collapse with local Stokes number.

Dissipation and enstrophy multiplier PDFs reach an asymptote at small scales.

Predicted RDFs agree with high-Reynolds-number inertial range theory.

Abstract

Turbulent flows preferentially concentrate inertial particles depending on their stopping time or Stokes number, which can lead to significant spatial variations in the particle concentration. Cascade models are one way to describe this process in statistical terms. Here, we use a direct numerical simulation (DNS) dataset of homogeneous, isotropic turbulence to determine probability distribution functions (PDFs) for cascade multipliers, which determine the ratio by which a property is partitioned into sub-volumes as an eddy is envisioned to decay into smaller eddies. We present a technique for correcting effects of small particle numbers in the statistics. We determine multiplier PDFs for particle number, flow dissipation, and enstrophy, all of which are shown to be scale dependent. However, the particle multiplier PDFs collapse when scaled with an appropriately defined local Stokes number. As anticipated from earlier works, dissipation and enstrophy multiplier PDFs reach an asymptote for sufficiently small spatial scales. From the DNS measurements, we derive a cascade model that is used it to make predictions for the radial distribution function (RDF) for arbitrarily high Reynolds numbers, $Re$, finding good agreement with the asymptotic, infinite $Re$ inertial range theory of Zaichik and Alipchenkov [New Journal of Physics 11, 103018 (2009)]. We discuss implications of these results for the statistical modeling of the turbulent clustering process in the inertial range for high Reynolds numbers inaccessible to numerical simulations.

Figures

Tables

Table 1. Table 1: RDF asymptotic values, g 0 subscript 𝑔 0 g_{0} , for S ​ t L − 3 / 2 ​ r / L ≪ 1 much-less-than 𝑆 superscript subscript 𝑡 𝐿 3 2 𝑟 𝐿 1 St_{L}^{-3/2}r/L\ll 1 and S ​ t L ≪ 1 much-less-than 𝑆 subscript 𝑡 𝐿 1 St_{L}\ll 1 for the different cascade simulation methods, and the  Zaichik and Alipchenkov ( 2009 ) model.

Simulation / Model		$g_0$
Cascade simulation – Method A		13
Cascade simulation – Method B		5.9
Cascade simulation – Method C		3.9
Zaichik and Alipchenkov (2009), figure 3		3.9

Equations

$$r/L=2^{-N/3}.$$

$$P_P(n;\bar{n})=\frac{{\bar{n}}^n{e}^{-\bar{n}}}{n!}$$

$$P_B(k;n,p)=\left(\genfrac{}{}{0ex}{}{n}{k}\right)p^n(1-p{)}^{n-k}.$$

$$P(m;\bar{n},p)=\sum _n{P}_P(n;\bar{n})P_B(k=mn;n,p).$$

$$f(m;\beta )={\left(m-m^2\right)}^{\beta -1}\frac{\Gamma (2\beta )}{2\Gamma (\beta )}$$

$$\sigma (\beta )=\sqrt{\frac{1}{4\left(2\beta +1\right)}}.$$

$$\bar{\beta }=\beta (\bar{\sigma })$$

$${\bar{\sigma }}^2=\frac{1}{M}\sum _{i=1}^M\sigma ({\beta }_i{)}^2$$

$$\text{VAR}(\ln \bar{\beta })=\frac{\text{VAR}(\bar{\beta })}{\bar{\beta }}\approx$$

$${\left(\frac{\partial \beta }{\partial {\sigma }^2}({\bar{\sigma }}^2)\right)}^2\text{VAR}({\bar{\sigma }}^2)=\frac{\text{VAR}({\bar{\sigma }}^2)}{64{\bar{\sigma }}^4},$$

$$\text{VAR}({\bar{\sigma }}^2)=\frac{1}{M}\sum _{i=1}^M{\left({\bar{\sigma }}^2-\sigma ({\beta }_i{)}^2\right)}^2.$$

$$S{t}_r\equiv \frac{{\tau }_s}{{\tau }_r}=St{\left(\frac{r}{\eta }\right)}^{-2/3},$$

$${\tau }_r={\tau }_{\eta }{\left(\frac{r}{\eta }\right)}^{2/3}$$

$$S{t}_r=2^{2N/9}S{t}_L,$$

$$S{t}_L\equiv \frac{{\tau }_s}{{\tau }_L}=St{\left(\frac{L}{\eta }\right)}^{-2/3}$$

$$ϵ=2\nu S_{ij}{S}_{ij}=\nu \left[({\partial }_i{u}_j)({\partial }_i{u}_j)+({\partial }_i{u}_j)({\partial }_j{u}_i)\right],$$

$$\mathcal{E}=|\nabla \times \vec{u}{|}^2=({\partial }_i{u}_j)({\partial }_i{u}_j)-({\partial }_i{u}_j)({\partial }_j{u}_i).$$

$$\beta (S{t}_r)\approx {\beta }_{\text{min}}\left({\left(\frac{S{t}_r}{a_1}\right)}^{b_1}+{\left(\frac{S{t}_r}{a_2}\right)}^{b_2}\right),$$

$$a_1=0.15,$$

$$a_2=0.45-0.25\exp \left(-\frac{2}{3}\ln {\left(20.5\frac{r}{L}\right)}^2\right),$$

$$b_1=-1.2,$$

$$b_2=0.85+0.35\left[1+\mathrm{erf}\left(-1.8\ln \left(29.3\frac{r}{L}\right)\right)\right],$$

$${\beta }_{\text{min}}=4+4\left[1+\mathrm{erf}\left(\ln \left(4\frac{r}{L}\right)\right)\right],$$

$$a_1=0.15,a_2=0.45,$$

$$b_1=-1.2,b_2=1.55,{\beta }_{\text{min}}=4.$$

$$P_P(n;C_{\infty }\bar{n})=\frac{(C_{\infty }\bar{n}{)}^n{e}^{-C_{\infty }\bar{n}}}{n!}.$$

$$\sum _{i,j,k}\frac{1}{8}{C}_{ijk}^{(N+3)}=C^{(N)}.$$

$$(C_{111}^{(N+3)}+C_{211}^{(N+3)})m_1^{(N+3)}$$

$$(C_{112}^{(N+3)}+C_{212}^{(N+3)})m_3^{(N+3)}$$

$$(C_{111}^{(N+3)}+C_{121}^{(N+3)})m_5^{(N+3)}$$

$$(C_{112}^{(N+3)}+C_{122}^{(N+3)})m_7^{(N+3)}$$

$$(C_{121}^{(N+3)}+C_{221}^{(N+3)})m_2^{(N+3)}$$

$$(C_{122}^{(N+3)}+C_{222}^{(N+3)})m_4^{(N+3)}$$

$$C_{121}^{(N+3)}+C_{112}^{(N+3)})m_6^{(N+3)}$$

$$\sum _{j,k}{C}_{1ik}^{(N+3)}=2C^N{m}_1^{(N+1)},$$

$$\sum _{i,j}{C}_{ij1}^{(N+3)}=2C^N{m}_3^{(N+1)},$$

$$\frac{C_{111}^{(N+3)}}{C^N}$$

Full text

I Background and Introduction

Clustering of inertial (finite-stopping-time) particles into dense zones in fluid turbulence has applications in many fields (for a general review see Toschi and Bodenschatz, 2009). A number of recent papers have focussed on understanding the basic mechanisms responsible for this effect; several of these (Zaichik and Alipchenkov, 2009; Bragg and Collins, 2014a, b; Gustavsson and Mehlig, 2014) provide thorough reviews and comparisons of previous studies dating back to the early work of Maxey (1987) and Squires and Eaton (1991) which we will only sketch briefly. The early work emphasized the role of centrifugation of finite-inertia particles out of vortical structures in turbulence. More recent evidence that clustering arises even in random, irrotational flows suggests that, while vorticity still plays a role, the dominant role is played by so-called “history effects”, in which inertial particle velocity dispersions at any location carry a memory of particle encounters with more remote flow regimes which have larger characteristic velocity differences (Pan and Padoan, 2010, 2013; Bragg et al., 2015a). These history effects lead to spatial gradients in particle random relative velocities, and these gradients in turn generate systematic flows or currents which can outweigh dispersive effects and produce zones of highly variable particle concentration (Zaichik and Alipchenkov, 2003, 2009; Bragg and Collins, 2014a, b).

To date, by far the most attention regarding particle clustering in turbulence has been devoted to very small spatial scales $r<\eta$ or even $r\ll\eta$, where $\eta$ is the Kolmogorov scale, partly because it is on these scales that particle collisions occur and partly because numerical simulations to date have produced only very limited inertial ranges, at best (see however Bec et al., 2010a; Ireland et al., 2015). Theories by Zaichik and Alipchenkov (2003) et seq., and Pan and Padoan (2010) et seq. have been shown to be promising in explaining the cause of particle clustering in terms of history effects, with helpful contributions from the traditional local centrifugation mechanism (Zaichik and Alipchenkov, 2009; Bragg and Collins, 2014b; Bragg et al., 2015a; Ireland et al., 2015). A thorough review of the effects of clustering and relative velocity effects on particle collisions, emphasizing the astronomical literature, can be found in Pan and Padoan (2014, 2015).

Our focus is on clustering at larger scales in the inertial range $\eta<r<L$, where $L$ is the integral scale. Inertial range clustering has important applications for remote sensing of terrestrial clouds (Shaw et al., 2002), the formation of primitive planetesimals (asteroids and comets) in the early solar nebula (Cuzzi et al., 2008, 2010; Chambers, 2010; Pan et al., 2011; Hopkins, 2014; Johansen et al., 2015), and even the structure of the interstellar medium (Hopkins and Lee, 2015). While little studied in the context of particle clustering, inertial range scaling is known to have different properties than seen in the dissipation range $r<\eta$ (Bec et al., 2007; Bragg et al., 2015b). Only limited predictions have been made of its scaling properties at very high Reynolds number $Re$ (Zaichik and Alipchenkov, 2003, 2009).

In the inertial range, so-called cascade models which reproduce the statistics of fluid behavior, even if not realistic flow structures, may be valuable for modeling high Reynolds number ($Re$) regimes too demanding for direct numerical simulations. Their application is quite general (Meneveau and Sreenivasan, 1987a, 1991; Chhabra et al., 1989; Chhabra and Sreenivasan, 1992; Sreenivasan and Stolovitzky, 1994, 1995, see Hogan and Cuzzi (2007) for more references). We and others have used cascades to model particle clustering in turbulence in the astronomical applications mentioned above.

Cascade models operate by simply applying a partition function or multiplier $0\leq m\leq 1$ to any property $\cal P$ in some given volume of the flow, thus determining the ratio by which the property (dissipation, particle density, etc.) is partitioned into sub-volumes as an eddy is envisioned to decay into smaller eddies. The most common treatment is a binary cascade, in which $\cal P$ is partitioned into two equal subvolumes; however the approach can be applied to arbitrary numbers of subvolumes (Sreenivasan and Stolovitzky, 1994). The binary cascade operates on each volume of space, partitioning $\cal P$ into two equal subvolumes by multipliers $m$ and $1-m$, with the multiplier $m$ at each bifurcation drawn from a probability distribution function (PDF) of multipliers $P(m)$. If $P(m)=\delta(m-0.5)$, where $\delta$ is the delta function, the cascade has no effect because the property $\cal P$ is evenly divided, and remains constant per unit volume. On the other hand, broad $P(m)$ functions generate highly intermittent spatial distributions in which $\cal P$ has a wide range of values, fluctuating dramatically on small scales such as seen in dissipation (Meneveau and Sreenivasan, 1991; Sreenivasan and Stolovitzky, 1995) (figure 1b).

The dissipation range, a range of small scales approaching the Kolmogorov scale $\eta$, is found where $r<20-30\eta$ (Tennekes and Lumley, 1972; Frisch, 1995); in this range, where viscosity is important, the equations of motion are no longer completely scale-free, and fluid scaling properties differ from those in the inertial range. The properties of particle clustering do seem to be scale independent in this region, however (Bec et al., 2007; Bragg et al., 2015b), and one expects this regime to be flow-independent for high $Re$. There is also a range of large scales near the integral scale $L$, over which deviation from scale invariance surely occurs, but this range has not been well studied and is surely flow-dependent. The application to planetesimal formation has become focussed on particle concentration at scales much larger than the Kolmogorov scale (Cuzzi et al., 2008, 2010) because large clumps are needed for sufficiently rapid gravitational collapse. In previous particle clustering cascade models, Hogan and Cuzzi (2007) determined the multiplier PDFs for particle concentration and fluid enstrophy at small spatial scales, not too far from $\eta$ (to obtain better statistics), and applied them across all scales ranging up to the integral scale (see section V.1 for more discussion). Realizing the risks in this, they performed tests which seemed to validate the approach. However, discrepancies between Hogan and Cuzzi (2007) and Pan et al. (2011) at the low probabilities of interest for the planetesimal problem (Cuzzi et al., 2010) have led us to explore the scale dependence of $P(m)$ in more detail, in order to improve the fidelity of the cascade models.

It is worth noting at this point that it is not a requirement of cascade models that the PDFs be scale-independent; it is merely the first and most obvious assumption. In this paper we present evidence that the multiplier PDFs for particle concentration are scale-dependent and present simple guidelines for how this scale-dependence can be included in cascades. Multiplier PDFs can also be conditioned on local properties (Sreenivasan and Stolovitzky, 1995), and indeed were treated this way in our previous work to allow for particle mass-loading on the process (Hogan and Cuzzi, 2007). Scale-dependence per se is, however, a different effect than local conditioning, and in this paper we do not address local conditioning.

Before describing our own work, we review some experimental results on high-Reynolds number atmospheric boundary layer turbulence, which provide a useful background in scale invariance and complement the more typical, but lower-$Re$, numerical simulations. Studies of the properties of turbulence in atmospheric boundary layer flows have been conducted by Kholmyanskiy (1973); van Atta and Yeh (1975) and Meneveau and Sreenivasan (1987b); further analysis of the Meneveau and Sreenivasan (1987b) data was done by Chhabra et al. (1989).111These results were cited and reanalyzed by Sreenivasan and Stolovitzky (1995), however the reference to the basic data given by Sreenivasan and Stolovitzky (1995) is confused with an interpretive article by Chhabra and Jensen (1989), who themselves cite Meneveau and Sreenivasan (1987a) and Chhabra et al. (1989, at the time unpublished) for discussions and analysis of the basic data. The best reference for the basic experimental data is Meneveau and Sreenivasan (1987b, see their Table 1), who conducted an experiment on boundary layer turbulence using a sensor mounted 2 m above the flat roof of a four story building. The Reynolds number for the flow is calculated using the free stream velocity $U$ = 6 m/s and the height $h$ = 2 m of the sensor above the roof: $Re=Uh/\nu=8\times 10^{5}$ where the kinematic viscosity is $1.5\times 10^{-5}~{}$m2/sec, consistent with tabulated values in Meneveau and Sreenivasan (1987b) of the Taylor scale Reynolds number $Re_{\lambda}$, and its characteristic lengthscale $\lambda$ and velocity $u^{\prime}$. Meneveau and Sreenivasan (1987b) give the Kolmogorov scale as $\eta=7\times 10^{-4}$ m (we retain their preferred units). Analyses of flow structures by Chhabra et al. (1989) (their figure 5 and our figure 2) show fairly well-behaved power law scaling of dissipation for weightings which suppress regions that are strongly anomalous (panels g and h), i.e. strongly differing from the mean, to almost $r\sim 1.8\times 10^{4}~{}\eta$ = 12.6 m $\gg h$, suggesting an extensive inertial range. The large or integral scale $L$, which contains the energy in this flow, is thus apparently much larger than the vertical distance of the sensor from the boundary ($h=2$ m) and plausibly the same as the longitudinal integral scale given as $L$ = 180 m (Meneveau and Sreenivasan, 1987b), see also Hunt and Morrison (2000), and thus $L/\eta\sim 2\times 10^{5}$. However, when the role of strongly anomalous regions is emphasized (panels i through l of figure 2) the scalable inertial range contracts.

In cascade applications, it may be more meaningful to assess the scale dependence of $P(m)$ at large scales not in terms of multiples of $\eta$ as in Sreenivasan and Stolovitzky (1995) and most other work (e.g., Bec et al., 2007), but in terms of fractions of $L$ which more closely connects to causality and energy flow. We will also express scale fractions $r/L$ in terms of cascade bifurcation levels $N$ needed to achieve cubes $r$ on a side:

[TABLE]

For example, Sreenivasan and Stolovitzky (1995) compared multiplier PDFs for dissipation in the atmospheric boundary layer over a range of scales (see figure 1b) and showed that $P(m)$ is highly scale independent over a wide range of scales: $372\eta-3072\eta$, or $372\eta$ to $L/86$  (see also Chhabra and Sreenivasan, 1992). Dissipation depends on higher-order moments of the velocity gradients, so we are drawn for guidance to the behavior seen in the higher-order moments (larger $|q|$) in figure 2 (Chhabra et al., 1989). The results of Sreenivasan and Stolovitzky (1995) are consistent with the generally power law behavior seen for $25\eta-10000\eta$ (roughly $25\eta$ to $L/26$) in figure 2. That is, one might infer from where the plots in Chhabra et al. (1989) deviate from power law behavior, that the scale-free behavior demonstrated by Sreenivasan and Stolovitzky (1995) (figure 1b) might carry on to larger sizes than they actually presented, possibly until $r\sim L/26$ or $10000\eta$, but deviate noticeably for scales larger than $r\sim L/15$ (and at the smaller end below $25\eta$). Moreover, we can conclude from these comparisons that the scale-free behavior seen by Sreenivasan and Stolovitzky (1995) was safely out of the viscous range, and continued through the inertial range at scales up to $L/86<L/10$.

The goal of this paper is to use DNS data to derive probability distribution functions for cascade multipliers and construct a cascade model that can be used for modeling higher $Re$-number flows not accessible to direct numerical simulations. The paper is organized as follows: section II describes the DNS dataset used in this study; section III describes the data analysis including a novel technique for correcting the effects of small particle number statistics, and presents results for the multiplier PDFs for particle concentration, dissipation and enstrophy; section IV presents predictions of the new cascade model and comparison with DNS data at two different $Re$; and section V discusses the results and their implications. A summary and conclusions are given in section VI.

II Dataset

In this paper, we use data from the direct numerical simulations of Bec et al. (2010a); see also Arnèodo et al. (2008) and Benzi et al. (2009). The simulation computes forced, homogenous and isotropic turbulence in an incompressible fluid, and the dynamics of inertial particles suspended in the flow. The fluid flow is solved on a $2048^{3}$ Cartesian grid with a grid spacing that is approximately the Kolmogorov length scale $\eta\approx\Delta x=\Delta y=\Delta z$. Tracer and inertial particles are introduced into the flow and their trajectories are tracked. Particles are considered point particles and are dragged with the flow by viscous forces only; there is no back-reaction on the flow. Particles of different Stokes numbers $St\equiv\tau_{s}/\tau_{\eta}$ are considered, where $\tau_{s}$ is the aerodynamic stopping time of the particle ($\tau_{s}=0$ for tracers) and $\tau_{\eta}$ is the Kolmogorov time.

Figure 1a shows the second order structure function for particle velocity for this numerical flow (taken along trajectories of different $St$ particles, from figure 2 of Bec et al., 2010b), as a function of normalized scale $r/\eta$. While the structure function seems to show an inertial range to several thousand $\eta$, in reality the integral eddy scale for this simulation seems to be about half the computational box size, $L\sim 1024\eta$, and for this flow $Re_{\lambda}\approx 4(L/\eta)^{2/3}\sim$ 400 (table 1 of Bec et al., 2010b). In blue-green below the horizontal axis we indicate the corresponding values for $r/L$, and the corresponding cascade level $N$. The blue dashed line indicates the expected inertial range for the atmospheric flow of Meneveau and Sreenivasan (1987b), with corresponding values of $r/L$ and $N$ also indicated in blue above the figure. Note that the range where $P(m)$ for dissipation was observed to be scale independent by Sreenivasan and Stolovitzky (1995) corresponds roughly to the scale range $L/860-L/86$, well below the expected integral scale for that flow and well above the viscous subrange.

Data from this simulation are available publicly online (Lanotte et al., 2011), and we have downloaded and analyzed all of the publicly available data in the present work. This data consists of the entire flow field sampled at 13 instances in time, and particle trajectories sampled at 4,720 equidistant times, both covering about 6 large-eddy time scales $\tau_{L}$. All flow components and their first derivatives are available at the particle locations. In total there are $3\times 64$ files of particle trajectories each containing 3,184 particles (a total of $N_{p}\approx 600\textrm{k}$ particles) for each $St=0,0.16,0.6$ and $1.0$, and 64 files containing 3,184 particles each (i.e., a total of $N_{p}\approx 200\textrm{k}$ particles) for each $St=2,3,5,10,20,30,40,50,$ and $70$.

III Analysis

Determining concentration multipliers amounts to counting particles in cubic sub-volumes of size $r^{3}$ and calculating the fraction of particles falling in each half of the sampling box. We bisect each cube in all three orthogonal directions $x$, $y$, and $z$ each yielding 2 multiplier values totaling 6 multiplier values for $r^{3}$ cube. The available trajectory data is highly resolved in time (4,720 instances of time over approximately 6 large eddy times $\tau_{L}$), much more than what is needed for this analysis. The number of snapshots required for good statistics depends on the scale of interest since structures at large scale evolve more slowly than structures at small scale (and contain more particles), and therefore can be sampled less often. We choose to sample the particle data with a temporal spacing of $\tau_{\textrm{sample}}\approx 0.55~{}\tau_{r}$ where $\tau_{r}$ is the characteristic eddy life time at spatial scale $r$ estimated using Kolomogorov 1941 arguments as $\tau_{r}=\tau_{L}(r/L)^{2/3}$. For the box sizes considered, $512\eta$, $256\eta$, $128\eta$, $64\eta$, $45\eta$, $32\eta$, $24\eta$, $16\eta$ and $12\eta$, this results in sampling intervals of $0.34\tau_{L}$, $0.22\tau_{L}$, $0.14\tau_{L}$, $0.086\tau_{L}$, $0.067\tau_{L}$, $0.055\tau_{L}$, $0.044\tau_{L}$, $0.034\tau_{L}$ and $0.028\tau_{L}$, respectively. We populate the sample volume using the positions of all particles from the high-resolution trajectory files at these various discrete times.

III.1 Tracer particles

First, we will look at the tracer particles (Stokes number $St=0$) which follow the flow exactly. Particles are initially homogeneously distributed, and since the flow is incompressible, will stay that way on the average, with particle multiplier distributions given by a delta function at $m=0.5$. However, the observed distributions will only approach this behavior for very large number of particles; otherwise, effects of small number statistics complicate matters. The simulation dataset at hand has $N_{p}\approx 6\times 10^{5}$ particles with $St=0$, and this is small enough to cause significant deviations from ideal behavior at small scales.

For $St=0$ particles, we can quantify these small-number effects analytically. First, the probability of finding $n$ particles in a given box of size $r$ is governed by a Poisson distribution

[TABLE]

with the expectation value $E(n)=\bar{n}(r)$, which is here the average number of particles in a box of size $r$, that is $\bar{n}(r)=N_{p}r^{3}/{\cal L}^{3}$, and where $\cal L$ is the size of the simulation domain. This probability is also the probability of finding a concentration $C=n/{\bar{n}}$. A comparison of the observed probability distribution functions with equation (3.1) is shown in figure 3a.

Then, for a given box with exactly $n$ particles, consider the particles one at a time and ask if they fall into one half of the box, say the left half, or the other, right half. For tracer particles, the probability to fall into the left side is the same as falling into the right side of the box, i.e. $p_{\textrm{left}}=p_{\textrm{right}}=p=0.5$, and the probability of having exactly $k$ particles fall into one side of the box is then given by a binomial distribution with the probability

[TABLE]

This is then also the probability of finding a multiplier $m=k/n$ in a box of $n$ particles.

By combining the probabilities (3.1) and (3.2), we see that the probability of finding a multiplier $m$ in the entire simulation domain is

[TABLE]

A comparison of this analytical relation with the multiplier PDFs computed from the tracer particle trajectories in the simulation is shown in figure 3b. It shows that equation (3.3) models the observed distributions very accurately, and also that the distributions become rather wide at small scales even though the underlying probability distributions are delta functions at $m=0.5$. This effect is a kind of “false intermittency” due to small-number statistics alone.

III.2 Particle multipliers for non-zero Stokes numbers

III.2.1 Correcting for finite particle numbers

As we have seen in the previous section, the number of particles in the dataset is small enough to significantly affect the observed multiplier distributions. In the following we will describe how we can account for these effects and estimate what the underlying PDFs would be, given infinite particle numbers. The goal is to separate the finite-particle number effects, which may be important in many applications, from the effects of the turbulent concentration process. The finite-particle-number effects in a specific application can always be added back into our model later (see, e.g., section IV.2).

For the following analysis, we will assume a shape for these PDFs. It has been suggested that, at least in the atmospheric context (Sreenivasan and Stolovitzky, 1995), symmetric beta-distributions provide a good approximation for multiplier distributions of dissipation. Such distributions have also been used in previous studies of particle concentrations (e.g., Hogan and Cuzzi, 2007) and are defined by

[TABLE]

with $\Gamma$ being the Gamma function. The parameter $\beta$ determines the width of the distribution with small values of $\beta$ corresponding to wide, i.e. more intermittent, distributions. The width of the $\beta$ distribution (its standard deviation) is given by

[TABLE]

However, similarly to the tracers, the observed distribution width, $\sigma_{0}(\beta)$, will not only depend on the underlying $\beta$ value but also on the number of particles in a given sample, $n$, and the number of samples, $N_{s}$, used to compute the distribution. In order to characterize this dependence, we have conducted Monte-Carlo experiments. They mimic the finite-particle-number effects in the DNS under the assumption that the underlying probability distributions are $\beta$ distributions. The procedure works in the following way: First, for a given value of $\beta$, we draw a random multiplier $m$ from the $\beta$ distribution. Given the number of particles $n$ in a sample volume (a given box), this corresponds to a partition into $n_{\textrm{left}}=mn$ and $n_{\textrm{right}}=(1-m)n$ particles for the two halves of the sample volume, where $n_{\textrm{left}}$ and $n_{\textrm{right}}$ are non-integers in general. Then, we draw a random particle number $k_{\textrm{left}}$ from the corresponding Poisson distribution with expectation value $E(k_{\textrm{left}})=n_{\textrm{left}}$. This value (and $k_{\textrm{right}}=n-k_{\textrm{left}}$) represent one random sample of the number of particles found in two halves of a box, and correspond to “observed” multiplier values $m_{\textrm{left}}=k_{\textrm{left}}/n$ and $m_{\textrm{right}}=k_{\textrm{right}}/n$. We repeat the procedure many times and compute the standard deviation $\sigma$ from all random samples $m_{\textrm{left}}$ and $m_{\textrm{right}}$ combined ($N_{s}$ total number of samples). The result is a random sample of the “observed” distribution width given $n$, $N_{s}$, and the true underlying value of $\beta$. Figure 4 shows the results from such experiments for selected parameters $n$, $N_{s}$. As one would expect, the scatter is large for a small number of samples, $N_{s}$, and becomes smaller with increasing $N_{s}$. Also, one can see that a small number of particles per box, $n$, causes the observed distribution width to be systematically larger than the value from equation (3.5) (red and green symbols), but approaches the exact value as the number of particles gets large (blue and orange symbols).

Now that we understand how we can model the small-number statistics effects, we can proceed to correct the observations. Given values of $n$ and $N_{s}$, we conduct the above described Monte-Carlo experiments for many random values of $\beta$, and find the value of $\beta$ that results in a distribution width closest to an observed width $\sigma_{0}$. This gives us an estimate of the true underlying $\beta$ value. In practice, we consider random values of $\beta$ between $1$ and $\infty$ by drawing samples with equal probability with respect to their true width $\sigma(\beta)$ (equation 3.5). Once we have accumulated at least $100$ samples resulting in distribution widths falling in a tolerance range between $\sigma_{0}(1+\Delta)^{-1}$ and $\sigma_{0}(1+\Delta)$, we compute the appropriate mean and standard deviation of those $\beta$ values. Initially, a large value of $\Delta$ is chosen, but as more and more samples are accumulated, $\Delta$ is reduced step by step by factors of $2$ until the standard deviation converges, i.e. does not change significantly if $\Delta$ is reduced further. We also require that the final $\Delta$ is no more than $0.01$ (the width is matched to within at least $1\%$ of the measured width). To speed up the process, once $\Delta\leq 0.1$, we restrict the range of $\beta$ samples drawn to values between $\bar{\beta}e^{-2\sqrt{\text{VAR}(\ln\bar{\beta})}}$ and $\bar{\beta}e^{2\sqrt{\text{VAR}(\ln\bar{\beta})}}$ where $\bar{\beta}$ and $\sqrt{\text{VAR}(\ln\beta)}$ are the mean value and standard deviation of $\beta$ in the tolerance range. Specifically, we define $\bar{\beta}$ as the $\beta$ value of the distribution whose variance (width squared) is the same as the arithmetic mean of the individual variances, i.e.

[TABLE]

with

[TABLE]

using equation (3.5) and its inverse $\beta(\sigma^{2})=(8\sigma^{2})^{-1}-\frac{1}{2}$, and where $\beta_{i}$ are the $\beta$ values of all the samples having widths within the tolerance range around $\sigma_{0}$. The uncertainty in $\bar{\beta}$ is measured by its variance:

[TABLE]

following non-linear uncertainty propagation truncating the series after the first order, and where the variance of $\bar{\sigma}^{2}$ is computed by:

[TABLE]

Once converged, the final values of $\bar{\beta}$ and $\text{VAR}(\ln\bar{\beta})$ are estimates of the true $\beta$ value corrected for finite-particle-number effects, and an estimate of its uncertainty.

III.2.2 Volume-averaged $\beta$-distributions

Now that we know how to remove the effects that a finite number of particles has on the multiplier distributions for a given box size or scale, we compute corrected volume averages of the $\beta$ values for any given spatial scale and Stokes number. Specifically, at any given scale $r$, we first compute the multiplier values for all the sampling boxes. We then compute the distribution widths $\sigma_{0,j}$ from each subset of boxes having the same number of particles $n_{j}$. This width is then used to derive a corrected $\beta$ value, which corresponds to a corrected width $\sigma_{j}$ (through equation 3.5). We then combine all results by summing the square of these corrected $\sigma_{j}$ values weighted with the fractional volume that the boxes with $n_{j}$ particles occupy. The final, combined $\beta$ value is then given by this average width through the inverse of (3.5).

Figure 5a shows the combined concentration multiplier $\beta$ values for all Stokes numbers and all spatial scales considered. We did not consider scales smaller than $12\eta$ since the number of particles in such small boxes is too small to make reliable statistical inferences, even with our correction procedure. From the results, it is very apparent that the concentration multiplier distributions are not only functions of $St$, as has been known, but are also very much scale dependent. It seems intuitive, however, that at any given scale, the multipliers may depend only on the ratio of the particle stopping time and the dynamical time at that spatial scale. An effect of this sort was seen in particle concentration PDFs by Bec et al. (2007). Along these lines, we construct a local Stokes number

[TABLE]

where we have assumed that the dynamical time at scale $r$ is given by

[TABLE]

following Kolmogorov (1962). By plotting the multiplier $\beta$ results against the rescaled Stokes number, the curves approximately collapse into one (figure 5b), at least for scales not too close to the integral scale. This means that when local scale and stopping time are accounted for, the scaled multiplier $\beta$ curves are cascade level independent for scales $r<L/10$ or so, and are thus highly amenable to cascade models at, in principle, arbitrarily large $Re$ at least within the inertial range. For some caveats about $Re$-dependence however, see section V.

Also note that $St_{r}$ can be written in terms of an integral-scale Stokes number $St_{L}$:

[TABLE]

where

[TABLE]

using the same Kolmogorov scaling as in equation (3.11). Equation (3.12) separates terms that depend only on cascade level (first term), and particle properties (second term), and indicates that particles in different flows behave the same (have the same statistical cascade) if they have the same integral-scale Stokes number $St_{L}$. We will use this fact later in section IV.2 when we compare the cascade with DNS results from two different simulations at different $Re$.

It should be noted our observed $r^{-2/3}$-scaling seems to contradict the scaling found by Bec et al. (2007) for “quasi-Lagrangian” probability distribution functions of the mass density. Following an idea by Maxey (1987), they approximated the dynamics of inertial particle by those of tracers in an appropriate synthetic compressible velocity field and derived a scaling for the rate at which an $r$-sized “blob” of particles contracts. They argued that the scaling of the contraction rate relates to the scaling of the pressure field, give the contraction rate as being proportional to $r^{-5/3}$ for the $Re$ of their simulation 222They argue that the scaling should change to $r^{-4/3}$ for very high Reynolds numbers $Re_{\lambda}\geq 600$, and find that their density PDFs collapse when scaled with this contraction rate.

III.2.3 Composite PDFs

Combining the statistics from boxes with different particle number to form a single $\beta$-distribution that represents the concentration multipliers assumes that multiplier distributions are the same for all concentrations. Since the particles are independent, that is, they do not feel the presence of each other, one would assume that there is no such concentration dependence. However, particles concentrate in particular regions of the flow, and therefore flow properties differ in regions of different particle density and so the multiplier distributions may also be different. We can relax the assumption of equal multiplier distributions and compute a composite multiplier distribution by first computing $\beta$-distributions for each particle concentration separately, and then summing these distributions weighted with the fractional volume that the boxes with the particular density occupy. These composite distributions are shown in figures 6, 7 and 8.

The following observations can be made: First, composite and mean $\beta$-distributions may differ in shape, although by construction they have identical widths (second moments). That is, in general, the composite $\beta$-distribution is itself not a $\beta$-distribution. At large spatial scales, flat or exponential tails are apparent in the composite PDFs; these do not have the same shape as any $\beta$-distribution. At small scales, however, differences disappear within the margin of accuracy. Also, at the largest scale ($512\eta$) there are enough particles such that finite-particle-number effects are small and the raw PDFs (shown in green) are essentially identical to the composite PDFs (shown in black). However, the importance of correcting for sampling effects becomes apparent as we look at smaller scales where fewer particles cause spurious widening of the raw distributions or “false intermittency.” Finally, we should mention that at the smallest scales, the number of particles is so small that we are only measuring multipliers in high-concentration regions, which causes a sampling bias since particles are known to avoid vorticity, and such regions may produce more intermittency or broader multiplier PDFs. For instance, the average number of particles in a $32\eta$ box, for the Stokes numbers for which we have a total of $6\times 10^{5}$ particles (see section II) is only $\left(32\eta/2048\eta\right)^{3}N_{p}\approx 0.44$. The situation is even worse for, say, $r=12\eta$ and a case with only $2\times 10^{5}$ total particles. The average is then only 0.04 particles per sampling box. Multipliers are measured only in regions with a particle concentration that is at least 100 times larger than the mean concentrations since we can only reasonably measure multipliers if we have at least several particles in a sampling box.

III.3 Dissipation and enstrophy multiplier distributions

We also calculated the multipliers for fluid dissipation and enstrophy. The rate of turbulent dissipation is given by

[TABLE]

where $S_{ij}$ and $u_{i}$ are the strain rate tensor and the components of the velocity field, respectively, and where we use the Einstein summation convention. The enstrophy on the other hand is defined as the square of the vorticity:

[TABLE]

All of these flow velocity derivatives are available in the dataset, both in the particle trajectory data files and in the flow snapshots. We have computed multiplier distribution $\beta$ values for $\epsilon$ and $\cal E$ ($\beta_{\epsilon}$ and $\beta_{\cal E}$) from the trajectory of tracer particles (as they sample the flow more homogeneously than non-zero Stokes number particles), and from the full resolution flow snapshots. The results are shown in figure 9, although the tracer data is only shown for the largest spatial scales since it also suffers from finite-particle effects (not corrected here).

Note the presence of asymptotes for $r\lesssim 20\eta$ (perhaps better thought of as $r\lesssim L/50$, see figure 14) for both, as anticipated (Sreenivasan and Stolovitzky, 1995). Enstrophy is shown to have wider multiplier distributions (smaller $\beta$ values) than dissipation, and is therefore more intermittent. This is consistent with the findings of Meneveau et al. (1990) in several flows including atmospheric flow, and in numerical simulations (e.g., Cao et al., 1996; Chen et al., 1997). For a review, see Sreenivasan and Antonia (1997). Also, we note that even for the smallest spatial scales considered, still well within the inertial range, the dissipation rate multiplier $\beta$ does not reach the atmospheric flow values of $\beta_{\epsilon}\sim 3$ (Sreenivasan and Stolovitzky, 1995). See section V for more discussion.

IV New cascade model with level-dependent multipliers

IV.1 Cascade simulations

From the collapsed $\beta(St_{r})$ curves (figure 5b), we can build an empirical model for the particle multiplier distributions. A sum of two power laws approximates the curves for fixed scale $r$ well:

[TABLE]

with parameters $a_{1}$, $a_{2}$, $b_{1}$, $b_{2}$ determining the slopes and positions of the exponentials, respectively, and $\beta_{\textrm{min}}$ setting the minimum $\beta$ value. From the figure it is evident that there is some residual scale dependence – the curves for different spatial scales don’t overlap exactly. The following parameterization approximates this residual dependence:

[TABLE]

where $\operatorname{erf}$ and $\ln$ are the error function and the natural logarithm. The parameters asymptote for small $r/L\lesssim 3\times 10^{-3}$ (or large cascade level $N\gtrsim 25$) to:

[TABLE]

The model is shown in figure 10 compared to the DNS results. An even simpler model could probably be constructed using a single average curve of all $\frac{r}{L}\leq\frac{1}{8}$.

With this model for the particle multiplier distributions, we can perform statistical cascade simulations to predict the probability distribution function for the particle concentration. We start at cascade level 0 with a single concentration value of $C^{(0)}=1.0$. At every cascade level $N$, we then draw a random multiplier value $m^{(N)}$ from the corresponding distribution with a $\beta$ value given by the model (equation 4.2), and split the concentration value from the previous level into two values $C^{(N)}_{\textrm{left}}=2C^{(N-1)}m^{(N)}$ and $C^{(N)}_{\textrm{right}}=2C^{(N-1)}(1-m^{(N)})$. The factor 2 here comes from the fact that the concentration is the ratio of the particle number in a half box and the mean number in a half box (which itself is one half of the mean number of particles in full box). Such a cascade produces $2^{N}$ random samples of concentrations values at each cascade level $N$ which we use to compute concentration PDFs. For good statistics, however, we need many more samples, and for the predictions shown in the following section we computed 50,000 such cascade simulations.

IV.2 Level-dependent cascade predictions compared to DNS at two different $Re$

In order to demonstrate and assess the cascade predictions, we compare the probability distribution functions of the concentration factor generated by the cascade model with those measured directly from DNS datasets. In order to do so, we need to account for the small-number effects present in the DNS results which, as we have seen before, can cause observed distributions to be significantly widened relative to ideal ones that the cascade produces. Instead of correcting the DNS PDFs as we have done before for the measured multiplier distributions, we will degrade the cascade PDFs for this purpose by introducing finite-particle-number effects into them.

For motivating the procedure, let us imagine a hypothetical simulation $\cal H$ with a number of particles so large that finite-particle-number effects are negligible, and let $\bar{n}_{\infty}(r)$ be the average number of particles in a sampling box at length scale $r$ in that simulation. Our present dataset, in comparison, has on average only $\bar{n}(r)=N_{p}r^{3}/{\cal L}^{3}$ particles in a box of scale $r$, where as before $N_{p}$ is the total number of particles in the dataset with a given Stokes number, and $\cal{L}$ is the linear extent of the simulation domain, respectively. One can think of our current dataset as a randomly selected subset of the hypothetical simulation $\cal H$, generated by retaining particles from $\cal H$ with a probability of $p(r)=\bar{n}(r)/\bar{n}_{\infty}(r)$. For brevity, we will suppress the $r$ below. Specifically, let’s say some sampling box in $\cal H$ has $n_{\infty}$ particles in it (and therefore a concentration factor $C_{\infty}$ = $n_{\infty}/\bar{n}_{\infty}$). From these, we select particles with a probability of $p$, retaining in total $n$ particles, where $n$ is an integer random number with an expectation value of $E(n)=C_{\infty}\bar{n}$. For $\bar{n}_{\infty}\rightarrow\infty$, this is a Poisson process and $n$ is a random number with a probability mass function

[TABLE]

Using this idea, the recipe for introducing finite-particle-number effects into the cascade PDFs is as follows: First, we draw a random sample $C_{\infty}$ from a cascade-derived concentration PDF. Second, we draw a random sample $n$ from a Poisson distribution with the corresponding expectation value $E(n)=C_{\infty}\bar{n}$, where $\bar{n}$ is again the average number of particles at the spatial scale of interest in the DNS data we want to compare. The (integer) particle number $n$ corresponds to a discrete concentration factor $C=n/\bar{n}$. By repeating the procedure $N_{s}$ times, we can build from the samples a discrete probability distribution function of $C$. It accounts for the finite-particle-number effects and can be directly compared to PDFs measured from the DNS dataset.

Figure 11 shows, for different Stokes numbers, a comparison between the PDFs predicted by the cascade model, and degraded in this way using $N_{s}=10^{7}$, with those calculated directly from the DNS dataset we analyzed in this paper (Lanotte et al., 2011).

Under the assumption that the collapsed multiplier $\beta$ curve is universal and does not depend on Reynolds number, we can use cascade simulation to model conditions at different Reynolds numbers. For caveats to this and a suggestion regarding plausible $Re$ dependence, see the discussion in section V. It follows from equation (3.12) that particles of different $Re$ flows behave the same and have the same cascade statistics, if they have identical integral-scale Stokes numbers $St_{L}$.

Here, we compare our cascade with Pan et al. (2011) who performed direct numerical simulation of a compressible flow with suspended initial particles. Their simulation is on a $512^{3}$ Cartesian grid with an estimated $L/\eta\sim 200$ (compared to $L/\eta\sim 1024$ in the dataset we used here), and contains $8.6\times 10^{6}$ particles per Stokes number. Initial comparison of their uncorrected multiplier PDFs with ours showed a clear disagreement but most of that disagreement disappears once we take into account finite-particle number effects. Figure 12 shows a comparison of cumulative PDFs of the concentration between Pan et al. (2011) and our cascade model showing reasonable agreement, bearing in mind that the inviscid simulations of Pan et al. (2011) leave a little uncertainty about the value of $St$.

IV.3 Model prediction for the radial distribution function

In many applications, e.g. in terrestrial clouds, particle collisions play an important role, and it is therefore of great interest to model this process. The rate of collisions depends on two statistical quantities: the radial relative velocity between particles, and the radial distribution function (RDF), $g(r)$, defined as the probability of finding two particles at a given separation normalized with respect to homogeneously distributed particles (Pan et al., 2011; Bragg and Collins, 2014a). Relative velocities are beyond the scope of the present model, but the cascade model can be used to make predictions for the RDF. For this purpose, we performed statistical simulations similar to section IV.1 but with an important difference: we are here interested in the spatial distribution of the concentration, which then can be used to compute the RDF.

Our cascade model, however, only describes the particle multiplier PDF as a function of cascade level, and does not explicitly contain information about spatial correlations. There is therefore some ambiguity in how to compute concentration fields from the cascade model. We have explored three different methods to assess the range of possible solutions. Starting at the largest scale, we divide a cube of space in half along each spatial direction. This results in eight sub-cubes with half the linear size. In order to solve for the concentrations in these sub-cubes uniquely, we need eight equations. The first method – method A – makes the following choice: one constraint it given by the fact that the average of the concentrations over all eight sub-cubes is equal to the mean concentration, and seven additional constraints are given by relating the concentration in seven sets of neighboring sub-cubes through multipliers chosen randomly from the cascade model. The specific choice of equations is given in appendix A. Solving for concentrations to ever smaller cubes until some small cutoff length scale $(r/L)_{\textrm{min}}$ then yields a statistical realization of the concentration field that can be used to compute the RDF. The probability of finding two particles at a given distance is simply the product of the concentrations at the two points in space a given distance apart, averaged over the whole domain. If we start with a unit concentration at the largest scale, the normalization, that is the probability for homogeneously distributed particles, is simply 1. In order to reduce the computational cost and storage requirements to trackable amounts, we do not follow all sub-cubes to ever smaller scale but only a random selection of them. One half of sub-cubes are followed at each cascade level. Two more methods are obtained by relating the concentrations in the two half-cubes, for each direction separately, through a random multiplier. This set of equations is underdetermined. For method B we pick one particular solution, while for method C we use a least-square solver to determine the minimum-norm solution. For specifics, again, we refer to appendix A. Conceptually, it is clear that the three methods allow for different amounts of spatial randomness. Method A clearly maximizes intermittency while method C leads to the least spatially intermittent solution.

Figure 13 shows RDFs predicted by our cascade model simulations down to cascade level 72 (spatial scales of $(r/L)_{\textrm{min}}\approx 6\times 10^{-8}$). All methods give qualitatively the same results, all in good agreement with the Zaichik and Alipchenkov (2009) asymptotic $Re=\infty$ analytical solution, especially regarding the shape and the active range of scales. In fact, the magnitudes are even close enough, within a factor of order unity. RDFs for different $St_{L}$ computed using method A are plotted in figure 13a as a function of the scaled distance $St_{L}^{-3/2}r/L$, and are shown to collapse for small $St_{L}\lesssim 10^{-2}$. Effectively, the scaling behavior of the multiplier PDFs (figure 5) is carried over to the RDF. At large scales, $St_{L}^{-3/2}r/L\gtrsim 10^{2}$, the RDF has a value of 1, that is, particles are homogeneously distributed, and over an active range of scales, $10^{-2}\gtrsim St_{L}^{-3/2}r/L\gtrsim 10^{2}$, it rises and reaches an asymptote, $g_{0}$. For larger $St_{L}$ in our sample, however, the active range is shortened by being too close to the integral scale $L$, and $g(r)$ asymptotes at smaller values that are $St_{L}$-dependent. The theory of Zaichik and Alipchenkov (2009) predicts a very similar behavior and their curve for infinite Reynolds number, effectively for infinitely small $St_{L}$, is also shown in figure 13.

In figure 13b, we compare RDFs from the different methods and the Zaichik and Alipchenkov (2009) model by scaling them with their respective $g_{0}$. The curves are nearly identical. Method A produces the highest RDF values, i.e. the most intermittent concentration distributions, and method C, as it is biased towards lower intermittency, results in the smallest asymptotic value (see table 1) for values of $g_{0}$. Interestingly, $g_{0}$ for method C is the same as for the Zaichik and Alipchenkov (2009) theory. Zaichik and Alipchenkov (2009) use various approximations in their derivation. Among them, they model the turbulence by a Gaussian process which would underestimate the tails of their probability distributions and therefore underestimate intermittency.

V Discussion

V.1 The nature of scale-dependent particle concentration

As described in section I, the traditional explanation of clustering in terms of centrifugation of particles from eddies has been replaced with a somewhat more complicated and nuanced combination of physical processes (Bragg and Collins, 2014a, b). We believe that our results (figures 5 (right) and 10), stripped of their minor variations, represent a kind of universal curve (“U-curve”) for $\beta(St_{r})$ that can be interpreted in terms of these different processes acting on a particle of some $St$ over a range of scales $r$.

In the $St_{r}<0.1$ regime, the effect is dominated by centrifugation, which weakens as $St_{r}$ decreases (Chun et al., 2005, and others); thus $\beta$ increases (the multiplier PDF narrows) with decreasing $St_{r}$ in this regime. As $St_{r}$ increases beyond 1, the concentration effect is weakened by the decreasing sensitivity of particles to perturbations of any kind by eddies with timescales much shorter than their stopping times, so $\beta$ again increases. This effect, which could be thought of as an inertia impedance mismatch, has been described in terms of response functions (Meek and Jones, 1973; Völk et al., 1980; Markiewicz et al., 1991; Cuzzi and Hogan, 2003; Ormel and Cuzzi, 2007), but see Bec et al. (2010b); Pan and Padoan (2013); Hubbard (2012) and Hopkins (2014) for other more recent and more sophisticated analyses.

The strongest clustering effect is produced (the multiplier PDF has the lowest $\beta$) across perhaps one or two decades of eddy scale $r$ for a given $St$, centered on the combined parameter $St_{r}\sim 0.3$, presumably the regime where history effects in particle velocities play the dominant role. While our results support the idea that concentration is generically due to “eddies on the scale of $\eta St_{\eta}^{3/2}$..” (Ireland et al., 2015; Bec et al., 2007), we think a more refined description one could infer from the U-curve is that clustering is the cumulative result of a history of interactions with the flow of energy as it cascades over eddies ranging over two decades in size, driving particles ever deeper into a concentration “attractor” even in the inertial range (Toschi and Bodenschatz, 2009; Bragg and Collins, 2014a, b). The other “universal curve” of Zaichik and Alipchenkov (2003) (their figure 1) and their improved model (Zaichik and Alipchenkov, 2009, their figure 3) reproduced in figure 13 also has this sense. That is, we see a parallel based on causality, between time-asymmetrical “history effects” on particles of some $St$ as they are affected by energy flowing down the cascade through eddies of different scale, and a trajectory down one side and then back up the other side of our U-curve. Such a picture would lead the particle concentration as a function of spatial binning scale to increase sharply over some particular range of scales $r/L$ or $r/\eta$ related only to $St$ (the active range), and then remain constant towards smaller scales where eddy perturbations are felt only weakly because of, essentially, the poor impedance match with the particle stopping time.

A natural prediction of this model is thus that, at infinite $Re$ where energy is available on all timescales, far from the dissipation range, and in the absence of complications such as gravitational settling or fragmentation limits, the maximum particle concentration should not only arise over a similar range of scales near $r\sim\eta St_{\eta}^{3/2}$ (Zaichik and Alipchenkov, 2003, 2009; Bec et al., 2007; Ireland et al., 2015), but also should have a “saturation” amplitude that is $St$-independent. Indeed this would seem to be the prediction of Zaichik and Alipchenkov (2009, their figure 3). In current simulations (Ireland et al., 2015; Zaichik and Alipchenkov, 2009), as well in simulations we have conducted using the cascade, the clustering of larger St particles ($St_{L}\gtrsim 0.03$) asymptotes at smaller values of the RDF than seen for smaller particles (figure 13). We expect this is because the scale at which the larger particles reach $St_{r}\sim 0.3$ is too close to the integral scale, so their potentially two-decade-wide range of interest, that our U-curve indicates is needed to reach a true asymptote, is truncated at large scales.

Within the dissipation range at $r<20-30\eta$, the energy spectrum of the flow changes as a result of the now-fixed eddy timescale $t_{r}=t_{\eta}$ (Bragg et al., 2015b). Particles of $St_{\eta}\sim 1$ are now unique in that they do not experience the usual impedance mismatch with faster eddy forcing going to smaller spatial scales, so can continue to increase in concentration going to smaller scales. As noted by Bragg et al. (2015b) the question remains as to whether there is any sort of rollover at $r\ll\eta$ in the RDF of $St_{\eta}$=1 particles, or whether an actual singularity would exist for point particles. In terrestrial applications, finite particle sizes comparable to $\eta$ preclude unlimited singular behavior; however, in protoplanetary nebula applications (Pan and Padoan, 2010, 2013; Johansen et al., 2015), particle sizes of interest (submillimeter to dm) are orders of magnitude smaller than the Kolmogorov scale (km) so this is a question of significant interest.

Bec et al. (2007) explicitly described the dissipation regime as characterized by an “attractor” having fractal properties, and indeed multifractal properties were demonstrated by Hogan et al (1999) for clustering in this regime. It is known that cascades lead to fractal and multifractal spatial distributions (Meneveau et al., 1990; Sreenivasan and Stolovitzky, 1995, and references therein). We now suspect that the cascade of Hogan and Cuzzi (2007), in which the multiplier distributions do seem to obey level-independent scaling, were effectively dissipation range cascades. Tests by Hogan and Cuzzi (2007) showed good agreement between their level-independent cascade model and DNS. However, the multiplier PDFs were determined at $3\eta$ and all their DNS results were for low-to-moderate $Re_{\lambda}<140$ such that the integral scales were 14$\eta$, 24$\eta$, 45$\eta$, and 86$\eta$. At least the first three of these runs lie mostly within the dissipation range, where scaling does support level independence (Bec et al., 2007; Bragg et al., 2015b). It might be worth exploring the use of dissipation range cascades further from the standpoint of modeling fractal structure or to study higher moments of the particle density PDF. Indeed Bec and Chétrite (2007) present what is, essentially, a cascade model that reproduces aspects of the particle concentration PDF.

Moreover, to our knowledge, while fractal/multifractal behavior has been shown for particle clustering within the dissipation regime, either at scales of a few to tens of $\eta$ (Hogan et al., 1999), or scales smaller than $\eta$ (Bec, 2003; Bec et al., 2006), no explicit study of this property in the inertial range has been done. It would be of interest to find whether the inertial range cascade as described by the U-function (figure 5), which is level-dependent but in a predictable way that is level-independent, would also produce such a distribution, when suitably scaled for St. This could be of use in modeling radiative transfer properties (Shaw et al., 2002).

V.2 Dissipation, enstrophy, and $Re$-dependence

As mentioned earlier in section III.3, our dissipation multiplier PDFs have larger $\beta$ (are narrower) than the expected $\beta_{\epsilon}\sim 3$ for all binning sizes we could usefully study, reaching an asymptote of $\beta_{\epsilon}\sim 8$ at $r\lesssim L/86$. Based on the very extensive inertial range manifested in figure 1 (Bec et al., 2010b), apparently extending up to $>2000\eta$, we had expected to find scale-free behavior in the dissipation multiplier PDF over most of this range. However, recalling figure 2, especially as selected by larger $|q|$, which more strongly weight the structures where most dissipation occurs, the properties of dissipation are not invariant over as wide a range as is the second order velocity structure function that defines the inertial range in Bec et al. (2010b).

Figure 14 summarizes the overall scale variation of $\beta$ for dissipation and enstrophy, showing the spatial scale both in terms of $\eta$ and $L$ following figure 1. At large scales, the PDF is narrow for both (large $\beta$) but widens with decreasing scale. At scales of $\sim 12\eta$ ($\sim L/86$) it asymptotes at a value which seems to remain scale-free to smaller scales. Our new asymptotic values do not agree with values ($\beta_{\cal E}\sim 10$ for enstrophy) found in Hogan and Cuzzi (2007), or ($\beta_{\epsilon}\sim 3$ for dissipation) in Sreenivasan and Stolovitzky (1995). It may be that the Hogan and Cuzzi (2007) $\beta_{\cal E}$ is more properly associated with the dissipation range, but the discrepancy in $\beta_{\epsilon}$ alone merits some discussion.

We hypothesize that at much higher $Re$ than we can study here, the scale-dependence of $\beta_{\epsilon}$ morphs in a fashion so as to be consistent with Chhabra et al. (1989) and Sreenivasan and Stolovitzky (1995); that is, has a scale-free $\beta_{\epsilon}\sim 3$ for all scales less than at least 3000$\eta$ (based on Chhabra et al. (1989)) and probably less than $L/16$ (based on Chhabra et al. (1989) and Chhabra and Jensen (1989)). At larger scales we expect $\beta$ must increase in some smooth fashion similar to ours, with an overall behavior schematically shown by the red dotted line in figure 14. Meanwhile, by the logic that enstrophy $\cal E$ is always more intermittent (has smaller $\beta$) than dissipation, we then also hypothesize that $\beta_{\cal E}$ varies as suggested by the blue dotted line.

We suspect that our observed $\beta_{\epsilon}$ and $\beta_{\cal E}$ asymptote (for $r\lesssim L/86$) at larger values than would be true for much higher $Re$, because the viscous or dissipation range, which bounds the inertial range on its small-scale end and extends to 20-30$\eta$ in general, here impinges on the small-scale end of the nascent inertial range, and may prevent the dissipation and enstrophy from ever fully realizing their high-$Re$ intermittency. In contrast, at high atmospheric $Re$, the large-scale onset of the inertial range, at $3000-10^{4}\eta$ based on Chhabra et al. (1989) and Sreenivasan and Stolovitzky (1995), is completely isolated from the viscous range, as indicated in figure 1 by the blue axis labels on the upper horizontal axis.

The scale-dependence and asymptotic value of $\beta_{\cal E}$ is important, because the process of particle concentration may track the properties of $\cal E$ rather than those of $\epsilon$ based on the physics involved (section V.1). While it may be coincidental, in our DNS results, the $\beta$ for inertial particles minimizes at a value very close to our value for $\beta_{\cal E}$, and considerably smaller (more intermittent) than our value of $\beta_{\epsilon}$. A secondary, related hypothesis is that the particle concentration multipliers may track the behavior of enstrophy (if velocities and accelerations are dominated by vorticity), and the minima seen in the collapsed curves of figure 5, which now never fall below 3.0, might drop to significantly lower values, making the particle concentration field more intermittent. For this reason, cascades developed using our current collapsed $\beta(St_{r})$ curves may underestimate the abundance of zones of high concentration at high $Re$ to some degree. A better understanding of this $Re$-dependence will be needed to put cascade modeling of particle concentration on quantitatively solid ground. The original dataset of Sreenivasan and Stolovitzky (1995) probably contains enough information to assess the validity of these hypotheses regarding dissipation, but not for enstrophy.

V.3 Speculations regarding the effect of higher $Re$-numbers

As described in section IV, a level-dependent cascade can be described which captures the inertial range behavior for arbitrary $Re$ (at least to the point where mass loading starts to affect the physics, e.g., see Hogan and Cuzzi (2007)). This might be of use in modeling particle concentration in rain clouds, in the protoplanetary nebula, or in other applications at very high $Re$ where numerical simulations are impractical (such as spatial variations in microwave opacity). There are several reasons to expect different behavior in the dissipation range.

Recalling, however, our numerical discrepancy with Sreenivasan and Stolovitzky (1995) and others regarding the asymptotic value of $\beta$ for dissipation, we think the possibility of $Re$-dependence of even our collapsed $\beta(St_{r})$ might imply that our results (and other inertial range results at low $Re$) underestimate effects of concentration, in the sense that at higher $Re$, the minimum in the universal curve would move to lower $\beta$ (more intermittency and higher concentrations). One might speculate that the minimum $\beta$ for particles should track the $\beta$ for enstrophy instead of for dissipation.

In the application to planetesimal formation proposed by Cuzzi et al. (2010), it is necessary to create a joint PDF of particle concentration and enstrophy. This motivates a better understanding of high $Re$ behavior of $\beta$ for both particles and enstrophy.

VI Conclusions

Our results indicate that the multiplier PDFs for particle concentration (section III.2), dissipation, and enstrophy (section III.3) vary with scale, at least over the largest decade of spatial scales. We also find that the multiplier PDFs for particle concentration have two components: a traditional “$\beta$-function” component, and an exponential-tail component (section III.2.3). We find that the concentration multiplier $\beta$ values collapse to a scale-independent universal curve when plotted against an appropriately scaled local Stokes number $St_{r}=St(r/\eta)^{-2/3}$ (section III.2.2, in particular equation (3.10) and figure 5), allowing the cascade model to be used for modeling higher $Re$ conditions not accessible to numerical simulation.

For dissipation, the “$\beta\sim 3$” asymptotic behavior of Sreenivasan and Stolovitzky (1995) in high $Re$ atmospheric flows probably appears at around $r\sim L/30$ or $L/40$, and remains constant to smaller scales, at least to $r\sim 20-30\eta$ where the dissipation range begins. In the present simulation, the integral scale and the dissipation range are not separated far enough for the dissipation $\beta_{\epsilon}$ to reach such low values, and instead it asymptotes for scales below $r\sim 20\eta$ to a value of $\beta_{\epsilon}\sim 8$. Enstrophy, believed to be always more intermittent (smaller $\beta$) than dissipation, asymptotes in the DNS to a value of $\beta_{\cal E}\sim 3$. In light of the connection between vorticity and the acceleration, centrifuging and concentrating of particles, it may not be surprising that this value coincides with the minimum $\beta$ value for particle concentration multipliers at the optimal $St_{r}$.

Given that dissipation, and presumably enstrophy, have not reached their scale-independent, asymptotic values seen in very high $Re$ atmospheric flows, it could be expected that our collapsed particle multiplier $\beta(St_{r})$ curve is also $Re$ dependent. Analyses of this sort for DNS of particle-laden flows at significantly higher Reynolds number are therefore highly desirable, as are measurements of enstrophy multipliers in very high Reynolds number flows such as atmospheric flows.

We have also found that the cascade model can be used to construct a spatial distribution of particle concentration, that can be carried to arbitrarily high Reynolds numbers, and has a very good resemblance to the analytical theory of Zaichik and Alipchenkov (2009). More work is needed to assess the asymptotic level of maximum concentration for particles of any size (which we find, as did Zaichik and Alipchenkov (2009), is size invariant, in the infinite $Re$ limit).

Acknowledgements: The authors would like to thank Nic Brummell, Karim Shariff, Katepalli Sreenivasan, Federico Toschi and Alan Wray for insightful discussions on the topic of this paper, and Federico Toschi and Enrico Calzavarini for help accessing and using the DNS data.

We would like to dedicate this paper to Bob Hogan, who passed away in February 2012. Bob was responsible for all the previous computational work done on turbulent concentration by our group from 1992-2010, including multifractal behavior and our first cascades. He would be very interested in these new results, which explain a discrepancy that arose too late for him to resolve.

Appendix A Three methods for computing the spatial distribution of concentrations

We here provide more details about the three methods for computing spatial concentration fields from multiplier distributions that we used in section IV.3 for computing radial distribution functions.

Let us consider a cube of size $r/L=2^{-N/3}$ (cascade level $N$) having a concentration $C^{(N)}$, initially we start with a cascade level 0 cube having a unit concentration $C^{(0)}=1$, and divide it in half along each spatial direction. This results in eight sub-cubes with half the linear size which correspond to a cascade level of $N+3$ (equation 1.1). We denote the concentrations in these sub-cubes as $C_{ijk}^{(N+3)}$, where $i,j,k\in[1,2]$ are indices denoting the left (1) or right (2) sub-cube in the three spatial directions. Since there are eight unknowns, we need eight equations to uniquely determine the concentrations.

The first method – method A – makes the following choice: One constraint is given by the fact that the average of the concentrations over all eight sub-cubes is equal to the mean concentration, that is

[TABLE]

We get seven more constraints by relating the concentration of neighboring sub-cubes to multipliers that are chosen randomly from the cascade model. A possible choice are the combinations

[TABLE]

where $m^{(N+3)}_{i}$ with $i\in[1,...,7]$ are seven random multiplier values at cascade level $N+3$. Equations (A.1) and (A.2) are linearly independent and can be solved directly. Applying this procedure recursively to ever smaller cubes until some small cutoff length scale $(r/L)_{\textrm{min}}$ yields one statistical realization of the concentration field.

The two other methods are obtained by relating the concentrations in the two half-cubes, for each direction separately, through a random multiplier, that is

[TABLE]

where $m^{(N+1)}_{i}$ with $i\in[1,...,3]$ are multipliers randomly drawn from the level $N+1$ cascade model. Since these are only three equations, the linear system is underdetermined.

For method B we choose one particular solution to equations (A.3), namely

[TABLE]

There is an ambiguity, of course, as to whether multiplier $m^{(N+1)}_{i}$ is applied to the “left” or “right” half of each box, but since the multipliers $m^{(N+1)}_{i}$ and $(1-m^{(N+1)}_{i})$ have equal probability, either choice would be paired with the opposite given enough random samples (and we do average many random samples to construct each cascade). The eight sub-cubes from each large box have the identical concentrations, just differently distributed, for each set of $m^{(N+1)}_{i},i\in[1,...,3]$ whether $m^{(N+1)}_{i}$ goes to the left or right box in each case.

For the third and final method, method C, we use a least-square solver to determine the minimum-norm solution of (A.3), that is the solution that is closest to equally distributed concentrations. Clearly, this causes a bias towards the least intermittent spatial distribution. For strongly intermittent multipliers, this method can even lead to negative concentrations in one of the sub-cubes. Such solutions have to be discarded, and this further biases method C towards minimal intermittency.