Description

On the onset of high-Reynolds-number grid-generated wind tunnel turbulence
By L. MYDLARSKI AND Z. WARHAFT
Using an active grid devised by Makita (1991), shearless decaying turbulence is studied for the Taylor-microscale Reynolds number, Ri, varying from 50 to 473 in a small (40 x 40 cm2 cross-section) wind tunnel. The turbulence generator consists of grid bars with triangular wings that rotate and flap in a random way. The value of Rj, is determined by the mean speed of the air (varied from 3 to 14 m s-’) as it passes the rotating grid, and to a lesser extent by the randomness and rotation rate of the grid bars. Our main findings are as follows. A weak, not particularly well-defined scaling range (i.e. a power-law dependence of both the longitudinal (u) and transverse (v) spectra, – Fll(kl) and F22(kl) respectively, on wavenumber – – k l ) first appears at R, 50, with a slope, n1, (for the u spectrum) of approximately 1.3. As Ri was increased, n1 increased rapidly until Ri. 200 where n 1.5. From there on the increase in n1 was slow, and even by Rl = 473 it was still significantly below the Kolmogorov value of 1.67. Over the entire range, 50 < Ri < 473, the data were well described by the empirical fit: n1 = i(1 – 3.15RT2’3). Using a modified form of the Kolmogorov similarity law: F1 1 (kl) = C1,~2/3k,5’3(kly)5/3-nl where E is the turbulence energy dissipation rate and y is the Kolmogorov microscale, we determined a linear dependence between nl and C,,: C1, = 4.5 – 2.4nl. Thus for n1 =- 5/3 (which extrapolation of our results suggests will occur in this flow for Ri. 104),C1. = 0.5, the accepted high-Reynolds-number value of the Kolmogorov constant. Analysis of the p.d.f. of velocity differences Au(r) and Au(r) where r is an inertial subrange interval, conditional dissipation, and other statistics showed that there was a qualitative difference between the turbulence for R, < 100 (which we call weak turbulence) and that for R, > 200 (strong turbulence). For the latter, the p.d.f.s of Au(r) and Av(r) had super Gaussian tails and the dissipation (both of the u and u components) conditioned on Au(r) and Au(r) was a strong function of the velocity difference. For Ri. < 100, p.d.f.s of Au(r) and Au(r) were Gaussian and conditional dissipation statistics were weak. Our results for R;. > 200 are consistent with the predictions of the Kolmogorov refined similarity hypothesis (and make a distinction between the dynamical and kinematical contributions to the conditional statistics).
They have much in common with similar statistics done in shear flows at much higher R)., with which they are compared.
1. Introduction
We address here, from an experimental viewpoint, the way in which the character- istics of turbulence, particularly its spectrum and probability density function, evolve
332 L. Mydlarski and 2. Warhaft
The central postulate of the Kolmogorov (1941a, b) phenomenological theory of turbulence is that if the Reynolds number is sufficiently high, there is a range of scales that are locally isotropic. Here there exists a scaling region, known as the
Fll(kl) = C 1 ~ ~ / ~ k , ~ ’ ~ ,
F22(kl) = F33(kl) = C2E 2/3k-5/3 1 . (2)
Here Fll(kl), F22(kl) and F33(kl) are respectively the streamwise power spectra of the longitudinal velocity fluctuations, u, and the transverse velocity fluctuations, u and w, kl is the longitudinal wavenumber, E is the turbulence energy dissipation rate per unit mass (defined below) and C1 and C2 are the Kolmogorov constants (Monin & Yaglom
1975). Since the Kolmogorov theory assumes local isotropy, it can be deduced from (1) and (2) on purely kinematical grounds that C2/C1 = 4/3.
For equations (1) and (2) to hold, the Reynolds number must be high. Yet it is well known that three-dimensional turbulence, exhibiting a multiplicity of scales, can exist at Reynolds numbers that are too low to have the -5/3 scaling region. Thus traditional grid-generated turbulence (e.g. Comte-Bellot & Corrsin 1971) with an RA of around 50 does not have a scaling region, but is sufficiently well developed so that the dissipation rate can be estimated by the scaling law
E = A ( u ~ ) ~ / ~ / G . (3)
Here, G is the turbulence integral scale (defined below), A is a constant close to unity and the angle brackets denote averaging. This law implies a cascade since the dissipation (a small-scale phenomenon) is determined from the large-scale energetics, and the viscosity does not appear explicitly. It has been shown that even in low- Reynolds-number grid turbulence (see for example Sirivat & Warhaft 1983, figure 5) E determined from (3) is close to the model-independent estimate (assuming only isotropy at the dissipation scales) given by
E = 15v((du/d~)~) (4)
where du/dx is the streamwise derivative of the longitudinal velocity fluctuations. Does such low-Reynolds-number turbulence differ qualitatively from that at high-
Reynolds-number ? And if so, what are discernible differences? Finally, at what R:, does the low-Reynolds-number turbulence, without a scaling region, undergo transition to high-Reynolds-number turbulence?
High-Reynolds-number grid-generated wind tunnel turbulence 333
In order to answer these questions, it appears that the Reynolds number should be varied from low to high values in a systematic way using the same type of initial and boundary conditions. Moreover, ideally the turbulence should be shearless and homogeneous, and preferably isotropic at all scales. If it is only locally isotropic (with large-scale anisotropy imposed by the boundary conditions), then below a certain Reynolds number the whole spectrum will be anisotropic (and hence it will be affected by its boundary conditions) since there will not be a broad enough wavenumber range to allow the development of local isotropy. The recent boundary layer measurements of Saddoughi – & Veeravalli (1994, referred to herein as S & V), indicate that not until R, 600, a very high value in terms of laboratory flows, does a significant locally isotropic region of the spectrum develop. Below this the spectrum is affected by the Reynolds shear stress up to dissipation wavenumbers and there is little, if any, region of local isotropy. While these shear flow measurements clearly confirm Kolmogorov’s hypothesis that at sufficiently high Reynolds number a locally isotropic region does exist, they are unable to examine the characteristics of isotropic turbulence at lower Reynolds numbers. To do this, grid-generated turbulence is ideally suited since it has no shear (and hence no Reynolds shear stress spectrum) and although the large-scale field is not quite isotropic, it is much closer to isotropy than that of shear flows (see 63).
The traditional method of generating (approximately) isotropic turbulence is by means of a grid: a rectangular array of (usually square sectioned) bars placed at the entrance of the test section (Batchelor 1953; Comte-Bellot & Corrsin 1966). As mentioned above, the Rj. is low, typically in the range 50 to 150, although there are some isolated results at high R;.t. The reason for this is that grids produce low-intensity turbulence; the ratio of the r.m.s. to mean velocity, u/U, is generally 2% or less in the downstream region where the flow has reached isotropy. On the other hand, high R, can be easily obtained in shear flows because of their high intensities (jets, for example, have u / U of around 25%, e.g. Townsend 1976). In order to produce high-R;. grid turbulence, the grid mesh must be large and so too must be the mean speed. The former requires extremely large tunnels since the tunnel must be many meshes wide to provide isotropy, the latter gives rise to probe resolution problems and to unwanted compressibility effects.
? Kistler & Vrebalovich (1966) achieved an R1 of over 500 in a massive wind tunnel built for aircraft research. Their results show a clear inertial subrange. However, there was no systematic
variation of R; and many of the statistics relevant to the present work were not computed. The tunnel was dismantled soon after the measurements were made.

334 L. Mydlarski and Z. Warhaft
trum showed well over one decade of inertial subrange, of slope close to -5/3. The only slight drawback was that the anisotropy, u/v, was 1.22. This is a lit- tle higher than that of conventional grid turbulence which has values around 1.1 (Comte-Bellot & Corrsin 1967). This aspect, which we will see does not appear to affect the inertial-range dynamics, does preclude accurate determination of the three-dimensional spectrum from the u and u components, assuming isotropy of the large-scales.
Briefly, Makita places triangular agitator wings on each grid mesh. A regular pulse generator rotates the grid bars (and hence the wings) and at the same time the motor is fed a random pulse which reverses the rotation of the grid bar. Thus the wing, always in rotational motion, reverses its direction randomly, providing a flapping motion. Each bar is separately controlled, providing random flapping between bars. The resultant turbulence has high intensity (up to 15%) and an integral scale larger than the mesh (rather than smaller produced by conventional grids). Both act to produce the high RA.
We decided to build a Makita-style active grid in order to study scalar mixing at high-Reynolds-numbers. Our preliminary experiments led us to explore the questions raised above concerning the evolution of grid turbulence with Reynolds number, and this is the subject of the present work. Results on scalar mixing are in preparation.
The outline of the work is as follows. After describing the apparatus in $2 and the grid performance in 53, we examine the nature of the velocity spectra over the Rl range 50 to nearly 500 in $4. Emphasis is placed on how the magnitude of the scaling exponent, n, and thus the Kolmogorov ‘constant’, varies with Ri. We then turn to an examination of the p.d.f. and conditional dissipation of velocity differences determined over inertial subrange intervals. These results are discussed in terms of the internal intermittency of the flow and are related to recent investigations in shear flows.
2. Apparatus
The experiments were conducted in our vertical wind tunnel which has a 40.65 x 40.65 cm2 cross-section and is 4.5 m, long (see Sirivat & Warhaft 1983 for a sketch). At the entrance of the test section, we replaced the conventional static grid with an active grid. Our design follows that of Makita (1991) although there are some variations. A sketch with dimensions is shown in figure l(a) and a photo of the grid in the tunnel test section is shown in figure l(b). The basic operation of the grid has been outlined in the Introduction. Here we provide its specifications. The mesh spacing, M , between the grid bars is 5 cm, providing a tunnel cross-section of 8 x 8 M. Each of the 14 grid bars is a 0.64 cm diameter aluminium rod with 0.38 mm thick aluminium wings, attached as shown in figure l(a). Every rod is independently driven by a Superior Electric 5 W D.C. synchronous stepping motor (shown in the photo, figure lb), with 200 steps per revolution. The rotation speed of the rods could be varied up to a maximum speed of about 4 revolutions per second
High-Reynolds-number grid-generared wind tunnel turhulence 335
(r.p.s.), but for nearly all of the measurements reported here it was set to 2 r.p.s. (see $3). A random signal was fed into the motor to alter its rotation direction. The average time between switching was of the order of a rotation period. Thus sometimes a bar would undergo more than a full cycle before switching while other times it would only undergo a fraction of a cycle. The appearance of the whole grid when operating in this random mode was of a shimmering nature since all of the 14 rods were acting independently. (Note that all the wings are in the same plane along a particular grid bar. This aspect could also be randomized but we have not done so.)
The tunnel speed was varied from 3 m s-‘ to 14 m s-l giving a variation in R, from 99 to 473. (Some supplementary measurements at Ri, = 50 and 100 were done using conventional static grids.) The rate of rotation of the grid bars only slightly
affected the turbulence level or R,. More significantly, it did produce a small spike in the spectrum at twice its rotation frequency (see $3).
To obtain good horizontal homogeneity in the velocity field, the grid had to be carefully tuned. Initially there was a large velocity deficit near the walls. Holes were drilled in the wings adjacent to the walls to reduce this. (These ‘half‘ wings were attached to the frame of the grid and did not rotate, see figure la.) Notice that in figure l(b) there appear to be holes also in the rotating wings. All of these (excepting those on the grid bars closest to the walls (figure la)) were taped over for the experiments to be described below. The mean velocity profiles will be discussed in the next section.

High-Reynolds-number grid-generated wind tunnel turbulence 337
N – 0.6
N ,-.
$ 0.4 2 0.2
N ,-. v s
0.2
0 10 20 30 40
Distance from wall (cm)
FIGURE 2. Mean (U), longitudinal and transverse r.m.s. profiles, (u2)’I2 (solid squares) and (u2)”* (open squares) respectively, and the cross-correlation coefficient between u and u,p,,. RJ, = 275, random mode.
3. The grid performance
Before describing the main results concerning the evolution of the spectrum and
p.d.f. with Reynolds number, we will outline the performance of the active grid, paying particular attention to its departure from ‘ideal’ grid-generated turbulence.
Cross-stream profiles of the mean speed, U , the root-mean-square longitudinal ve- locity and transverse velocity, (u2)’/* and (v2) l/* respectively, and the cross-correlation coefficient between u and v, puo(- ( u ~ ) / [ ( u ~ ) ~ / * ( v * ) ~ / ~ ] where (uv) is the kinematic
Reynolds stress) are shown in figure 2 for the grid run in random mode at 6.9 m s-l (& = 275). There is good homogeneity in U for nearly 20 cm in the transverse direction and there is no Reynolds stress for even a larger core than this. The r.m.s. u and v profiles curve in opposite directions giving quite strong anisotropy as the walls are approached. In the central core region ((u~)/(v~))’/~ is about 1.21. This result, which is consistent with that of Makita (19911, is higher than that of static grid experiments which yield values of 1.1 or slightly less (e.g. Sirivat & Warhaft 1983). We will return to this departure from isotropy in a moment. The profiles in figure 2 are typical of those done at other Reynolds numbers.
Figure 3 shows the decay of the longitudinal velocity variance and third moment as a function of x / M . Because of the large mesh (M = 5 cm) and also because the test section had to be modified to accommodate the active grid, the test section was only 80M in extent. Nevertheless, the form of the variance decay law is similar to that observed in conventional grid-generated turbulence (e.g. Sirivat & Warhaft
1983). The decay law for U = 6.7 m s-I was (u2)/U2 = 1,23(x/M)-‘,*’. Nearly all of 338 L. Mydlarski and 2. Warhaft
X l M
FIGURE 3. A typical longitudinal variance decay (circles) and third moment decay (squares). U = 6.7 m s-‘, R?, = 319 (at x / M = 68). For these measurements, done in the random mode, the bar rotation speed was 1 r.p.s. This produced a slightly higher RA (for the same mean speed) than did the higher rotation speed of 2 r.p.s., used for all the other measurements reported in this paper.
P.d.f.s of the u and u fluctuating components of velocity are shown in figure 5, for U = 7.1 m s-‘ (R?b = 262, random mode). For conventional grid turbulence, these

340 L. Mydlarski and Z. Warhaft
–
–
h h
&-
I I I l l l l l l I 1 1 1 1 1 I l 1 I I IIIlld
10-1 100 10’ lo2 103 104
k, (m-9
FIGURE 4. Longitudinal velocity spectrum, Fil(kl), solid line, and transverse spectrum, F22(kl), dashed line. RA = 262, random mode. The insert shows klFll(kl) (solid line) and klF22(kl) (dashed line); kl is the longitudinal wavenumber = 27c:f/U.
are Gaussian with the skewness S,(= (u3)/( ( u ~ ) ) ~ ‘ ~ ) and kurtosis K,(= (u4)/(u2)*) having values of 0 and 3 respectively. Here there are small departures: S, = 0.19 and
K , = 3.17. For u, the skewness and kurtosis, S, and K,, were respectively 0.03 and 3.41. The relatively high turbulence intensity (table 1) and its rapid decay produces a small divergence of the turbulence kinetic energy, a/axj( $ (uiuiuj)). Although we determined this to be less than 5% of the turbulence energy dissipation rate E, the existence of the triple moment induces a slight skewness in u. The experimentally derived law for the third moment is also plotted in figure 3. It is (u3)/U3 = 2 . 5 8 ( ~ / M ) – ~ . ~ ~ for U = 6.7 m s-‘. This in turn possibly produces the small departure from the Gaussian state for u via higher-order terms in the equation for v4. Of course the skewness in u must be zero by symmetry and our results confirm this to a high degree. Although the large-scale velocity field is not strictly Gaussian we will show that its departure from the Gaussian state is very small compared to the highly non-Gaussian statistics of velocity differences and derivatives determined in the inertial subrange.
We now return to the departure from isotropy observed in the r.m.s. u and u profiles of figure 2. As noted in the Introduction, the anisotropic effects of shear flows are evident in the inertial subrange as well as at large-scales. Thus it is important to determine whether the anisotropy of this flow (which is much less than that observed in shear flows) affects the inertial subrange. Figure 6 shows the coherence between u and u in the laboratory coordinates and in transformed coordinates, with a 45″ rotation. Here u’ = (u + u)/$ and u’ = (u – v)/$. In the laboratory coordinates, there is no coherence between u and u while with the – 45″ rotation there is coherence at low wavenumbers but this drops to zero by klq 3 x (The cross-correlation coefficient between u’ and u’ was 0.21.) Also shown in figure 6 is the u spectrum compensated to produce a plateau in the scaling region (i.e. the inertial subrange, see 54.1). It is clear that the coherence between u and u is essentially zero by the time the inertial subrange begins. We note that the 45″ rotation produced the largest

346 L. Mydlarski and Z . Warhaft
where nl and n2 are the slopes of the scaling region for the u and u spectra respectively and C1, and C2* are now Kolmogorov variables: both C, and n are functions of Rl (and as Rl -+ m,nl,n2 –+ 5/3 and C, + C, equations (1) and (2)). In figure 9(b) we have plotted C1, = Fll(k)~-~/~kl’r”~-~/~ vs. klq for the four spectra of figure 9(a). The value of nl, which varied from 1.40(R~ = 99) to 1.58(Rl = 448) was determined by trial and error such that the scaling region would be horizontal. Note that C1, decreases as Rl increases but even for the high-Rl case its value is approximately 0.7, well above the accepted high-Reynolds-number estimate of approximately 0.5 (Monin & Yaglom 1975; Sreenivasan 1995). Before we look at the evolution of C, we will examine the way the slope of the spectra varies with Rl.
Figure 10(a) is a summary of the best-fit scaling exponent, nl, for all the u spectra we measured over the range 50 < Rl < 473. For low RA, we have used some of our conventional (static) grid data. Both the synchronous and random modes of running the grid have been included. In figure 10(b) we have plotted the same data as lO(a) as a function of &. We have also included the low-Reynolds-number static grid data of Jayesh, Tong & Warhaft (1994, figure 7). Finally in figure 1O(c) we have plotted the slope of the transverse spectra, n2, as a function of RA. – –
The transverse, u, spectra (figure 1Oc) follow a similar trend to that of u but the slope at a particular R:, is less and it corresponded to a shorter scaling region. This is consistent with the spectra of S -& V (e.g. their figure 12). Note, here too, the rapid change of spectral slope up to Rl 200.
In figures 10(a) and 1O(c) we have fitted curves (which were empirically described well by -2/3 power laws) to the u and v spectral exponents, n1 and 112. Defining p 1 = 5/3 – nl and p2 = 5/3 – n2, we find
pi = 5.25RT 213 (15)
and p2 = 7.51RT 2/3 . (16)
These -2/3 power laws suggest that a 5/3 scaling region (to within a measurement – error of 0.01) will not occur until Rn lo4. We will show below that this is consistent with recent high-& measurements done in large wind tunnels and in the atmosphere.
We now return to the value of the Kolmogorov constant. Figure 11 shows a plot of C1, and C2* as a function of p . (The values of C, were determined from all of the measured spectra in the same manner as for the four spectra in figure 9(b).) The best fit line to C1, is
C1, = 0.51 + 2.39~1. (17)
Thus, when p1 = O(n1 = 5/3),C1, = C1 = 0.51. The generally accepted value of the three-dimensional Kolmogorov constant C is 1.5 (e.g. Monin & Yaglom 1975) and the one- and three-dimensional constants are related by C1 = 18C/55. Thus

348 L. Mydlarski and 2. Warhaft
2
c2*
r 1
CIS
1
0 0.2 0.4 0.6
p=5/3-n
FIGURE 11. C1. (open circles) and Cz. (closed circles) (equations (13) and (14)) plotted as a function of p = 5/3 – n where n is the slope of the respective spectrum. For Cl., the line is that of best fit. For C2*, the best fit line has been forced through Cl.(p = 0) x 4/3 at p = 0.
C1 = 0.49. Our extrapolated value of 0.51 is remarkably consistent with this value. We emphasize that equation (17) is a best fit.
The ratio of the Kolmogorov constants C2/C1 must be 4/3 if n1 = n2 = 5/3. In order to determine the best fit line for C2. we have used the value C2 = 4/3 x 0.51 for p2 = 0 yielding
C2. = 0.68 + 3.07~2. (18)
Finally, the dependence between C, and RAc an be determined by substituting equa- tions (15) and (16) into (17) and (18):
C1, = 0.51 + 12.6RT2I3, (19)
C2. = 0.68 + 23.1RY2l3. (20)
Figure 12 shows a plot of C1. and C2* as a function of RL with the fitted curves (equations (19) and (20)). We have also included C1, determined from the data of S & V. These data will be discussed in $5.
We now turn to the issue of isotropy. Assuming our flow is axisymmetric (the w statistics were determined to be the same as the v statistics), the isotropy is determined by the ratio of the u and v statistics. Figure 13 shows F22(kl)/F11(kl), i.e. the ratio of the v to u spectra, for RA = 50,262 and 377. At small wavenumbers, this ratio should be 0.5 if the energy-containing scales are isotropic. Our ratios are smaller than this, reflecting the large-scale anisotropy discussed in $3. In the region for which our spectra exhibit power-law scaling (signified by arrows in figure 13) there is, at the higher Reynolds numbers, a significant wavenumber region for which the ratio of the spectra is approximately 4/3, the value that must occur if both the u and u spectra have slopes of -5/3 and the turbulence is locally isotropic (this is the ratio C2/C1, equation (1)). The width of the 4/3 region is approximately one decade for the highest RA. For higher wavenumbers (the dissipation range), the ratio F22(kl)/Fll(kl) must increase as the spectra roll off exponentially and our spectra show this. (Noise – appears to be limiting the upper and lower spectra at klq 0.5.)

High-Reynolds-number grid-generated wind I , , , / , r –tunnel turbulence 349
I 10′
10′ 1 100
r – – – – _ _ _ _ _ – – – . – – – – – — – – _ _ _ – — – – – – -_ c2*
c,* 100 y :lo-‘
-.. . -3- L — -2 — – – –
10 , , I I , I , , I I I , I I I , I I I I I , , , .
Further insight into local isotropy in the inertial subrange can be gained by examining the ratio C2.;/C1*. From equations (13) and (14) this is
Our experiments (figures 10(a,c)) show that n2 # nl. A plot of nl – n2 is shown in figure 14. Since y12 is always less than nl, equation (21) indicates that F22(kl)/F11(kl) should increase with wavenumber and this is apparent for the lowest curve in figure 13. (We emphasize that this discussion holds for the scaling range only.) As RL increases, the slope of F22(kl)/FIl(kl) decreases (figure 13) since the difference between nl and nz diminishes. If n1 = 112 = n then it can be shown from the isotropic relation between Fll(k1) and Fzz(k1) (S. B. Pope, private communication), that
cz*/c1* = i(l+ n) (22)
(yielding 4/3 for – n = 5/3). For our high-Reynolds-number experiments n1 – 112 – 1.6 so Cp/Cl. 1.3 which cannot be experimentally differentiated from 1.33 (figure 13).
The ratio C2*/C1* (from the data of figure 12) is shown as a function of RI. in figure 15. Notice that C2./C1, (figure 15) is not equal to F22(kl)/Fll(kl) (figure 131, even for the high-Reynolds-number cases. Even here the small difference between nl and n2 is significant and the relationship between C2./C1.: and F22(k1)/F11(k1) must be determined from equation (21). Figure – 15 suggests that the Kolmogorov value of
4/3 will be reached by Rl 2200. S & V comment that only for their highest RA case (of RA = 1450) do they see appreciable local isotropy. Our extrapolation is broadly consistent with their results.
In spite of the lack of complete similarity even at our highest Ri, we have shown that all of our results are in agreement with the high-Reynolds-number scaling relationship / / q K Rr3/4 over the entire Reynolds number range (figure – 7b). As stated in the Introduction, this law assumes that E ( u ~ ) ~ / ~ / [ , i.e. that there is a cascade

352 L. Mydlarski and 2. Warhaft
0 200 400
R,
FIGURE 17. The value of I,, the vertical ordinate (figure 16) where the energy and dissipation curves cross. The best fit curve is I, = 1.2 R,”2.
A plot of I, us. RA is shown in figure 17. As the two examples of figure 16 have already shown, its value decreases with Rl, from about 0.11 at – R, N 100 to 0.06 for R, 450. A -1/2 power law fitted the data well, yielding
– I, = 1.20R~’/*. (23)
We note that above RA 300 the dependence of I, on R, is very weak. Here too, the variation in n was small (figure 10a,c). This is the region of strong or fully developed turbulence, to be more fully explored in the next section.
We conclude this section with a discussion of the third-order structure function,
( A u ( ~ ) ~ ) where Au(r) is the longitudinal velocity difference in the longitudinal direction. Unlike the second-order structure function (or spectrum), the form of ( A u ( ~ ) ~ ) is model independent. For high-Reynolds-number turbulence it follows directly from
High-Reynolds-number grid-generated wind tunnel turbulence 353
1
100 2
101 102 103
R,
FIGURE 18. The width of the Fll(kl) scaling region as a function of R,. The line increases as R;’*.
the Navier-Stokes equations (Kolmogorov 19416) that ( A u ( ~ ) ~ ) = -(4/5)~r. In figure 19 we have plotted – ( A ~ ( r ) ~ ) / ( r & ) as a function of r / q for various R,. At high Ri there should be a plateau at 4/5. Our highest R;. come close to this value although the plateau is not very wide. The quality of these structure functions is similar to that observed by Anselmet et al. (1984) for a jet at RA = 536. Our plots clearly show that for Ri < 150 there is significant departure from the four-fifths law and this is consistent with our observations of the nature of the spectra in this region, a region which we have referred to as weak turbulence.
4.2. The p.d$s, conditional statistics and jine-scale structure
While approximately Gaussian at the large-scales, turbulence at the small scales is strongly intermittent and non-Gaussian. If the Reynolds number is high enough, the intermittent structure should be observable not only in the dissipation range, but also in the inertial subrange (K62). It is the main purpose of this section to show how the intermittent structure evolves with Reynolds number by studying the p.d.f.s and conditional statistics of inertial subrange quantities. But first we will document some of the derivative (dissipation) statistics.
Figure 20 shows the p.d.f.s of du/dt and du/dt at Ri = 262. Both are super Gaussian. Their kurtoses (Ka,pt and K8u,81 defined as ( ( d ~ / d t ) ~ ) / ( ( d u / d t ) ~ ) * and similarly for dv/at) are 7.3 and 9.7 respectively. The du/dt p.d.f. must have zero skewness (Sarpt = { ( d ~ / d t ) ~ ) / ( ( d v / d t ) ~ ) ~ ’ ~ and similarly defined for duldt) by symmetry. (Its measured value was 0.05.) However, must have skewness from purely kinematical considerations (e.g. Wyngaard & Tennekes 1970). Its value for the p.d.f. of figure 20 is 0.49. Figures 21(a) and 21(h) show plots of SaUlat, Kauldr and Kaupt as a function of Ri. We have also included some high-Reynolds-number shear flow data taken from the compilation of Champagne (1978) as well as some low-Reynolds-number data from Tong & Warhaft (1994). The limited Reynolds number range and scatter of the present data are apparent from the graph, but the trend of the data (particularly for the kurtosis) as well as their magnitude is consistent with the shear flow data.

356 L. Mydlarski and 2. Warhaft
k, r
FIGURE 22. The dissipation spectra, Rn = 377, random mode. The solid and dashed lines are the u
and v dissipation spectra respectively.
RA. The scatter was around &-5% about an average value of 0.23 in the peak of the u dissipation spectrum over the range 100 < RL < 473.
We now turn to the p.d.f. and conditional statistics of the velocity difference
Au(r) = u(x + r ) – u(x), (24)
Av(r) = v(x + r ) – v(x), (25)
where r is an inertial range distance. The behaviour of the statistics of Au and Av are central to understanding models of internal intermittency and thus form the basis for modifications to K41. Beginning with K62, there has been a vast literature on this subject, known as the Kolmogorov refined similarity hypothesis (KRSH). The subject is still in ferment. Reviews are given by Kraichnan (1991) and Nelkin (1994). More recent work will be referred to where appropriate.
For our highest and lowest values of RA (473 and 50 respectively) we varied r (equations (24) and (25)) across the whole of the scaling range. We will show that the statistics of Au(r) (as well as Av(r)) are quite different for those two extreme Reynolds numbers. But in order to examine the trend of Au(r) and Av(r) as a function of Reynolds number, we chose a single value of r (denoted as r,) corresponding to the wavenumber k,(r = 271./kl) that is halfway between the beginning and end of the scaling range (on a linear plot) for each spectrum. Thus as the scaling range dilates with Reynolds number, r, remains in the same relative position within the scaling region. Figure 23 shows two spectra (Rl = 50 and 473) with arrows pointing to the beginning and end of the scaling range as well as to the wavenumber corresponding to r,. (For the low-Reynolds-number case, the scaling range is very slight indeed.) On some occasions we will also look at the statistics of Au(r) and Av(r) for larger r, chosen so that its value is midway between the beginning and end of the scaling region on the logarithmic plot. This is denoted as rb in figure 23.
Figure 24 shows p.d.f.s of Au and Av, normalized by their respective r.m.s. values for the r, separation. For the two lowest values of RL the p.d.f.s are close to Gaussian

358 L. Mydlarski and 2. Warhaft
:
0
0 i
0 I 1 II1llll I I 1 1 1 1 1 1 1 I I Illllll I I IIIIIII I I I , I , I ,
10-4 10-2 100
(y42 n) -Me- 7)
FIGURE 25. The kurtosis of the velocity difference as a function of r for Ri. = 100 (conventional grid) (circles) and 473 (active grid, random mode) (diamonds). The r normalization facilitates proper comparison of the two data sets.
Of course varying r changes the shape of the velocity-difference p.d.f.s: as r increases they should become more Gaussian in character, while for smaller r, as dissipation distances are approached, the p.d.f.s should look more like those of the derivatives. Thus as r increases, the kurtosis of Au,KA~, should decrease. The variation of Kau with r varied across the full extent of the inertial subrange is shown in figure 25 for RL = 100 and 473. For both cases the kurtosis approaches 3 as r approaches e. For intermediate values of r, K A ~ is lower for the low-Reynolds-number case and this is
consistent with figure 24.
We now study the conditional statistics of Au and Av. K62 hypothesised that for high RA, the energy dissipation rate averaged over a radius r, E,, (where r is within the inertial subrange) is related to Au(r) by
Au(r) = V ( ~ E , ) ‘ / ~ (26) where V is a stochastic variable independent of r and 6,. Thus the quantities Au(r) and E~ must be statistically dependent. In figure 26 we have plotted ( r , ~ : ~ ‘ ) ~ ‘ ~ and (r&?)1/3 conditioned on Au(r,). Here E,!~~(= 15~U-~((i?u/dt)~)) and 7 . 5U-2((Lb/8t)2)) ~ are (one-dimensional surrogates for the total dissipation) determined over a record of length ra, from which we also obtained Au(ra) from velocity difference between the start and end of the record. The data have been normalized by ( r , ( ~ ” ) ) ‘ / ~ (or (r,(~~’))’/~) and are plotted as a function of Au(ra) normalized by the r.m.s. value of Au(r,). For high RA, figure 26 shows that both ~ f f and E;: are statistically dependent on Au(r,):

3 60 L. Mydlarski and 2. Warhaji
.
-d
A
(r, €J”3ku (4Y3
FIGURE 27. The expectation of IAu(ra)l conditioned on E;: and E::. Lower curves (left-hand axis):
( IAu(ra) I l(r,&1)1’3) / IAUJ,~~. Upper curves (right-hand axis) : ( IAu(ra)l J(r,q!; )’I3) / IAUI,~~. The symbols are the same as for figure 26.
Recently Chen et al. (1995) and Thoroddsen (1995) have pointed out that a statistical dependence between Au(Y) and rq!’ must occur, even if KRHS does not hold. A correlation will occur between 8:’ and Au(r) on purely kinematic grounds. For a given Au(r), there exists a minimum possible value of the dissipation E,!’ which corresponds to a linear variation in u over the distance r ( ~ F ‘ l ~ i ~ = 15v(A~(r)/r)~). On the other hand, the existence of a statistical dependence between 8;’ and Au(r) suggests a dynamical contribution. For our high-R:, experiments (figure 26a, b) both 8:’ and 8;’ conditioned on Au(r) show statistical dependence on Au(r), giving strong support for KRHS.
In figure 27 we have plotted IAu(ru)l conditioned on €Fa1 and on This type of conditional statistic has been previously determined experimentally (Stolovitzky,
Kailasnath & Sreenivasan 1992; Praskovsky 1992; Zhu, Antonia & Hosokawa 1995; Thoroddsen 1995) and computationally (Chen et al. 1993; Wang et al. 1996) but has not been systematically studied over a wide range of Reynolds numbers. Note that there is not a simple correspondence between these and the conditional statistics of figure 26. Our results (figure 27) show that IAu(ru)I conditioned on E,!~’ shows similar statistical dependence on cia1 for high and low Reynolds numbers suggesting that kinematics are largely responsible. However, IAu(r,)l conditioned on E;:, while showing statistical dependence for high Rl, is nearly flat for the low-RJ. cases. Here the dynamical effects of the intermittency are evidently playing a role at high RA. This is consistent with figure 26(b). Thus the combination of our high- and low-Reynolds- number cases and the conditional statistics of both 8:’ and 8;’ on Au(r) enable us to separate the kinematical from the dynamical contributions to the statistical dependence: at low Ri. the dependence appears to be kinematical while at high Ri a significant dynamical effect is observed. We note that recent data of Thoroddsen (1995) at Ri = 208 shows independence of $'(r) and IAu(r)l. From this single observation, he casts doubt on KRSH. Our findings suggest his R:, was probably not

362 L. Mydlarski and 2. Warhaft
0.8 F
5 0.4
Q
0
yf rl
FIGURE 29. The cross-correlations (rE,)~/; as a function of r/q for RJ. = 473 (solid symbols) and R). = 100 (open symbols). Squares: p,Au(r)l, (rEf1)1/3 ; circles: pIAu+),, (rF;l)1/3.
r, and rb the cross-correlation p l A u ~ r ~ l , ~ r ~ ~ l ) l / 3 begins to become significant at – Ri – 100, reaching values of 0.35 and 0.25 for r, and r b respectively at RA 200. The value of ~~Au(,)~,(,~~1)l/; is close to zero for RA < 100 suggesting that the small correlation between (Au(r)( and cll for low Rl is kinematical only. (Notice that the cross- correlation of IAu(r)l and c2* is shifted down with respect to the cross-correlation of IAu(r)l and E” figure 28(a, b).) It is clear from these correlations that there is a significant change in the turbulence characteristics for 100 < Rn < 200. Over this interval there is a transition from weak turbulence (with almost Gaussian p.d.f.s of Au(r), figure 24, and weak statistical dependence between Au(r) and E, figures 26 and 27) to strong turbulence which exhibits non-Gaussian p.d.f.s and strong statistical dependence between Au(r) and E,. We observed in the section on spectra that there was also a distinct change in the nature of the spectrum over the same RA range.
The variation of plAU(r)l,(re,)l/; with r / q is shown in figure 29, for Ri = 100 and 473. Notice that for IAu(r)I correlated with there is a rapid decrease (to zero) as r/q increases. This variation has also been observed by Praskovsky (1992) and by Zhu et al. (1995). On the other hand, for the high-&, case, the cross-correlation between IAu(r)( and E~~ declines relatively slowly, remaining almost constant (with a value of around 0.3) within the inertial subrange. Praskovsky (1992) points out (contrary to his own findings) that p, should not depend on r within the inertial subrange for high-Reynolds-numbers. Our results show the importance of not using the same velocity component for Au(r) and E,.
High-Reynolds-number grid-generated wind tunnel turbulence 363
Finally, we have determined the intermittency exponent (K62) from the autocorre- lation of E (Praskovsky & Oncley 1994; Nelkin 1994):
Figure 30(a) shows pEE(r) for the one-dimensional surrogates E” and E ~ plotted as ~ , a function of r / q for RA = 100 and 377. From these (and estimates at other RJ,p was determined by compensation, i.e. by multiplying pEE(r) by r p and adjusting p so the curve became horizontal in the inertial range. These estimates of – – ,u are plotted as a function of R, in figure 30(b). Below R; 1 0 0 , ~ 0. This result is consistent with the nature of the velocity difference p.d.f.s (figure 24) and the conditional dissipation (figures 26 and 27) in this region. There is then a steep rise to a value of around 0.11 at R, NN 450. Measurements in very high-Reynolds-number flows show that p is approximately constant, with a value of approximately 0.2 (e.g. Praskovsky & Oncley 1994; Chambers & Antonia 1983). Evidently, we are not yet at a sufficiently high RA
to attain this value. This is consistent with the spectra ($4.1) which are still evolving – at R, 500.
We point out that the curves in figure 30(a) appear to asymptote to a non-zero value at r / q + GO (particularly at the lower R?. where the inertial range ends at a lower value of r/q). The asymptotic value can be calculated by decomposition of the dissipation into mean (E) and fluctuating ( E ’ ) components as follows:
((c(x) + e’(x))(~(x + r ) + E’(X + r)))
1-5lim PcAr) = ((E + e Y )
— (E(X)E(X + r)) + (E’(x)E’(x + r ) )
(2)+ (E’2)
The asymptotic limits of pCE(r), assuming E(x) = E(x + r ) (homogeneity) and limr+m, (F’(x)E’(x + r ) ) = 0 are indicated on figure 30(a).
The intermittency exponent provides the correction to K41, modifying the 5/3 scaling exponent to 5/3 + p / m where m is model dependent but is probably close to 9 (Nelkin 1994). Thus the slope of the spectrum is only slightly increased at high-Reynolds-numbers. Such a correction has not yet been observed. We will argue in the next section that this is because the RI. must be greater than lo4, a value even
larger than can be obtained from atmospheric measurements.
5. Discussion –
Our experiments have shown that by RA 200 the turbulence spectrum has a broad scaling range with a slope close to the Kolmogorov value of 5/3 and that the internal intermittency is well developed, showing interaction between the inertial and dissipation scales. As the Reynolds number is further increased the variation of the slope of the scaling exponent is only weakly dependent on RA.
How do these findings relate to other experiments done at much higher Reynolds numbers? An extrapolation of the scaling exponent, nl, (figure 10 and equation – (15)) suggests that a true 5/3 exponent will not occur until R?, lo4. At this value of Ri we estimate C, (equation (1)) to be 0.51. Recently Praskovsky & Oncley (1994) have examined data from the atmosphere and from large wind tunnels in the Ri range

High-Reynolds-number grid-generated wind tunnel turbulence 365
(near the dissipation range). We suggest that the scaling exponent has not quite reached 5/3. We have calculated C1, from their data (by replotting their spectra in the manner of our figure 9b). Its value is plotted in figure 12 and is consistent with the trend in our own data.
It appears, then, that an RA of approximately lo4 is required before the prediction of K41 is properly fulfilled in terms of the spectrum (equation (1)). Presumably, the value of R; must be significantly greater than this before the Kolmogorov (1962) intermittency correction (which steepens the slope even further) will be observable. We note that an Ri of lo5 in air implies (assuming an r.m.s. velocity of 1 m s-l) an integral length scale of 10 km! However, we believe that insight into the value – of intermittency will be gained from controlled experiments at R?. 1000. We are
currently building a large active grid to achieve this.
While the comparatively wide Reynolds number range of this study has been facilitated by new developments in grid design, it is curious that there appear to be no systematic studies of the variation of the turbulence spectrum with Ri. for other flows such as jets and boundary layers. These are flows for which a wide range of Reynolds numbers can be achieved with existing technology. It could be anticipated that the evolution of the spectrum for these flows will be different to that described here. For instance, in oscillating grid turbulence (a flow in which a conventional grid oscillates in a quiescent fluid) the 5/3 slope tends to occur at low Ri, and as soon as a scaling region is observed. This is possibly because the turbulence undergoes more straining by the time it reaches the measuring point (Hunt & Vassilicos 1991). Similarly, for a passive scalar in conventional grid turbulence, a scaling region occurs earlier than for the velocity field (Jayesh et al. 1994). The nature of turbulence statistics appears to be sensitive to their mode of generation even to quite high-Reynolds-numbers. It still remains to be seen whether true universality is ever achieved between flows of different origin, although the results presented here show at least for second-order moments there appears to be a consistency between grid turbulence and shear flows at high R;..
Our results have implications for the modelling of turbulence at moderate Reynolds numbers and – are therefore of practical importance. We have shown that even by R;. 500,C1, is 0.7, a value significantly different from the asymptotic value of 0.5. A value of R, = 500 is relatively high in terms of industrial fluid mechanics. For example, the Ri of turbulence in a commercial combustor is of this order. Large- eddy simulations usually model the inertial subrange assuming that its behaviour is independent of R;. We suggest that errors will accrue unless the Reynolds number dependence is accounted for in these and other modelling procedures.
6. Conclusions
Using the active grid design of Makita (1991) we have explored the evolution of grid turbulence with Reynolds number. We have shown that as Ri is varied from 50 to nearly 500 there is a qualitative change in the nature of – the turbulence. For R;, < 100 it has a weak scaling range (none at all below Rj, 50) and there is little or no effect of the intermittency on the inertial range. Nevertheless here, as for high – – R;.,e (u’)~/’// implying a fully developed cascade in which the dissipation can be inferred from large-scale quantities. We call this weak turbulence. Above Ri 200 the scaling region is well developed with an exponent close to 5/3. Here the internal intermittency is reflected in the inertial subrange. We call this strong turbulence and have shown that it has similar characteristics to the high-Reynolds-number turbulence 366 L. Mydlarski and Z. Warhap
studied in shear flows and in the atmosphere. We find that the change from weak to strong turbulence occurs in the range 100 < R), < 200. Here there is rapid change in the value of the scaling exponent as well as in the conditional and other statistics (figures 10, 24, 26 and 28).
Our experiment appears to be the first to bridge the gap between the low-Reynolds- – number (Rl 50) laboratory studies of grid-generated turbulence which have pro- & Corrsin 1971) and the high-Reynolds-number experiments (R). – lo3) done in the vided so much insight into its kinematics and dynamics (Batchelor 1953; Comte-Bellot
atmosphere and the oceans that have provided broad confirmation of K41 as well as estimates of the scaling and intermittency exponents, and the Kolmogorov constant (see Champagne 1978; and Nelkin 1994 for references to more recent experiments). We note that many turbulent flows of practical importance, such as in industrial mixers or machines, occur at the intermediate Reynolds numbers studied here.
Our results suggests that much can be learned about the behaviour of turbulence at high-Reynolds-numbers using a small wind tunnel. We note that the insights gathered from the atmospheric experiments as well as from the canonical experiments of Saddoughi & Veeravalli (1994) are costly and varying flow conditions is difficult. (Indeed control is impossible in the atmosphere.) Nevertheless, the RE. of 500 achieved in the present experiment is still too low to explore the detailed structure of high- Rl turbulence. We have shown that at this Rl there is still a significant evolution in the value of C1, and Cz* (figure 12), the scaling exponent is still less than the Kolmogorov value of 5/3 (figure 10) and p, the intermittency exponent, is still less than the observed high-Ra value of 0.2 (figure 30). We are presently constructing a larger Makita style grid with the aim of achieving an Ri of around 1000, to study statistics higher than those of second order, to which we have confined ourselves
here.
We once again thank Mr E. P. Jordan who expertly constructed the grid. Dr Chen- ning Tong carried out some preliminary experiments and Professors S. B. Pope and E. D. Siggia provided advice as well as encouragement. We sincerely thank them. The work was supported by the Department of Energy (Basic Energy Sciences) whom we also thank.
Note added in proof: Recently we have constructed a larger version of the active grid, obtaining an Rj, of 800. The new results are entirely consistent with the trends – reported here; for example, the slope of the u spectrum is 1.6, at Ra 800 and this is consistent with the trend in figure 10(a). These, and results of passive scalar mixing experiments in the same flow, are presently being prepared for publication.
REFERENCES
ANSELMET,F., GAGNE, Y., HOPFINGER,E. J. & ANTONIA, R. A. 1984 High-order velocity structure functions in turbulent shear flows. J. Fluid Mech. 140, 63-89.
BROWNEL,. W. B., ANTONIA, R. A. & CHUAL, . P. 1989 Calibration of X-probes for turbulent flow measurements. Exps. Fluids 7, 201-208.
CHAMBERS, A. J. & ANTONIA, R. A. 1984 Atmospheric estimates of power law exponents p and pe. Boundary Layer Met. 28, 343-352.
CHAMPAGNE, F. H. 1978 The fine-scale structure of the turbulent velocity field. J. Fluid Mech. 86, 67-108.
CHEN, S., DOOLEN, S. D., KRAICHNARN., H. & SHE, Z.-S. 1993 On statistical correlations between
Hr,~lz-Relnolds-numher grid-generated wind tunnel turbulence 367
velocity increments and locally avcraged dissipation in homogeneous turbulence. Phys. Fluids A 5, 458-463.
CHEN, S.: DOOLEN, S. D., KRAICHAKANR,. H. & WANG, L.-P. 1995 Is the Komogorov refined similarity relation dynamic or kinematic? Phys. Reu. Lett. 74, 1755-1 758.
COMTE-BFLLOTG. . & CORKSIN, S. 1966 The use of a contraction to improve the isotropy of grid-generated turbulence. J. Fluid Mech. 25, 657-682.
COMTE-BELLOT, G. & CoRRsiN, S. 1971 Siniple Eulerian time correlation of full and narrow-band velocity signals in grid-generated ‘isotropic’ turbulence. J. Fluid Mech. 48, 273-337.
GAD-EL-HAK, M. & COKRSIN, S. 1974 Measurements of the nearly isotropic turbulence behind a uniform jet grid. d. Fluid Mech. 62, 11 5-143.
HUNT, J. C. R. & VASSILICOS, J. C. 1991 Kolmogorov’s contributions to the physical and geometrical understanding of small-scale turbulence and recent developments. Proc. R. Soc. Lond. A 434, 183-210.
JAYESH & WARHAFT, Z. 1992 Probability distribution, conditional dissipation, and transport of passive temperature fluctuations in grid-generated turbulence. Phys. Fluids A 4, 2292-2307.
JAYESHT, ONG, C. & WARHAFT, 2. 1994 On temperature spectra in grid turbulence. Phys. Fluids 6, 306- 3 12.
KISTLER, A. L. & VREBALOVICH, T. 1966 Grid turbulence at large Reynolds numbers. J. Fluid Mech. 26, 37-47.
KOLMOGOROV, A. N. 1941a The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Dokl. Akad. Nauk. SSSR 30, 301-305 (referred to herein as K41).
KOLMOGOROV,A. N. 1941b Dissipation of energy in locally isotropic turbulence. Dokl. Akad. Nauk. SSSR 32, 16-18 (referred to herein as K41).
KOLMOGOROAV., N. 1962 A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at very high Reynolds numbers. J. Fluid Mech. 13, 82-85 (referred to herein as K62).
KRAICHNAN,R. 1991 Turbulent cascade and intermittency growth. J. Fluid Mech. 434, 65-78. LESIEUR, M. 1990 Turbulence in Fluids, 2nd Edn. Kluwer.
LING, S. C. & WAN, C. A. 1972 Decay of isotropic turbulence generated by a mechanically agitated grid. Phys. Fluids 15, 1363-1369.
MAKITA, H. 1991 Realization of a large-scale turbulence field in a small wind tunnel. Fluid Dyn. R ~ s .8, 53-64.
MAKITAH. . & SASSA, K. 1991 Active turbulence generation in a laboratory wind tunnel. Aduances in Turbulence 3 (ed. A. V. Johansson & P. H. Alfredsson). Springer.
MAKITA, H., SASSA. K., IWASAKI, T. & IIDU, A. 1987 Trans. Japan SOC. Mech. Engng B 53, 3180. (In
Japanese).
MONINA, . S. & YAGLOM, A. M. 1975 Statistical Fluid Mechanics, Volume 2. MIT Press.
NELKINM, . 1994 Universality and scaling in fully developed turbulence. Adu. Phys. 43, 143-181.
PRASKOVSKY, A. A. 1992 Experimental verification of the Kolmogorov refined similarity hypothesis. Phvs. Fluids A 4, 2589-2591.
PRASKOVSKY, A. A. & ONCLEY, S. 1994 Measurements of the Kolmogorov constant and intermittency exponent at very high Reynolds numbers. Phys. Fluids 6, 2886-2888.
SADDOUGHI,S. G. & VEERAVALLI, S. V. 1994 Local isotropy in turbulent boundary layers at high Reynolds number. J. Fluid Mech. 268, 333-372 (referred to herein as S & V).
SIRIVAT, A. & WARHAFT Z. 1983 The effect of a passive cross-stream temperature gradient on the evolution of temperature variance and heat flux in grid turbulence. J. Fluid Mech. 128, 323-346.
SREENivAsAN, K. R. 1991 On local isotropy of passive scalars in turbulent shear flows. Proc. R. Soc. Land. A 434, 165-182.
SREENIVASAKN. , R. 1995 On the universality of the Kolmogorov constant. Phys. Fluids 7,2778-2784.
STOLOVITZKY, G. KAILASNATH, P. & SREENIVASAKN., R. 1992 Kolmogorov’s refined similarity hypotheses. Phys. Rev. Lett. 69, 1178-1181.
THORODDSEN,S. T. 1995 Reevaluation of the experimental support for the Kolmogorov refined similarity hypothesis. Phys. Fluids 7, 691-693.
TONG,C . & WARHAFT, Z. 1994 On passive scalar derivative statistics in grid turbulence. Phys. Fluids
6. 2165-2176.