A Cascade Model for Particle Concentration and Enstrophy in Fully Developed
Turbulence with Mass Loading Feedback
R. C. Hogan
Bay Area Environmental Research Institute; MS 245-3 Moffett Field, CA 94035-1000
J. N. Cuzzi
NASA Ames Research Center; MS 245-3 Moffett Field, CA 94035-1000
(Dated: November 4, 2018)
A cascade model is described based on multiplier distributions determined from 3D direct numer-
ical simulations (DNS) of turbulent particle laden flows, which include two-way coupling between
the phases at global mass loadings equal to unity. The governing Eulerian equations are solved
using psuedo-spectral methods on up to 5123 computional grid points. DNS results for particle
concentration and enstrophy at Taylor microscale Reynolds numbers in the range 34 - 170 were
used to directly determine multiplier distributions on spatial scales 3 times the Kolmogorov length
scale. The multiplier probability distribution functions (PDFs) are well characterized by the β dis-
tribution function. The width of the PDFs, which is a measure of intermittency, decreases with
increasing mass loading within the local region where the multipliers are measured. The functional
form of this dependence is not sensitive to Reynolds numbers in the range considered. A partition
correlation probability is included in the cascade model to account for the observed spatial anti-
correlation between particle concentration and enstrophy. Joint probability distribution functions
of concentration and enstrophy generated using the cascade model are shown to be in excellent
agreement with those derived directly from our 3D simulations. Probabilities predicted by the cas-
cade model are presented at Reynolds numbers well beyond what is achievable by direct simulation.
These results clearly indicate that particle mass loading significantly reduces the probabilities of
high particle concentration and enstrophy relative to those resulting from unloaded runs. Particle
mass density appears to reach a limit at around 100 times the gas density. This approach has
promise for significant computational savings in certain applications.
PACS numbers: 47.61.Jd, 47.27.E-, 47.27.eb
Keywords: Turbulence, Multiphase Flows, Statistical Distributions
I. INTRODUCTION
The study of turbulent flows incorporating heavy par-
ticles in suspension (particles with finite stopping times)
is an important endeavor that has both fundamental and
practical relevance to many scientific and engineering
problems. Such flows have been investigated mainly in
numerical simulations where detailed statistical analysis
of the flow fields is possible [1, 2, 3, 4] These simulations,
limited to relatively low Taylor microscale Reynolds num-
bers Reλ (∼ 40), demonstrated that particles whose fluid
response times are comparable to the lifetime of the
smallest turbulent eddies produce a highly nonuniform
field with intense regions of concentration. Preliminary
indications were that the feedback from such concentra-
tions of particles could locally damp turbulence - how-
ever, the role of this “mass loading” effect in determining
the statistical distributions of particle density and vari-
ous fluid scalars has not been thoroughly studied. Ex-
perimental investigations of turbulence modification by
particles have demonstrated that the degree of turbu-
lence damping increases with particle mass loading and
∗Electronic address: hogan@cosmic.arc.nasa.gov
†Electronic address: jcuzzi@mail.arc.nasa.gov
concentration [4].
The phenomenon known as intermittency can be de-
scribed as intense fluctuations, on small spatial and tem-
poral scales in the turbulent field, that contribute to
the exponential tails of probability distribution functions
(PDFs) of scalars such as velocity increments and gradi-
ents [5, 6, 7], dissipation [8], pressure [9, 10], enstrophy
[11, 12] and velocity circulation [13]. Intermittency in
the density field of preferentially concentrated particles
has also been observed and studied [14, 15].
Although intermittency in turbulence still lacks a com-
plete theoretical understanding, progress has been made
with phenomenological models that capture intermit-
tency in a cascade process. Richardson [16] and later Kol-
mogorov [17] suggested that such models might be used
to explain the process of eddy fragmentation initiated
by unstable large scale structures in a turbulent fluid.
Intermittency in the context of fragmentation though a
cascading process has been studied for large-scale gravi-
tating masses [18] and velocity increments in turbulence
[19]. Simple cascade models were explored by Meneveau
and Sreenivasan [20] and were reviewed by Sreenivasan
and Stolovitzky [21] The scale similarity of random fields
was explored by Novikov [22, 23], with a focus on the
energy dissipation cascade. In Novikov’s work, the ratio
of dissipation averaged over two spheres, one embedded
within the other, served as a measure of enstrophy par-
http://arxiv.org/abs/0704.1810v2
mailto:hogan@cosmic.arc.nasa.gov
mailto:jcuzzi@mail.arc.nasa.gov
titioning between larger and smaller scales. The prob-
ability distribution of these ratios, known as multipliers
or breakdown coefficients, was shown to relate to multi-
fractal and statistical measures (moments) of the velocity
and dissipation fields. A recent review of intermittency
in multiplicative cascades stresses that this theory is a
kinematic description and its connection with the real
dynamics remains unclear [24].
Our previous numerical study of particle concentration
in turbulent flows showed that the particle density field
is a multifractal on scales comparable to the Kolmogorov
length scale [14]. This result suggests that a deeper de-
scription of the statistical properties of the particle con-
centration field, based on multiplier PDFs, may also be
possible. Analytical efforts have suggested that dissipa-
tion and vorticity in the fluid phase should be locally
linked with particle concentration [25]. Numerical work
in this regard has demonstrated that preferential concen-
tration is statistically anticorrelated with low vorticity:
particles tend to concentrate in regions where enstrophy
is relatively weak [26, 27].
In this paper we present a cascade model in the spirit
of Novikov [22, 23] that follows the partitioning of pos-
itive definite scalars associated with both the fluid and
the particles. Multipliers controlling the partitioning of
enstrophy and particle density at each step in the cascade
are drawn from probability distribution functions (PDFs)
which are determined empirically from direct numerical
simulations (DNS). Moreover, the multiplier PDFs are
dependent on, or conditioned by, the particle mass den-
sity or mass loading. The cascade model then generates
joint PDFs for particle concentration and enstrophy at
arbitrary cascade levels. A partitioning correlation prob-
ability is also applied at each cascade level to account for
the observed spatial anticorrelation between enstrophy
and particle concentration [26, 28].
In Section II we describe the cascade model and its
parameters, which are empirically determined from DNS
calculations. Details of the DNS equations, and our nu-
merical methods, are discussed in the Appendix. Results
are shown in section III, including comparisons of joint
PDFs of enstrophy and particle concentration as pre-
dicted by the cascade model with those obtained directly
from the DNS results. Cascade model PDF predictions
at Reynolds numbers well beyond the DNS values are
also presented. In section IV, we summarize our results
and discuss their implications.
II. CASCADE MODEL
A turbulent cascade can be envisioned as an hierar-
chical breakdown of larger eddies into smaller ones that
halts when the fluid viscosity alone can dissipate eddy ki-
netic energy. Eddies or similar turbulent structures such
as vortex tubes are bundles of energy containing vorticity
and dissipation. These structures start with a size com-
parable to the integral scale Λ of the flow, and break down
in steps to a size comparable to the Kolmogorov scale
η before being dissipated away by viscosity. The fluid
vorticity and dissipation exhibit spatial fluctuations that
increase in intensity as the spatial scale decreases. This
phenonemon is known as intermittency and has been ob-
served in a variety of processes with strong nonlinear in-
teractions.
In previous numerical and experimental studies, locally
averaged intermittent dissipation fields with scale at or
near η were used to quantify the statistical properties
of multiplier distributions [21]. Multipliers are random
variables that govern the partitioning of a positive defi-
nite scalar as turbulent structures break down along the
cascade. In these studies the statistical distribution of
multipliers (their PDF) were shown to be invariant over
spatial scales that fall within the turbulent inertial range.
Multifractal properties of the cascading field are deriv-
able from such multiplier distributions [23], and cascade
models based on the iterative application of multipliers
to a cascading variable have been shown to mimic inter-
mittency.
While invariant with level in the inertial range of a
cascade, multiplier PDFs might depend on local proper-
ties of the environment. For instance, Sreenivasan and
Stolovitzky [21] showed that the degree of intermittency
in dissipation increases with the degree of local strain
rate, and constructed multiplier distributions for local
energy dissipation conditioned on the local strain rate.
The physical mechanism behind this effect is believed to
be related to vortex stretching dynamics creating intense
bursts of dissipation.
All the multiplier PDFs measured by Sreenivasan and
Stolovitzky [21], whether conditioned or unconditioned
by local properties, are well characterized by the β dis-
tribution function,
p(m) =
Γ(2β)
Γ(β)2
mβ−1(1−m)β−1 (1)
where m is the multiplier variable and β is a shape con-
trolling parameter. A large β produces a narrow, delta-
function-like curve centered at m = 0.5, whereas β = 1
produces a flat distribution between m = 0 and 1. These
limits for β correspond to uniform and highly intermit-
tent processes respectively. In conditioned multipliers,
the value of β varies with some local property of the
fluid.
Concentration of particles in turbulence is a result of
the active dynamics of eddies on all scales. The process
depends on the scale of the eddies and the corresponding
particle response to those eddies. Intense particle den-
sity fluctuations, akin to intermittency, were observed
in a previous numerical study where it was also shown
that nonuniform particle concentrations have multifrac-
tal scaling properties [14]. These results strongly suggest
that a phenomenological cascade model based on multi-
pliers may adequately describe the particle density field.
Simulations that have included particle feedback on the
fluid through the mass loading effect show that damp-
ing of local turbulence occurs [2, 29]. The latter have
shown that vorticity dynamics is affected locally by par-
ticle feedback. This interplay between the phases could
attenuate vortex stretching and, thereby, diminish local
turbulent intermittency. Multiplier distributions condi-
tioned on local mass loading should therefore be an inte-
gral part of a realistic fluid-particle cascade model.
A. Two-Phase Cascade model
Below we describe a two-phase cascade model that in-
corporates simultaneous multiplier processes for parti-
cle concentration C and fluid enstrophy S, in addition
to a process that models their spatial anticorrelation.
The multiplier distributions are conditioned by the local
particle concentration, as determined empirically from
DNS fields equilibrated to Reλ = 34, 60, 107, and 170.
The spatial anticorrelation was also quantified from these
fields. Local measures of particle concentration (C) and
enstrophy (S) used are defined in the Appendix.
A schematic illustration of our two-phase partition-
ing process is shown in FIG. 1. The cascading vector
(S,C) has components representing enstrophy and parti-
cle concentration. Initially the components are assigned
the value unity and are associated with a common cell
having a volume of unity. Each component is partitioned
into two parts; (mSS, (1−mS)S) and (mCC, (1−mC)C),
respectively, where mS ,mC are multipliers for S and C
whose values are between zero and one inclusive and are
random members of the corresponding multiplier distri-
butions. The parts are associated with two daughter cells
each containing half the volume of the starting cell. In
the example shown in FIG. 1, mS and mC are assumed
to be greater than 0.5. The largest parts of S and C are
placed in the same daughter cell with probability Γ (and
in different cells with probability 1− Γ). This partition-
ing process is repeated for each daughter cell down the
cascade until the ratio of the daughter cell size to the
initial cell size equals a specified cutoff. When this cutoff
is set to the ratio of the turbulent lengthscales Λ and η,
the cascade corresponds to turbulence characterized by
Reλ ∼ (Λ/η)2/3 [30].
B. Conditioned Multipliers
The parameters of the cascade model are empirically
derived from the particle density and enstrophy fields C
and S as calculated by DNS (see Appendix). The simu-
lation parameters for four DNS runs representing Reλ =
36, 60, 104, and 170 are shown in Table I. The turbu-
lence kinetic energy q, the volume averaged dissipation
ǫ, and Λ are calculated from the 3-D turbulent energy
spectrum E(k) and kinematic viscosity ν,
E(k)dk (2)
mC mC *C
1−Γ = .7
Γ = .3
mS *S
mC *C (1− )*CmC
(1− )*S
(1− )*SmS *S
FIG. 1: Figure depicting the breakdown of a parcel of en-
strophy (S) and particle concentration (C) into two parcels
each with half the volume of the parent. The corresponding
multipliers mS and mC are assumed to be greater than 0.5 in
this figure. These measures are broken down and distributed
between the two parcels in one of two ways - the larger por-
tions are partitioned together with probability Γ= 0.3 (upper
figure), or in opposite directions with probability 1− Γ= 0.7
(lower figure).
ǫ = 2ν
E(k)k2dk (3)
dk (4)
where k is wavenumber. kmax =
times the number
of computational nodes per side is the maximum effective
wavenumber. Thus kmaxη > 1 indicates an adequate
resolution of the Kolmogorov scale.
Parameter Case I Case II Case III Case IV
Nodes/side 64 128 256 512
ν .01 .003 .0007 .0002
Reλ 34 60. 104 170
q 1.5 .65 .28 .14
23. 22.8 22.4 23
kmaxη 1.4 1.5 1.45 1.56
14.1 23.3 45.8 86.2
Γ .31 .29 .27 .32
D .0001 .00003 .000007 .000002
νp .001 .0003. .00007 .00002
TABLE I: Case Parameters for DNS runs. The quantities D
and νp are defined in the Appendix. Other quantities above
are defined in Section II.
The 3-D DNS computational box is uniformly subdi-
vided into spatial cells 3η on a side, and the average value
of C and S is determined for each cell ( see Appendix ).
The cells are divided into groups associated with disjoint
ranges of C. Each cell is then divided into two parts of
equal volume and averages for C and S are determined
for each part. The C and S multipliers for each cell are
evaluated as the ratio of these averages to the averages
in the parent cell. A conditional multiplier distribution
p(m) is then determined for each binned value of C from
the corresponding set of cell multipliers. Plots of p(m) for
three values of C are shown in FIG. 2. The points repre-
sent distributions derived from all DNS runs and the solid
lines are least squares fits to the β distribution function
(Eq. 1). For the lower values of C, Reλ-independence
is apparent; only the Reλ = 170 case provided data for
the largest C range. The plots clearly indicate that the
intermittency in C is reduced (multiplier PDFs narrow)
as C is increased. Derived values of βC(C) and βS(C)
are shown as a function of C in FIG. 3. Least squares
fits to the functional form p1 exp(p2C
p3) are drawn as
solid lines and the best fit parameter values for this func-
tion are tabulated in Table II. Bounding curves (dashed
lines) are defined by setting p2 and p3 to their 2σ lim-
its, to establish a plausible range of uncertainty in the
predictions.
Scalar p1 p2 p3
C 2.7 .045 1.02
S 9. .03 1.06
TABLE II: β model parameters
It is certainly of interest that such large solid/gas mass
loadings as C = 100 appear in the DNS runs at all, given
published reports that particle mass loading significantly
dampens turbulent intensity even for mass loadings on
the order of unity [1, 4]. These diverse results might be
reconciled since the particles we study herein are all far
smaller than the Kolmogorov scale and also have only a
very small lag velocity relative to the gas. Recall that
we force the turbulence, as might be the case if it were
FIG. 2: Empirically determined conditional multiplier distri-
butions p(m|C) for particle concentration at three different
mass loading values, C = 1, 20 and 50. The distributions are
obtained from bifurcations of cells with a spatial scale equal
to 3η. Results at Reλ = 34 ( square ), 60 (triangle), 107 (cir-
cle) and 170 ( cross ) are overlain. Only the simulation with
Reλ = 170 provided results for C = 50. At each mass loading
the p(m) at all Reynolds numbers are very well approximated
with the β distribution function ( solid line ). The distribu-
tion widths narrow as the mass loading increases, indicating
a decrease in the intermittency.
being constantly forced by energetic sources operating on
larger scales than our computational volume. However,
FIG. 3 strongly suggests an upper limit for C ( ∼ 100 )
for both βS and βC .
The cascade anticorrelation parameter Γ was deter-
mined by counting the number of parent cells within
which the larger partitions of C and S were found to
share the same daughter cell. This number divided by
the total number of parent cells defines Γ. The derived
Γ value is approximately constant across the DNS cases,
as indicated in Table I. Operationally, the Γ used in
the cascade model was determined by taking a simple
average of the Γ values in Table I.
Overall, the invariance of Γ and the βC(C) and βS(C)
functions across our range of Reλ justifies their treat-
ment as level independent parameters in the two-phase
cascade model. One caveat remains, which would be of
interest to address in future work. While it has been
shown that multiplier distributions leading to βC and βS
are level-invariant over a range of scales within an inertial
range [21], our simulations were numerically restricted
to values of Re in which the inertial range has not yet
become fully developed. Our reliance on the smallest
available scales of 3η to 1.5η (those providing the largest
available intermittency) might lead to some concern that
FIG. 3: The β parameters as functions of local mass loading
C for enstrophy and particle concentration at 3η. Results for
all DNS cases are indicated as described in FIG. 2. A least
squares fit of an exponential function to the points over the
entire mass loading range is shown ( solid line ). Dashed lines
correspond to the upper and lower limits of the function, and
are derived using the 2σ errors of p2 and p3.
they were already sampling the dissipation range of our
calculations, and thus may not be appropriate for a cas-
cade code. We tested this possibility by calculating mul-
tipliers for the next largest level bifurcation (6η to 3η)
for the Reλ = 170 case. The β values for those multi-
plier distributions are slightly larger in value, but consis-
tent with the C-dependence shown in FIG. 2 (6η scales
don’t provide good distribution functions beyond C ∼
15). Thus we believe that for the purpose of demonstrat-
ing this technique, and for the purpose of estimating the
occurrence statistics of C under particle mass loading,
our results are satisfactory. For applications requiring
quantitatively detailed and/or more accurate P (S,C), it
would certainly be of interest to extend the DNS calcu-
lations to larger Re, at which a true inertial range might
be found.
III. MODEL RESULTS
The 2D joint probability distribution function or PDF
of concentration and enstrophy, a fractional volume mea-
sure, was generated from the cascade model and com-
pared with results derived directly from numerical DNS
simulations. The basic probability density P (S,C) gives
the fractional volume occupied by cells having enstrophy
S and concentration C, per unit S and C; thus the frac-
tional volume having C and S in some range ∆S,∆C
is P (S,C)∆S∆C. For quantities varying over orders
of magnitude, it is convenient to adopt ∆S = S and
∆C = C, and we will present the results in the form
P (S,C)SC.
We started by binning results at spatial scale 3η, ob-
tained from the semi-final level of a cascade model run,
into a uniform logarithmic grid of S,C bins each having
width ∆(logS) = ∆(logC) = δ, with corresponding val-
ues of ∆S and ∆C. The number of 3η cells accumulated
in each bin was normalized by the total number of such
cells in the sample to convert it to a fractional volume
∆V (S,C) = P (S,C)∆S∆C. Then
∆V (S,C)
P (S,C)∆S∆C
∆(logS)∆(logC)
→ P (S,C)SC as δ → 0.
In practice of course, the binning ranges δ are not van-
ishingly small.
The plots in FIGs. 4 5 and 7 then, show the PDF
as the volume fraction P (S,C)SC. Cascade levels 9, 12,
15, and 18 correspond approximately to the Reλ of the
four simulation cases shown in Table I. These levels
were determined from the ratio of Λ and η for each case:
level = 3log
(Λ/η). The factor 3 accounts for cascade
bifurcations of 3D cells, because it takes three partition-
ings, along three orthogonal planes, to generate eight
subvolumes of linear dimension one-half that of the par-
ent volume. That is, 2level is equal to the number of
η cells within a 3D volume having linear dimension Λ
and (2level/3)2/3 is the corresponding Reλ. The number
of cascade realizations is, in turn, equal to the product
of the number of Λ-size volumes in the computational
box and the number of simulation snapshots processed.
In general it is difficult to generate DNS results with a
ratio of Λ and η that is an exact power of two. In or-
der to correctly compare DNS simulations with the cas-
cade model it was necessary to interpolate between two
cascade generated P (S,C)SC computed at scale ratios
(levels) that bracketed the ratios that were actually sim-
ulated. In FIG. 4 we compare iso-probability contours of
P (S,C)SC predicted by cascade models representing the
four DNS cases with the same contours derived directly
from the simulated S and C fields. The agreement is very
good.
A. Predictions at higher Reynolds number
The cascade model was used to generate PDFs at
deeper levels in order to assess the effect of mass loading
on the probabilities of high C and S. We generated 256
realizations of a level 24 cascade, 20 realizations of a level
30 cascade, and one realization of a level 36 cascade.
FIG. 5(a) shows the average of 256 realizations of a 24
level cascade, taken to lower probability values. The pro-
nounced crowding of the contours at the top of the figure
indicates the effect of particle mass loading on reducing
the intermittency of C at high values of C. For compar-
FIG. 4: Comparisons of cascade model predictions of
P (S,C)SC with DNS results at Reλ = 34 (a), 60 (b), 107
(c) , and 170 (d). Contours indicate probabilities .001, .01, .1
and .3. Dashed contours are cascade model predictions and
solid ones are DNS results.
ison, FIG. 5(b) shows a control run of a 24 level cascade
with all conditioning turned off. In this control case, the
exponential tails characterizing intermittent fluctuations
are seen at both low and high C.
In order to evaluate the effect of the uncertainties in the
extrapolations of the β curves for C and S on the PDF,
two cascade runs to level 24 were generated using the pa-
rameters for the upper and lower dotted curves in FIG. 3.
In FIG. 6 we show cross-sections of the PDFs produced
by these runs along the C axis through the distribution
modes to compare with the same cross-section for a run
using the nominal parameters in Table II. Both models
diverge from the mean model beyond C > 40, with the
upper (lower) curve corresponding to the outside (inside)
βC(C) and βS(C) bounds in FIG. 3. Figure 6 indicates
that the sensitivity of the PDF to the β model parame-
ters at the 2σ level is only apparent at large C, and all
models show a sharp dropoff in the probability for C >
A crowding effect similar to the one seen in FIG. 5(a)
is shown in FIG. 7 for iso-probability contours equal to
5× 10−4, for cascade levels 6, 12, 18, 24, 30 and 36.
Figures 8(a) and 8(b) compare 1D cuts through the
modes of the PDFs for cascades of 18 - 36 levels, indi-
cating that going to deeper levels (higher Reλ) results in
larger intermittency at the low-C end (as expected), re-
taining the exponential tail characteristic of intermittent
processes, but the highest particle concentration end of
the distribution is extended more slowly. Certainly at
the order of magnitude level, a particle mass loading ra-
tio of 100 times the gas density appears to be as high as
preferential concentration can produce. This result could
be inferred directly from inspection of the conditioned β
distributions of FIG. 3.
FIG. 5: (a) Cascade model predictions for a 24 level case,
taken to lower probability levels, using 256 realizations of the
cascade. Contours are labeled by log(P (S,C)SC). Note the
crowding of contours at high C values, indicating the high-C
limit of the process under conditions of mass loading.(b) A
control cascade to level 24, as in FIG. 5(a), with conditioning
turned off. The difference between (a) and (b) clearly shows
the “choking” effects of particle mass loading on intermittency
in C.
IV. SUMMARY
A two-phase cascade model for enstrophy and parti-
cle concentration in 3-D, isotropic, fully developed tur-
bulence with particle loading feedback has been devel-
oped and tested. Multiplier distributions for enstrophy
and particle concentration were empirically determined
from direct numerical simulation fields at Taylor scale
Reynolds numbers between 34 and 170. These simula-
tions included ‘two-way’ coupling between the phases at
global particle/gas mass loadings equal to unity. The
shape of all multiplier distributions is well characterized
by the β distribution function, with a value of β that
depends systematically on the local degree of mass load-
ing. The values of β increase monotonically with mass
loading and begin to rapidly increase at mass loadings
FIG. 6: 1D cuts through the mode of the PDF of FIG. 5(a)
parallel to the C axis, showing the effects of uncertainty in
the conditioning curve βC(C). The solid curve is the nominal
model and the dashed curves are obtained by allowing the
parameters p2 and p3 to take their 2σ extreme values.
FIG. 7: Cascade model predictions for P (S,C)SC = 5×10−4
for levels 6, 12, 18, 24, 30, and 36. Contour labels indicate
the cascade levels.
greater than 100.
The C-dependent multiplier distributions were used as
input to a cascade model that simulates the breakdown,
or cascade, of enstrophy S and particle concentration C
from large to small spatial scales. The spatial anticorre-
lation between enstrophy and particle concentration was
empirically determined from 3D DNS models and shown
to be constant with Reλ. This constant was used as
a correlation probability governing the relative spatial
distribution of S and C at each bifurcation step in the
cascade model.
The cascade model we have developed clearly repro-
duces the statistical distributions and spatial correlations
observed in our DNS calculations. The cascade parame-
ter values we have derived appear to be universal within
FIG. 8: (a) 1D global cuts through the cascade model PDFs
P (S,C)SC for runs with 18, 24, 30, and 36 levels. (b) closeup
of 1-D cuts through high-C regime.
the range of Reλ of our simulations. We thus specu-
late that they can be used to predict approximate joint
probabilities of enstrophy and particle concentration at
higher Reynolds numbers, at great savings in computer
time. For example, a typical DNS run to Reλ = 170
takes about 170 cpu hours on an Origins 3000 machine,
while a cascade model to an equivalent level takes 0.1 cpu
hours.
We have presented joint probabilites of S and C de-
rived from cascade runs up to level 36. The contours
shown in FIG. 5(a) and FIG. 6 clearly show the effects
of particle mass loading on the probability distribution
functions of C in the regimes where C is large. It appears
that particle mass loadings greater than 100 are rare in
turbulent flows.
The properties of the cascade rest on the physics of our
DNS simulations, and we speculate that two separate ef-
fects are involved. First, particle mass loading dampens
fluid motions of all types, decreasing vorticity stretching
and all other forms of ongoing eddy bifurcation which are
needed to produce intermittency. Second, as a byproduct
of this, particle mass loading may alter the Kolmogorov
timescale locally and shift the most effectively concen-
trated particle Stokes number St to a larger value than
that characterizing particles already lying in the local
volume, reducing the probability of preferentially con-
centrating the local particles any further.
Caveats and Future Work:
As described in section II, our multiplier distributions
were taken from the most numerous cells, with the largest
intermittency, which are at the smallest scales possible
(furthest from the forcing scale). At Reynolds numbers
accessible to DNS, a true inertial range is only beginning
to appear, and while, sampling at the smallest spatial
scales possible, we are as closely approaching the asymp-
totic values within the true inertial range as possible,
where level-independence has been demonstrated in the
past [21], it is possible that our values are subject to
inaccuracy by virtue of being sampled too close to the
dissipation scale. Any such inaccuracy will affect our
cascade results quantitatively but not qualitatively. As
computer power increases, it would be a sensible thing
to continue experiments like these at higher Reλ.
A more general model that treats enstrophy and strain
as independent cascading scalarsmight allow for a higher-
fidelity particle concentration cascade, since C is known
to be linked to the difference between these two scalars
[25] (the so-called second invariant tensor II). However,
such an effort would introduce further complexity of its
own, as II is no longer positive definite. We consider the
development of such a model a suitable task for future
work.
APPENDIX
We used an Eulerian scheme developed by Dr. Alan
Wray to solve the coupled set of fluid/particle equations
used in this study. This was done to maximize the
computational efficiency of the calculations and, more
importantly, to accurately evaluate multipliers over the
wide range of particle concentrations and enstrophies ex-
pected. In this study the effects of particle collisions and
external forces on the particles (e.g., gravity) are not con-
sidered. The turbulence is spectrally forced at k =
such that moments of the Fourier coefficients of the force
satisfy isotropy up to the fourth order. The instanta-
neous Navier-Stokes equations describing the conserva-
tion of mass and momentum for an incompressible fluid
∇ ·U = 0 (A.1)
+(U · ∇)U = −
+ν∇2U−α
(U −V) (A.2)
where U is fluid velocity, V is particle velocity, ρf and
ρp are the fluid and particle mass densities, ν is fluid
viscosity, P is pressure, and α is the inverse of the particle
gas drag stopping time τp.
The compressible equations for the particles are
+∇(ρpV) = D∇2ρp (A.3)
∂(ρpV)
+∇(ρpVV) = νp∇2(ρpV)+αρp(U−V) (A.4)
where νp is a “particle viscosity”, and D is a “particle
diffusivity”. The particle diffusivity and viscosity terms
numerically smooth out particle mass and momentum,
alleviating the formation of steep gradients of ρp that
can lead to numerical instabilities eg. [31].
The right hand sides of Eqs. A.2 and A.4 contain phase
coupling terms which are linearly dependent on (U−V).
The linear form of the coupling follows from the assump-
tions that the particle size is much less than η, and that
the material density of the particles is much greater than
ρf [2]. Additional contributions to the particle-gas cou-
plings involving pressure, viscous and Basset forces [29]
have not been added since they are expected to be weak
in our size regime of interest. The particle field is intro-
duced with a constant mass density and an initial veloc-
ity given by the local gas velocity in a field of statisti-
cally stationary turbulence. All runs are continued until
the particle statistics (RMS of conentration distribution)
have equilibrated.
The particle Stokes number St is defined relative to
the Kolmogorov time scale τη as St = τp/τη, and Φ =
Mp/Mf is the global mass loading, whereMp and Mf are
the total mass of particles and fluid respectively. In this
study ρf , St, and Φ are set to unity, D/ν = 0.01, and
νp/ν = 0.1. Explicitly setting St = 1 guarantees that
the particles are preferentially concentrated. When Φ is
unity, ρp is a surrogate for the local mass loading or local
concentration factor C. The values of νp and D minimize
the diluting effects of numerical particle diffusion while
preventing numerical blowups; their values were deter-
mined from a set of DNS runs in which their values were
decreased systematically until numerical instabilities set
Eqs. A.1 - A.4 are solved using psuedo-spectral meth-
ods commonly used to solve Naviers-Stokes equations for
a turbulent fluid. The Fast Fourier Transform (FFT) al-
gorithm is used to efficiently evaluate the dynamical vari-
ables U, V and ρp on a 3D uniform grid of computional
nodes with periodic boundary conditions. The computa-
tional algorithm is parallelized using MPI and is written
in Fortran 90. All runs for this study were executed on
SGI Origins supercomputers with up to 1024 processors.
Enstrophy is defined as
(∂iUj − ∂jUi)2 (A.5)
where i, j are summed over the three coordinate dimen-
sions of U.
The local spatial average of a scalar over a sample vol-
ume is estimated as,
Fidv (A.6)
where Fi is the scalar’s value on computational node i
centered within a cube of volume dv and the sum is over
all n nodes covering the sample volume. We normalized
this average by the global average value to get a quan-
tity that measures the scalar’s local value relative to its
mean. In this paper C and S will denote normalized spa-
tial averages of particle concentration and enstrophy over
cubes 3η on a side.
ACKNOWLEDGMENTS
We are very grateful to Dr. Alan Wray for providing
the 3-D code and for useful comments on its use. We
thank Robert Last for parallelizing the cascade code on
the SGI Origins 3000. We also would like to thank the
consultants and support staff at the NAS facility for pro-
viding invaluable assistance, and the Science Mission Di-
rectorate of NASA for generous grants of computer time.
We thank Prof. K. Sreenivasan for several helpful con-
versations in the preliminary stages of this project and
the internal reviewers Drs. Alan Wray and Denis Richard
for their suggestions for improving the manuscript. This
research has been made possible by a grant from NASA’s
Planetary Geology and Geophysics program.
[1] K. D. Squires and J. K. Eaton, Phys. Fluids A 2, 1191
(1990).
[2] K. D. Squires and J. K. Eaton, Tech. Rep. MD-55, Stan-
ford University (1990).
[3] K. D. Squires and J. K. Eaton, Phys. Fluids. A 3, 1159
(1990).
[4] J. D. Kulick, J. R. Fessler, and J. K. Eaton, J. Fluid
Mech. 227, 109 (1994).
[5] B. Castaing, Y. Gagne, and E. J. Hopfinger, Physica D
46, 177 (1990).
[6] S. P. G. Dinavahi, K. S. Breuer, and L. Sirovich, Phys.
Fluids 7, 1122 (1995).
[7] P. Kailasnath, K. R. Sreenivasan, and G. Stolovitzky,
Phys. Rev. Lett. 68, 2766 (1992).
[8] A. Vincent and M. Meneguzzi, J. Fluid Mech. 225, 1
(1991).
[9] A. Pumir, Phys. Fluids 6, 2071 (1994).
[10] E. Lamballais, M. Lesieur, and O. Métais, Phys. Rev. E
56, 6761 (1997).
[11] J. Jiménez, A. A. Wray, P. G. Saffman, and R. S. Rogallo,
J. Fluid Mech. 255, 65 (1993).
[12] G. He, S. Chen, R. H. Kraichnan, R. Zhang, and Y. Zhou,
Phys. Rev. Lett. 81, 4636 (1998).
[13] N. Cao, S. Chen, and K. R. Sreenivasan, Phys. Rev. Lett.
76, 616 (1996).
[14] R. C. Hogan, J. N. Cuzzi, and A. R. Dobrovolskis, Phys.
Rev. E 60, 1674 (1999).
[15] E. Balkovsky, G. Falkovich, and A. Fouxon, Phys. Rev.
Lett. 86, 2790 (2001).
[16] L. F. Richardson, Weather Prediction by Numerical Pro-
cess. (Cambridge University Press, Cambridge U.K.,
1922).
[17] A. N. Komolgorov, J. Fluid Mech. 13, 82 (1962).
[18] T. Chiueh, Chin. J. Phys. 32, 319 (1994).
[19] M. Gorokhovski, Tech. Rep., Center for Turbulence Re-
search, Annual Research Briefs (2003).
[20] C. Meneveau and K. R. Sreenivasan, Phys. Rev. Lett.
59, 1424 (1987).
[21] K. R. Sreenivasan and G. Stolovitzky, J. Fluid Mech.
379, 105 (1995).
[22] E. A. Novikov, Phys. Fluids A 2, 814 (1990).
[23] E. A. Novikov, Phys. Rev. E 50, R3303 (1994).
[24] J. Jiménez, J. Fluid Mech. 409, 99 (2000).
[25] M. R. Maxey, Phys. Fluids 30, 1915 (1987).
[26] K. D. Squires and J. K. Eaton, J. Fluid Mech. 226, 1
(1991).
[27] A. M. Ahmed and S. Elghobashi, Phys. Fluids 13, 3346
(2001).
[28] J. K. Eaton and J. R. Fessler, Int. J. Multiphase Flow
20, Suppl., 169 (1994).
[29] S. Elghobashi and G. C. Truesdell, Phys. Fluids A 5,
1790 (1993).
[30] U. Frisch, Turbulence (Cambridge University Press,
Cambridge, U.K., 1995), chap. 8.
[31] A. Johansen, A. C. Anderson, and A. Brandenburg, As-
tron. Astrophys. (2004).

  A cascade model is described based on multiplier distributions determined
from 3D direct numerical simulations (DNS) of turbulent particle laden flows,
which include two-way coupling between the phases at global mass loadings equal
to unity. The governing Eulerian equations are solved using pseudo-spectral
methods on up to 512**3 computional grid points. DNS results for particle
concentration and enstrophy at Taylor microscale Reynolds numbers in the range
34 - 170 were used to directly determine multiplier distributions (PDFs) on
spatial scales 3 times the Kolmogorov length scale. The width of the PDFs,
which is a measure of intermittency, decreases with increasing mass loading
within the local region where the multipliers are measured. The functional form
of this dependence is not sensitive to Reynolds numbers in the range
considered. A partition correlation probability is included in the cascade
model to account for the observed spatial anticorrelation between particle
concentration and enstrophy. Joint probability distribution functions of
concentration and enstrophy generated using the cascade model are shown to be
in excellent agreement with those derived directly from our 3D simulations.
Probabilities predicted by the cascade model are presented at Reynolds numbers
well beyond what is achievable by direct simulation. These results clearly
indicate that particle mass loading significantly reduces the probabilities of
high particle concentration and enstrophy relative to those resulting from
unloaded runs. Particle mass density appears to reach a limit at around 100
times the gas density. This approach has promise for significant computational
savings in certain applications.

A Cascade Model for Particle Concentration and Enstrophy in Fully Developed
Turbulence with Mass Loading Feedback
R. C. Hogan
Bay Area Environmental Research Institute; MS 245-3 Moffett Field, CA 94035-1000
J. N. Cuzzi
NASA Ames Research Center; MS 245-3 Moffett Field, CA 94035-1000
(Dated: November 4, 2018)
A cascade model is described based on multiplier distributions determined from 3D direct numer-
ical simulations (DNS) of turbulent particle laden flows, which include two-way coupling between
the phases at global mass loadings equal to unity. The governing Eulerian equations are solved
using psuedo-spectral methods on up to 5123 computional grid points. DNS results for particle
concentration and enstrophy at Taylor microscale Reynolds numbers in the range 34 - 170 were
used to directly determine multiplier distributions on spatial scales 3 times the Kolmogorov length
scale. The multiplier probability distribution functions (PDFs) are well characterized by the β dis-
tribution function. The width of the PDFs, which is a measure of intermittency, decreases with
increasing mass loading within the local region where the multipliers are measured. The functional
form of this dependence is not sensitive to Reynolds numbers in the range considered. A partition
correlation probability is included in the cascade model to account for the observed spatial anti-
correlation between particle concentration and enstrophy. Joint probability distribution functions
of concentration and enstrophy generated using the cascade model are shown to be in excellent
agreement with those derived directly from our 3D simulations. Probabilities predicted by the cas-
cade model are presented at Reynolds numbers well beyond what is achievable by direct simulation.
These results clearly indicate that particle mass loading significantly reduces the probabilities of
high particle concentration and enstrophy relative to those resulting from unloaded runs. Particle
mass density appears to reach a limit at around 100 times the gas density. This approach has
promise for significant computational savings in certain applications.
PACS numbers: 47.61.Jd, 47.27.E-, 47.27.eb
Keywords: Turbulence, Multiphase Flows, Statistical Distributions
I. INTRODUCTION
The study of turbulent flows incorporating heavy par-
ticles in suspension (particles with finite stopping times)
is an important endeavor that has both fundamental and
practical relevance to many scientific and engineering
problems. Such flows have been investigated mainly in
numerical simulations where detailed statistical analysis
of the flow fields is possible [1, 2, 3, 4] These simulations,
limited to relatively low Taylor microscale Reynolds num-
bers Reλ (∼ 40), demonstrated that particles whose fluid
response times are comparable to the lifetime of the
smallest turbulent eddies produce a highly nonuniform
field with intense regions of concentration. Preliminary
indications were that the feedback from such concentra-
tions of particles could locally damp turbulence - how-
ever, the role of this “mass loading” effect in determining
the statistical distributions of particle density and vari-
ous fluid scalars has not been thoroughly studied. Ex-
perimental investigations of turbulence modification by
particles have demonstrated that the degree of turbu-
lence damping increases with particle mass loading and
∗Electronic address: hogan@cosmic.arc.nasa.gov
†Electronic address: jcuzzi@mail.arc.nasa.gov
concentration [4].
The phenomenon known as intermittency can be de-
scribed as intense fluctuations, on small spatial and tem-
poral scales in the turbulent field, that contribute to
the exponential tails of probability distribution functions
(PDFs) of scalars such as velocity increments and gradi-
ents [5, 6, 7], dissipation [8], pressure [9, 10], enstrophy
[11, 12] and velocity circulation [13]. Intermittency in
the density field of preferentially concentrated particles
has also been observed and studied [14, 15].
Although intermittency in turbulence still lacks a com-
plete theoretical understanding, progress has been made
with phenomenological models that capture intermit-
tency in a cascade process. Richardson [16] and later Kol-
mogorov [17] suggested that such models might be used
to explain the process of eddy fragmentation initiated
by unstable large scale structures in a turbulent fluid.
Intermittency in the context of fragmentation though a
cascading process has been studied for large-scale gravi-
tating masses [18] and velocity increments in turbulence
[19]. Simple cascade models were explored by Meneveau
and Sreenivasan [20] and were reviewed by Sreenivasan
and Stolovitzky [21] The scale similarity of random fields
was explored by Novikov [22, 23], with a focus on the
energy dissipation cascade. In Novikov’s work, the ratio
of dissipation averaged over two spheres, one embedded
within the other, served as a measure of enstrophy par-
http://arxiv.org/abs/0704.1810v2
mailto:hogan@cosmic.arc.nasa.gov
mailto:jcuzzi@mail.arc.nasa.gov
titioning between larger and smaller scales. The prob-
ability distribution of these ratios, known as multipliers
or breakdown coefficients, was shown to relate to multi-
fractal and statistical measures (moments) of the velocity
and dissipation fields. A recent review of intermittency
in multiplicative cascades stresses that this theory is a
kinematic description and its connection with the real
dynamics remains unclear [24].
Our previous numerical study of particle concentration
in turbulent flows showed that the particle density field
is a multifractal on scales comparable to the Kolmogorov
length scale [14]. This result suggests that a deeper de-
scription of the statistical properties of the particle con-
centration field, based on multiplier PDFs, may also be
possible. Analytical efforts have suggested that dissipa-
tion and vorticity in the fluid phase should be locally
linked with particle concentration [25]. Numerical work
in this regard has demonstrated that preferential concen-
tration is statistically anticorrelated with low vorticity:
particles tend to concentrate in regions where enstrophy
is relatively weak [26, 27].
In this paper we present a cascade model in the spirit
of Novikov [22, 23] that follows the partitioning of pos-
itive definite scalars associated with both the fluid and
the particles. Multipliers controlling the partitioning of
enstrophy and particle density at each step in the cascade
are drawn from probability distribution functions (PDFs)
which are determined empirically from direct numerical
simulations (DNS). Moreover, the multiplier PDFs are
dependent on, or conditioned by, the particle mass den-
sity or mass loading. The cascade model then generates
joint PDFs for particle concentration and enstrophy at
arbitrary cascade levels. A partitioning correlation prob-
ability is also applied at each cascade level to account for
the observed spatial anticorrelation between enstrophy
and particle concentration [26, 28].
In Section II we describe the cascade model and its
parameters, which are empirically determined from DNS
calculations. Details of the DNS equations, and our nu-
merical methods, are discussed in the Appendix. Results
are shown in section III, including comparisons of joint
PDFs of enstrophy and particle concentration as pre-
dicted by the cascade model with those obtained directly
from the DNS results. Cascade model PDF predictions
at Reynolds numbers well beyond the DNS values are
also presented. In section IV, we summarize our results
and discuss their implications.
II. CASCADE MODEL
A turbulent cascade can be envisioned as an hierar-
chical breakdown of larger eddies into smaller ones that
halts when the fluid viscosity alone can dissipate eddy ki-
netic energy. Eddies or similar turbulent structures such
as vortex tubes are bundles of energy containing vorticity
and dissipation. These structures start with a size com-
parable to the integral scale Λ of the flow, and break down
in steps to a size comparable to the Kolmogorov scale
η before being dissipated away by viscosity. The fluid
vorticity and dissipation exhibit spatial fluctuations that
increase in intensity as the spatial scale decreases. This
phenonemon is known as intermittency and has been ob-
served in a variety of processes with strong nonlinear in-
teractions.
In previous numerical and experimental studies, locally
averaged intermittent dissipation fields with scale at or
near η were used to quantify the statistical properties
of multiplier distributions [21]. Multipliers are random
variables that govern the partitioning of a positive defi-
nite scalar as turbulent structures break down along the
cascade. In these studies the statistical distribution of
multipliers (their PDF) were shown to be invariant over
spatial scales that fall within the turbulent inertial range.
Multifractal properties of the cascading field are deriv-
able from such multiplier distributions [23], and cascade
models based on the iterative application of multipliers
to a cascading variable have been shown to mimic inter-
mittency.
While invariant with level in the inertial range of a
cascade, multiplier PDFs might depend on local proper-
ties of the environment. For instance, Sreenivasan and
Stolovitzky [21] showed that the degree of intermittency
in dissipation increases with the degree of local strain
rate, and constructed multiplier distributions for local
energy dissipation conditioned on the local strain rate.
The physical mechanism behind this effect is believed to
be related to vortex stretching dynamics creating intense
bursts of dissipation.
All the multiplier PDFs measured by Sreenivasan and
Stolovitzky [21], whether conditioned or unconditioned
by local properties, are well characterized by the β dis-
tribution function,
p(m) =
Γ(2β)
Γ(β)2
mβ−1(1−m)β−1 (1)
where m is the multiplier variable and β is a shape con-
trolling parameter. A large β produces a narrow, delta-
function-like curve centered at m = 0.5, whereas β = 1
produces a flat distribution between m = 0 and 1. These
limits for β correspond to uniform and highly intermit-
tent processes respectively. In conditioned multipliers,
the value of β varies with some local property of the
fluid.
Concentration of particles in turbulence is a result of
the active dynamics of eddies on all scales. The process
depends on the scale of the eddies and the corresponding
particle response to those eddies. Intense particle den-
sity fluctuations, akin to intermittency, were observed
in a previous numerical study where it was also shown
that nonuniform particle concentrations have multifrac-
tal scaling properties [14]. These results strongly suggest
that a phenomenological cascade model based on multi-
pliers may adequately describe the particle density field.
Simulations that have included particle feedback on the
fluid through the mass loading effect show that damp-
ing of local turbulence occurs [2, 29]. The latter have
shown that vorticity dynamics is affected locally by par-
ticle feedback. This interplay between the phases could
attenuate vortex stretching and, thereby, diminish local
turbulent intermittency. Multiplier distributions condi-
tioned on local mass loading should therefore be an inte-
gral part of a realistic fluid-particle cascade model.
A. Two-Phase Cascade model
Below we describe a two-phase cascade model that in-
corporates simultaneous multiplier processes for parti-
cle concentration C and fluid enstrophy S, in addition
to a process that models their spatial anticorrelation.
The multiplier distributions are conditioned by the local
particle concentration, as determined empirically from
DNS fields equilibrated to Reλ = 34, 60, 107, and 170.
The spatial anticorrelation was also quantified from these
fields. Local measures of particle concentration (C) and
enstrophy (S) used are defined in the Appendix.
A schematic illustration of our two-phase partition-
ing process is shown in FIG. 1. The cascading vector
(S,C) has components representing enstrophy and parti-
cle concentration. Initially the components are assigned
the value unity and are associated with a common cell
having a volume of unity. Each component is partitioned
into two parts; (mSS, (1−mS)S) and (mCC, (1−mC)C),
respectively, where mS ,mC are multipliers for S and C
whose values are between zero and one inclusive and are
random members of the corresponding multiplier distri-
butions. The parts are associated with two daughter cells
each containing half the volume of the starting cell. In
the example shown in FIG. 1, mS and mC are assumed
to be greater than 0.5. The largest parts of S and C are
placed in the same daughter cell with probability Γ (and
in different cells with probability 1− Γ). This partition-
ing process is repeated for each daughter cell down the
cascade until the ratio of the daughter cell size to the
initial cell size equals a specified cutoff. When this cutoff
is set to the ratio of the turbulent lengthscales Λ and η,
the cascade corresponds to turbulence characterized by
Reλ ∼ (Λ/η)2/3 [30].
B. Conditioned Multipliers
The parameters of the cascade model are empirically
derived from the particle density and enstrophy fields C
and S as calculated by DNS (see Appendix). The simu-
lation parameters for four DNS runs representing Reλ =
36, 60, 104, and 170 are shown in Table I. The turbu-
lence kinetic energy q, the volume averaged dissipation
ǫ, and Λ are calculated from the 3-D turbulent energy
spectrum E(k) and kinematic viscosity ν,
E(k)dk (2)
mC mC *C
1−Γ = .7
Γ = .3
mS *S
mC *C (1− )*CmC
(1− )*S
(1− )*SmS *S
FIG. 1: Figure depicting the breakdown of a parcel of en-
strophy (S) and particle concentration (C) into two parcels
each with half the volume of the parent. The corresponding
multipliers mS and mC are assumed to be greater than 0.5 in
this figure. These measures are broken down and distributed
between the two parcels in one of two ways - the larger por-
tions are partitioned together with probability Γ= 0.3 (upper
figure), or in opposite directions with probability 1− Γ= 0.7
(lower figure).
ǫ = 2ν
E(k)k2dk (3)
dk (4)
where k is wavenumber. kmax =
times the number
of computational nodes per side is the maximum effective
wavenumber. Thus kmaxη > 1 indicates an adequate
resolution of the Kolmogorov scale.
Parameter Case I Case II Case III Case IV
Nodes/side 64 128 256 512
ν .01 .003 .0007 .0002
Reλ 34 60. 104 170
q 1.5 .65 .28 .14
23. 22.8 22.4 23
kmaxη 1.4 1.5 1.45 1.56
14.1 23.3 45.8 86.2
Γ .31 .29 .27 .32
D .0001 .00003 .000007 .000002
νp .001 .0003. .00007 .00002
TABLE I: Case Parameters for DNS runs. The quantities D
and νp are defined in the Appendix. Other quantities above
are defined in Section II.
The 3-D DNS computational box is uniformly subdi-
vided into spatial cells 3η on a side, and the average value
of C and S is determined for each cell ( see Appendix ).
The cells are divided into groups associated with disjoint
ranges of C. Each cell is then divided into two parts of
equal volume and averages for C and S are determined
for each part. The C and S multipliers for each cell are
evaluated as the ratio of these averages to the averages
in the parent cell. A conditional multiplier distribution
p(m) is then determined for each binned value of C from
the corresponding set of cell multipliers. Plots of p(m) for
three values of C are shown in FIG. 2. The points repre-
sent distributions derived from all DNS runs and the solid
lines are least squares fits to the β distribution function
(Eq. 1). For the lower values of C, Reλ-independence
is apparent; only the Reλ = 170 case provided data for
the largest C range. The plots clearly indicate that the
intermittency in C is reduced (multiplier PDFs narrow)
as C is increased. Derived values of βC(C) and βS(C)
are shown as a function of C in FIG. 3. Least squares
fits to the functional form p1 exp(p2C
p3) are drawn as
solid lines and the best fit parameter values for this func-
tion are tabulated in Table II. Bounding curves (dashed
lines) are defined by setting p2 and p3 to their 2σ lim-
its, to establish a plausible range of uncertainty in the
predictions.
Scalar p1 p2 p3
C 2.7 .045 1.02
S 9. .03 1.06
TABLE II: β model parameters
It is certainly of interest that such large solid/gas mass
loadings as C = 100 appear in the DNS runs at all, given
published reports that particle mass loading significantly
dampens turbulent intensity even for mass loadings on
the order of unity [1, 4]. These diverse results might be
reconciled since the particles we study herein are all far
smaller than the Kolmogorov scale and also have only a
very small lag velocity relative to the gas. Recall that
we force the turbulence, as might be the case if it were
FIG. 2: Empirically determined conditional multiplier distri-
butions p(m|C) for particle concentration at three different
mass loading values, C = 1, 20 and 50. The distributions are
obtained from bifurcations of cells with a spatial scale equal
to 3η. Results at Reλ = 34 ( square ), 60 (triangle), 107 (cir-
cle) and 170 ( cross ) are overlain. Only the simulation with
Reλ = 170 provided results for C = 50. At each mass loading
the p(m) at all Reynolds numbers are very well approximated
with the β distribution function ( solid line ). The distribu-
tion widths narrow as the mass loading increases, indicating
a decrease in the intermittency.
being constantly forced by energetic sources operating on
larger scales than our computational volume. However,
FIG. 3 strongly suggests an upper limit for C ( ∼ 100 )
for both βS and βC .
The cascade anticorrelation parameter Γ was deter-
mined by counting the number of parent cells within
which the larger partitions of C and S were found to
share the same daughter cell. This number divided by
the total number of parent cells defines Γ. The derived
Γ value is approximately constant across the DNS cases,
as indicated in Table I. Operationally, the Γ used in
the cascade model was determined by taking a simple
average of the Γ values in Table I.
Overall, the invariance of Γ and the βC(C) and βS(C)
functions across our range of Reλ justifies their treat-
ment as level independent parameters in the two-phase
cascade model. One caveat remains, which would be of
interest to address in future work. While it has been
shown that multiplier distributions leading to βC and βS
are level-invariant over a range of scales within an inertial
range [21], our simulations were numerically restricted
to values of Re in which the inertial range has not yet
become fully developed. Our reliance on the smallest
available scales of 3η to 1.5η (those providing the largest
available intermittency) might lead to some concern that
FIG. 3: The β parameters as functions of local mass loading
C for enstrophy and particle concentration at 3η. Results for
all DNS cases are indicated as described in FIG. 2. A least
squares fit of an exponential function to the points over the
entire mass loading range is shown ( solid line ). Dashed lines
correspond to the upper and lower limits of the function, and
are derived using the 2σ errors of p2 and p3.
they were already sampling the dissipation range of our
calculations, and thus may not be appropriate for a cas-
cade code. We tested this possibility by calculating mul-
tipliers for the next largest level bifurcation (6η to 3η)
for the Reλ = 170 case. The β values for those multi-
plier distributions are slightly larger in value, but consis-
tent with the C-dependence shown in FIG. 2 (6η scales
don’t provide good distribution functions beyond C ∼
15). Thus we believe that for the purpose of demonstrat-
ing this technique, and for the purpose of estimating the
occurrence statistics of C under particle mass loading,
our results are satisfactory. For applications requiring
quantitatively detailed and/or more accurate P (S,C), it
would certainly be of interest to extend the DNS calcu-
lations to larger Re, at which a true inertial range might
be found.
III. MODEL RESULTS
The 2D joint probability distribution function or PDF
of concentration and enstrophy, a fractional volume mea-
sure, was generated from the cascade model and com-
pared with results derived directly from numerical DNS
simulations. The basic probability density P (S,C) gives
the fractional volume occupied by cells having enstrophy
S and concentration C, per unit S and C; thus the frac-
tional volume having C and S in some range ∆S,∆C
is P (S,C)∆S∆C. For quantities varying over orders
of magnitude, it is convenient to adopt ∆S = S and
∆C = C, and we will present the results in the form
P (S,C)SC.
We started by binning results at spatial scale 3η, ob-
tained from the semi-final level of a cascade model run,
into a uniform logarithmic grid of S,C bins each having
width ∆(logS) = ∆(logC) = δ, with corresponding val-
ues of ∆S and ∆C. The number of 3η cells accumulated
in each bin was normalized by the total number of such
cells in the sample to convert it to a fractional volume
∆V (S,C) = P (S,C)∆S∆C. Then
∆V (S,C)
P (S,C)∆S∆C
∆(logS)∆(logC)
→ P (S,C)SC as δ → 0.
In practice of course, the binning ranges δ are not van-
ishingly small.
The plots in FIGs. 4 5 and 7 then, show the PDF
as the volume fraction P (S,C)SC. Cascade levels 9, 12,
15, and 18 correspond approximately to the Reλ of the
four simulation cases shown in Table I. These levels
were determined from the ratio of Λ and η for each case:
level = 3log
(Λ/η). The factor 3 accounts for cascade
bifurcations of 3D cells, because it takes three partition-
ings, along three orthogonal planes, to generate eight
subvolumes of linear dimension one-half that of the par-
ent volume. That is, 2level is equal to the number of
η cells within a 3D volume having linear dimension Λ
and (2level/3)2/3 is the corresponding Reλ. The number
of cascade realizations is, in turn, equal to the product
of the number of Λ-size volumes in the computational
box and the number of simulation snapshots processed.
In general it is difficult to generate DNS results with a
ratio of Λ and η that is an exact power of two. In or-
der to correctly compare DNS simulations with the cas-
cade model it was necessary to interpolate between two
cascade generated P (S,C)SC computed at scale ratios
(levels) that bracketed the ratios that were actually sim-
ulated. In FIG. 4 we compare iso-probability contours of
P (S,C)SC predicted by cascade models representing the
four DNS cases with the same contours derived directly
from the simulated S and C fields. The agreement is very
good.
A. Predictions at higher Reynolds number
The cascade model was used to generate PDFs at
deeper levels in order to assess the effect of mass loading
on the probabilities of high C and S. We generated 256
realizations of a level 24 cascade, 20 realizations of a level
30 cascade, and one realization of a level 36 cascade.
FIG. 5(a) shows the average of 256 realizations of a 24
level cascade, taken to lower probability values. The pro-
nounced crowding of the contours at the top of the figure
indicates the effect of particle mass loading on reducing
the intermittency of C at high values of C. For compar-
FIG. 4: Comparisons of cascade model predictions of
P (S,C)SC with DNS results at Reλ = 34 (a), 60 (b), 107
(c) , and 170 (d). Contours indicate probabilities .001, .01, .1
and .3. Dashed contours are cascade model predictions and
solid ones are DNS results.
ison, FIG. 5(b) shows a control run of a 24 level cascade
with all conditioning turned off. In this control case, the
exponential tails characterizing intermittent fluctuations
are seen at both low and high C.
In order to evaluate the effect of the uncertainties in the
extrapolations of the β curves for C and S on the PDF,
two cascade runs to level 24 were generated using the pa-
rameters for the upper and lower dotted curves in FIG. 3.
In FIG. 6 we show cross-sections of the PDFs produced
by these runs along the C axis through the distribution
modes to compare with the same cross-section for a run
using the nominal parameters in Table II. Both models
diverge from the mean model beyond C > 40, with the
upper (lower) curve corresponding to the outside (inside)
βC(C) and βS(C) bounds in FIG. 3. Figure 6 indicates
that the sensitivity of the PDF to the β model parame-
ters at the 2σ level is only apparent at large C, and all
models show a sharp dropoff in the probability for C >
A crowding effect similar to the one seen in FIG. 5(a)
is shown in FIG. 7 for iso-probability contours equal to
5× 10−4, for cascade levels 6, 12, 18, 24, 30 and 36.
Figures 8(a) and 8(b) compare 1D cuts through the
modes of the PDFs for cascades of 18 - 36 levels, indi-
cating that going to deeper levels (higher Reλ) results in
larger intermittency at the low-C end (as expected), re-
taining the exponential tail characteristic of intermittent
processes, but the highest particle concentration end of
the distribution is extended more slowly. Certainly at
the order of magnitude level, a particle mass loading ra-
tio of 100 times the gas density appears to be as high as
preferential concentration can produce. This result could
be inferred directly from inspection of the conditioned β
distributions of FIG. 3.
FIG. 5: (a) Cascade model predictions for a 24 level case,
taken to lower probability levels, using 256 realizations of the
cascade. Contours are labeled by log(P (S,C)SC). Note the
crowding of contours at high C values, indicating the high-C
limit of the process under conditions of mass loading.(b) A
control cascade to level 24, as in FIG. 5(a), with conditioning
turned off. The difference between (a) and (b) clearly shows
the “choking” effects of particle mass loading on intermittency
in C.
IV. SUMMARY
A two-phase cascade model for enstrophy and parti-
cle concentration in 3-D, isotropic, fully developed tur-
bulence with particle loading feedback has been devel-
oped and tested. Multiplier distributions for enstrophy
and particle concentration were empirically determined
from direct numerical simulation fields at Taylor scale
Reynolds numbers between 34 and 170. These simula-
tions included ‘two-way’ coupling between the phases at
global particle/gas mass loadings equal to unity. The
shape of all multiplier distributions is well characterized
by the β distribution function, with a value of β that
depends systematically on the local degree of mass load-
ing. The values of β increase monotonically with mass
loading and begin to rapidly increase at mass loadings
FIG. 6: 1D cuts through the mode of the PDF of FIG. 5(a)
parallel to the C axis, showing the effects of uncertainty in
the conditioning curve βC(C). The solid curve is the nominal
model and the dashed curves are obtained by allowing the
parameters p2 and p3 to take their 2σ extreme values.
FIG. 7: Cascade model predictions for P (S,C)SC = 5×10−4
for levels 6, 12, 18, 24, 30, and 36. Contour labels indicate
the cascade levels.
greater than 100.
The C-dependent multiplier distributions were used as
input to a cascade model that simulates the breakdown,
or cascade, of enstrophy S and particle concentration C
from large to small spatial scales. The spatial anticorre-
lation between enstrophy and particle concentration was
empirically determined from 3D DNS models and shown
to be constant with Reλ. This constant was used as
a correlation probability governing the relative spatial
distribution of S and C at each bifurcation step in the
cascade model.
The cascade model we have developed clearly repro-
duces the statistical distributions and spatial correlations
observed in our DNS calculations. The cascade parame-
ter values we have derived appear to be universal within
FIG. 8: (a) 1D global cuts through the cascade model PDFs
P (S,C)SC for runs with 18, 24, 30, and 36 levels. (b) closeup
of 1-D cuts through high-C regime.
the range of Reλ of our simulations. We thus specu-
late that they can be used to predict approximate joint
probabilities of enstrophy and particle concentration at
higher Reynolds numbers, at great savings in computer
time. For example, a typical DNS run to Reλ = 170
takes about 170 cpu hours on an Origins 3000 machine,
while a cascade model to an equivalent level takes 0.1 cpu
hours.
We have presented joint probabilites of S and C de-
rived from cascade runs up to level 36. The contours
shown in FIG. 5(a) and FIG. 6 clearly show the effects
of particle mass loading on the probability distribution
functions of C in the regimes where C is large. It appears
that particle mass loadings greater than 100 are rare in
turbulent flows.
The properties of the cascade rest on the physics of our
DNS simulations, and we speculate that two separate ef-
fects are involved. First, particle mass loading dampens
fluid motions of all types, decreasing vorticity stretching
and all other forms of ongoing eddy bifurcation which are
needed to produce intermittency. Second, as a byproduct
of this, particle mass loading may alter the Kolmogorov
timescale locally and shift the most effectively concen-
trated particle Stokes number St to a larger value than
that characterizing particles already lying in the local
volume, reducing the probability of preferentially con-
centrating the local particles any further.
Caveats and Future Work:
As described in section II, our multiplier distributions
were taken from the most numerous cells, with the largest
intermittency, which are at the smallest scales possible
(furthest from the forcing scale). At Reynolds numbers
accessible to DNS, a true inertial range is only beginning
to appear, and while, sampling at the smallest spatial
scales possible, we are as closely approaching the asymp-
totic values within the true inertial range as possible,
where level-independence has been demonstrated in the
past [21], it is possible that our values are subject to
inaccuracy by virtue of being sampled too close to the
dissipation scale. Any such inaccuracy will affect our
cascade results quantitatively but not qualitatively. As
computer power increases, it would be a sensible thing
to continue experiments like these at higher Reλ.
A more general model that treats enstrophy and strain
as independent cascading scalarsmight allow for a higher-
fidelity particle concentration cascade, since C is known
to be linked to the difference between these two scalars
[25] (the so-called second invariant tensor II). However,
such an effort would introduce further complexity of its
own, as II is no longer positive definite. We consider the
development of such a model a suitable task for future
work.
APPENDIX
We used an Eulerian scheme developed by Dr. Alan
Wray to solve the coupled set of fluid/particle equations
used in this study. This was done to maximize the
computational efficiency of the calculations and, more
importantly, to accurately evaluate multipliers over the
wide range of particle concentrations and enstrophies ex-
pected. In this study the effects of particle collisions and
external forces on the particles (e.g., gravity) are not con-
sidered. The turbulence is spectrally forced at k =
such that moments of the Fourier coefficients of the force
satisfy isotropy up to the fourth order. The instanta-
neous Navier-Stokes equations describing the conserva-
tion of mass and momentum for an incompressible fluid
∇ ·U = 0 (A.1)
+(U · ∇)U = −
+ν∇2U−α
(U −V) (A.2)
where U is fluid velocity, V is particle velocity, ρf and
ρp are the fluid and particle mass densities, ν is fluid
viscosity, P is pressure, and α is the inverse of the particle
gas drag stopping time τp.
The compressible equations for the particles are
+∇(ρpV) = D∇2ρp (A.3)
∂(ρpV)
+∇(ρpVV) = νp∇2(ρpV)+αρp(U−V) (A.4)
where νp is a “particle viscosity”, and D is a “particle
diffusivity”. The particle diffusivity and viscosity terms
numerically smooth out particle mass and momentum,
alleviating the formation of steep gradients of ρp that
can lead to numerical instabilities eg. [31].
The right hand sides of Eqs. A.2 and A.4 contain phase
coupling terms which are linearly dependent on (U−V).
The linear form of the coupling follows from the assump-
tions that the particle size is much less than η, and that
the material density of the particles is much greater than
ρf [2]. Additional contributions to the particle-gas cou-
plings involving pressure, viscous and Basset forces [29]
have not been added since they are expected to be weak
in our size regime of interest. The particle field is intro-
duced with a constant mass density and an initial veloc-
ity given by the local gas velocity in a field of statisti-
cally stationary turbulence. All runs are continued until
the particle statistics (RMS of conentration distribution)
have equilibrated.
The particle Stokes number St is defined relative to
the Kolmogorov time scale τη as St = τp/τη, and Φ =
Mp/Mf is the global mass loading, whereMp and Mf are
the total mass of particles and fluid respectively. In this
study ρf , St, and Φ are set to unity, D/ν = 0.01, and
νp/ν = 0.1. Explicitly setting St = 1 guarantees that
the particles are preferentially concentrated. When Φ is
unity, ρp is a surrogate for the local mass loading or local
concentration factor C. The values of νp and D minimize
the diluting effects of numerical particle diffusion while
preventing numerical blowups; their values were deter-
mined from a set of DNS runs in which their values were
decreased systematically until numerical instabilities set
Eqs. A.1 - A.4 are solved using psuedo-spectral meth-
ods commonly used to solve Naviers-Stokes equations for
a turbulent fluid. The Fast Fourier Transform (FFT) al-
gorithm is used to efficiently evaluate the dynamical vari-
ables U, V and ρp on a 3D uniform grid of computional
nodes with periodic boundary conditions. The computa-
tional algorithm is parallelized using MPI and is written
in Fortran 90. All runs for this study were executed on
SGI Origins supercomputers with up to 1024 processors.
Enstrophy is defined as
(∂iUj − ∂jUi)2 (A.5)
where i, j are summed over the three coordinate dimen-
sions of U.
The local spatial average of a scalar over a sample vol-
ume is estimated as,
Fidv (A.6)
where Fi is the scalar’s value on computational node i
centered within a cube of volume dv and the sum is over
all n nodes covering the sample volume. We normalized
this average by the global average value to get a quan-
tity that measures the scalar’s local value relative to its
mean. In this paper C and S will denote normalized spa-
tial averages of particle concentration and enstrophy over
cubes 3η on a side.
ACKNOWLEDGMENTS
We are very grateful to Dr. Alan Wray for providing
the 3-D code and for useful comments on its use. We
thank Robert Last for parallelizing the cascade code on
the SGI Origins 3000. We also would like to thank the
consultants and support staff at the NAS facility for pro-
viding invaluable assistance, and the Science Mission Di-
rectorate of NASA for generous grants of computer time.
We thank Prof. K. Sreenivasan for several helpful con-
versations in the preliminary stages of this project and
the internal reviewers Drs. Alan Wray and Denis Richard
for their suggestions for improving the manuscript. This
research has been made possible by a grant from NASA’s
Planetary Geology and Geophysics program.
[1] K. D. Squires and J. K. Eaton, Phys. Fluids A 2, 1191
(1990).
[2] K. D. Squires and J. K. Eaton, Tech. Rep. MD-55, Stan-
ford University (1990).
[3] K. D. Squires and J. K. Eaton, Phys. Fluids. A 3, 1159
(1990).
[4] J. D. Kulick, J. R. Fessler, and J. K. Eaton, J. Fluid
Mech. 227, 109 (1994).
[5] B. Castaing, Y. Gagne, and E. J. Hopfinger, Physica D
46, 177 (1990).
[6] S. P. G. Dinavahi, K. S. Breuer, and L. Sirovich, Phys.
Fluids 7, 1122 (1995).
[7] P. Kailasnath, K. R. Sreenivasan, and G. Stolovitzky,
Phys. Rev. Lett. 68, 2766 (1992).
[8] A. Vincent and M. Meneguzzi, J. Fluid Mech. 225, 1
(1991).
[9] A. Pumir, Phys. Fluids 6, 2071 (1994).
[10] E. Lamballais, M. Lesieur, and O. Métais, Phys. Rev. E
56, 6761 (1997).
[11] J. Jiménez, A. A. Wray, P. G. Saffman, and R. S. Rogallo,
J. Fluid Mech. 255, 65 (1993).
[12] G. He, S. Chen, R. H. Kraichnan, R. Zhang, and Y. Zhou,
Phys. Rev. Lett. 81, 4636 (1998).
[13] N. Cao, S. Chen, and K. R. Sreenivasan, Phys. Rev. Lett.
76, 616 (1996).
[14] R. C. Hogan, J. N. Cuzzi, and A. R. Dobrovolskis, Phys.
Rev. E 60, 1674 (1999).
[15] E. Balkovsky, G. Falkovich, and A. Fouxon, Phys. Rev.
Lett. 86, 2790 (2001).
[16] L. F. Richardson, Weather Prediction by Numerical Pro-
cess. (Cambridge University Press, Cambridge U.K.,
1922).
[17] A. N. Komolgorov, J. Fluid Mech. 13, 82 (1962).
[18] T. Chiueh, Chin. J. Phys. 32, 319 (1994).
[19] M. Gorokhovski, Tech. Rep., Center for Turbulence Re-
search, Annual Research Briefs (2003).
[20] C. Meneveau and K. R. Sreenivasan, Phys. Rev. Lett.
59, 1424 (1987).
[21] K. R. Sreenivasan and G. Stolovitzky, J. Fluid Mech.
379, 105 (1995).
[22] E. A. Novikov, Phys. Fluids A 2, 814 (1990).
[23] E. A. Novikov, Phys. Rev. E 50, R3303 (1994).
[24] J. Jiménez, J. Fluid Mech. 409, 99 (2000).
[25] M. R. Maxey, Phys. Fluids 30, 1915 (1987).
[26] K. D. Squires and J. K. Eaton, J. Fluid Mech. 226, 1
(1991).
[27] A. M. Ahmed and S. Elghobashi, Phys. Fluids 13, 3346
(2001).
[28] J. K. Eaton and J. R. Fessler, Int. J. Multiphase Flow
20, Suppl., 169 (1994).
[29] S. Elghobashi and G. C. Truesdell, Phys. Fluids A 5,
1790 (1993).
[30] U. Frisch, Turbulence (Cambridge University Press,
Cambridge, U.K., 1995), chap. 8.
[31] A. Johansen, A. C. Anderson, and A. Brandenburg, As-
tron. Astrophys. (2004).
