## Spectral analysis of structure functions and their scaling exponents in forced isotropic turbulence

W. D. McComb, 1, ∗ S. R. Yoffe, 2 M. F. Linkmann, 1 and A. Berera 1

1

SUPA, School of Physics and Astronomy, University of Edinburgh, James Clerk Maxwell Building, The King's Buildings, Edinburgh EH9 3JZ, UK 2 SUPA, Department of Physics, University of Strathclyde, John Anderson Building, 107 Rottenrow East, Glasgow G4 0NG, UK

## Abstract

The pseudospectral method, in conjunction with a new technique for obtaining scaling exponents ζ n from the structure functions S n ( r ), is presented as an alternative to the extended selfsimilarity (ESS) method and the use of generalized structure functions. We propose plotting the ratio | S n ( r /S ) 3 ( r ) | against the separation r in accordance with a standard technique for analyzing experimental data. This method differs from the ESS technique, which plots S n ( r ) against S 3 ( r ), with the assumption S 3 ( r ) ∼ r . Using our method for the particular case of S 2 ( r ) we obtain the new result that the exponent ζ 2 decreases as the Taylor-Reynolds number increases, with ζ 2 → 0 679 . ± 0 013 as . R λ → ∞ . This supports the idea of finite-viscosity corrections to the K41 prediction for S 2 , and is the opposite of the result obtained by ESS. The pseudospectral method also permits the forcing to be taken into account exactly through the calculation of the energy input in real space from the work spectrum of the stirring forces.

PACS numbers: 47.11.Kb, 47.27.Ak, 47.27.er, 47.27.Gs

∗ wdm@ph.ed.ac.uk

## I. INTRODUCTION

In this paper we revisit an old, but unresolved, issue in turbulence: the controversy that continues to surround the Kolmogorov theory (or K41) [1, 2]. This controversy began with the publication in 1962 of Kolmogorov's 'refinement of previous hypotheses' (K62), which gave a role to the intermittency of the dissipation rate [3]. From this beginning, the search for 'intermittency corrections' has grown into a veritable industry over the years: for a general discussion, see the book by Frisch [4] and the review by Boffetta, Mazzino and Vulpiani [5]. The term 'intermittency corrections' is rather tendentious, as no relationship has ever been demonstrated between intermittency, which is a property of a single realization, and the ensemble-averaged energy fluxes which underlie K41, and it is now increasingly replaced by 'anomalous exponents'. It has also been observed by Kraichnan [6], Saffman [7], Sreenivasan [8] and Qian [9] that the title of K62 is misleading. It in fact represents a profoundly different view of the underlying physics of turbulence, as compared to K41. For this reason alone it is important to resolve this controversy.

While this search has been a dominant theme in turbulence for many decades, at the same time there has been a small but significant number of theoretical papers exploring the effect of finite Reynolds numbers on the Kolmogorov exponents, such as the work by Effinger and Grossmann [10], Barenblatt and Chorin [11], Qian [9], Gamard and George [12] and Lundgren [13]. All of these papers have something to say; but the last one is perhaps the most compelling, as it appears to offer a rigorous proof of the validity of K41 in the limit of infinite Reynolds number. This is reinforced by the author's comparison with the experimental results of Mydlarski and Warhaft [14].

The controversy surrounding K41 basically amounts to: 'intermittency corrections' versus 'finite Reynolds number effects'. The former are expected to increase with increasing Reynolds number, the latter to decrease. In time, direct numerical simulation (DNS) should establish the nature of high-Reynolds-number asymptotics, and so decide between the two. In the meantime, one would like to find some way of extracting the 'signature' of this information from current simulations.

As is well known, one way of doing this is by Extended Self Similarity (ESS). The study of turbulence structure functions ( e.g. see van Atta and Chen [15] and Anselmet, Gagne, Hopfinger and Antonia [16]) was transformed in the mid-1990s by the introduction of ESS

by Benzi and co-workers [17, 18]. Their method of plotting results for S n ( r ) against S 3 ( r ), rather than against the separation r , showed extended regions of apparent scaling behavior even at low-to-moderate values of the Reynolds number, and was widely taken up by others e.g. Fukayama, Oyamada, Nakano, Gotoh and Yamamoto [19], Stolovitzky, Sreenivasan and Juneja [20], Meneveau [21], Grossmann, Lohse and Rech [22] and Sain, Manu and Pandit [23]. A key feature of this work was the implication that corrections to the exponents of structure functions increase with increasing Reynolds number, which suggests that intermittency is the dominant effect. In the next section, we will explain ESS in more detail, in terms of how it relates to other ways of obtaining exponents.

## II. STRUCTURE FUNCTIONS, EXPONENTS AND ESS

In this section we discuss the various ways in which exponents are defined and measured. The longitudinal structure functions are defined as

$$S _ { n } ( r ) = \langle \delta u _ { L } ^ { n } ( r ) \rangle$$

where the (longitudinal) velocity increment is given by

$$\delta u _ { L } ( r ) = \left [ u ( x + r, t ) - u ( x, t ) \right ] \cdot \hat { r } \.$$

Integration of the K´rm´n-Howarth equation (KHE) leads, a a in the limit of infinite Reynolds number , to the Kolmogorov '4/5' law, S 3 ( r ) = -(4 / 5) εr . If the S n , for n ≥ 4, exhibit a range of power-law behavior, then; in general, and solely on dimensional grounds, the structure functions of order n are expected to take the form

$$S _ { n } ( r ) = C _ { n } ( \varepsilon r ) ^ { n / 3 } \.$$

Measurement of the structure functions has repeatedly found a deviation from the above dimensional prediction for the exponents. If the structure functions are taken to scale with exponents ζ n , thus:

$$S _ { n } ( r ) & \sim r ^ { \zeta _ { n } }$$

then it has been found [16, 18] that the difference ∆ n = | n/ 3 -ζ n | is non-zero and increases with order n . Exponents ζ n which differ from n/ 3 are often referred to as anomalous exponents [18].

FIG. 1. Comparison of our values of local-slope exponents ζ n ( r ) (symbols) given in (5) with our values of the ESS exponents Σ n ( r ) (lines) as calculated from equation (8), divided by n/ 3 to show their relationship to K41 exponents. Results are given for n = 2 4 6. Both sets of exponents were , , calculated from the real-space velocity field and are presented here for R λ = 177. The separation, r , has been scaled on the dissipation scale, η = ( ν 3 0 /ε ) 1 4 / .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-29e896e1-42d7-4ac2-82e5-5350a767fcb1.jpg)

In order to study the behavior of the exponents ζ n , it is usual to make a log-log plot of S n against r , and measure the local slope:

$$\zeta _ { n } ( r ) = \frac { d \log S _ { n } ( r ) } { d \log r } \.$$

Following Fukayama et al. [19], the presence of a plateau when any ζ n ( r ) is plotted against r indicates a constant exponent, and hence a scaling region. Yet, it is not until comparatively high Reynolds numbers are attained that such a plateau is found. Instead, as seen in Fig. 1 (symbols), even for the relatively large value of Reynolds number, R λ = 177, a scaling region cannot be identified. (We note that Grossmann et al. [22] have argued that a minimum value of R λ ∼ 500 is needed for satisfactory direct measurement of local scaling exponents.)

The introduction of ESS relied on the fact that S 3 scales with ζ 3 = 1 in the inertial range.

Benzi et al. [17] argued that if

$$S _ { n } ( r ) \sim [ S _ { 3 } ( r ) ] ^ { \zeta _ { n } ^ { * } } \, \quad \text{with} \quad \zeta _ { n } ^ { * } = \zeta _ { n } / \zeta _ { 3 } \.$$

ζ ∗ n should then be equivalent to ζ n in the scaling region.

A practical difficulty led to a further step. The statistical convergence of odd-order structure functions is significantly slower than that for even-orders, due to the delicate balance of positive and negative values involved in the former [19]. To overcome this, generalized structure functions, where the velocity difference is replaced by its modulus, have been introduced [17] (see also [19, 20]),

$$G _ { n } ( r ) = \langle | \delta u _ { L } ( r ) | ^ { n } \rangle \sim r ^ { \zeta _ { n } ^ { \prime } },$$

with scaling exponents ζ ′ n . The fact that S 3 ∼ r in the inertial range does not rigorously imply that G 3 ∼ r in the same range. But, by plotting G 3 ( r ) against | S 3 ( r ) , Benzi | et al. [18] showed that, for R λ = 225 - 800, the third-order exponents satisfied ζ ′ 3 glyph[similarequal] 1 006 . ζ 3 . Hence, it is now generally assumed that ζ ′ n and ζ n are equal (although, Fig. 2 in Belin, Tabeling and Willaime [24] implies some discrepancy at the largest length scales, and the authors note that the exponents ζ ′ n and ζ n need not be the same). Thus, by extension, G 3 with ζ ′ 3 = 1, leads to

$$G _ { n } ( r ) \sim [ G _ { 3 } ( r ) ] ^ { \Sigma _ { n } }, \quad \text{with} \quad \Sigma _ { n } = \zeta _ { n } ^ { \prime } / \zeta _ { 3 } ^ { \prime } \.$$

Benzi et al. [17] found that plotting their results on this basis gave a larger scaling region. This extended well into the dissipative lengthscales and allowed exponents to be more easily extracted from the data. Also, Grossmann et al. [22] state that the use of generalized structure functions is essential to take full advantage of ESS.

There is, however, an alternative to the use of generalized structure functions. This is the pseudospectral method . In using this for some of the present work, we followed the example of Qian [25, 26], Tchoufag, Sagaut and Cambon [27] and Bos, Chevillard, Scott and Rubinstein [28], who obtained S 2 and S 3 from the energy and energy transfer spectra, respectively, by means of exact quadratures.

The organization of our own work in this paper is now as follows. We begin with a description of our DNS before illustrating ESS, using results from our own simulations, in section IV, where we show that our results for ESS agree closely with those of other investigations [18, 19]. These particular results were obtained in the usual way by direct

convolution sums, using a statistical ensemble, and the generalized structure functions. In section V we describe the theoretical basis for using the pseudospectral method [25-28] which includes a rigorous derivation of the forcing in the real-space energy balance equation. This is followed by the introduction of a new scaling exponent in section VI and a presentation of our numerical results for seven Taylor-scale Reynolds numbers spanning the range 101 3 . ≤ R λ ≤ 435 2. .

## III. NUMERICAL METHOD

We used a pseudospectral DNS, with full dealiasing implemented by truncation of the velocity field according to the two-thirds rule [29]. Time advancement for the viscous term was performed exactly using an integrating factor, while the non-linear term was stepped forward in time using Heun's method [30], which is a second-order predictor-corrector routine. Each simulation was started from a Gaussian-distributed random field with a specified energy spectrum, which followed k 4 for the lowk modes. Measurements were taken after the simulations had reached a stationary state. The system was forced by negative damping, with the Fourier transform of the force f given by

$$f ( k, t ) & = ( \varepsilon _ { W } / 2 E _ { f } ) u ( k, t ) \quad \text{for} \quad 0 < | k | < k _ { f } ; \\ & = 0 \quad \text{otherwise},$$

where u k ( , t ) is the instantaneous velocity field (in wavenumber space). The highest forced wavenumber, k f , was chosen to be k f = 2 5 . k min , where k min = 2 π/L box = 1 is the lowest resolved wavenumber. As E f was the total energy contained in the forcing band, this ensured that the energy injection rate was ε W = constant. It is worth noting that any method of energy injection employed in the numerical simulation of isotropic turbulence is not experimentally realizable. The present method of negative damping has also been used in other investigations [31-34], albeit not necessarily such that ε W is maintained constant (although note the theoretical analysis of this type of forcing by Doering and Petrov [35]). Also, note that the correlation between the force and the velocity is restricted to the very lowest wavenumbers.

For each Reynolds number studied, we used the same initial spectrum and input rate ε W . The only initial condition changed was the value assigned to the (kinematic) viscosity.

Once the initial transient had passed, the velocity field was sampled every half a large-eddy turnover time, τ = L/U , where L denotes the average integral scale and U the rms velocity. The ensemble populated with these sampled realizations was used, in conjunction with the usual shell averaging, to calculate statistics. Simulations were run using lattices of size 128 3 , 256 3 , 512 3 , 1024 3 and 2048 , with corresponding Reynolds numbers ranging from 3 R λ = 41 8 up to 435 2. The smallest wavenumber was . . k min = 2 π/L box = 1 in all simulations, while the maximum wavenumber satisfied k max η glyph[greaterorequalslant] 1 30 for all runs except one which satisfied . k max η glyph[greaterorequalslant] 1 01, where . η is the Kolmogorov dissipation lengthscale. The integral scale, L , was found to lie between 0 23 . L box and 0 17 . L box . It can be seen in Figure 2 of McComb, Hunter and Johnston [36] that a small-scale resolution of k max η &gt; 1 6 is desirable in order to cap-. ture the relevant dissipative physics. Evidently, this would restrict the attainable Reynolds number of the simulated flow, and the reference suggests that k max η glyph[greaterorequalslant] 1 3 would still be . acceptable (containing ∼ 99 5% of dissipative dynamics [37]). . In contrast, at k max η glyph[similarequal] 1 a non-negligible part of dissipation is not taken into account. Most high resolution DNSs of isotropic turbulence try to attain Reynolds numbers as high possible and thus opt for minimal resolution requirements. In this paper the simulations have been conducted following a more conservative approach, where the emphasis has been put on higher resolution, thus necessarily compromising to some extent on Reynolds number. Large-scale resolution has only relatively recently received attention in the literature. As mentioned above, the largest scales of the flow are smaller than a quarter of the simulation box size. Details of the individual runs are summarized in Table I.

Our simulations have been well validated by means of extensive and detailed comparison with the results of other investigations. Further details of the performance of our code including verification of isotropy may be found in the thesis by Yoffe [37], along with values for the Kolmogorov constant and velocity-derivative skewness; and a direct comparison with the freely-available pseudospectral code hit3d [38, 39]. Furthermore our data reproduces the characteristic behavior for the plot of the dimensionless dissipation rate C ε against R λ [40], and agree closely with other representative results in the literature, such as the work by Wang, Chen, Brasseur and Wyngaard [41], Cao, Chen and Doolen [42], Gotoh, Fukayama and Nakano [43], Kaneda, Ishihara, Yokokawa, Itakura and Uno [33], Donzis, Sreenivasan and Yeung [44] and Yeung, Donzis and Sreenivasan [45], although there are some differences in forcing methods.

|   R λ |     ν 0 |    N |     ε |     U |   L/L box |   k max η | M   |
|-------|---------|------|-------|-------|-----------|-----------|-----|
|  42.5 | 0.01    |  128 | 0.094 | 0.581 |      0.23 |      2.34 | 101 |
|  64.2 | 0.005   |  128 | 0.099 | 0.607 |      0.21 |      1.37 | 101 |
| 101.3 | 0.002   |  256 | 0.099 | 0.607 |      0.19 |      1.41 | 101 |
| 113.3 | 0.0018  |  256 | 0.1   | 0.626 |      0.2  |      1.31 |     |
| 176.9 | 0.00072 |  512 | 0.102 | 0.626 |      0.19 |      1.31 | 15  |
| 203.7 | 0.0005  |  512 | 0.099 | 0.608 |      0.18 |      1.01 |     |
| 217   | 0.0005  | 1024 | 0.1   | 0.63  |      0.19 |      2.02 |     |
| 276.2 | 0.0003  | 1024 | 0.1   | 0.626 |      0.18 |      1.38 |     |
| 335.2 | 0.0002  | 1024 | 0.102 | 0.626 |      0.18 |      1.01 |     |
| 435.2 | 0.00011 | 2048 | 0.102 | 0.614 |      0.17 |      1.3  |     |

TABLE I. A summary of the main parameters of our numerical simulations. The values quoted for the dissipation rate ε , the rms velocity U and the integral scale L are ensemble- and shell-averaged mean values, where the ensembles have been populated with snapshots of the steady-state velocity field taken every half large eddy turnover time. M counts individual realizations used to calulate ESS exponents for those runs for which the ESS method has been performed, not the size of the ensemble used for calculating statistics.

## IV. REAL-SPACE CALCULATION OF THE STRUCTURE FUNCTIONS AND ESS

In order to calculate the structure functions in real space, we calculate the longitudinal correlation of one lattice site with all other sites

$$S _ { n } ( r ) = \frac { 1 } { 3 N ^ { 3 } } \sum _ { x } \left [ \left ( u _ { x } ( x + r e _ { x } ) - u _ { x } ( x ) \right ) ^ { n } + \left ( u _ { y } ( x + r e _ { y } ) - u _ { y } ( x ) \right ) ^ { n } + \left ( u _ { z } ( x + r e _ { z } ) - u _ { z } ( x ) \right ) ^ { n } \right ],$$

for each realization. The results are subsequently ensemble-averaged over many realizations. Figure 2 shows the calculated standard local slopes for four different Reynolds numbers R λ = 42, 64, 101 and 177. A plateau would indicate a constant exponent, that is a scaling region, but the figure does not show the formation of plateaux for these Reynolds numbers, implying that there is no scaling region. If, in contrast the generalized structure functions

FIG. 2. (Color online) The local slopes ζ n ( r ) of the structure functions S n ( r ) for n = 2 , · · · , 8 from our DNS. The slopes are calculated using the structure functions evaluated in real space for Reynolds numbers R λ = 42, 64, 101 and 177. Note the absence of a scaling region in all cases.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-b0fb5c5f-bf00-4364-b380-6359ea275a82.jpg)

are used to evaluate the local slopes Σ n as defined in (8), scaling regions are obtained and thus the ESS scaling exponents can be measured. As shown in Fig. 3, extended plateaux for each order of generalized structure function are observed, which reach well into the dissipation scales where we do not expect power-law behavior. Since S n ∼ r n as r → 0, we note that Σ n ( r ) → n/ 3 as r → 0. It should be borne in mind that this K41-type behavior is an artefact of ESS, as has been pointed out before [18, 46, 47]. ESS exponents obtained from our calculations are shown to be consistent with relevant results from the literature in Table II.

## V. SPECTRAL METHODS

The use of spectral methods to calculate the structure functions is known in the literature through the work by Bos et al. [28], Tchoufag et al. [27] and Qian [25, 26]. We will now

FIG. 3. (Color online) The same results as in Fig. 2 evaluated using the ESS method for structure functions calculated in real space. The vertical lines indicate the scaling region used to evaluate the scaling exponent, Σ n .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-98d4acc0-af93-417b-abce-c65b1dbe0d42.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-02ae7b61-d444-4196-9c43-a44c780dffe6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-4a2246e0-5b3a-4c71-b940-c6a8a2b69f81.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-2aa6f887-16a6-489e-a607-237596d4b3e3.jpg)

briefly explain the reason for doing it this way.

Since the calculation of correlation and structure functions in real space requires a convolution in which the correlation of each site with every other (longitudinal) site needs to be measured, moving to large lattice sizes in order to reach larger Reynolds numbers results in a significant increase in computational workload. In addition, the number of realizations required to generate the ensemble takes both longer to produce and occupies substantially more storage space. In order to reach higher Reynolds numbers we have therefore used spectral expressions for the correlation and structure functions. Real-space quantities are then calculated by Fourier-transforming the appropriate spectral density. For example, the two-point correlation tensor may be found using

$$C _ { \alpha \beta } ( r ) = \, \int d ^ { 3 } k \, \left \langle u _ { \alpha } ( k ) u _ { \beta } ^ { * } ( k ) \right \rangle \ e ^ { i k \cdot x } \.$$

The assumption of isotropy then allows the transformation for the calculation of the isotropic

TABLE II. Measurement of the scaling exponents from our DNS data using ESS. For purposes of comparison, we also show the values predicted from dimensional analysis (referred to as 'K41 theory'), along with results from Fukayama et al. [19], Gotoh et al. [43] and Benzi et al. [18].

| R λ   |   Σ 2 |   Σ 4 |   Σ 5 |   Σ 6 |   Σ 7 |   Σ 8 | Source         |
|-------|-------|-------|-------|-------|-------|-------|----------------|
|       | 0.667 | 1.333 | 1.667 | 2     | 2.333 | 2.667 | K41 theory     |
| 42.5  | 0.69  | 1.287 | 1.551 | 1.796 | 2.023 | 2.234 | Our DNS, Σ n   |
| 64.2  | 0.692 | 1.284 | 1.544 | 1.783 | 2.002 | 2.201 | Our DNS, Σ n   |
| 101.3 | 0.692 | 1.283 | 1.544 | 1.785 | 2.008 | 2.215 | Our DNS, Σ n   |
| 176.9 | 0.694 | 1.279 | 1.533 | 1.762 | 1.967 | 2.151 | Our DNS, Σ n   |
| 70    | 0.69  | 1.288 | 1.555 | 1.804 | 2.037 | 2.254 | DNS [19], Σ n  |
| 125   | 0.692 | 1.284 | 1.546 | 1.788 | 2.011 | 2.217 | DNS [19], Σ n  |
| 381   | 0.709 | 1.3   | 1.56  | 1.79  | 1.99  | 2.18  | DNS [43], ζ n  |
| 460   | 0.701 | 1.29  | 1.54  | 1.77  | 1.98  | 2.17  | DNS [43], ζ n  |
| 800   | 0.7   | 1.28  | 1.54  | 1.78  | 2     | 2.23  | Exp. [18], Σ n |

correlation function to be reduced to

$$C ( r ) = \frac { 1 } { 2 } C _ { \alpha \alpha } ( r ) = \int d k \ E ( k ) \frac { \sin k r } { k r } \.$$

The real space correlation tensor C αβ ( r ), when written as a spatial average instead of an ensemble average (assuming ergodicity, as is the usual practice), is a convolution

$$C _ { \alpha \beta } ( r ) = \lim _ { V \to \infty } \frac { 1 } { V } \int _ { V } d x \ u _ { \alpha } ( x ) u _ { \beta } ( x + r )$$

of the velocity components u α and u β . In a similar manner to the standard pseudospectral DNS technique, that is switching to real space in order to calculate the non-linear term and to avoid the convolution in Fourier space [37], this approach replaces the convolution in real space with a Fourier transform of local Fourier-space spectra. Furthermore, the shellaveraged spectra require substantially less in the way of storage and processing capabilities than real-space ensembles.

The derivations of the spectral representation of the second- and third-order structure

functions are given in Appendix A, leading to

$$S _ { 2 } ( r ) = 4 \int _ { 0 } ^ { \infty } d k \ E ( k ) a ( k r )$$

and where

$$a ( x ) = \frac { 1 } { 3 } - \frac { \sin x - x \cos x } { x ^ { 3 } }$$

in agreement with Bos et al. [28]. The local slopes ζ n can now be found by taking derivatives of the spectral forms for the structure functions , as shown in further detail in Appendix A. 1

The spectral approach has the consequence that we are now evaluating the conventional structure functions S n ( r ), rather than the generalized structure functions, G n ( r ), as commonly used (including by us) for ESS. Before proceeding to the calculation of the scaling exponents, we will now discuss the effects of finite forcing and the calculation of viscous corrections to the second- and third-order structure functions, and subsequently validate the pseudospectral approach by comparing results for the spectrally obtained second- and third-order structure functions to corresponding real-space results.

## A. Effects of finite forcing on the structure functions

The pseudospectral method can also be used to calculate the corrections due to finite forcing on the structure functions. The second- and third-order structure functions are related by energy conservation, that is in real space by the K´rm´n-Howarth equation a a

$$- \, \frac { 3 } { 2 } \frac { \partial U ^ { 2 } } { \partial t } = - \frac { 3 } { 4 } \frac { \partial S _ { 2 } } { \partial t } - \frac { 1 } { 4 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } S _ { 3 } \right ) + \frac { 3 \nu _ { 0 } } { 2 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } \frac { \partial S _ { 2 } } { \partial r } \right )$$

where (3 / 2) ∂U /∂t 2 = ∂E/∂t = -ε D , with ε D denoting the decay rate. For decaying turbulence the decay rate equals the dissipation rate ε , that is ε D = . For stationary turbulence ε ε D = 0, and the dissipation rate then has to be equal to the energy input rate, that is ε = ε W , where ε W denotes the energy input rate, and the time-dependent S 2 term vanishes. There is, however, a further complication. Replacing the dissipation rate ε with the energy input rate ε W , we see that the dissipation rate is acting as an input. In the K41 theory, we

1 Currently this is only possible for S 2 ( r ) and S 3 ( r ), as the spectral expressions for the higher order structure functions have not been derived yet.

$$S _ { 3 } ( r ) = 6 C _ { L L, L } ( r ) = 1 2 \int _ { 0 } ^ { \infty } d k \ \frac { T ( k ) } { k ^ { 2 } } \ \frac { \partial a ( k r ) } { \partial r }$$

have an equivalence between the inertial transfer, dissipation, and input rates, in the infinite Reynolds number limit. However, in using ε (= ε W ) as the input, it has been implicitly assumed that the forcing does not depend on the scale r . This is in general not the case, and as such some method for accounting for the effects of (finite) forcing ( i.e. Re &lt; ∞ ) must be introduced. Various proposals have been put forward, mainly through the inclusion of corrections to the KHE. These are discussed in Appendix B. An alternative treatment containing the exact energy input is given in the following section.

## 1. The KHE derived from the Fourier space energy balance

Instead of including a correction term in the KHE for forced turbulence, the effect of finite forcing can be calculated exactly using the spectral approach. We begin with the energy balance equation in spectral space (nowadays referred to as the Lin equation [40, 48])

$$\frac { \partial E ( k, t ) } { \partial t } = T ( k, t ) - 2 \nu _ { 0 } k ^ { 2 } E ( k, t ) + W ( k, t )$$

where W k,t ( ) is the work spectrum of the stirring forces and thus contains the relevant information about the forcing. In order to obtain the energy balance equation in real space including the effects of (finite) forcing, we assume isotropy, take the Fourier transform of the Lin equation and use the definitions of the structure functions S 2 and S 3 in order to obtain the energy balance equation in real space (the KHE) relating S 2 and S 3

$$\frac { \partial U ^ { 2 } } { \partial t } - \frac { 1 } { 2 } \frac { \partial S _ { 2 } ( r ) } { \partial t } = \frac { 1 } { 6 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } S _ { 3 } ( r ) \right ) - \frac { \nu _ { 0 } } { r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } \frac { \partial S _ { 2 } ( r ) } { \partial r } \right ) + \frac { 2 } { 3 } I ( r, t )$$

where the input term I ( r, t ) is defined as

$$I ( r, t ) = \frac { 3 } { r ^ { 3 } } \int _ { 0 } ^ { r } d y \ y ^ { 2 } W ( y, t )$$

and W y,t ( ) is the (three-dimensional) Fourier transform of the work spectrum W k,t ( ). A detailed derivation of this equation from the Lin equation (18) can be found in Appendix C.

In order to write (19) in terms of energy loss, we multiply by -3 / 2 on both sides and obtain

$$\varepsilon _ { D } = - \frac { 3 } { 4 } \frac { \partial S _ { 2 } ( r ) } { \partial t } - \frac { 1 } { 4 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } S _ { 3 } ( r ) \right ) + \frac { 3 \nu _ { 0 } } { 2 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } \frac { \partial S _ { 2 } ( r ) } { \partial r } \right ) - I ( r, t ) \.$$

This is now the general form of the KHE including the effects of general, unspecified forcing. For the case of free decay, ε D = ε and the input term I ( r, t ) = 0, which leads to the wellknown KHE for free decay. On the other hand, for the stationary case ε D = 0, ∂S /∂t 2 = 0

and the input term is independent of time I ( r, t ) = I ( r ), thus

$$I ( r ) = - \frac { 1 } { 4 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } S _ { 3 } ( r ) \right ) + \frac { 3 \nu _ { 0 } } { 2 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } \frac { \partial S _ { 2 } ( r ) } { \partial r } \right ) \.$$

This derivation shows that the dissipation rate should not be present in the KHE for stationary turbulence with finite forcing. Only for the limit of δ ( k )-forcing do we obtain I ( r ) = ε W = ε and hence recover the dissipation rate as an input term.

In contrast to previous attempts to include the effects of forcing into the KHE, we do not approximate the work term. Instead, we use full information of this term as supplied by the work spectrum. In this way an explicit form for the actual energy input by the stirring forces can be calculated, as we shall now show. By integrating (22) with respect to r one obtains

$$S _ { 3 } ( r ) = X ( r ) + 6 \nu _ { 0 } \frac { \partial S _ { 2 } } { \partial r }$$

with the input due to finite forcing

$$X ( r ) = - \frac { 4 } { r ^ { 4 } } \int _ { 0 } ^ { r } d y \ y ^ { 4 } I ( y )$$

evaluated using the spectral method

$$X ( r ) = - 1 2 r \int _ { 0 } ^ { \infty } d k \ W ( k ) \ \left [ \frac { 3 \sin k r - 3 k r \cos k r - ( k r ) ^ { 2 } \sin k r } { ( k r ) ^ { 5 } } \right ] \. \quad \ \ ( 2 5 )$$

Note that X r ( ) is not a correction to K41, as used in previous studies. Instead, it replaces the erroneous use of the dissipation rate and contains all the information of the energy input at all scales. In the limit of δ ( k )-forcing, I ( y ) = ε W = ε , such that X r ( ) = -4 εr/ 5, giving K41 in the infinite Reynolds number limit.

## B. Real-space versus spectral structure functions

Results from pseudospectral calculations of the structure functions were compared to real space ensemble averaged results for R λ = 101 and R λ = 177. As can be seen in Fig. 4, agreement for S 2 ( r ) is very good for all r . For S 3 ( r ) we observe good agreement for small r , but the curves diverge at large r . This could be due to DNS data being periodic in L box = 2 π . Since S 3 ( r ) is an odd function of r , it must go to zero in the center of the domain. The pseudospectral method, however, involves a (weighted) superposition of damped oscillating functions which does not necessarily require that S 3 ( L box / 2) = 0. Figure

5 shows the compensated second- and third-order structure function calculated from realspace ensembles for the Taylor-scale Reynolds number range 43 glyph[lessorequalslant] R λ glyph[lessorequalslant] 177. This may be compared to figure 3(b) in the paper by Ishihara, Gotoh and Kaneda [49].

## C. Spectral calculation from DNS

The pseudospectral method was used to calculate structure functions for R λ = 101, 113 and also for the higher Reynolds number range 177 glyph[lessorequalslant] R λ glyph[lessorequalslant] 435, which can be seen in Fig. 6, where the arrows indicate the direction of increasing Reynolds number. The lower horizontal dotted line in the picture indicates K41 scaling for the third-order structure function, while the upper horizontal line indicates our measured value C K,S 2 = 2 07 for the . prefactor C K,S 2 of the second-order structure function. The prefactor C K,S 2 is related to the prefactor C K,E LL of the longitudinal energy spectrum E LL ( k ) by C K,S 2 = 4 02 . C K,E LL [50], which has been measured by Sreenivasan [51] to be C K,E LL = 0 53 . ± 0 055. . This results in C K,S 2 = 2 13 . ± 0 22, compared to our measured value of . C K,S 2 = 2 07. Comparison of Figs. 5 . and 6 shows how the spectrally calculated results for the second- and third-order structure functions extend the real-space results to higher Reynolds numbers.

The viscous correction to the 'four-fifths'-law and to the exact input contribution X r ( ) calculated by the spectral method are shown in Fig. 7 and Fig. 8, respectively, with the thirdorder structure function plotted for comparison. This is presented for our highest resolved simulation at R λ = 435 2. . Together with the viscous correction, the input X r ( ) given in (25) can be seen to account for differences between the third-order structure function and the 'four-fifths'-law at all scales, as can be seen in Fig. 8. In contrast, Fig. 7 shows that the viscous correction alone only accounts for the difference between DNS data for the third-order structure function and K41 at the small scales, as expected, since at scales much smaller than the forcing scale the system becomes insensitive to the details of the (large-scale) forcing.

## VI. A NEW SCALING EXPONENT

We now arrive at our proposal to introduce a new local-scaling exponent Γ n , which can be used to determine the ζ n . We work with S n ( r ) and consider the quantity | S n ( r /S ) 3 ( r ) . |

FIG. 4. (Color online) Comparison of second- and third-order (scaled) structure functions calculated from the real-space ensembles and from energy and transfer spectra for two different Reynolds numbers. The agreement for S 2 is good at all scales, while for S 3 the results from real-space and spectral space calculations diverge at the large scales.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-ae555ea3-1f27-4e47-80a0-d0db216ad35a.jpg)

FIG. 5. Second- and third-order (scaled) structure functions, calculated from the real-space ensembles. This figure should be compared to Fig. 6. The comparison shows how the real-space results can be extended to higher Reynolds numbers using the spectral method.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-b4307228-56b4-4b57-927a-27f65c357d1f.jpg)

In this procedure, the exponent Γ n is defined by

$$\left | \frac { S _ { n } ( r ) } { S _ { 3 } ( r ) } \right | & \sim r ^ { \Gamma _ { n } } \, \quad \text{ where } \quad \Gamma _ { n } = \zeta _ { n } - \zeta _ { 3 } \.$$

The definition of Γ n is motivated by a long-established technique in experimental physics, where the effective experimental error can be reduced by plotting the ratio of two dependent variables: see e.g. Chapter 3 in the well-known book by Bevington and Robinson on data analysis [52]. Of course this does not work in all cases, but only where the quantities are positively correlated, and we have verified that this is the case for S 2 and | S 3 | .

The error in the measurement of the n th -order structure function can be expressed as

$$S _ { n } ( r ) = \left [ 1 + \epsilon _ { n } ( r ) + O ( \epsilon _ { n } ^ { 2 } ) \right ] \overline { S } _ { n } ( r )$$

where S n ( r ) is the 'true' value and glyph[epsilon1] n ( r ) a measurement of the systematic error and considered small. Hence if glyph[epsilon1] n ( r ) ∼ glyph[epsilon1] ( r ), then the ratio ∣ ∣ ∣ S n ( r ) S 3 ( r ) ∣ ∣ ∣ has an error proportional to the

10

FIG. 6. Compensated second- and third-order (scaled) structure functions, calculated from energy and transfer spectra. The horizontal lines indicate K41 scaling, where the upper horizontal line indicates our measured value C K,S 2 = 2 07 for the prefactor . C K,S 2 of the second-order structure function. The arrows indicate the direction of increasing Reynolds number, and one observes the formation of longer plateaux for both S 2 and S 3 with increasing Reynolds number indicating an approach to K41 scaling.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-4f85b60d-c01c-4d0f-8823-922622fbaa19.jpg)

second-order of a small quantity,

$$\frac { S _ { n } ( r ) } { S _ { 3 } ( r ) } = \frac { \overline { S } _ { n } ( r ) } { \overline { S } _ { 3 } ( r ) } \left [ 1 - \epsilon ^ { 2 } ( r ) + O ( \epsilon ^ { 3 } ) \right ] \,.$$

For illustration purposes we assumed here that S n and S 3 were perfectly correlated, note that for imperfect correlation there is still a reduction in error. The local slope now is found by considering Γ ( 2 r ) = ζ 2 ( r ) -ζ 3 ( r ). By once again assuming that ζ 3 ( r ) = 1, the local slope for the second-order structure function is found as Γ ( 2 r ) + 1 = ζ 2 ( r ).

## A. Scaling exponents from spectra

The new scaling exponent Γ n , based on the conventional (as opposed to the generalized ) structure functions S n ( r ), is compatible with spectral methods and has been tested for the

FIG. 7. Evaluation of the viscous correction to the 'four-fifths'-law for the third-order structure function using the spectral method, with R λ = 435 2. . The correction can be seen to account for the difference between the 'four-fifths'-law and DNS data at small scales.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-964ee23b-02f4-41b8-b4c1-4de320aa861d.jpg)

case n = 2 in Fig. 9. The dimensionless quantity U S | 2 ( r /S ) 3 ( r ) , where | U is the rms velocity, is plotted against r/η , for three values of R λ . Note that, since K41 predicts Γ 2 = -1 / 3, we have plotted a compensated form, in which we multiply the ratio by ( r/η ) 1 3 / , such that K41 scaling would correspond to a plateau. From the figure, we can see a trend towards K41 scaling as the Reynolds number is increased. Note that this figure also illustrates the ranges used to find values for our new exponent Γ , for the following cases. 2 Γ 2 was fitted to the ranges λ &lt; r &lt; cλ , with c = 2 0 2 5 2 6, and 2.7 for . , . , . R λ = 101 217 276 and 435, , , respectively.

Figure 10 summarizes the comparison between our results for our new method of determining the second-order exponent and those based on ESS (our own and others [18, 19]) or on direct measurement [43], in terms of their overall dependence on the Taylor-Reynolds number. In order to establish the form of the dependence of the exponents on Reynolds number, we fitted curves to the data points using the nonlinear-least-squares MarquardtLevenberg algorithm, with the error quoted being one standard error. Using the data ob-

FIG. 8. Evaluation of the viscous correction to the input X r ( ) (in place of the 'four-fifths'-law) for the third-order structure function. Calculations were performed using the spectral method, with R λ = 435 2. The contributions can be seen to account for the differences between the 'four-fifths'-. law and DNS data at all scales.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-6b07f56e-916f-4fe9-966a-fc5b4714416a.jpg)

tained from our new method, we fitted a curve Γ 2 +1 = A + BR p λ , to find the asymptotic value A = 0 679 . ± 0 013, . which is consistent with the deviations from K41 scaling being finite Reynolds number effects. In order to compare this result with results obtained by the ESS method, we fitted the curve Σ 2 = C + DR q λ to our own data plus that of Fukayama et al. [19]. Evidently the two fitted curves show very different trends, with results for Σ 2 increasing with increasing Reynolds number, whereas ζ 2 = Γ 2 +1 decreases and approaches 2/3 (within one standard error) as R λ increases.

It should be emphasized that with both methods, that is ESS and our new method, it is necessary to take ζ 3 = 1 in the inertial range, in order to obtain the inertial-range value of either Σ 2 = ζ 2 (by ESS) or Γ 2 = ζ 2 -1 (our new method). For this reason, we plot Γ 2 +1, rather than Γ 2 in Fig. 10. An obvious difference between our proposed method and ESS is apparent as r → 0. This is readily understood in terms of the regularity condition for the velocity field, which leads to S n ( r ) ∼ r n as r → 0 [20, 53]. This yields Γ n ( r ) + 1 → n -2,

FIG. 9. The compensated ratio ( r/η ) 1 3 / U S | 2 ( r /S ) 3 ( r ) | plotted against r and scaled on the dissipation scale, η . K41 scaling would correspond to a plateau. Arrows indicate λ/η , while the vertical dotted lines show the region used to fit each exponent. Note that the measured compensated slopes become flatter with increasing Reynolds number, pointing towards the deviation from Kolmogorov scaling being a finite Reynolds number effect.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-0f19ad4c-e775-4605-91a4-640fa8af15a6.jpg)

whereas ESS gives Σ n ( r ) → n/ 3.

In this context it may be of interest to briefly discuss the experimental results of Mydlarski and Warhaft [14], who measured the exponent of the longitudinal energy spectrum for a range of Taylor-Reynolds numbers from R λ = 50 to R λ = 473, which is similar to the range of Taylor-Reynolds numbers studied in the present paper. The authors found the inertial-range exponent α of the longitudinal energy spectrum E LL ( k ) to depend on R λ in the following way

$$\alpha = - \frac { 5 } { 3 } \left ( 1 + 3. 1 5 R _ { \lambda } ^ { - 2 / 3 } \right )$$

hence α → -5 3 in the limit of infinite Reynolds number. / Thus the results of Mydlarski and Warhaft support our result for the exponent of S 2 , since α → -5 3 implies / ζ 2 → 2 3 / as R λ →∞ .

FIG. 10. Dependence of our new exponent ζ 2 = Γ +1 on Reynolds number, compared to measured 2 ESS exponents denoted by Σ 2 from Fukayama et al. [19], Benzi et al. [18] and our DNS, as well as results for ζ 2 by Gotoh et al. [43]. The horizontal line indicates the K41 value of 2/3. The dash-dot line indicates the fit to Γ 2 +1, while the dashed line shows a fit to Σ 2 using our results and those of Fukayama et al. [19].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0539v2-figure-7ac831bf-a70a-478d-9214-8407fe148d06.jpg)

## VII. CONCLUSIONS

As we have said in the introduction, the point at issue is essentially 'intermittency corrections versus finite Reynolds number effects'. The former has received much more attention; but, in recent years, there has been a growing interest in studying finite Reynolds number effects, experimentally and by DNS, for the case of S 3 : see [27, 43, 54] and references therein. (Although we note that in this case the emphasis is on the prefactor rather than the exponent.)

Our new result that ζ 2 = Γ 2 +1 → 2 3 is an indication that anomalous values of / ζ 2 are due to finite Reynolds number effects, consistent with the experimental results of Mydlarski et al. which point in the same direction. Previously it had been suggested by Barenblatt et al. that ESS could be interpreted in this way [46], but this was disputed by Benzi et al. [47].

There is much remaining to be understood about these matters and we suggest that our new method of analyzing data can help. It should, of course, be noted that our use of S 3 (as evaluated by pseudospectral methods) rather than G 3 (as used with ESS), may also be a factor in our result. As a matter of interest, we conclude by noting that our analysis could provide a stimulus for further study of ESS and may lead to an understanding of the relationship between the two methods. It is also the case that the pseudospectral method could be used for the general study of higher-order structure functions, but this awaits the derivation of the requisite Fourier transforms.

## ACKNOWLEDGMENTS

The authors would like to thank Matthew Salewski, who read a first draft and made many helpful suggestions. This work has made use of the resources provided by HECToR [55], made available through the Edinburgh Compute and Data Facility (ECDF) [56]. A. B. is supported by STFC, S. R. Y. and M. F. L. are funded by EPSRC.

## Appendix A: Derivation of the spectral representation of the structure functions

The second order structure function can be derived from the relationship between the isotropic correlation function C r ( ) and the energy spectrum E k ( ) where

$$C ( r ) = \int _ { 0 } ^ { \infty } d k \ E ( k ) \ \frac { \sin k r } { k r }$$

and

$$E ( k ) = \frac { 2 } { \pi } \int _ { 0 } ^ { \infty } d r \ C ( r ) \ k r \sin k r$$

as e.g. in Batchelor [57] equation (3.4.15) p49. Equation (A1) can be used to derive a spectral expression for the longitudinal correlation function C LL ( r ), since C r ( ) and C LL ( r ) are related through

$$C ( r ) = \frac { 1 } { 2 r ^ { 2 } } \frac { \partial } { \partial r } \left [ r ^ { 3 } C _ { L L } ( r ) \right ] \,,$$

which can be integrated to give

$$C _ { L L } ( r ) = \frac { 2 } { r ^ { 3 } } \int _ { 0 } ^ { r } d y \ y ^ { 2 } C ( y ) \.$$

Now the integral over y is done analytically to obtain

$$C _ { L L } ( r ) = 2 \int _ { 0 } ^ { \infty } d k \ E ( k ) \ \left [ \frac { \sin k r - k r \cos k r } { ( k r ) ^ { 3 } } \right ]$$

as in Monin and Yaglom [50] vol. 2, equation (12.75). The spectral expression for the second-order structure function is now readily seen to be

$$S _ { 2 } ( r ) = 2 U ^ { 2 } - 2 C _ { L L } ( r ) = 2 \cdot \underbrace { \frac { 2 } { 3 } \int _ { 0 } ^ { \infty } d k \ E ( k ) } _ { U ^ { 2 } } - 4 \int _ { 0 } ^ { \infty } d k \ E ( k ) \ \left [ \frac { \sin k r - k r \cos k r } { ( k r ) ^ { 3 } } \right ] \, \ ( A 6 )$$

which can be written in a more concise form

$$S _ { 2 } ( r ) = 4 \int _ { 0 } ^ { \infty } d k \ E ( k ) a ( k r )$$

using

Note that a (0) = 0 since

$$\lim _ { x \to 0 } \frac { \sin x - x \cos x } { x ^ { 3 } } = \frac { 1 } { 3 } \.$$

Similarly, the spectral expression for the third-order structure function can be derived from the relationship between the isotropic third-order correlation function and the transfer spectrum

$$\frac { 1 } { 2 } \left ( 3 + r \frac { \partial } { \partial r } \right ) \left ( \frac { \partial } { \partial r } + \frac { 4 } { r } \right ) C _ { L L, L } ( r ) = \frac { 1 } { 2 r ^ { 2 } } \frac { \partial } { \partial r } \left [ \frac { 1 } { r } \frac { \partial } { \partial r } \left ( r ^ { 4 } C _ { L L, L } ( r ) \right ) \right ] = \int _ { 0 } ^ { \infty } d k \ T ( k ) \ \frac { \sin k r } { k r }$$

as in e.g. Batchelor [57] equation (5.5.14) p. 101, where Batchelor's K r ( ) corresponds to (1 /r 4 )( ∂/∂r r C ) 4 LL,L ( r ) in our notation. After integrating by parts (with respect to r ) one obtains

$$C _ { L L, L } ( r ) & = 2 r \int _ { 0 } ^ { \infty } d k \ T ( k ) \ \left [ \frac { 3 \sin k r - 3 k r \cos k r - ( k r ) ^ { 2 } \sin k r } { ( k r ) ^ { 5 } } \right ] \quad \text{(A11)} \\ & = 2 \int _ { 0 } ^ { \infty } d k \ \frac { T ( k ) } { k ^ { 2 } } \ \frac { \partial a ( k r ) } { \partial r }$$

$$\bar { \ } J _ { 0 } \quad ^ { \omega \nu } k ^ { 2 } \quad \partial r$$

see e.g. Monin and Yaglom [50] vol. 2, equation (12.141 ′′′ ). The spectral expression for the third-order structure function follows directly

$$S _ { 3 } ( r ) = 6 C _ { L L, L } ( r ) = 1 2 \int _ { 0 } ^ { \infty } d k \ \frac { T ( k ) } { k ^ { 2 } } \ \frac { \partial a ( k r ) } { \partial r } \.$$

$$a ( x ) = \frac { 1 } { 3 } - \frac { \sin x - x \cos x } { x ^ { 3 } } \.$$

The local slopes ζ 1 and ζ 2 can now be found by taking derivatives of the spectral forms for the structure functions

$$\text{for the structure functions} \\ \zeta _ { 2 } ( r ) & = \frac { r } { S _ { 2 } ( r ) } \frac { \partial S _ { 2 } ( r ) } { \partial r } = \frac { 4 r } { S _ { 2 } ( r ) } \int _ { 0 } ^ { \infty } d k \ E ( k ) \ \frac { \partial a ( k r ) } { \partial r } \\ & = \frac { 4 } { S _ { 2 } ( r ) } \int _ { 0 } ^ { \infty } d k \ E ( k ) \ \left [ \frac { 3 \sin k r - 3 k r \cos k r - ( k r ) ^ { 2 } \sin k r } { ( k r ) ^ { 3 } } \right ] \ ; \\ \zeta _ { 3 } ( r ) & = \frac { r } { S _ { 3 } ( r ) } \frac { \partial S _ { 3 } ( r ) } { \partial r } = \frac { 1 2 r } { S _ { 3 } ( r ) } \int _ { 0 } ^ { \infty } d k \ \frac { T ( k ) } { k ^ { 2 } } \ \frac { \partial ^ { 2 } a ( k r ) } { \partial r ^ { 2 } } \\ & = \frac { 1 2 r } { S _ { 3 } ( r ) } \int _ { 0 } ^ { \infty } d k \ T ( k ) \ \left [ \frac { ( 5 ( k r ) ^ { 2 } - 1 2 ) \sin k r - \left ( ( k r ) ^ { 2 } - 1 2 ) k r \cos k r } { ( k r ) ^ { 5 } } \right ] \.$$

## Appendix B: Previous attempts to include the effects of forcing in studies of the structure functions

Gotoh et al. [43] studied a 'generalized' KHE equation defined through their equation (27), which is rewritten here as

$$\varepsilon = - \frac { 1 } { 4 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } S _ { 3 } \right ) + \frac { 3 \nu _ { 0 } } { 2 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } \frac { \partial S _ { 2 } } { \partial r } \right ) + \frac { 3 } { 4 r } I _ { G } ( r )$$

with the input defined from their equation (28) in terms of our input term, I ( r ), as

$$I _ { G } ( r ) = 4 r \int _ { 0 } ^ { \infty } d k \ W ( k ) \ \left [ \frac { 1 } { 3 } - \frac { \sin k r - k r \cos k r } { ( k r ) ^ { 3 } } \right ] = \frac { 4 r } { 3 } \left ( \varepsilon _ { W } - I ( r ) \right ) \,.$$

Gotoh et al. retain the dissipation rate in the KHE, despite its origin as ∂U /∂t 2 = 0 for forced turbulence. They then find a correction in order to compensate for the retained dissipation rate, which also includes the work done by the stirring forces. Note that ε = ε W cancels on both sides of the equation. The correction is then approximated, since the forcing is confined to low wavenumbers

$$I _ { G } ( r ) \simeq \frac { 2 } { 1 5 } \varepsilon _ { W } K ^ { 2 } r ^ { 3 }$$

where

$$K ^ { 2 } = \frac { \int _ { 0 } ^ { \infty } d k \ k ^ { 2 } W ( k ) } { \int _ { 0 } ^ { \infty } d k \ W ( k ) } = \frac { 1 } { \varepsilon _ { W } } \int _ { 0 } ^ { \infty } d k \ k ^ { 2 } W ( k ) \.$$

This approximation is plotted in figure 13 of [43], and can be seen to give good agreement to DNS data for most scales. It is also used to express S 3 ( r ) in terms of S 2

$$S _ { 3 } ( r ) = - \frac { 4 \varepsilon r } { 5 } + 6 \nu _ { 0 } \frac { \partial S _ { 2 } } { \partial r } + \frac { 2 } { 3 5 } \varepsilon _ { W } K ^ { 2 } r ^ { 3 }$$

as shown in (their) figure 11. This approach has also been discussed in Kaneda et al. Sirovich et al. [53] use

[58].

$$S _ { 3 } ( r ) = - \frac { 4 \varepsilon r } { 5 } + 6 \nu _ { 0 } \frac { \partial S _ { 2 } } { \partial r } + \frac { 6 } { r ^ { 4 } } \int _ { 0 } ^ { r } d y \ y ^ { 4 } \langle \delta u _ { L } ( y ) \delta f _ { L } ( y ) \rangle$$

where the longitudinal force increment is defined (in a similar manner to δu L ( r )) as

$$\delta f _ { L } ( r ) = \left [ f ( x + r ) - f ( x ) \right ] \cdot \hat { r } \.$$

This approach also retains the dissipation rate alongside a correction term. The integral in (B6) is approximated to give

$$S _ { 3 } ( r ) \simeq - \frac { 4 \varepsilon r } { 5 } + 6 \nu _ { 0 } \frac { \partial S _ { 2 } } { \partial r } + \frac { 2 } { 7 } \varepsilon k _ { 0 } ^ { 2 } r ^ { 3 }$$

with k 0 the forcing wavenumber. This can be compared to (B5) for the result obtained by Gotoh et al. .

## Appendix C: Derivation of the KHE from the Fourier space energy balance

In order to obtain the energy balance equation in real space including the effects of (finite) forcing, we assume isotropy and take the Fourier transform of the Lin equation (18)

$$\lambda & \lambda \lambda & \lambda \lambda & \lambda \lambda \\ \frac { \partial } { \partial t } \int _ { 0 } ^ { \infty } & d k \ E ( k, t ) \ \frac { \sin k r } { k r } \\ & = \int _ { 0 } ^ { \infty } d k \ T ( k, t ) \ \frac { \sin k r } { k r } - 2 \nu _ { 0 } \int _ { 0 } ^ { \infty } d k \ k ^ { 2 } E ( k, t ) \ \frac { \sin k r } { k r } + \int _ { 0 } ^ { \infty } d k \ W ( k, t ) \ \frac { \sin k r } { k r } \.$$

For the work term we obtain

$$W ( r, t ) = \int _ { 0 } ^ { \infty } d k \ W ( k, t ) \ \frac { \sin k r } { k r } = \frac { 1 } { 2 } \left [ \langle u _ { \alpha } ( x ) f _ { \alpha } ( x + r ) \rangle + \langle f _ { \alpha } ( x ) u _ { \alpha } ( x + r ) \rangle \right ] \,. \quad \text{(C2)}$$

For the dissipation term, note that

$$\left ( \frac { 2 } { r } + \frac { \partial } { \partial r } \right ) \frac { \partial } { \partial r } \int _ { 0 } ^ { \infty } d k \ E ( k, t ) \ \frac { \sin k r } { k r } = - \int _ { 0 } ^ { \infty } d k \ k ^ { 2 } E ( k, t ) \ \frac { \sin k r } { k r }$$

thus, using (A10) and (A1), we obtain for the energy balance equation in real space

$$\frac { \partial C ( r ) } { \partial t } = \frac { 1 } { 2 r ^ { 2 } } \frac { \partial } { \partial r } \left [ \frac { 1 } { r } \frac { \partial } { \partial r } \left ( r ^ { 4 } C _ { L L, L, L } ( r ) \right ) \right ] + \frac { 2 \nu _ { 0 } } { r ^ { 2 } } \frac { \partial } { \partial r } \left ( r ^ { 2 } \frac { \partial C ( r ) } { \partial r } \right ) + W ( r, t ) \.$$

By using the relation (A3) between C r ( ) and C LL ( r ), multiplying by 2 r 2 integrating once over r and finally dividing by r 3 we obtain the KHE in terms of the longitudinal correlation functions

$$\frac { \partial C _ { L L } ( r ) } { \partial t } = \frac { 1 } { r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } C _ { L L, L } ( r ) \right ) + \frac { 2 \nu _ { 0 } } { r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } \frac { \partial C _ { L L } ( r ) } { \partial r } \right ) + \frac { 2 } { r ^ { 3 } } \int _ { 0 } ^ { r } d y \ y ^ { 2 } W ( y, t ) \.$$

Now we can insert the definitions of the structure functions and arrive at the KHE relating S 2 and S 3

$$\frac { \partial U ^ { 2 } } { \partial t } - \frac { 1 } { 2 } \frac { \partial S _ { 2 } ( r ) } { \partial t } = \frac { 1 } { 6 r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } S _ { 3 } ( r ) \right ) - \frac { \nu _ { 0 } } { r ^ { 4 } } \frac { \partial } { \partial r } \left ( r ^ { 4 } \frac { \partial S _ { 2 } ( r ) } { \partial r } \right ) + \frac { 2 } { 3 } I ( r, t )$$

where the input term I ( r, t ) is defined as

$$I ( r, t ) = \frac { 3 } { r ^ { 3 } } \int _ { 0 } ^ { r } d y \ y ^ { 2 } W ( y, t ) \.$$

- [1] A. N. Kolmogorov. The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. C. R. Acad. Sci. URSS , 30:301, 1941.
- [2] A. N. Kolmogorov. Dissipation of energy in locally isotropic turbulence. C. R. Acad. Sci. URSS , 32:16, 1941.
- [3] A. N. Kolmogorov. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. J. Fluid Mech. , 13:8285, 1962.
- [4] U. Frisch. Turbulence: the legacy of A. N. Kolmogorov . Cambridge University Press, 1995.
- [5] G. Boffetta, A. Mazzino, and A. Vulpiani. Twenty-five years of multifractals in fully developed turbulence: a tribute to Giovanni Paladin. J. Phys. A: Math. Theor. , 41:363001, 2008.
- [6] R. H. Kraichnan. On Kolmogorov's inertial-range theories. J. Fluid Mech. , 62:305, 1974.
- [7] P. G. Saffman. Problems and progress in the theory of turbulence. In H. Fiedler, editor, Structure and Mechanisms of Turbulence II , volume 76 of Lecture Notes in Physics , pages 273-306. Springer-Verlag, 1977.
- [8] K. R. Sreenivasan. Fluid Turbulence. Rev. Mod. Phys. , 71:S383, 1999.
- [9] J. Qian. Closure Approach to High-Order Structure Functions of Turbulence. Physical Review Letters , 84(4):646-649, 2000.

- [10] H. Effinger and S. Grossmann. Static Structure Function of Turbulent Flow from the NavierStokes Equations. Z. Phys. B , 66:289-304, 1987.
- [11] G. I. Barenblatt and A. J. Chorin. New perspectives in turbulence: scaling laws, asymptotics and intermittency. SIAM Rev. , 40:265-291, 1998.
- [12] S. Gamard and W. K. George. Reynolds number dependence of energy spectra in the overlap region of isotropic turbulence. Flow, turbulence and combustion , 63:443-477, 1999.
- [13] T. S. Lundgren. Kolmogorov two-thirds law by matched asymptotic expansion. Phys. Fluids , 14:638, 2002.
- [14] L. Mydlarski and Z. Warhaft. On the onset of high-Reynolds-number grid-generated wind tunnel turbulence. J. Fluid Mech. , 320:331-368, 1996.
- [15] C. W. van Atta and W. Y. Chen. Structure functions of turbulence in the atmospheric boundary layer over the ocean. J. Fluid Mech. , 44:145, 1970.
- [16] F. Anselmet, Y. Gagne, E. J. Hopfinger, and R. A. Antonia. High-order velocity structure functions in turbulent shear flows. J. Fluid Mech. , 140:63, 1984.
- [17] R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Massaioli, and S. Succi. Extended selfsimilarity in turbulent flows. Phys. Rev. E , 48:R29-R32, 1993.
- [18] R. Benzi, S. Ciliberto, C. Baudet, and G. R. Chavarria. On the scaling of three-dimensional homogeneous and isotropic turbulence. Physica D: Nonlinear Phenomena , 80(4):385-398, 1995.
- [19] D. Fukayama, T. Oyamada, T. Nakano, T. Gotoh, and K. Yamamoto. Longitudinal structure functions in decaying and forced turbulence. J. Phys. Soc. Japan , 69:701, 2000.
- [20] G. Stolovitzky, K. R. Sreenivasan, and A. Juneja. Scaling functions and scaling exponents in turbulence. Phys. Rev. E , 48:R3217, 1993.
- [21] C. Meneveau. Transition between viscous and inertial-range scaling of turbulence structure functions. Phys. Rev. E , 54:3657, 1996.
- [22] S. Grossmann, D. Lohse, and A. Reeh. Application of extended self-similarity in turbulence. Phys. Rev. E , 56:5473, 1997.
- [23] A. Sain, Manu, and R. Pandit. Turbulence and multiscaling in the randomly forced NavierStokes Equation. Phys. Rev. Lett. , 81:4377, 1998.
- [24] F. Belin, P. Tabeling, and H. Willaime. Exponents of the structure functions in a low temperature helium experiment. Physica D , 93:52-63, 1996.

- [25] J. Qian. Inertial range and the finite Reynolds number effect of turbulence. Physical Review E , 55:337, 1997.
- [26] J. Qian. Slow decay of the finite Reynolds number effect of turbulence. Physical Review E , 60(3):3409-3412, 1999.
- [27] J. Tchoufag, P. Sagaut, and C. Cambon. Spectral approach to finite Reynolds number effects on Kolmogorov's 4/5 law in isotropic turbulence. Phys. Fluids , 24:015107, 2012.
- [28] W. J. T. Bos, L. Chevillard, J. F. Scott, and R. Rubinstein. Reynolds number effect on the velocity increment skewness in isotropic turbulence. Phys. Fluids , 24:015108, 2012.
- [29] S. A. Orszag. On the Elimination of Aliasing in Finite-Difference Schemes by Filtering HighWavenumber Components. J. Atmos. Sc. , 28:1074, 1971.
- [30] K. Heun. Neue Methoden zur approximativen Integration der Differentialgleichungen einer unabh¨ angigen Ver¨nderlichen. a Z. Math. Phys. , 45:23-38, 1900.
- [31] J. Jim´nez, A. A. Wray, P. G. Saffman, and R. S. Rogallo. e The structure of intense vorticity in isotropic turbulence. J. Fluid Mech. , 255:65, 1993.
- [32] Y. Yamazaki, T. Ishihara, and Y. Kaneda. Effects of Wavenumber Truncation on HighResolution Direct Numerical Simulation of Turbulence. J. Phys. Soc. Jap. , 71:777-781, 2002.
- [33] Y. Kaneda, T. Ishihara, M. Yokokawa, K. Itakura, and A. Uno. Energy dissipation and energy spectrum in high resolution direct numerical simulations of turbulence in a periodic box. Phys. Fluids , 15:L21, 2003.
- [34] Y. Kaneda and T. Ishihara. High-resolution direct numerical simulation of turbulence. Journal of Turbulence , 7:1-17, 2006.
- [35] C. R. Doering and N. P. Petrov. Low-wavenumber forcing and turbulent energy dissipation. Progress in Turbulence , 101(1):11-18, 2005.
- [36] W. D. McComb, A. Hunter, and C. Johnston. Conditional mode-elimination and the subgridmodelling problem for isotropic turbulence. Phys. Fluids , 13:2030, 2001.
- [37] S. R. Yoffe. Investigation of the transfer and dissipation of energy in isotropic turbulence . PhD thesis, University of Edinburgh, 2012. http://arxiv.org/pdf/1306.3408v1.pdf .
- [38] S. G. Chumakov. Scaling properties of subgrid-scale energy dissipation. Phys. Fluids , 19:058104, 2007.
- [39] S. G. Chumakov. A priori study of subgrid-scale flux of a passive scalar in turbulence. Phys. Rev. E , 78:036313, 2008.

- [40] W. D. McComb. Homogeneous, Isotropic Turbulence: Phenomenology, Renormalization and Statistical Closures . Oxford University Press, 2014.
- [41] L.-P. Wang, S. Chen, J. G. Brasseur, and J. C. Wyngaard. Examination of hypotheses in the Kolmogorov refined turbulence theory through high-resolution simulations. Part 1. Velocity field. J. Fluid Mech. , 309:113, 1996.
- [42] N. Cao, S. Chen, and G. D. Doolen. Statistics and structures of pressure in isotropic turbulence. Phys. Fluids , 11:2235-2250, 1999.
- [43] T. Gotoh, D. Fukayama, and T. Nakano. Velocity field statistics in homogeneous steady turbulence obtained using a high-resolution direct numerical simulation. Phys. Fluids , 14:1065, 2002.
- [44] D. A. Donzis, K. R. Sreenivasan, and P. K. Yeung. Scalar dissipation rate and dissipative anomaly in isotropic turbulence. J. Fluid Mech. , 532:199-216, 2005.
- [45] P. K. Yeung, D. A. Donzis, and K. R. Sreenivasan. Dissipation, enstrophy and pressure statistics in turbulence simulations at high Reynolds numbers. J. Fluid Mech. , 700:5-15, 2012.
- [46] G. I. Barenblatt, A. J. Chorin, and V. M. Prostokishin. Comment on the paper 'On the scaling of three-dimensional homogeneous and isotropic turbulence' by Benzi et al. Physica D , 127:105-110, 1999.
- [47] R. Benzi, S. Ciliberto, C. Baudet, and G. Ruiz-Chavarria. Reply to the comment of Barenblatt et al. Physica D , 127:111-112, 1999.
- [48] P. Sagaut and C. Cambon. Homogeneous Turbulence Dynamics . Cambridge University Press, Cambridge, 2008.
- [49] T. Ishihara, T. Gotoh, and Y. Kaneda. Study of high-Reynolds number isotropic turbulence by direct numerical simulation. Ann. Rev. Fluid Mech. , 41:165, 2009.
- [50] A. S. Monin and A. M. Yaglom. Statistical Fluid Mechanics . MIT Press, 1975.
- [51] K. R. Sreenivasan. On the universality of the Kolmogorov constant. Phys. Fluids , 7:2778, 1995.
- [52] P. R. Bevington and D. K. Robinson. Data Reduction and Error Analysis for the Physical Sciences . McGraw-Hill, third edition, 2003.
- [53] L. Sirovich, L. Smith, and V. Yakhot. Energy spectrum of homogeneous and isotropic turbulence in far dissipation range. Phys. Rev. Lett. , 72:344, 1994.

- [54] R. A. Antonia and P. Burattini. Approach to the 4/5 law in homogeneous isotropic turbulence. J. Fluid Mech. , 550:175, 2006.
- [55] Until 2014 HECToR was the UK's national HPC resource. More information can be found at http://www.hector.ac.uk/ .
- [56] More information on ECDF can be obtained at http://www.ecdf.ed.ac.uk/ .
- [57] G. K. Batchelor. The theory of homogeneous turbulence . Cambridge University Press, Cambridge, 1st edition, 1953.
- [58] Y. Kaneda, J. Yoshino, and T. Ishihara. Examination of Kolmogorov's 4/5 Law by HighResolution Direct Numerical Simulation Data of Turbulence. J. Phys. Soc. Japan , 77:064401, 2008.
W. D. McComb S. R. Yoffe M. F. Linkmann A. Berera
Spectral analysis of structure functions and their scaling exponents in forced isotropic turbulence
The pseudospectral method, in conjunction with a new technique for obtaining
scaling exponents $\zeta_n$ from the structure functions $S_n(r)$, is presented
as an alternative to the extended self-similarity (ESS) method and the use of
generalized structure functions. We propose plotting the ratio
$|S_n(r)/S_3(r)|$ against the separation $r$ in accordance with a standard
technique for analysing experimental data. This method differs from the ESS
technique, which plots $S_n(r)$ against $S_3(r)$, with the assumption $S_3(r)
\sim r$. Using our method for the particular case of $S_2(r)$ we obtain the new
result that the exponent $\zeta_2$ decreases as the Taylor-Reynolds number
increases, with $\zeta_2 \to 0.679 \pm 0.013$ as $R_{\lambda} \to \infty$. This
supports the idea of finite-viscosity corrections to the K41 prediction for
$S_2$, and is the opposite of the result obtained by ESS. The pseudospectral
method also permits the forcing to be taken into account exactly through the
calculation of the energy input in real space from the work spectrum of the
stirring forces.