Scale Dependence of Multiplier Distributions for Particle Concentration, Enstrophy and Dissipation in the Inertial Range of Homogeneous Turbulence
SScale Dependence of Multiplier Distributions for Particle Concentration, Enstrophyand Dissipation in the Inertial Range of Homogeneous Turbulence
Thomas Hartlep
Bay Area Environmental Research Institute, Petaluma, CA 94952, USA ∗ Jeffrey N. Cuzzi
NASA Ames Research Center, Moffett Field, CA 94035, USA
Brian Weston
University of California Davis, Mechanical & Aerospace Engineering, Davis, CA 95616, USA (Dated: 12 December 2016 (Submitted to
Physical Review E ), 28 February 2017 (Accepted))Turbulent flows preferentially concentrate inertial particles depending on their stopping timeor 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 directnumerical simulation (DNS) dataset of homogeneous, isotropic turbulence to determine probabilitydistribution functions (PDFs) for cascade multipliers , which determine the ratio by which a propertyis partitioned into sub-volumes as an eddy is envisioned to decay into smaller eddies. We present atechnique for correcting effects of small particle numbers in the statistics. We determine multiplierPDFs for particle number, flow dissipation, and enstrophy, all of which are shown to be scaledependent. However, the particle multiplier PDFs collapse when scaled with an appropriatelydefined local
Stokes number. As anticipated from earlier works, dissipation and enstrophy multiplierPDFs reach an asymptote for sufficiently small spatial scales. From the DNS measurements, wederive 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 & Alipchenkov [New Journal of Physics 11, 103018 (2009)]. Wediscuss implications of these results for the statistical modeling of the turbulent clustering processin the inertial range for high Reynolds numbers inaccessible to numerical simulations. I. BACKGROUND AND INTRODUCTION
Clustering of inertial (finite-stopping-time) particlesinto dense zones in fluid turbulence has applications inmany fields [for a general review see 1]. A number ofrecent papers have focussed on understanding the basicmechanisms responsible for this effect; several of these [2–5] provide thorough reviews and comparisons of previousstudies dating back to the early work of Maxey [6] andSquires and Eaton [7] which we will only sketch briefly.The early work emphasized the role of centrifugation offinite-inertia particles out of vortical structures in turbu-lence. More recent evidence that clustering arises evenin random, irrotational flows suggests that, while vortic-ity still plays a role, the dominant role is played by so-called “history effects”, in which inertial particle velocitydispersions at any location carry a memory of particleencounters with more remote flow regimes which havelarger characteristic velocity differences [8–10]. Thesehistory effects lead to spatial gradients in particle randomrelative velocities, and these gradients in turn generatesystematic flows or currents which can outweigh disper-sive effects and produce zones of highly variable particleconcentration [2–4, 11]. ∗ Also, NASA Ames Research Center, Moffett Field, CA 94035,USA; [email protected]
To date, by far the most attention regarding particleclustering in turbulence has been devoted to very smallspatial scales r < η or even r (cid:28) η , where η is theKolmogorov scale, partly because it is on these scalesthat particle collisions occur and partly because numeri-cal simulations to date have produced only very limitedinertial ranges, at best [see however 12, 13]. Theoriesby Zaichik and Alipchenkov [11] et seq. , and Pan andPadoan [8] et seq. have been shown to be promising inexplaining the cause of particle clustering in terms ofhistory effects, with helpful contributions from the tra-ditional local centrifugation mechanism [2, 4, 10, 13]. Athorough review of the effects of clustering and relativevelocity effects on particle collisions, emphasizing the as-tronomical literature, can be found in Pan and Padoan[14, 15].Our focus is on clustering at larger scales in the in-ertial range η < r < L , where L is the integral scale.Inertial range clustering has important applications forremote sensing of terrestrial clouds [16], the formationof primitive planetesimals (asteroids and comets) in theearly solar nebula [17–22], and even the structure of theinterstellar medium [23]. While little studied in the con-text of particle clustering, inertial range scaling is knownto have different properties than seen in the dissipationrange r < η [24, 25]. Only limited predictions have beenmade of its scaling properties at very high Reynolds num-ber Re [2, 11]. a r X i v : . [ phy s i c s . f l u - dyn ] M a r In the inertial range, so-called cascade models whichreproduce the statistics of fluid behavior, even if not re-alistic flow structures, may be valuable for modeling highReynolds number ( Re ) regimes too demanding for directnumerical simulations. Their application is quite gen-eral [26–31, see [32] for more references]. We and othershave used cascades to model particle clustering in turbu-lence in the astronomical applications mentioned above.Cascade models operate by simply applying a parti-tion function or multiplier ≤ m ≤ P in some given volume of the flow, thus determining theratio by which the property (dissipation, particle density,etc.) is partitioned into sub-volumes as an eddy is envi-sioned to decay into smaller eddies. The most commontreatment is a binary cascade, in which P is partitionedinto two equal subvolumes; however the approach can beapplied to arbitrary numbers of subvolumes [30]. Thebinary cascade operates on each volume of space, parti-tioning P into two equal subvolumes by multipliers m and1 − m , with the multiplier m at each bifurcation drawnfrom a probability distribution function (PDF) of multi-pliers P ( m ). If P ( m ) = δ ( m − . δ is the deltafunction, the cascade has no effect because the property P is evenly divided, and remains constant per unit vol-ume. On the other hand, broad P ( m ) functions generatehighly intermittent spatial distributions in which P hasa wide range of values, fluctuating dramatically on smallscales such as seen in dissipation [27, 31] (figure 1b).The dissipation range , a range of small scales ap-proaching the Kolmogorov scale η , is found where r < − η [33, 34]; in this range, where viscosity is impor-tant, the equations of motion are no longer completelyscale-free, and fluid scaling properties differ from thosein the inertial range. The properties of particle clus-tering do seem to be scale independent in this region,however [24, 25], 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 fromscale invariance surely occurs, but this range has not beenwell studied and is surely flow-dependent. The applica-tion to planetesimal formation has become focussed onparticle concentration at scales much larger than the Kol-mogorov scale [17, 18] because large clumps are neededfor sufficiently rapid gravitational collapse. In previousparticle clustering cascade models, Hogan and Cuzzi [32]determined the multiplier PDFs for particle concentra-tion and fluid enstrophy at small spatial scales, not toofar from η (to obtain better statistics), and applied themacross all scales ranging up to the integral scale (seesection V.1 for more discussion). Realizing the risks inthis, they performed tests which seemed to validate theapproach. However, discrepancies between Hogan andCuzzi [32] and Pan et al. [20] at the low probabilities ofinterest for the planetesimal problem [18] have led us toexplore the scale dependence of P ( m ) in more detail, inorder to improve the fidelity of the cascade models. It is worth noting at this point that it is not a re-quirement of cascade models that the PDFs be scale-independent; it is merely the first and most obviousassumption. In this paper we present evidence thatthe multiplier PDFs for particle concentration are scale- dependent and present simple guidelines for how thisscale-dependence can be included in cascades. MultiplierPDFs can also be conditioned on local properties [31],and indeed were treated this way in our previous work toallow for particle mass-loading on the process [32]. Scale-dependence per se is, however, a different effect than localconditioning, and in this paper we do not address localconditioning.Before describing our own work, we review some ex-perimental results on high-Reynolds number atmosphericboundary layer turbulence, which provide a useful back-ground in scale invariance and complement the more typ-ical, but lower- Re , numerical simulations. Studies of theproperties of turbulence in atmospheric boundary layerflows have been conducted by Kholmyanskiy [35], vanAtta and Yeh [36] and Meneveau and Sreenivasan [37];further analysis of the Meneveau and Sreenivasan [37]data was done by Chhabra et al. [28].[38] The best ref-erence for the basic experimental data is Meneveau andSreenivasan [37, see their Table 1], who conducted anexperiment on boundary layer turbulence using a sensormounted 2 m above the flat roof of a four story building.The Reynolds number for the flow is calculated using thefree stream velocity U = 6 m/s and the height h = 2 mof the sensor above the roof: Re = U h/ν = 8 × where the kinematic viscosity is 1 . × − m /sec, con-sistent with tabulated values in Meneveau and Sreeni-vasan [37] of the Taylor scale Reynolds number Re λ , andits characteristic lengthscale λ and velocity u (cid:48) . Mene-veau and Sreenivasan [37] give the Kolmogorov scale as η = 7 × − m (we retain their preferred units). Analy-ses of flow structures by Chhabra et al. [28] (their figure5 and our figure 2) show fairly well-behaved power lawscaling of dissipation for weightings which suppress re-gions that are strongly anomalous (panels g and h), i.e.strongly differing from the mean, to almost r ∼ . × η = 12.6 m (cid:29) h , suggesting an extensive inertial range.The large or integral scale L , which contains the energy inthis flow, is thus apparently much larger than the verticaldistance of the sensor from the boundary ( h = 2 m) andplausibly the same as the longitudinal integral scale givenas L = 180 m [37], see also Hunt and Morrison [39], andthus L/η ∼ × . However, when the role of stronglyanomalous regions is emphasized (panels i through l offigure 2) the scalable inertial range contracts.In cascade applications, it may be more meaningful toassess the scale dependence of P ( m ) at large scales not interms of multiples of η as in Sreenivasan and Stolovitzky[31] and most other work [e.g., 24], but in terms of frac-tions of L which more closely connects to causality andenergy flow. We will also express scale fractions r/L in FIG. 1. (a) The second order structure function for particle velocity for the DNS data we analyze. Statistics are obtainedover the trajectories of particles of different St , from figure 2 of Bec et al. [12]. In particular circles and downward facingtriangles refer to tracer-like particles that more faithfully follow the fluid flow. Dotted power laws labeled “2” and “2/3” arethe theoretical expectations for the flow velocity structure function in the viscous range and the inertial range, respectively.The bottom axis is labeled as in Bec et al. [12] by r/η , and also by us (in green) with our estimated values for r/L , where L is the integral scale, and also by the cascade level N that gives equidimensional volumes of side r . Extending out the topof the plot is an offset blue short-dashed line indicating the expected structure function for the high- Re atmospheric flowof Meneveau and Sreenivasan [37], analyzed by Chhabra et al. [28] and Sreenivasan and Stolovitzky [31], with integral scalewhich we estimate as L ∼ × η . We have attempted to place a corresponding scale of r/L and N on the top axis (in blue).(b) Scale-independent multiplier PDFs for dissipation (cid:15) as determined by Sreenivasan and Stolovitzky [31] in the atmosphericboundary layer, connected to the left panel by dotted arrows indicating where those measurements lie on the structure function(300-3000 η , well away from both L and η ). It is our expectation that the scale-independent β -distribution with β ∼ η . Our analyses of the DNS dataset are binned on scales between r = 12 η − η , corresponding to r ∼ L/ L/
86 on the lower scale as discussed in in section I. terms of cascade bifurcation levels N needed to achievecubes r on a side: r/L = 2 − N/ . (1.1)For example, Sreenivasan and Stolovitzky [31] com-pared multiplier PDFs for dissipation in the atmosphericboundary layer over a range of scales (see figure 1b) andshowed that P ( m ) is highly scale independent over a widerange of scales: 372 η − η , or 372 η to L/
86 [see also29]. Dissipation depends on higher-order moments of thevelocity gradients, so we are drawn for guidance to thebehavior seen in the higher-order moments (larger | q | ) infigure 2 [28]. The results of Sreenivasan and Stolovitzky[31] are consistent with the generally power law behav-ior seen for 25 η − η (roughly 25 η to L/
26) in fig-ure 2. That is, one might infer from where the plotsin Chhabra et al. [28] deviate from power law behav-ior, that the scale-free behavior demonstrated by Sreeni-vasan and Stolovitzky [31] (figure 1b) might carry on tolarger sizes than they actually presented, possibly until r ∼ L/
26 or 10000 η , but deviate noticeably for scaleslarger than r ∼ L/
15 (and at the smaller end below25 η ). Moreover, we can conclude from these compar-isons that the scale-free behavior seen by Sreenivasanand Stolovitzky [31] was safely out of the viscous range, and continued through the inertial range at scales up to L/ < L/ Re -number flows not accessible to di-rect numerical simulations. The paper is organized asfollows: section II describes the DNS dataset used in thisstudy; section III describes the data analysis including anovel technique for correcting the effects of small particlenumber statistics, and presents results for the multiplierPDFs for particle concentration, dissipation and enstro-phy; section IV presents predictions of the new cascademodel and comparison with DNS data at two different Re ; and section V discusses the results and their impli-cations. A summary and conclusions are given in sec-tion VI. II. DATASET
In this paper, we use data from the direct numericalsimulations of Bec et al. [12]; see also Arn`eodo et al. [40]and [41]. The simulation computes forced, homogenousand isotropic turbulence in an incompressible fluid, and
FIG. 2. Plot of a family of normalized q -th moments of thedissipation E r in an atmospheric boundary layer, as averagedover binning lengthscale r , plotted against r/η (taken fromfigure 5b of Chhabra et al. [28]). The quantity µ i ( q, l ) ≡ ( E r ) qi / Σ j ( E r ) qj , where ( E r ) i is the dissipation in the i th binof size r . Smaller | q | values suppress the effect of stronglyanomalous regions, while large negative values of q select forregions of anomalously low turbulent dissipation. Abruptchanges in the slope of these plots (most obvious for large | q | ) might indicate departure from the true scale-free inertialrange, both in the dissipation range at < − η , and at verylarge scales where vortex stretching has yet to become effec-tive. Vertical arrows on horizontal axis are the authors [28]estimate of the inertial scaling range, but the scaling range isnarrower for larger | q | . the dynamics of inertial particles suspended in the flow.The fluid flow is solved on a 2048 Cartesian grid witha grid spacing that is approximately the Kolmogorovlength scale η ≈ ∆ x = ∆ y = ∆ z . Tracer and inertialparticles are introduced into the flow and their trajecto-ries are tracked. Particles are considered point particlesand are dragged with the flow by viscous forces only;there is no back-reaction on the flow. Particles of differ-ent Stokes numbers St ≡ τ s /τ η are considered, where τ s is the aerodynamic stopping time of the particle ( τ s = 0for tracers) and τ η is the Kolmogorov time.Figure 1a shows the second order structure functionfor particle velocity for this numerical flow [taken alongtrajectories of different St particles, from figure 2 of 42],as a function of normalized scale r/η . While the struc-ture function seems to show an inertial range to severalthousand η , in reality the integral eddy scale for thissimulation seems to be about half the computational boxsize, L ∼ η , and for this flow Re λ ≈ L/η ) / ∼
400 [table 1 of 42]. In blue-green below the horizontalaxis we indicate the corresponding values for r/L , and the corresponding cascade level N . The blue dashed lineindicates the expected inertial range for the atmosphericflow of Meneveau and Sreenivasan [37], with correspond-ing values of r/L and N also indicated in blue above thefigure. Note that the range where P ( m ) for dissipationwas observed to be scale independent by Sreenivasan andStolovitzky [31] corresponds roughly to the scale range L/ − L/
86, well below the expected integral scale forthat flow and well above the viscous subrange.Data from this simulation are available publicly on-line [43], and we have downloaded and analyzed all ofthe publicly available data in the present work. Thisdata consists of the entire flow field sampled at 13 in-stances in time, and particle trajectories sampled at 4,720equidistant times, both covering about 6 large-eddy timescales τ L . All flow components and their first derivativesare available at the particle locations. In total thereare 3 ×
64 files of particle trajectories each containing3,184 particles (a total of N p ≈ St = 0 , . , . .
0, and 64 files containing 3,184particles each (i.e., a total of N p ≈ St = 2 , , , , , , , , and 70. III. ANALYSIS
Determining concentration multipliers amounts tocounting particles in cubic sub-volumes of size r and cal-culating the fraction of particles falling in each half of thesampling box. We bisect each cube in all three orthogo-nal directions x , y , and z each yielding 2 multiplier valuestotaling 6 multiplier values for r cube. The available tra-jectory data is highly resolved in time (4,720 instances oftime over approximately 6 large eddy times τ L ), muchmore than what is needed for this analysis. The num-ber of snapshots required for good statistics depends onthe scale of interest since structures at large scale evolvemore slowly than structures at small scale (and containmore particles), and therefore can be sampled less often.We choose to sample the particle data with a tempo-ral spacing of τ sample ≈ . τ r where τ r is the charac-teristic eddy life time at spatial scale r estimated usingKolomogorov 1941 arguments as τ r = τ L ( r/L ) / . Forthe box sizes considered, 512 η , 256 η , 128 η , 64 η , 45 η ,32 η , 24 η , 16 η and 12 η , this results in sampling inter-vals of 0 . τ L , 0 . τ L , 0 . τ L , 0 . τ L , 0 . τ L , 0 . τ L ,0 . τ L , 0 . τ L and 0 . τ L , respectively. We popu-late the sample volume using the positions of all particlesfrom the high-resolution trajectory files at these variousdiscrete times. FIG. 3. (a) Probability distribution functions (PDFs) of the particle concentration itself for tracer particles ( St = 0) from thesimulation dataset, compared to Poisson distributions, for box sizes of r = 512 η , 256 η , 128 η , and 64 η . (b) PDFs of the particleconcentration multipliers for St = 0 particles computed from the simulation dataset, and the analytical Poisson/Binomialmodel from equation (3.3). III.1. Tracer particles
First, we will look at the tracer particles (Stokes num-ber St = 0) which follow the flow exactly. Particles areinitially homogeneously distributed, and since the flow isincompressible, will stay that way on the average, withparticle multiplier distributions given by a delta func-tion at m = 0 .
5. However, the observed distributionswill only approach this behavior for very large numberof particles; otherwise, effects of small number statisticscomplicate matters. The simulation dataset at hand has N p ≈ × particles with St = 0, and this is smallenough to cause significant deviations from ideal behav-ior at small scales.For St = 0 particles, we can quantify these small-number effects analytically. First, the probability of find-ing n particles in a given box of size r is governed by aPoisson distribution P P ( n ; ¯ n ) = ¯ n n e − ¯ n n ! (3.1)with the expectation value E ( n ) = ¯ n ( r ), which is herethe average number of particles in a box of size r , that is¯ n ( r ) = N p r / L , and where L is the size of the simulationdomain. This probability is also the probability of findinga concentration C = n/ ¯ n . A comparison of the observedprobability distribution functions with equation (3.1) isshown in figure 3a.Then, for a given box with exactly n particles, considerthe particles one at a time and ask if they fall into onehalf of the box, say the left half, or the other, right half.For tracer particles, the probability to fall into the leftside is the same as falling into the right side of the box,i.e. p left = p right = p = 0 .
5, and the probability of havingexactly k particles fall into one side of the box is then given by a binomial distribution with the probability P B ( k ; n, p ) = (cid:18) nk (cid:19) p n (1 − p ) n − k . (3.2)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 seethat the probability of finding a multiplier m in the entiresimulation domain is P ( m ; ¯ n, p ) = (cid:88) n P P ( n ; ¯ n ) P B ( k = mn ; n, p ) . (3.3)A comparison of this analytical relation with the multi-plier PDFs computed from the tracer particle trajecto-ries in the simulation is shown in figure 3b. It shows thatequation (3.3) models the observed distributions very ac-curately, and also that the distributions become ratherwide at small scales even though the underlying proba-bility distributions are delta functions at m = 0 .
5. Thiseffect is a kind of “false intermittency” due to small-number statistics alone.
III.2. Particle multipliers for non-zero Stokesnumbers
Correcting for finite particle numbers
As we have seen in the previous section, the numberof particles in the dataset is small enough to significantlyaffect the observed multiplier distributions. In the follow-ing we will describe how we can account for these effectsand estimate what the underlying PDFs would be, giveninfinite particle numbers. The goal is to separate thefinite-particle number effects, which may be importantin many applications, from the effects of the turbulentconcentration process. The finite-particle-number effectsin a specific application can always be added back intoour model later (see, e.g., section IV.2).For the following analysis, we will assume a shape forthese PDFs. It has been suggested that, at least in theatmospheric context [31], symmetric beta-distributionsprovide a good approximation for multiplier distributionsof dissipation. Such distributions have also been used inprevious studies of particle concentrations [e.g., 32] andare defined by f ( m ; β ) = (cid:0) m − m (cid:1) β − Γ(2 β )2Γ( β ) (3.4)with Γ being the Gamma function. The parameter β de-termines the width of the distribution with small valuesof β corresponding to wide, i.e. more intermittent, dis-tributions. The width of the β distribution (its standarddeviation) is given by σ ( β ) = (cid:115)
14 (2 β + 1) . (3.5)However, similarly to the tracers, the observed distribu-tion width, σ ( β ), will not only depend on the under-lying β value but also on the number of particles in agiven sample, n , and the number of samples, N s , usedto compute the distribution. In order to characterizethis dependence, we have conducted Monte-Carlo exper-iments. They mimic the finite-particle-number effects inthe DNS under the assumption that the underlying prob-ability distributions are β distributions. The procedureworks in the following way: First, for a given value of β , we draw a random multiplier m from the β distribu-tion. Given the number of particles n in a sample vol-ume (a given box), this corresponds to a partition into n left = mn and n right = (1 − m ) n particles for the twohalves of the sample volume, where n left and n right arenon-integers in general. Then, we draw a random par-ticle number k left from the corresponding Poisson dis-tribution with expectation value E ( k left ) = n left . Thisvalue (and k right = n − k left ) represent one random sam-ple of the number of particles found in two halves ofa box, and correspond to “observed” multiplier values m left = k left /n and m right = k right /n . We repeat the pro-cedure many times and compute the standard deviation σ from all random samples m left and m right combined ( N s total number of samples). The result is a random sam-ple of the “observed” distribution width given n , N s , andthe true underlying value of β . Figure 4 shows the resultsfrom such experiments for selected parameters n , N s . Asone would expect, the scatter is large for a small num-ber of samples, N s , and becomes smaller with increasing N s . Also, one can see that a small number of particlesper box, n , causes the observed distribution width to besystematically larger than the value from equation (3.5)(red and green symbols), but approaches the exact value FIG. 4.
Symbols:
Randomly sampled distribution widths for4 sets of parameters N s , n as a function of β . For each pa-rameter combination, 50 random samples are plotted. Smoothcurve: β ( σ ) given by equation (3.5). as the number of particles gets large (blue and orangesymbols).Now that we understand how we can model the small-number statistics effects, we can proceed to correct theobservations. Given values of n and N s , we conduct theabove described Monte-Carlo experiments for many ran-dom values of β , and find the value of β that resultsin a distribution width closest to an observed width σ .This gives us an estimate of the true underlying β value.In practice, we consider random values of β between 1and ∞ by drawing samples with equal probability withrespect to their true width σ ( β ) (equation 3.5). Oncewe have accumulated at least 100 samples resulting indistribution widths falling in a tolerance range between σ (1 + ∆) − and σ (1 + ∆), we compute the appropri-ate mean and standard deviation of those β values. Ini-tially, a large value of ∆ is chosen, but as more andmore samples are accumulated, ∆ is reduced step bystep by factors of 2 until the standard deviation con-verges, i.e. does not change significantly if ∆ is reducedfurther. We also require that the final ∆ is no morethan 0 .
01 (the width is matched to within at least 1%of the measured width). To speed up the process, once∆ ≤ .
1, we restrict the range of β samples drawn to val-ues between ¯ βe − √ VAR(ln ¯ β ) and ¯ βe √ VAR(ln ¯ β ) where ¯ β and (cid:112) VAR(ln β ) are the mean value and standard devia-tion of β in the tolerance range. Specifically, we define ¯ β as the β value of the distribution whose variance (widthsquared) is the same as the arithmetic mean of the indi-vidual variances, i.e. ¯ β = β (¯ σ ) (3.6)with ¯ σ = 1 M M (cid:88) i =1 σ ( β i ) (3.7)using equation (3.5) and its inverse β ( σ ) = (8 σ ) − − ,and where β i are the β values of all the samples hav-ing widths within the tolerance range around σ . Theuncertainty in ¯ β is measured by its variance:VAR(ln ¯ β ) = VAR( ¯ β )¯ β ≈ (cid:18) ∂β∂σ (¯ σ ) (cid:19) VAR(¯ σ ) = VAR(¯ σ )64¯ σ , (3.8)following non-linear uncertainty propagation truncatingthe series after the first order, and where the variance of¯ σ is computed by:VAR(¯ σ ) = 1 M M (cid:88) i =1 (cid:0) ¯ σ − σ ( β i ) (cid:1) . (3.9)Once converged, the final values of ¯ β and VAR(ln ¯ β ) areestimates of the true β value corrected for finite-particle-number effects, and an estimate of its uncertainty. Volume-averaged β -distributions Now that we know how to remove the effects that afinite number of particles has on the multiplier distribu-tions for a given box size or scale, we compute correctedvolume averages of the β values for any given spatialscale and Stokes number. Specifically, at any given scale r , we first compute the multiplier values for all the sam-pling boxes. We then compute the distribution widths σ ,j from each subset of boxes having the same numberof particles n j . This width is then used to derive a cor-rected β value, which corresponds to a corrected width σ j (through equation 3.5). We then combine all results bysumming the square of these corrected σ j values weightedwith the fractional volume that the boxes with n j par-ticles occupy. The final, combined β value is then givenby this average width through the inverse of (3.5).Figure 5a shows the combined concentration multiplier β values for all Stokes numbers and all spatial scales con-sidered. We did not consider scales smaller than 12 η sincethe number of particles in such small boxes is too smallto make reliable statistical inferences, even with our cor-rection procedure. From the results, it is very apparentthat the concentration multiplier distributions are notonly functions of St , as has been known, but are alsovery much scale dependent. It seems intuitive, however,that at any given scale, the multipliers may depend onlyon the ratio of the particle stopping time and the dynam-ical time at that spatial scale. An effect of this sort wasseen in particle concentration PDFs by Bec et al. [24].Along these lines, we construct a local Stokes number St r ≡ τ s τ r = St (cid:18) rη (cid:19) − / , (3.10) where we have assumed that the dynamical time at scale r is given by τ r = τ η (cid:18) rη (cid:19) / (3.11)following Kolmogorov [44]. By plotting the multiplier β results against the rescaled Stokes number, the curvesapproximately collapse into one (figure 5b), at least forscales not too close to the integral scale. This means thatwhen local scale and stopping time are accounted for, the scaled multiplier β curves are cascade level independent for scales r < L/
10 or so, and are thus highly amenableto cascade models at, in principle, arbitrarily large Re atleast 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 : St r = 2 N/ St L , (3.12)where St L ≡ τ s τ L = St (cid:18) Lη (cid:19) − / (3.13)using the same Kolmogorov scaling as in equation (3.11).Equation (3.12) separates terms that depend only on cas-cade level ( first term ), and particle properties ( secondterm ), and indicates that particles in different flows be-have the same (have the same statistical cascade) if theyhave the same integral-scale Stokes number St L . We willuse this fact later in section IV.2 when we compare thecascade with DNS results from two different simulationsat different Re .It should be noted our observed r − / -scaling seemsto contradict the scaling found by Bec et al. [24] for“quasi-Lagrangian” probability distribution functions ofthe mass density. Following an idea by Maxey [6], theyapproximated the dynamics of inertial particle by thoseof tracers in an appropriate synthetic compressible veloc-ity field and derived a scaling for the rate at which an r -sized “blob” of particles contracts. They argued thatthe scaling of the contraction rate relates to the scalingof the pressure field, give the contraction rate as beingproportional to r − / for the Re of their simulation [45],and find that their density PDFs collapse when scaledwith this contraction rate. Composite PDFs
Combining the statistics from boxes with different par-ticle number to form a single β -distribution that repre-sents the concentration multipliers assumes that multi-plier distributions are the same for all concentrations.Since the particles are independent, that is, they do not FIG. 5. Concentration multiplier distribution β values corrected for finite-particle number sampling effect as a function of theordinary Stokes number St (a), and in terms of the local Stokes number St r (b). Solid lines connect results for the same binningscale for guiding the eye. feel the presence of each other, one would assume thatthere is no such concentration dependence. However,particles concentrate in particular regions of the flow,and therefore flow properties differ in regions of differentparticle density and so the multiplier distributions mayalso be different. We can relax the assumption of equalmultiplier distributions and compute a composite multi-plier distribution by first computing β -distributions foreach particle concentration separately, and then summingthese distributions weighted with the fractional volumethat the boxes with the particular density occupy. Thesecomposite distributions are shown in figures 6, 7 and 8.The following observations can be made: First, com-posite and mean β -distributions may differ in shape, al-though by construction they have identical widths (sec-ond moments). That is, in general, the composite β -distribution is itself not a β -distribution. At large spa-tial scales, flat or exponential tails are apparent in thecomposite PDFs; these do not have the same shape as any β -distribution. At small scales, however, differencesdisappear within the margin of accuracy. Also, at thelargest scale (512 η ) there are enough particles such thatfinite-particle-number effects are small and the raw PDFs(shown in green) are essentially identical to the compos-ite PDFs (shown in black). However, the importance ofcorrecting for sampling effects becomes apparent as welook at smaller scales where fewer particles cause spuri-ous widening of the raw distributions or “false intermit-tency.” Finally, we should mention that at the smallestscales, the number of particles is so small that we areonly measuring multipliers in high-concentration regions,which causes a sampling bias since particles are known toavoid vorticity, and such regions may produce more inter-mittency or broader multiplier PDFs. For instance, theaverage number of particles in a 32 η box, for the Stokes numbers for which we have a total of 6 × particles(see section II) is only (32 η/ η ) N p ≈ .
44. The situ-ation is even worse for, say, r = 12 η and a case with only2 × total particles. The average is then only 0.04 par-ticles per sampling box. Multipliers are measured only inregions with a particle concentration that is at least 100times larger than the mean concentrations since we canonly reasonably measure multipliers if we have at leastseveral particles in a sampling box. III.3. Dissipation and enstrophy multiplierdistributions
We also calculated the multipliers for fluid dissipationand enstrophy. The rate of turbulent dissipation is givenby (cid:15) = 2 νS ij S ij = ν [( ∂ i u j )( ∂ i u j ) + ( ∂ i u j )( ∂ j u i )] , (3.14)where S ij and u i are the strain rate tensor and the com-ponents of the velocity field, respectively, and where weuse the Einstein summation convention. The enstrophyon the other hand is defined as the square of the vorticity: E = |∇ × (cid:126)u | = ( ∂ i u j )( ∂ i u j ) − ( ∂ i u j )( ∂ j u i ) . (3.15)All of these flow velocity derivatives are available in thedataset, both in the particle trajectory data files and inthe flow snapshots. We have computed multiplier distri-bution β values for (cid:15) and E ( β (cid:15) and β E ) from the tra-jectory of tracer particles (as they sample the flow morehomogeneously than non-zero Stokes number particles),and from the full resolution flow snapshots. The resultsare shown in figure 9, although the tracer data is onlyshown for the largest spatial scales since it also suffersfrom finite-particle effects (not corrected here). (a) (b) (c)(d) (e) (f)(g) (h) (i)(j) (k) (l)(m) (n) (o)(p) (q) (r)(s) (t) (u)(v) (w) (x)FIG. 6. Probability distribution functions (PDFs) of particle multipliers measured from DNS particle trajectories for di ↵ erentStokes numbers and a box size of r = 512 ⌘ (a through l) and 256 ⌘ (m through x). Shown are the volume-averaged, mean -distributions ( red dashed curves ) as described in section III.2, the composite PDFs ( black dashed curves ) as defined insection III.2, and raw multiplier histograms uncorrected for finite-particle-number e ↵ ects ( green curves ). Corresponding shadingin red and grey shows the uncertainty in the measured PDFs (plus/minus one standard deviation). FIG. 6. Probability distribution functions (PDFs) of particle multipliers measured from DNS particle trajectories for differentStokes numbers and a box size of r = 512 η (a through l) and 256 η (m through x). Shown are the volume-averaged, mean β -distributions ( red dashed curves ) as described in section III.2, the composite PDFs ( black dashed curves ) as defined insection III.2, and raw multiplier histograms uncorrected for finite-particle-number effects ( green curves ). Corresponding shadingin red and grey shows the uncertainty in the measured PDFs (plus/minus one standard deviation). (a) (b) (c)(d) (e) (f)(g) (h) (i)(j) (k) (l)(m) (n) (o)(p) (q) (r)(s) (t) (u)(v) (w) (x)FIG. 7. Same as figure 6 but for box sizes r = 128 ⌘ (a through l) and 64 ⌘ (m through x). FIG. 7. Same as figure 6 but for box sizes r = 128 η (a through l) and 64 η (m through x). (a) (b) (c)(d) (e) (f)(g) (h) (i)(j) (k) (l)(m) (n) (o)(p) (q) (r)(s) (t) (u)(v) (w) (x)FIG. 8. Same as figure 6 but for box sizes r = 45 ⌘ (a through l) and 32 ⌘ (m through x). FIG. 8. Same as figure 6 but for box sizes r = 45 η (a through l) and 32 η (m through x). FIG. 9. Dissipation and enstrophy multiplier distribution β -values as a function of spatial scale computed from the flowdata and from tracer particle trajectories ( St = 0). Resultsfrom tracer particles are only shown for large spatial scalessince at smaller scales not enough particles are available for areasonable estimation of the multiplier distributions. Note the presence of asymptotes for r (cid:46) η (perhapsbetter thought of as r (cid:46) L/
50, see figure 14) for both,as anticipated [31]. Enstrophy is shown to have widermultiplier distributions (smaller β values) than dissipa-tion, and is therefore more intermittent. This is con-sistent with the findings of Meneveau et al. [46] in sev-eral flows including atmospheric flow, and in numericalsimulations [e.g., 47, 48]. For a review, see Sreenivasanand Antonia [49]. Also, we note that even for the small-est spatial scales considered, still well within the inertialrange, the dissipation rate multiplier β does not reachthe atmospheric flow values of β (cid:15) ∼ IV. NEW CASCADE MODEL WITHLEVEL-DEPENDENT MULTIPLIERSIV.1. Cascade simulations
From the collapsed β ( St r ) curves (figure 5b), we canbuild an empirical model for the particle multiplier dis-tributions. A sum of two power laws approximates thecurves for fixed scale r well: β ( St r ) ≈ β min (cid:32)(cid:18) St r a (cid:19) b + (cid:18) St r a (cid:19) b (cid:33) , (4.1)with parameters a , a , b , b determining the slopesand positions of the exponentials, respectively, and β min setting the minimum β value. From the figure it is evi-dent that there is some residual scale dependence – thecurves for different spatial scales don’t overlap exactly.The following parameterization approximates this resid- FIG. 10. Beta values of the cascade model following equa-tions (4.1) and (4.2) ( dashed lines ). Results from the presentDNS data are shown for comparison ( solid lines ), they arethe same curves as in figure 5 except that the symbols havebeen suppressed here for legibility. Spatial scales, r , are dif-ferentiated by color (see figure legend). ual dependence: a = 0 . ,a = 0 . − .
25 exp (cid:18) −
23 ln (cid:16) . rL (cid:17) (cid:19) ,b = − . ,b = 0 .
85 + 0 . (cid:104) (cid:16) − . (cid:16) . rL (cid:17)(cid:17)(cid:105) ,β min = 4 + 4 (cid:104) (cid:16) ln (cid:16) rL (cid:17)(cid:17)(cid:105) , (4.2)where erf and ln are the error function and the naturallogarithm. The parameters asymptote for small r/L (cid:46) × − (or large cascade level N (cid:38)
25) to: a = 0 . , a = 0 . ,b = − . , b = 1 . , β min = 4 . (4.3)The model is shown in figure 10 compared to the DNSresults. An even simpler model could probably be con-structed using a single average curve of all rL ≤ .With this model for the particle multiplier distribu-tions, we can perform statistical cascade simulations topredict the probability distribution function for the par-ticle concentration. We start at cascade level 0 with asingle concentration value of C (0) = 1 .
0. At every cas-cade level N , we then draw a random multiplier value m ( N ) from the corresponding distribution with a β valuegiven by the model (equation 4.2), and split the con-centration value from the previous level into two values C ( N )left = 2 C ( N − m ( N ) and C ( N )right = 2 C ( N − (1 − m ( N ) ).The factor 2 here comes from the fact that the concen-tration is the ratio of the particle number in a half boxand the mean number in a half box (which itself is one3half of the mean number of particles in full box). Such acascade produces 2 N random samples of concentrationsvalues at each cascade level N which we use to computeconcentration PDFs. For good statistics, however, weneed many more samples, and for the predictions shownin the following section we computed 50,000 such cascadesimulations. IV.2. Level-dependent cascade predictionscompared to DNS at two different Re In order to demonstrate and assess the cascade pre-dictions, we compare the probability distribution func-tions of the concentration factor generated by the cascademodel 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 wehave seen before, can cause observed distributions to besignificantly widened relative to ideal ones that the cas-cade produces. Instead of correcting the DNS PDFs aswe have done before for the measured multiplier distribu-tions, we will degrade the cascade PDFs for this purposeby introducing finite-particle-number effects into them.For motivating the procedure, let us imagine a hy-pothetical simulation H with a number of particles solarge that finite-particle-number effects are negligible,and let ¯ n ∞ ( r ) be the average number of particles ina sampling box at length scale r in that simulation.Our present dataset, in comparison, has on average only¯ n ( r ) = N p r / L particles in a box of scale r , where asbefore N p is the total number of particles in the datasetwith a given Stokes number, and L is the linear extentof the simulation domain, respectively. One can think ofour current dataset as a randomly selected subset of thehypothetical simulation H , generated by retaining parti-cles from H with a probability of p ( r ) = ¯ n ( r ) / ¯ n ∞ ( r ). Forbrevity, we will suppress the r below. Specifically, let’ssay some sampling box in H has n ∞ particles in it (andtherefore a concentration factor C ∞ = n ∞ / ¯ n ∞ ). Fromthese, we select particles with a probability of p , retain-ing in total n particles, where n is an integer randomnumber with an expectation value of E ( n ) = C ∞ ¯ n . For¯ n ∞ → ∞ , this is a Poisson process and n is a randomnumber with a probability mass function P P ( n ; C ∞ ¯ n ) = ( C ∞ ¯ n ) n e − C ∞ ¯ n n ! . (4.4)Using this idea, the recipe for introducing finite-particle-number effects into the cascade PDFs is as fol-lows: First, we draw a random sample C ∞ from acascade-derived concentration PDF. Second, we draw arandom sample n from a Poisson distribution with thecorresponding expectation value E ( n ) = C ∞ ¯ n , where ¯ n is again the average number of particles at the spatialscale of interest in the DNS data we want to compare. The (integer) particle number n corresponds to a dis-crete concentration factor C = n/ ¯ n . By repeating theprocedure N s times, we can build from the samples a dis-crete probability distribution function of C . It accountsfor the finite-particle-number effects and can be directlycompared to PDFs measured from the DNS dataset.Figure 11 shows, for different Stokes numbers, a com-parison between the PDFs predicted by the cascademodel, and degraded in this way using N s = 10 , withthose calculated directly from the DNS dataset we ana-lyzed in this paper [43].Under the assumption that the collapsed multiplier β curve is universal and does not depend on Reynolds num-ber, we can use cascade simulation to model conditionsat different Reynolds numbers. For caveats to this anda suggestion regarding plausible Re dependence, see thediscussion in section V. It follows from equation (3.12)that particles of different Re flows behave the same andhave the same cascade statistics, if they have identicalintegral-scale Stokes numbers St L .Here, we compare our cascade with Pan et al. [20] whoperformed direct numerical simulation of a compressibleflow with suspended initial particles. Their simulation ison a 512 Cartesian grid with an estimated
L/η ∼ L/η ∼ . × particles per Stokes number. Initialcomparison of their uncorrected multiplier PDFs withours showed a clear disagreement but most of that dis-agreement disappears once we take into account finite-particle number effects. Figure 12 shows a comparisonof cumulative PDFs of the concentration between Pan et al. [20] and our cascade model showing reasonableagreement, bearing in mind that the inviscid simulationsof Pan et al. [20] leave a little uncertainty about the valueof St . IV.3. Model prediction for the radial distributionfunction
In many applications, e.g. in terrestrial clouds, par-ticle collisions play an important role, and it is there-fore of great interest to model this process. The rateof collisions depends on two statistical quantities: theradial relative velocity between particles, and the radialdistribution function (RDF), g ( r ), defined as the proba-bility of finding two particles at a given separation nor-malized with respect to homogeneously distributed par-ticles [3, 20]. Relative velocities are beyond the scope ofthe present model, but the cascade model can be used tomake predictions for the RDF. For this purpose, we per-formed statistical simulations similar to section IV.1 butwith an important difference: we are here interested inthe spatial distribution of the concentration, which thencan be used to compute the RDF.Our cascade model, however, only describes the par-4 (a)(b) (c)FIG. 11. Cumulative probability distribution functions for the concentration factor C , for di ↵ erent Stokes numbers, comparingthe DNS measurements ( solid curves ) with the cascade predictions, degraded to account for finite-particle-number e ↵ ects ( dottedcurves ). with an important di ↵ erence: we are here interested inthe spatial distribution of the concentration, which thencan be used to compute the RDF.Our cascade model, however, only describes the par-ticle multiplier PDF as a function of cascade level, anddoes not explicitly contain information about spatial cor-relations . There is therefore some ambiguity in how tocompute concentration fields from the cascade model.We have explored three di ↵ erent methods to assess therange of possible solutions. Starting at the largest scale,we divide a cube of space in half along each spatial direc-tion. This results in eight sub-cubes with half the linearsize. In order to solve for the concentrations in thesesub-cubes uniquely, we need eight equations. The firstmethod – method A – makes the following choice: oneconstraint it given by the fact that the average of the con-centrations over all eight sub-cubes is equal to the meanconcentration, and seven additional constraints are givenby relating the concentration in seven sets of neighboringsub-cubes through multipliers chosen randomly from thecascade model. The specific choice of equations is givenin appendix A. Solving for concentrations to ever smallercubes until some small cuto ↵ length scale ( r/L ) min then yields a statistical realization of the concentration fieldthat can be used to compute the RDF. The probabilityof finding two particles at a given distance is simply theproduct of the concentrations at the two points in space agiven distance apart, averaged over the whole domain. Ifwe start with a unit concentration at the largest scale, thenormalization, that is the probability for homogeneouslydistributed particles, is simply 1. In order to reduce thecomputational cost and storage requirements to trackableamounts, we do not follow all sub-cubes to ever smallerscale but only a random selection of them. One half ofsub-cubes are followed at each cascade level. Two moremethods are obtained by relating the concentrations inthe two half-cubes, for each direction separately, througha random multiplier. This set of equations is underde-termined. For method B we pick one particular solution,while for method C we use a least-square solver to deter-mine the minimum-norm solution. For specifics, again,we refer to appendix A. Conceptually, it is clear thatthe three methods allow for di ↵ erent amounts of spatialrandomness. Method A clearly maximizes intermittencywhile method C leads to the least spatially intermittentsolution.
FIG. 11. Cumulative probability distribution functions for the concentration factor C , for different Stokes numbers, comparingthe DNS measurements ( solid curves ) with the cascade predictions, degraded to account for finite-particle-number effects ( dottedcurves ). ticle multiplier PDF as a function of cascade level, anddoes not explicitly contain information about spatial cor-relations . There is therefore some ambiguity in how tocompute concentration fields from the cascade model.We have explored three different methods to assess therange of possible solutions. Starting at the largest scale,we divide a cube of space in half along each spatial direc-tion. This results in eight sub-cubes with half the linearsize. In order to solve for the concentrations in thesesub-cubes uniquely, we need eight equations. The firstmethod – method A – makes the following choice: oneconstraint it given by the fact that the average of the con-centrations over all eight sub-cubes is equal to the meanconcentration, and seven additional constraints are givenby relating the concentration in seven sets of neighboringsub-cubes through multipliers chosen randomly from thecascade model. The specific choice of equations is givenin appendix A. Solving for concentrations to ever smallercubes until some small cutoff length scale ( r/L ) min thenyields a statistical realization of the concentration fieldthat can be used to compute the RDF. The probabilityof finding two particles at a given distance is simply theproduct of the concentrations at the two points in space agiven distance apart, averaged over the whole domain. Ifwe start with a unit concentration at the largest scale, the normalization, that is the probability for homogeneouslydistributed particles, is simply 1. In order to reduce thecomputational cost and storage requirements to trackableamounts, we do not follow all sub-cubes to ever smallerscale but only a random selection of them. One half ofsub-cubes are followed at each cascade level. Two moremethods are obtained by relating the concentrations inthe two half-cubes, for each direction separately, througha random multiplier. This set of equations is underde-termined. For method B we pick one particular solution,while for method C we use a least-square solver to deter-mine the minimum-norm solution. For specifics, again,we refer to appendix A. Conceptually, it is clear thatthe three methods allow for different amounts of spatialrandomness. Method A clearly maximizes intermittencywhile method C leads to the least spatially intermittentsolution.Figure 13 shows RDFs predicted by our cascade modelsimulations down to cascade level 72 (spatial scales of( r/L ) min ≈ × − ). All methods give qualitatively thesame results, all in good agreement with the Zaichik andAlipchenkov [2] asymptotic Re = ∞ analytical solution,especially regarding the shape and the active range ofscales. In fact, the magnitudes are even close enough,within a factor of order unity. RDFs for different St L FIG. 12. Cumulative probability distribution functions forthe concentration factor C , for different cascade levels com-paring DNS results from figure 8 of Pan et al. [20] with anestimated Stokes number St ≈ . St L ≈ . × − ( solidcurves ) with the prediction of our cascade model degradedto account for the finite-particle-number effects in Pan’s DNSdataset ( dotted curves ). For our cascade, we have used our St = 5 model with a St L ≈ × − close to the value of Pan. computed using method A are plotted in figure 13a as afunction of the scaled distance St − / L r/L , and are shownto collapse for small St L (cid:46) − . Effectively, the scal-ing behavior of the multiplier PDFs (figure 5) is carriedover to the RDF. At large scales, St − / L r/L (cid:38) , theRDF has a value of 1, that is, particles are homoge-neously distributed, and over an active range of scales,10 − (cid:38) St − / L r/L (cid:38) , it rises and reaches an asymp-tote, g . For larger St L in our sample, however, the ac-tive range is shortened by being too close to the integralscale L , and g ( r ) asymptotes at smaller values that are St L -dependent. The theory of Zaichik and Alipchenkov[2] predicts a very similar behavior and their curve forinfinite Reynolds number, effectively for infinitely small St L , is also shown in figure 13.In figure 13b, we compare RDFs from the differentmethods and the Zaichik and Alipchenkov [2] model byscaling them with their respective g . The curves arenearly identical. Method A produces the highest RDFvalues, i.e. the most intermittent concentration distribu-tions, and method C, as it is biased towards lower in-termittency, results in the smallest asymptotic value (seetable I) for values of g . Interestingly, g for method C isthe same as for the Zaichik and Alipchenkov [2] theory.Zaichik and Alipchenkov [2] use various approximationsin their derivation. Among them, they model the turbu- Simulation / Model g Cascade simulation – Method A 13Cascade simulation – Method B 5.9Cascade simulation – Method C 3.9Zaichik and Alipchenkov [2], figure 3 3.9TABLE I. RDF asymptotic values, g , for St − / L r/L (cid:28) St L (cid:28) lence by a Gaussian process which would underestimatethe tails of their probability distributions and thereforeunderestimate intermittency. V. DISCUSSIONV.1. The nature of scale-dependent particleconcentration
As described in section I, the traditional explanationof clustering in terms of centrifugation of particles fromeddies has been replaced with a somewhat more com-plicated and nuanced combination of physical processes[3, 4]. We believe that our results (figures 5 (right)and 10), stripped of their minor variations, represent akind of universal curve (“U-curve”) for β ( St r ) that canbe interpreted in terms of these different processes actingon a particle of some St over a range of scales r .In the St r < . St r decreases [50, andothers]; thus β increases (the multiplier PDF narrows)with decreasing St r in this regime. As St r increases be-yond 1, the concentration effect is weakened by the de-creasing sensitivity of particles to perturbations of anykind by eddies with timescales much shorter than theirstopping times, so β again increases. This effect, whichcould be thought of as an inertia impedance mismatch,has been described in terms of response functions [51–55], but see Pan and Padoan [9], Bec et al. [42], Hubbard[56] and Hopkins [21] for other more recent and moresophisticated analyses.The strongest clustering effect is produced (the multi-plier PDF has the lowest β ) across perhaps one or twodecades of eddy scale r for a given St , centered on thecombined parameter St r ∼ .
3, presumably the regimewhere history effects in particle velocities play the dom-inant role. While our results support the idea that con-centration is generically due to “eddies on the scale of ηSt / η ..” [13, 24], we think a more refined descriptionone could infer from the U-curve is that clustering is the cumulative result of a history of interactions with theflow of energy as it cascades over eddies ranging over two6 FIG. 13. (a) Radial distribution function (RDF) g ( r ) for various values of St L ( coloured solid lines ) determined from cascademodel simulations using method A with the asymptotic RDF value, g , indicated by an arrow. For comparison, the modelof Zaichik and Alipchenkov [2] from their figure 3 for infinite Reynolds number ( black dashed line ) is also shown. (b) RDFsfor St L = 3 × − for the different cascade simulation methods, and the Zaichik and Alipchenkov [2] model, scaled by theirrespective asymptotic value, g , given in table I. decades in size, driving particles ever deeper into a con-centration “attractor” even in the inertial range [1, 3, 4].The other “universal curve” of Zaichik and Alipchenkov[11] (their figure 1) and their improved model [2, their fig-ure 3] reproduced in figure 13 also has this sense. Thatis, 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 cascadethrough eddies of different scale, and a trajectory downone side and then back up the other side of our U-curve.Such a picture would lead the particle concentration as afunction of spatial binning scale to increase sharply oversome particular range of scales r/L or r/η related onlyto St (the active range), and then remain constant to-wards smaller scales where eddy perturbations are feltonly weakly because of, essentially, the poor impedancematch with the particle stopping time.A natural prediction of this model is thus that, at in-finite Re where energy is available on all timescales, farfrom the dissipation range, and in the absence of com-plications such as gravitational settling or fragmentationlimits, the maximum particle concentration should notonly arise over a similar range of scales near r ∼ ηSt / η [2, 11, 13, 24], but also should have a “saturation” am-plitude that is St -independent. Indeed this would seemto be the prediction of Zaichik and Alipchenkov [2, theirfigure 3]. In current simulations [2, 13], as well in simula-tions we have conducted using the cascade, the clusteringof larger St particles ( St L (cid:38) .
03) asymptotes at smallervalues of the RDF than seen for smaller particles (fig-ure 13). We expect this is because the scale at whichthe larger particles reach St r ∼ . r < − η , the en-ergy spectrum of the flow changes as a result of the now-fixed eddy timescale t r = t η [25]. Particles of St η ∼ et al. [25] the question remains as to whether there is anysort of rollover at r (cid:28) η in the RDF of St η =1 particles,or whether an actual singularity would exist for pointparticles. In terrestrial applications, finite particle sizescomparable to η preclude unlimited singular behavior;however, in protoplanetary nebula applications [8, 9, 22],particle sizes of interest (submillimeter to dm) are ordersof magnitude smaller than the Kolmogorov scale (km) sothis is a question of significant interest.Bec et al. [24] explicitly described the dissipationregime as characterized by an “attractor” having frac-tal properties, and indeed multifractal properties weredemonstrated by Hogan et al (1999) for clustering in thisregime. It is known that cascades lead to fractal andmultifractal spatial distributions [31, 46, and referencestherein]. We now suspect that the cascade of Hogan andCuzzi [32], in which the multiplier distributions do seemto obey level- independent scaling, were effectively dissi-pation range cascades . Tests by Hogan and Cuzzi [32]showed good agreement between their level-independentcascade model and DNS. However, the multiplier PDFswere determined at 3 η and all their DNS results were forlow-to-moderate Re λ <
140 such that the integral scaleswere 14 η , 24 η , 45 η , and 86 η . At least the first three ofthese runs lie mostly within the dissipation range, wherescaling does support level independence [24, 25]. It mightbe worth exploring the use of dissipation range cascadesfurther from the standpoint of modeling fractal structure7or to study higher moments of the particle density PDF.Indeed Bec and Ch´etrite [57] present what is, essentially,a cascade model that reproduces aspects of the particleconcentration PDF.Moreover, to our knowledge, while fractal/multifractalbehavior has been shown for particle clustering withinthe dissipation regime, either at scales of a few to tensof η [58], or scales smaller than η [59, 60], no explicitstudy of this property in the inertial range has beendone. It would be of interest to find whether the inertialrange cascade as described by the U-function (figure 5),which is level- dependent but in a predictable way thatis level- independent , would also produce such a distribu-tion, when suitably scaled for St. This could be of use inmodeling radiative transfer properties [16]. V.2. Dissipation, enstrophy, and Re -dependence As mentioned earlier in section III.3, our dissipationmultiplier PDFs have larger β (are narrower) than theexpected β (cid:15) ∼ β (cid:15) ∼ r (cid:46) L/ > η , wehad expected to find scale-free behavior in the dissipationmultiplier PDF over most of this range. However, recall-ing figure 2, especially as selected by larger | q | , whichmore strongly weight the structures where most dissipa-tion occurs, the properties of dissipation are not invariantover as wide a range as is the second order velocity struc-ture function that defines the inertial range in Bec et al. [42].Figure 14 summarizes the overall scale variation of β for dissipation and enstrophy, showing the spatial scaleboth in terms of η and L following figure 1. At largescales, the PDF is narrow for both (large β ) but widenswith decreasing scale. At scales of ∼ η ( ∼ L/
86) itasymptotes at a value which seems to remain scale-free tosmaller scales. Our new asymptotic values do not agreewith values ( β E ∼
10 for enstrophy) found in Hogan andCuzzi [32], or ( β (cid:15) ∼ β E is more properly associated with the dissipation range,but the discrepancy in β (cid:15) alone merits some discussion.We hypothesize that at much higher Re than we canstudy here, the scale-dependence of β (cid:15) morphs in afashion so as to be consistent with Chhabra et al. [28]and Sreenivasan and Stolovitzky [31]; that is, has a scale-free β (cid:15) ∼ η (basedon Chhabra et al. [28]) and probably less than L/ et al. [28] and Chhabra and Jensen[61]). At larger scales we expect β must increase in somesmooth fashion similar to ours, with an overall behaviorschematically shown by the red dotted line in figure 14.Meanwhile, by the logic that enstrophy E is always more intermittent (has smaller β ) than dissipation, we thenalso hypothesize that β E varies as suggested by the bluedotted line.We suspect that our observed β (cid:15) and β E asymptote(for r (cid:46) L/
86) at larger values than would be true formuch higher Re , because the viscous or dissipation range,which bounds the inertial range on its small-scale end andextends to 20-30 η in general, here impinges on the small-scale end of the nascent inertial range, and may preventthe dissipation and enstrophy from ever fully realizingtheir high- Re intermittency. In contrast, at high atmo-spheric Re , the large-scale onset of the inertial range, at3000 − η based on Chhabra et al. [28] and Sreenivasanand Stolovitzky [31], is completely isolated from the vis-cous range, as indicated in figure 1 by the blue axis labelson the upper horizontal axis.The scale-dependence and asymptotic value of β E isimportant, because the process of particle concentrationmay track the properties of E rather than those of (cid:15) basedon the physics involved (section V.1). While it may becoincidental, in our DNS results, the β for inertial parti-cles minimizes at a value very close to our value for β E ,and considerably smaller (more intermittent) than ourvalue of β (cid:15) . A secondary, related hypothesis is that theparticle concentration multipliers may track the behav-ior of enstrophy (if velocities and accelerations are domi-nated by vorticity), and the minima seen in the collapsedcurves of figure 5, which now never fall below 3.0, mightdrop to significantly lower values, making the particleconcentration field more intermittent. For this reason,cascades developed using our current collapsed β ( St r )curves may underestimate the abundance of zones of highconcentration at high Re to some degree. A better un-derstanding of this Re -dependence will be needed to putcascade modeling of particle concentration on quantita-tively solid ground. The original dataset of Sreenivasanand Stolovitzky [31] probably contains enough informa-tion to assess the validity of these hypotheses regardingdissipation, but not for enstrophy. V.3. Speculations regarding the effect of higher Re -numbers As described in section IV, a level-dependent cascadecan be described which captures the inertial range be-havior for arbitrary Re (at least to the point where massloading starts to affect the physics, e.g., see Hogan andCuzzi [32]). This might be of use in modeling particleconcentration in rain clouds, in the protoplanetary neb-ula, or in other applications at very high Re where nu-merical simulations are impractical (such as spatial vari-ations in microwave opacity). There are several reasonsto expect different behavior in the dissipation range.Recalling, however, our numerical discrepancywith Sreenivasan and Stolovitzky [31] and others re-8 FIG. 14. Hypothetical scaling behavior of dissipation andenstrophy: The solid red and blue symbols, as in figure 9,are those we calculate from the numerical flow at Re λ = 400.Dotted lines are hypothetical values at very high Re . High- Re atmospheric dissipation [31] is scale-free at β ≈ L/ − L/ β values for dissipation and enstropy donot seem to vary through the viscous (dissipation) subrange r/η <
20. It is plausible and expected that enstrophy willalways be more intermittent than dissipation (have smaller β values). The atmospheric (high- Re ) value for the scale-free asymptote for dissipation is roughly β =3 (black dottedline). We hypothesize that β values at high Re follow trajec-tories similar to the red and blue dotted lines for dissipationand enstrophy respectively. That is, the observed behaviors(symbols) is the effect of an incompletely developed inertialrange. garding the asymptotic value of β for dissipation, wethink the possibility of Re -dependence of even ourcollapsed β ( St r ) might imply that our results (and otherinertial range results at low Re ) underestimate effectsof concentration, in the sense that at higher Re , theminimum in the universal curve would move to lower β (more intermittency and higher concentrations). Onemight speculate that the minimum β for particles shouldtrack the β for enstrophy instead of for dissipation.In the application to planetesimal formation proposedby Cuzzi et al. [18], it is necessary to create a joint PDFof particle concentration and enstrophy. This motivatesa better understanding of high Re behavior of β for bothparticles and enstrophy. VI. CONCLUSIONS
Our results indicate that the multiplier PDFs for par-ticle concentration (section III.2), dissipation, and en- strophy (section III.3) vary with scale, at least over thelargest decade of spatial scales. We also find that themultiplier PDFs for particle concentration have two com-ponents: a traditional “ β -function” component, and anexponential-tail component (section III.2). We find thatthe concentration multiplier β values collapse to a scale -independent universal curve when plotted against an ap-propriately scaled local Stokes number St r = St ( r/η ) − / (section III.2, in particular equation (3.10) and figure 5),allowing the cascade model to be used for modelinghigher Re conditions not accessible to numerical simu-lation.For dissipation, the “ β ∼
3” asymptotic behaviorof Sreenivasan and Stolovitzky [31] in high Re atmo-spheric flows probably appears at around r ∼ L/
30 or L/
40, and remains constant to smaller scales, at least to r ∼ − η where the dissipation range begins. In thepresent simulation, the integral scale and the dissipationrange are not separated far enough for the dissipation β (cid:15) to reach such low values, and instead it asymptotes forscales below r ∼ η to a value of β (cid:15) ∼
8. Enstrophy,believed to be always more intermittent (smaller β ) thandissipation, asymptotes in the DNS to a value of β E ∼ β value for particle concentration multipliersat the optimal St r .Given that dissipation, and presumably enstrophy,have not reached their scale-independent, asymptotic val-ues seen in very high Re atmospheric flows, it could be ex-pected that our collapsed particle multiplier β ( St r ) curveis also Re dependent. Analyses of this sort for DNS ofparticle-laden flows at significantly higher Reynolds num-ber are therefore highly desirable, as are measurementsof enstrophy multipliers in very high Reynolds numberflows such as atmospheric flows.We have also found that the cascade model can be usedto construct a spatial distribution of particle concentra-tion, that can be carried to arbitrarily high Reynoldsnumbers, and has a very good resemblance to the ana-lytical theory of Zaichik and Alipchenkov [2]. More workis needed to assess the asymptotic level of maximum con-centration for particles of any size (which we find, as didZaichik and Alipchenkov [2], is size invariant, in the infi-nite Re limit). Acknowledgements:
The authors would like tothank Nic Brummell, Karim Shariff, Katepalli Sreeni-vasan, Federico Toschi and Alan Wray for insightful dis-cussions on the topic of this paper, and Federico Toschiand Enrico Calzavarini for help accessing and using theDNS data.We would like to dedicate this paper to Bob Hogan,who passed away in February 2012. Bob was respon-9sible for all the previous computational work done onturbulent concentration by our group from 1992-2010, in-cluding multifractal behavior and our first cascades. Hewould be very interested in these new results, which ex-plain a discrepancy that arose too late for him to resolve.
Appendix A: Three methods for computing thespatial distribution of concentrations
We here provide more details about the three methodsfor computing spatial concentration fields from multiplierdistributions that we used in section IV.3 for computingradial distribution functions.Let us consider a cube of size r/L = 2 − N/ (cascadelevel N ) having a concentration C ( N ) , initially we startwith a cascade level 0 cube having a unit concentration C (0) = 1, and divide it in half along each spatial direc-tion. This results in eight sub-cubes with half the linearsize which correspond to a cascade level of N + 3 (equa-tion 1.1). We denote the concentrations in these sub-cubes as C ( N +3) ijk , where i, j, k ∈ [1 ,
2] are indices denot-ing the left (1) or right (2) sub-cube in the three spatialdirections. Since there are eight unknowns, we need eightequations to uniquely determine the concentrations.The first method – method A – makes the followingchoice: One constraint is given by the fact that the aver-age of the concentrations over all eight sub-cubes is equalto the mean concentration, that is (cid:88) i,j,k C ( N +3) ijk = C ( N ) . (A.1)We get seven more constraints by relating the concen-tration of neighboring sub-cubes to multipliers that arechosen randomly from the cascade model. A possiblechoice are the combinations( C ( N +3)111 + C ( N +3)211 ) m ( N +3)1 = C ( N +3)111 , ( C ( N +3)112 + C ( N +3)212 ) m ( N +3)3 = C ( N +3)112 , ( C ( N +3)111 + C ( N +3)121 ) m ( N +3)5 = C ( N +3)111 , ( C ( N +3)112 + C ( N +3)122 ) m ( N +3)7 = C ( N +3)112 , ( C ( N +3)121 + C ( N +3)221 ) m ( N +3)2 = C ( N +3)221 , ( C ( N +3)122 + C ( N +3)222 ) m ( N +3)4 = C ( N +3)122 ,C ( N +3)121 + C ( N +3)112 ) m ( N +3)6 = C ( N +3)121 , (A.2)where m ( N +3) i with i ∈ [1 , ...,
7] are seven random multi-plier 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 cubesuntil some small cutoff length scale ( r/L ) min yields onestatistical realization of the concentration field.The two other methods are obtained by relating theconcentrations in the two half-cubes, for each direction separately, through a random multiplier, that is (cid:88) j,k C ( N +3)1 ik = 2 C N m ( N +1)1 , (cid:88) i,k C ( N +3) i k = 2 C N m ( N +1)2 , (cid:88) i,j C ( N +3) ij = 2 C N m ( N +1)3 , (A.3)where m ( N +1) i with i ∈ [1 , ...,
3] are multipliers randomlydrawn from the level N + 1 cascade model. Since theseare only three equations, the linear system is underdeter-mined.For method B we choose one particular solution toequations (A.3), namely C ( N +3)111 C N = (cid:16) m ( N +1)1 (cid:17) (cid:16) m ( N +1)2 (cid:17) (cid:16) m ( N +1)3 (cid:17) ,C ( N +3)211 C N = (cid:16) − m ( N +1)1 (cid:17) (cid:16) m ( N +1)2 (cid:17) (cid:16) m ( N +1)3 (cid:17) ,C ( N +3)121 C N = (cid:16) m ( N +1)1 (cid:17) (cid:16) − m ( N +1)2 (cid:17) (cid:16) m ( N +1)3 (cid:17) ,C ( N +3)221 C N = (cid:16) − m ( N +1)1 (cid:17) (cid:16) − m ( N +1)2 (cid:17) (cid:16) m ( N +1)3 (cid:17) ,C ( N +3)112 C N = (cid:16) m ( N +1)1 (cid:17) (cid:16) m ( N +1)2 (cid:17) (cid:16) − m ( N +1)3 (cid:17) ,C ( N +3)212 C N = (cid:16) − m ( N +1)1 (cid:17) (cid:16) m ( N +1)2 (cid:17) (cid:16) − m ( N +1)3 (cid:17) ,C ( N +3)122 C N = (cid:16) m ( N +1)1 (cid:17) (cid:16) − m ( N +1)2 (cid:17) (cid:16) − m ( N +1)3 (cid:17) ,C ( N +3)222 C N = (cid:16) − m ( N +1)1 (cid:17) (cid:16) − m ( N +1)2 (cid:17) (cid:16) − m ( N +1)3 (cid:17) . (A.4)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 ) haveequal probability, either choice would be paired with theopposite given enough random samples (and we do av-erage many random samples to construct each cascade).The eight sub-cubes from each large box have the iden-tical concentrations, just differently distributed, for eachset of m ( N +1) i , i ∈ [1 , ...,
3] whether m ( N +1) i goes to theleft or right box in each case.For the third and final method, method C , we use aleast-square solver to determine the minimum-norm solu-tion of (A.3), that is the solution that is closest to equallydistributed concentrations. Clearly, this causes a biastowards the least intermittent spatial distribution. Forstrongly intermittent multipliers, this method can evenlead to negative concentrations in one of the sub-cubes.Such solutions have to be discarded, and this further bi-ases method C towards minimal intermittency.0 [1] F. Toschi and E. Bodenschatz, Annual Review of FluidMechanics , 375 (2009).[2] L. I. Zaichik and V. M. Alipchenkov, New Journal ofPhysics , 103018 (2009).[3] A. D. Bragg and L. R. Collins, New Journal of Physics , 055013 (2014).[4] A. D. Bragg and L. R. Collins, New Journal of Physics , 055014 (2014).[5] K. Gustavsson and B. Mehlig, ArXiv e-prints; submit-ted to Ann. Rev. Flu. Mech. (2014), arXiv:1412.4374[physics.flu-dyn].[6] M. R. Maxey, Journal of Fluid Mechanics , 441(1987).[7] K. D. Squires and J. K. Eaton, Physics of Fluids , 1169(1991).[8] L. Pan and P. Padoan, Journal of Fluid Mechanics ,73 (2010), arXiv:1005.2419 [physics.flu-dyn].[9] L. Pan and P. Padoan, Astrophys. J. , 12 (2013),arXiv:1305.0307 [astro-ph.EP].[10] A. D. Bragg, P. J. Ireland, and L. R. Collins, Journalof Fluid Mechanics , 327 (2015), arXiv:1501.03842[physics.flu-dyn].[11] L. I. Zaichik and V. M. Alipchenkov, Physics of Fluids , 1776 (2003).[12] J. Bec, L. Biferale, M. Cencini, A. S. Lanotte, andF. Toschi, Journal of Fluid Mechanics , 527 (2010).[13] P. J. Ireland, A. D. Bragg, and L. R. Collins, ArXive-prints (2015), arXiv:1507.07026 [physics.flu-dyn].[14] L. Pan and P. Padoan, Astrophys. J. , 101 (2014),arXiv:1410.1989 [astro-ph.EP].[15] L. Pan and P. Padoan, Astrophys. J. , 10 (2015).[16] R. A. Shaw, A. B. Kostinski, and D. D. Lanterman, J.Quant. Spectr. Rad. Trans. , 13 (2002).[17] J. N. Cuzzi, R. C. Hogan, and K. Shariff, Astrophys. J. , 1432 (2008), arXiv:0804.3526.[18] J. N. Cuzzi, R. C. Hogan, and W. F. Bottke, Icarus ,518 (2010), arXiv:1004.0270 [astro-ph.EP].[19] J. E. Chambers, Icarus , 505 (2010).[20] L. Pan, P. Padoan, J. Scalo, A. G. Kritsuk, and M. L.Norman, The Astrophysical Journal , 6 (2011).[21] P. F. Hopkins, Astrophys. J. , 59 (2014),arXiv:1406.5509.[22] A. Johansen, E. Jacquet, J. N. Cuzzi, A. Morbidelli, andM. Gounelle, “New Paradigms for Asteroid Formation,”in Asteroids IV , edited by P. Michel, F. E. DeMeo, andW. F. Bottke (2015) pp. 471–492.[23] P. F. Hopkins and H. Lee, ArXiv e-prints (2015),arXiv:1510.02477.[24] J. Bec, L. Biferale, M. Cencini, A. Lanotte, S. Musac-chio, and F. Toschi, Physical Review Letters , 084502(2007), nlin/0608045.[25] A. D. Bragg, P. J. Ireland, and L. R. Collins, Phys. Rev.E , 023029 (2015), arXiv:1411.7422 [physics.flu-dyn].[26] C. Meneveau and K. R. Sreenivasan, Physical ReviewLetters , 1424 (1987).[27] C. Meneveau and K. R. Sreenivasan, Journal of FluidMechanics , 429 (1991).[28] A. B. Chhabra, C. Meneveau, R. V. Jensen, and K. R.Sreenivasan, Phys. Rev. A , 5284 (1989).[29] A. B. Chhabra and K. R. Sreenivasan, Physical ReviewLetters , 2762 (1992). [30] K. Sreenivasan and G. Stolovitzky, Acta Mechanica(Suppl.) , 113 (1994).[31] K. R. Sreenivasan and G. Stolovitzky, Journal of Statis-tical Physics , 311 (1995).[32] R. C. Hogan and J. N. Cuzzi, Phys. Rev. E , 056305(2007), arXiv:0704.1810.[33] H. Tennekes and J. L. Lumley, First Course in Turbu-lence (MIT Press, Cambridge, 1972).[34] U. Frisch,
Turbulence. The legacy of A. N. Kolmogorov. (Cambridge University Press, Cambridge (UK), 1995).[35] M. Z. Kholmyanskiy, Izv. Atmos. Ocean. Phys. , 801(1973).[36] C. W. van Atta and T. T. Yeh, Journal of Fluid Mechan-ics , 417 (1975).[37] C. Meneveau and K. R. Sreenivasan, Nuclear Physics BProceedings Supplements , 49 (1987).[38] These results were cited and reanalyzed by Sreenivasanand Stolovitzky [31], however the reference to the basicdata given by Sreenivasan and Stolovitzky [31] is con-fused with an interpretive article by Chhabra and Jensen[61], who themselves cite Meneveau and Sreenivasan [26]and Chhabra et al. [28, at the time unpublished] for dis-cussions and analysis of the basic data.[39] J. C. R. Hunt and J. F. Morrison, Eur. J. Mech. B -Fluids , 673 (2000).[40] A. Arn`eodo, R. Benzi, J. Berg, L. Biferale, E. Bo-denschatz, A. Busse, E. Calzavarini, B. Castaing,M. Cencini, L. Chevillard, R. T. Fisher, R. Grauer,H. Homann, D. Lamb, A. S. Lanotte, E. L´ev`eque,B. L¨uthi, J. Mann, N. Mordant, W.-C. M¨uller, S. Ott,N. T. Ouellette, J.-F. Pinton, S. B. Pope, S. G. Roux,F. Toschi, H. Xu, and P. K. Yeung, Physical ReviewLetters , 254504 (2008), arXiv:0802.3776 [nlin.CD].[41] R. Benzi, L. Biferale, E. Calzavarini, D. Lohse,and F. Toschi, Phys. Rev. E , 066318 (2009),arXiv:0806.4762 [physics.flu-dyn].[42] J. Bec, L. Biferale, A. S. Lanotte, A. Scagliarini, andF. Toschi, Journal of Fluid Mechanics , 497 (2010),arXiv:0904.2314 [physics.flu-dyn].[43] A. S. Lanotte, E. Calzavarini, F. Toschi,J. Bec, L. Biferale, and M. Cencini,(2011), http://data.3tu.nl/repository/uuid:f7cd7b9d-ae4e-498e-92b4-7efe2d350d86 [Online;accessed 16-Dec-2013 and 24-Jan-2014].[44] A. N. Kolmogorov, Journal of Fluid Mechanics , 82(1962).[45] They argue that the scaling should change to r − / forvery high Reynolds numbers Re λ ≥ , 894 (1990).[47] N. Cao, S. Chen, and K. R. Sreenivasan, Physical ReviewLetters , 3799 (1996).[48] S. Chen, K. R. Sreenivasan, and M. Nelkin, PhysicalReview Letters , 1253 (1997).[49] K. R. Sreenivasan and R. A. Antonia, Annual Review ofFluid Mechanics , 435 (1997).[50] J. Chun, D. L. Koch, S. L. Rani, A. Ahluwalia, and L. R.Collins, Journal of Fluid Mechanics , 219 (2005).[51] C. C. Meek and B. G. Jones, Journal of the AtmosphericSciences , Journal of the Atmospheric Sciences , 239(1973).[52] H. J. V¨olk, F. C. Jones, G. E. Morfill, and S. Roeser,Astronomy & Astrophysics , 316 (1980). [53] W. J. Markiewicz, H. Mizuno, and H. J. Voelk, Astron-omy & Astrophysics , 286 (1991).[54] J. N. Cuzzi and R. C. Hogan, Icarus , 127 (2003).[55] C. W. Ormel and J. N. Cuzzi, Astronomy & Astrophysics , 413 (2007), astro-ph/0702303.[56] A. Hubbard, Monthly Notices of the Royal AstronomicalSociety , 784 (2012), arXiv:1207.5365 [astro-ph.EP].[57] J. Bec and R. Ch´etrite, New Journal of Physics , 77(2007), nlin/0701033.[58] R. C. Hogan, J. N. Cuzzi, and A. R. Dobrovolskis, Phys. Rev. E , 1674 (1999).[59] J. Bec, Physics of Fluids , L81 (2003), nlin/0306049.[60] J. Bec, L. Biferale, G. Boffetta, M. Cencini, S. Musac-chio, and F. Toschi, Physics of Fluids , 091702 (2006),nlin/0606024.[61] A. Chhabra and R. V. Jensen, Physical Review Letters62