New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
On the existence and scaling of structure functions in turbulence according to the data

Contributed by Alexandre J. Chorin, January 23, 2006
Abstract
We sample a velocity field that has an inertial spectrum and a skewness that matches experimental data. In particular, we compute a selfconsistent correction to the Kolmogorov exponent and find that for our model it is zero. We find that the higherorder structure functions diverge for orders larger than a certain threshold, as theorized in some recent work. The significance of the results for the statistical theory of homogeneous turbulence is reviewed.
In 1941, Kolmogorov (1) formulated his famous scaling theory of the inertial range in turbulence, according to which the secondorder structure function, i.e., the function S _{2}(r) = 〈(u(x + r) − u(x))^{2}〉, scales as (εr)^{2/3}, where x, x + r are points in a turbulent flow field, u is the component of the velocity in the direction of r, ε is the mean rate of energy dissipation, r is the length r of r, and the angle brackets denote an average. A Fourier transform yields the Kolmogorov–Obukhov inertial range spectrum E(k) = Cε^{2/3} k ^{−5/3}, where C is a constant, and k is the wave number (2). These results are the lynchpins of turbulence theory, yet uncertainty lingers as to their general validity and the details of the derivation. Subsequently (see, e.g., ref. 3) it was claimed that this scaling result generalizes to structure functions of any order, i.e., to S_{p} (r) = 〈(u(x + r) − u(x)) ^{p} 〉 ∼ (εr)^{p/3} for any p > 0, but note that Kolmogorov and Obukhov themselves never went that far; Kolmogorov in ref. 4 gave only the further result for p = 3.
Shortly after the publication of this theory, Landau and Lifshitz (5) challenged its derivation on the ground that the rate of energy dissipation is intermittent, i.e., spatially inhomogenous, and cannot be treated as a constant. This observation fits in with the experimental observation that for p > 3, the exponents are smaller than the values given by the extended Kolmogorov theory, and it also is often claimed for the p = 2 the experimental value of the Kolmogorov–Obukhov exponent is larger than the prediction of this scaling theory. Obukhov (6) and Kolmogorov (7) themselves produced a “corrected” theory, and the current belief is that the theory has to be supplemented by “intermittency corrections”; much effort has been expended on the calculation of these corrections (see e.g., refs. 3 and 8).
A different analysis of possible corrections to the Kolmogorov exponents has been offered by Barenblatt and Chorin (9, 10). The structure functions depend on the Reynolds number R (for a definition of R suitable for the present case, see ref. 11) a bulk length scale L, and the mean rate of energy dissipation ε. Dimensional analysis yields S_{p} = (εr)^{p/3}Φ _{p} (r/L, R), where Φ _{p} is an unknown dimensionless function of two large arguments. If one makes a complete similarity assumption (see ref. 12) one finds Φ _{p} ∼ C for large arguments, and one recovers the Kolmogorov exponents above, but other assumptions are possible and indeed natural. In particular, an analogy with a successful scaling theory for turbulent boundary layers (13, 14) leads to the assumption Φ _{p} = C_{p} (R)(r/L)^{(α p/lnR)} where α _{p} is a constant and C_{p} (R) is a function of the Reynolds number R. This assumption yields structure functions S_{p} proportional to r ^{p/3+αp/lnR}, i.e., exponents that depend on R, and maybe even more important, a proportionality constant that depends on R and may well diverge as R → ∞ for p large enough. As pointed out in ref. 13, the dependence on 1/lnR is an assumption, and in reality may be weaker yet; such dependence could be hard to detect in experimental data.
This analysis is quite compatible with the Landau observation regarding intermittency, as was already pointed out in ref. 14: Turbulent flow is intrinsically intermittent, and indeed in the limit of infinite R, extremely intermittent, in the sense that the bulk of the vorticity occupies a negligible fraction of the available volume (15); this fact is a dominant fact about turbulence, and the notion of “intermittency corrections” is no more meaningful than a notion of “turbulence corrections” to turbulent flow. However, at finite R viscosity reduces the intermittency, and the scaling has to be corrected for “intermittency reductions” that depend on R. In related work (16), Barenblatt and Chorin made numerical estimates on the basis of an approximate theory of turbulence and came to the conclusion that, as R increases, the structure functions for p ≥ 4 diverge; similar conclusions had been reached by P. Constantin (personal communication) through mathematical considerations and by Mandelbrot (17) by a geometrical similarity argument. Thus, indeed C_{p} (R) would diverge as R increases for p ≥ 4.
In the absence of solutions of the Navier–Stokes equations, these various theories can only be checked by data from experiments, and this endeavor seems to be difficult at present. It occurred to us to interrogate further the available data about the inertial range by building a stochastic computer model of the inertial range, compatible with well accepted and reproducible data about their skewness (indicating the level of intermittency) and then examine the resulting structure functions. The field we construct is obviously not unique; however, any model that satisfies constraints imposed by the data is highly instructive and indicates what possibilities are worthy of further study. The results are striking enough to open new vistas for theory, as we discuss in Conclusions.
We thus construct on the computer a onedimensional (1D) Gaussian velocity field with a powerlaw, then modify it so that the velocity differences assume the skewnesses presented in Batchelor's book (18) on the basis of the data in ref. 19; the power law will be either the Kolmogorov–Obukhov spectrum or a modification that satisfies other constraints (see below). We use a 1D model so that we can perform the calculations with sufficient accuracy and think that this model is sufficient to carry our conclusions. We then calculate, when possible, the higherorder structure functions. The data in refs. 18 and 19 are tabulated there for various Reynolds numbers; we take the data that correspond to the largest R = 42,200; our model has a power law for all wave numbers and is thus inviscid. Note that by modern standards this value of R is not very large; we thus assume that the skewness does not vary drastically as R increases further. We first present the technical aspects of the construction, then the results, then provide a discussion.
Building a Model
We first explain how to sample effectively a homogenous Gaussian velocity field with a given spectrum. Our basic tool is a construction from Elliott et al. (20, 21). A homogeneous Gaussian random field is determined by its secondorder means The spectral representation of the Gaussian random field u is where w is a Wiener process. Expand dw in a complete orthonormal series φ _{m} where the γ _{m} form a sequence of independent random Gaussian variables. By substituting Eq. 3 into Eq. 2 we find where the coefficients c_{m} (x) are and ℱ is the Fourier transform. Elliott et al. (20, 21) proposed to use as basis of the decomposition, Fourier transforms of wavelets based on the Meyer wavelet; these wavelets are generated from the Meyer mother wavelet ϕ(k) by the wavelet relation where the index m refers to different scales (octaves) and n to different dilations. Using this prescription, Eq. 4 can be represented as where E _{m} ^{1/2} = 2^{−m/2} E ^{1/2}(2^{−m} k). In the particular case of the Kolmogorov–Obukhov energy spectrum, E(k) ∼ k ^{−5/3}, the Fourier transform of the c_{mn} coefficients represented in Eq. 7 becomes Similar formulas are obtained for every powerlaw spectrum. Defining f_{m} (2^{−m} x − n) = ℱ^{−1}(2^{−m/2}2^{−m/2} k^{−5/6}ψ _{mn} (k), the representation of the field given in Eq. 7 is The use of the summation (Eq. 9 ) requires a truncation in m and n that preserves accuracy. The truncation uses the spatial decay of the functions f _{m} (see ref. 20). It is convenient to center the summation around the term with the smallest (in magnitude) value of the argument (2^{−m} x − n). This centering can be done by defining for each term m of the outer sum the index n̄_{m} = 2^{−m} x and shifting indexes in the inner sum to n′ = n − n̄_{m} . The parameter m governs the scale; the truncation in m defines a frequency range, which can be modified by rescaling the distance x if necessary, to obtain a numerically convenient range −M ≤ m ≤ 0. The truncation in n (the number of dilations) defines the region where the support of f_{m} is concentrated; we suppose that this range lies between a pair of integers −N and N. We can then write More information can be found in refs. 20 and 22. A technical issue concerning the generation of Gaussian random variables should be pointed out. The number of random variables needed to sample the field scales as 2^{M+1} x, which imposes severe demands on the storage in the computer. To overcome this problem, a reversible pseudorandom number generation is necessary to sample γ _{mn} . Simple linear congruential generators with moduli around 2^{31} can exhaust their period in few minutes in a conventional PC, and the resulting poor distribution of the samples can dramatically bias simulation results for sample sizes much smaller than the period length. To overcome this problem, we used reversible multiple linear congruential generators with many long streams and substreams (23) that provide periods of ≈2^{191}.
We now explain how to modify the construction of the preceding section to impose on the sampled velocity field the nonGaussian characteristics observed in the experimental data.
We introduce nonGaussianity into the numerical experiment through the coefficients γ _{m} . In the previous paragraphs, the γ′_{m} s were Gaussian variables; now we construct a distribution for the γ′_{m} s that preserves the mean and variance of the original Gaussian distribution but with a skewness different from 0. The skewness of a random variable η is 〈η^{3}〉/〈η^{2}〉^{3/2}; it is zero for a Gaussian variable. We describe a transformation that maps a Gaussian variable on a nonGaussian variable with the same mean and variance but with a prescribed, nonzero skewness controlled by a single parameter. A Gaussian variable with mean zero and unit variance has the probability density
The following change of variables yields a skewed probability density function, with a negative skewness as in the cases we consider:
The resulting probability density function g(y) is obtained by calculating g(y) = p(x)/
Experimental results for the velocity field were obtained by Stewart in 1951 (19) and are also presented in ref. 18 (see Fig. 2). The skewness varies with scale, and this variation is easy to incorporate into our wavelet representation because the various wavelets describe motion on different scales. The parameter a(r) has been obtained by inversion of the formula (Eq. 13 ) for values of c _{3}(r) corresponding to the experiment for the highest value of R reported in ref. 18.
Results
Having constructed these fields, i.e., written a computer program that samples them, we proceeded to calculate their properties by Monte Carlo; we now list some of the results.
(i) First, we noticed that if the secondorder structure has the Kolmogorov form, and if the field is Gaussian, then all of the structure functions obey the Kolmogorov scaling. This result is quite natural and obvious and was well known to Kolmogorov (personal communication with G. I. Barenblatt), but we had never seen it stated in print. The obvious converses are not necessarily true.
(ii) It is “almost” a theorem that the thirdorder structure function scales like r (see ref. 4 as well as refs. 8 and 15); all one has to assume is an extremely likely bound on the rate of blowup of the dissipation as the Reynolds number grows. It is natural to ask then what form must the secondorder structure function take if the third order scales as r and the skewness is as observed. Let the secondorder structure function scale as r
^{2/3+α} and the thirdorder structure scale as r
^{ξ3}; in Fig. 3 we plot the computed values of ξ_{3} as a function of α. The relationship between them can be approximated by ξ_{3} = 1 −
(iii) The next question is maybe the most significant: Suppose the secondorder structure function has the Kolmogorov–Obukhov form and the skewness is as observed: what can one say about the higherorder structure functions? When we tried to sample these higherorder functions as above, we found that the variance of the Monte Carlo results was very large and did not go down as rapidly as one may expect as the number of trials went up, leading us to suspect that some moments of the velocity field diverged.
We therefore plotted the data for the variable x = S _{1}(r) = Δu at fixed values of r in log–log scale (see Fig. 4). The data for the nonGaussian field present a powerlaw tail that clearly contrasts with the exponential decay of the Gaussian field. To estimate the form of the decay, we construct the cumulative density function (CDF) P(x > X) of the variable x = S _{1}(r). We use CDFs rather than the density functions themselves because they fluctuate less. The behavior of the CDFs reveals the convergence, or lack thereof, of the structure function moments. In the limit as the number of samples n → ∞, 〈S_{p} 〉 _{r} = ∫(Δu) ^{p}f(Δu)d(Δu), where f is the probability density function, the derivative of the CDF. The CDF for the Gaussian field presents an exponential decay compatible with the finiteness of moments of any order. However, the CDF for the velocity differences in the nonGaussian case presents a powerlaw tail whose exponent determines the order of the moments that are finite. The tail of the CDF of S _{1} scales for small values as a power law with exponent −2.2 ± 0.01; however, the final part of the CDF, which contains the relevant information about the decay, is well described by an exponent −3.3 ± 0.05 within our experimental error. This result means that the CDF will decay approximately as a power law with exponent −4.3 ± 0.05 and then produces a threshold for the higherorder structure functions to converge at order p ≈ 3.3. We therefore conclude that, within our model, the structure functions of order ≥3.3 do not exist. In Fig. 5, we plotted the CDFs for several values of r, to show that the behavior we describe is at least approximately independent of r when r is small.
This finding does not of course show that the moments of the true velocity field in turbulence fail to exist, but only that the data contain a strong enough deviation from Gaussianity for these moments not to exist.
Finally, we provide a short discussion of sampling errors. All of the data presented have been obtained by the leastsquares method, except the data for the final part of the tail in Figs. 4 and 5, where we have also applied the maximum likelihood method to reduce the estimation error of the exponent (24), within a framework suggested in ref. 25. With 95% confidence, the errors are comparable with the thickness of the lines drawn in all cases. The regions of the sampling where aliasing problems appear in the fast Fourier transform have been discarded in the analysis.
Conclusions
Despite the nonuniqueness of the fields we have constructed, some conclusions can be drawn from the computations above, in particular the following.
(i) The original Kolmogorov–Obukhov −5/3 spectrum is consistent with the “exact” scaling of the thirdorder structure function even in the presence of intermittency; intermittency does not necessarily require that the spectrum be modified. As long as one considers only the lowerorder structure functions originally considered by Kolmogorov and Obukhov, our model provides no reason to modify their original scaling.
(ii) The Kolmogorov scaling of the structure functions and its extensions is exact at all orders for Gaussian fields; however, the data for real flows reveal enough departure from Gaussianity for structure functions of order higher than a threshold larger than three to diverge as the Reynolds number grows; if this happens, the higherorder exponents become strongly dependent on the Reynolds number, their scaling laws depend on the Reynolds number, and the proposals of Barenblatt, Chorin, and Prostokishin become eminently reasonable, as well as those of Mandelbrot and Constantin. Thus, within the limitations of our model, the complete similarity assumption, on which the extension of the Kolmogorov scaling rests, fails above the threshold but not below it; when it fails it should be replaced by an incomplete similarity assumption as in the papers cited above.
Acknowledgments
We thank Prof. G. I. Barenblatt for inspiration, discussions, and helpful comments and Prof. C. Stone for help with analysis of the data. This work was supported in part by National Science Foundation Grant DMS0410110 and by the Office of Science, Office of Advanced Scientific Computing Research, Mathematical, Information, and Computational Sciences Division, Applied Mathematical Sciences Subprogram of the U.S. Department of Energy, under Contract DEAC0376SF00098. A.A. also was supported by Direccion General de Enseñanza Superior of the Spanish Government Grant BFM200308258.
Footnotes
 ^{‡}To whom correspondence should be addressed. Email: chorin{at}math.berkeley.edu

Author contributions: A.A. and A.J.C. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.

Conflict of interest statement: No conflicts declared.
 Abbreviation:
 CDF,
 cumulative density function
 © 2006 by The National Academy of Sciences of the USA
References

↵
 Kolmogorov A. N.

↵
 Obukhov A. M.

↵
 Frisch U.

↵
 Kolmogorov A. N.

↵
 Landau L. D. ,
 Lifshitz E. M.
 ↵
 ↵

↵
 Monin A. S. ,
 Yaglom A.

↵
 Barenblatt G. I. ,
 Chorin A. J.
 ↵
 ↵

↵
 Barenblatt G. I.

↵
 Barenblatt G. I. ,
 Chorin A. J.
 ↵

↵
 Chorin A. J.

↵
 Barenblatt G. I. ,
 Chorin A. J.
 Giovannini A. ,
 Cottet G.H. ,
 Gagnon Y. ,
 Ghoniem A. F. ,
 Meiburg E.
 ↵

↵
 Batchelor G.

↵
 Stewart R.
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Breiman L. ,
 Kooperberg C. ,
 Stone C.
Citation Manager Formats
More Articles of This Classification
Physical Sciences
Related Content
 No related articles found.
Cited by...
 No citing articles found.