Generalized self-similar spectrum and the effect of large-scale in decaying homogeneous isotropic turbulence

In statistically stationary conditions, the turbulent energy spectrum in a high Reynolds number flow exhibits a k−5/3 (Kolmogorov) regime, with a faster decay at larger, dissipative wavenumbers. Here, we investigate how the energy spectrum of a turbulent flow evolves in time when turbulence decays freely, in the absence of forcing. Results from direct numerical simulation of decaying turbulence in a periodic box with several different initial conditions suggest a generalized self-similar spectrum, depending on k s = k × η ( t ) and k l = k × L ( t ) , where η(t) and L(t) are, respectively, the small (Kolmogorov) and large scales of the flow. A closure method allows us to obtain an explicit form of the spectrum, which reproduces the deviations from the Kolmogorov spectrum at small k observed numerically. The solution can also be used to determine the second and third order structure functions in decaying turbulent flows, and to study their scaling behavior. Our results compare favorably with high-Reynolds number wind tunnel data. Moreover, our theoretical results provide support to the interesting empirical observation by Pearson et al (2002 Phys. Fluids 14 1288–90) that, independent of the large scale structure of the flow, the dimensionless energy dissipation rate is a universal constant when scaled in terms of the turbulent kinetic energy of the flow, and of the length scale corresponding to the peak of the compensated energy spectrum.


Introduction
Fluid turbulence results from stirring a fluid, typically at large scales. The injected energy is dissipated by viscosity, which acts at much smaller scales. In a statistically stationary state, there is a constant net flux of energy from the energy injection scale down to the dissipation scale. Over this range of scales, energy 'cascades', i.e., is transferred through scales by nonlinear interactions without loss [1], as first clearly stated by Richardson [2]. Kolmogorov, using the equation describing the energy balance in a turbulent flow [3], later understood from first principles the meaning of this cascade, and the importance of the flux of energy through scales, ε [4]. It is worth pointing out that the ideas used to describe fluid turbulence also apply to a wide class of systems involving a broad range of length-scales, in the general framework of wave-turbulence [5], with applications to magnetohydrodynamics [6], plasma physics [7], oceanography [8], to name a few examples.
For stationary turbulence, the energy spectrum in the inertial range, i.e., the distribution of turbulence kinetic energy per scale, between the energy injection scale and the dissipation scale, can be obtained from the balance between energy cascade (transfer) and dissipation, resulting in the celebrated Kolmogorov's −5/3-law: e µ -( ) E k k 2 3 5 3 , as shown in many textbooks [9][10][11][12][13]. Note that in the steady state, ε is equal to the energy dissipation rate. Viscous effects prevail at scales smaller than the Kolmogorov scale, h n e º ( ) complicated than in stationary flows. In fact, the simple balance between energy transfer and dissipation, which is crucial to describe steady-state situations, breaks down [18][19][20]. We focus here on decaying turbulence, in the absence of energy injection. This is the most widely studied case of non-stationary turbulence, partly because of its apparent simplicity, and of its relevance to practical situations, such as turbulence behind a grid in a wind tunnel. As a result, many theoretical [21][22][23][24][25][26][27][28][29][30][31][32], experimental [21,[33][34][35][36], and numerical [37][38][39] studies have been devoted to decaying turbulence. The subject has been extensively reviewed in classical textbooks [10,13,40,41]. One of the lessons learned over the years is that energy decay does depend on the initial conditions and on the geometric constraints on the flow. For theoretical understanding, an enticing assumption is that some kind of self-similarity or self-preservation occurs [3,14,16,17], so the turbulence can be characterized by the evolution of a length scale and a velocity scale. Even with this simplifying assumption, there is no guarantee that the entire spectrum can be simply rescaled by a single length scale [9]. In decaying turbulence, the equation describing the energy balance [3,4] involves an unsteady term, in addition to the energy transfer and the energy dissipation terms. It has been shown that, behind a grid, all three terms of the energy balance equation need to be explicitly taken into account, as energy dissipation and energy transfer never properly balance [18][19][20]. This observation makes it very unlikely that a single length (and velocity) scale is sufficient to describe the energy spectrum or other statistics, since in principle the balance of any two physical mechanisms would fix a length scale. From direct numerical simulations (DNS) of decaying turbulent flows at moderate spatial resolution (N=768 grid points), starting from a spatially wellresolved turbulent flow in a steady state at R λ ≈230, we observe that in decaying turbulence, the velocity spectrum deviates from the standard Kolmogorov k −5/3 spectrum at small values of k. The difference is related to the energy flux, which varies as a function of the wavenumber in spectral space. We find, however, that the spectra at different times can be superposed, by scaling the wavenumber k by L(t), the instantaneous integral scale of the flow.
To capture this nontrivial structure, we propose here a generalized self-similar form of the energy spectra for decaying turbulence, which involves two length scales, namely L(t) and η(t). Our theoretical approach to determine the evolution and the functional form of the spectra is based on the energy balance equation in Fourier space, the so-called Lin equation [17,42]. The problem is amenable to elementary analytic solutions by using a closure that relates the energy transfer with the spectrum [43]. The resulting solution describes semiquantitatively not only the results of our own DNS of decaying turbulence, but also several important features from recent measurements in the variable density turbulence tunnel (VDTT) in Göttingen [36,44,45]. We note that corrections of the energy spectrum at low k have been already proposed theoretically [13,31,46], but with structures different from that predicted by our approach.
We begin, in section 2, by presenting the results of DNS of decaying turbulent flows, which exhibit clear deviations with respect to the steady-state solution in the presence of a large-scale forcing. In section 3, we explain the assumptions used to simplify the equations of motion. The solution of the problem is presented in section 4. We then return to the DNS results to demonstrate the quality of the model prediction, see section 5. We finally explore, see section 6, the implications of our models for the second and third order structure functions. Last, section 7 contains a brief summary of our results.

DNS of decaying turbulent flows
In this section, we briefly explain how we numerically simulated an ensemble of decaying flows, and present our main observations from the statistical analysis of the results.

Direct numerical simulations (DNS)
We numerically generated a flow by solving the Navier-Stokes equations: in a periodic box of size (2π) 3 . To this end, we used standard pseudo-spectral methods, with 768 3 Fourier modes. The code is dealiased using the standard 2/3-rule [47], which means that the largest faithfully simulated wavenumber is k max =256. A time-dependent forcing fis applied at the smallest wavelength ( < | | k 1.5) by injecting a constant amount of kinetic energy per unit time, ε i , as described in [48]. At the steady state, the amount of dissipated energy, ε, equals the rate of energy injection: ε=ε i . Here we present results obtained from turbulent flows that are initially at a Reynolds number R λ ≈230, where º being the mean kinetic energy per unit mass of the flow. With our choice of parameters, the value of the Kolmogorov length scale at the initial state is such that k max ×η≈1.7, which ensures that spatial resolution of the flow is not an issue since the Kolmogorov scale increases as the turbulence decays. As a consequence of applying the forcing at the smallest wavenumbers in the simulation, the integral scale of the flow  ( » 1.3) is comparable to the system size (2π). Whereas this could be seen as a shortcoming of our simulations [49], we note that it is known that the decay of turbulence notoriously depends on the specificities of flow at the largest scales [13,37], a property that will be discussed more thoroughly in section 3.4.
Starting with a DNS of steady-state turbulent flow at R λ ≈230, we suddenly turned off the energy injection rate by letting ε i =0 at a time, which we refer to in the following as t=0, and we monitored the decay of the flow. The procedure was repeated 8 times with different initial steady flow fields at R λ ≈230, and we averaged over the different flow realizations to obtain a description independent of the specific characters of a particular initial condition. In what follows, we refer to the energy and energy dissipation, averaged over the 8 different runs as˜( ) E t and ẽ( ) t . To characterize the variations ofẼ and ẽ over different runs, we also estimated the standard deviations of the fluctuations from run-to-run of energy and dissipation, δE(t) and δε(t).
This work starts from the dynamics in Fourier space. For this reason, a quantity of prime importance is the energy spectrum, E(k, t), defined in such a way that E(k, t) dk is the amount of kinetic energy at time t, contained in a shell in wavenumber space, of magnitude between k and k+dk. Throughout this text, E(k, t) refers to an ensemble average over many initial conditions. In practice, we determined this quantity by averaging over the eight different runs.

Evolution of energy and dissipation
To characterize the evolution of the turbulence, we use the eddy-turnover time of the flow at the initial state, i.e., before the energy injection was turned off, where L 0 is the velocity correlation length at the initial state, i.e., L 0 =L(t=0) and L(t) is defined as The integral scale, ( ) t , conventionally defined based on the longitudinal velocity correlation, is related to L by  p = ( )L 3 4 . Figure 1 shows the evolution of the average energy and dissipation˜( ) E t and ẽ( ) t with time normalized by T 0 . The energy, normalized by its mean value at the steady state with forcing,˜( ) E 0 , is as expected essentially constant for t<0 and decays thereafter, see figure 1(a). At the end of the simulation, at t≈18T 0 , the energy has decreased by a factor of approximately 35. The Reynolds number at the end of the simulation is R λ ≈80. The variation of energy over different initial conditions, δE(t), as shown in figure 1(b), fluctuates over time at t<0 (when the flow is forced), and decreases monotonically when normalized by˜( ) E 0 . Figure 1(b) shows that in the steady state regime (for t<0), the run-to-run variation d˜( ) E E 0 fluctuates within ≈15%, which provides an estimate of the uncertainty when determining δE. Whereas δE(t) decreases for t>0, d ( )˜( ) E t E t slightly increases during the decay phase, to reach ≈23% at the end of the runs (not shown). Figure 1(c) shows the dissipation of kinetic energy, ẽ( ) t . Contrary to energy, which decays linearly with time starting at time t=0, dissipation decays much slower at first: the derivative ẽ t d d is very close to 0 at t=0. The run-to-run fluctuations of ε, δε (t)/ε i , (not shown) starts at ≈10% at the steady state, and decreases to ≈2% at the end of the runs. As it was the case for the fluctuations of δE, the accuracy of our estimates of δε is approximately 20%.
To characterize more quantitatively the turbulence decay, figure 1(d) shows the evolution ofẼ as a function of t in a log-log representation. The dashed line suggests an asymptotic µ -E t 2 regime. This can be understood theoretically [13] by noticing that in our flow, the integral length scale, ( ) t , initially comparable to the system size, cannot grow significantly as turbulence decays. This is supported by the plot of ( ) t in the inset of figure 1(d). Other decay regimes, reviewed in section 3.4, are possible.

Spectral dynamics
In a homogeneous, isotropic turbulent (HIT) flow, it is convenient to re-write the Navier-Stokes equations, equations (1) and (2) in spectral form. In the absence of any forcing (f=0 in equation (1)), the equation of evolution for the energy spectrum E(k, t) defined in section 2, reads [17,42]: where ∂ k Π(k, t) represents the nonlinear term in equation (1), written as the derivative of a flux term, Π(k, t), from small to large wavenumbers (or equivalently, from large to small scales). One of the main difficulties in studying fluid turbulence comes from the need to find a relation between Π(k, t) and E(k, t), i.e., a closure. It is therefore interesting to consider the evolution of E(k, t) and Π(k, t), determined numerically. We now present our results for the spectral functions E(k, t) and the flux term Π(k, t), averaged over the eight different runs with different initial conditions, as discussed in section 2.2. The compensated average spectra, e -( ) k E k t , i 2 3 5 3 , are shown in figure 2(a), at several times during the decay, in which the time interval between two successive spectra is D » t T 2 0 . The uppermost curve, corresponding to the initial condition, is approximately flat for about a decade of wavenumbers, consistent with a steady state when the energy injection and dissipation are balanced. The decrease of the spectra at t>0 implies the decay of energy. The cutoff wavenumber where the spectrum begins to fall faster than algebraically also decreases with time. In addition, at small values of k the compensated spectrum deviates from being constant, contrary to what is observed in the statistically steady state. On the basis of phenomenological considerations inspired by Kolmogorov theory, it can be expected that in a solution slowly decaying in time, the spectrum at large k is very close to the steady-state spectrum corresponding to the instantaneous energy dissipation rate ẽ( ) t . In particular, one can define at each time t a value of the Kolmogorov length, which we use to compensate the various spectra, as done in a statistically steady state. The result, shown in figure 2(b), shows a remarkably good collapse of the spectra for kη(t)0.3, corresponding to the dissipative range in steady-state situation. On the contrary, the spectra in the small-wavenumber (large-scale) range show a systematic change with time. Namely, in the normalizations in figure 2(b), the compensated spectra at small k still decrease with time, reflecting the reduction of energy transfer with time.
The energy transfer, Π(k, t) in equation (5), is shown in figure 3 for the same times as in figure 2. In figure 3(a), Π(k, t) is made dimensionless by ε i . At t=0, i.e., just before the energy injection is turned off, Π(k, t) is equal to ε i over a decade of values of k, corresponding to a constant flux of energy from large to small scales. As time goes on, the energy flux rapidly decreases, consistent with the already documented decay of the kinetic energy dissipation ẽ( ) t , shown in figure 1(c). As a result, after » t T 10 0 , the values ofΠ(k, t) are very small. To compare with energy transfer in a statistically steady state, figure 3(b) shows Π(k, t) made dimensionless by the instantaneous rate of kinetic energy dissipation, ẽ( ) t , and plotted as a function of kη(t). Consistent with the observations of figure 2(b), all curves superpose very well for kη(t)0.3, which corresponds to the dissipative range of scales in stationary steady flows. At smaller wavenumbers, or larger scales, the energy transfer also manifestly decreases, as shown in figure 3(b). The very good superposition of E(k, t) and Π(k, t) at large values of k, observed in figures 2(b) and 3(b), when using only ẽ( ) t (equivalently, η(t)) to scale the variables, points to an effective 'equilibration' of the small-scale motion in decaying turbulent flows. This is consistent with the phenomenological notion that the characteristic time scales of eddies decrease rapidly with their size, so the small-scale eddies rapidly adjust, in a statistical sense, to energy flux changes and are effectively in a quasi-steady state. On the other hand, the larger-scale eddies do not instantaneously relax to the steady state. The form of in figure 3(b) shows that energy transfer is  reduced at small wavenumbers, the more so as time goes on. This is certainly related to the deviation of the spectrum E(k, t) from its steady-state e k 2 3 5 3 form, which also becomes more pronounced with time. These observations provide a justification for the approximate theoretical treatment presented in the following section.

Theoretical description: simplifying approximations
Whereas solving equation (5) for decaying turbulence is a challenging task, the following simplifying assumptions, inspired in part by observations above, provide some analytic insights on the nature of the solution. We discuss these assumptions in turn.

A generalized self-similarity assumption
The existence of distinct ranges of wavenumbers, clearly deduced from the observation that the small-scales reach a steady state, while eddies of larger scales keep evolving in a nontrivial manner, is manifestly a consequence of the multiplicity of time scales in the problem. Using a theoretical treatment rather standard in nonlinear dynamics, we introduce two variables for the wavenumber, namely a 'slow variable', characterizing the large scales, º ( ) k kL t l , and a 'fast variable', characterizing the small scales, where L(t) and η(t) are defined by equations (4) and (6), respectively. The standard use of 'fast' and 'slow' variables comes from the expression for the partial derivative: l s , which makes it clear that variations in k l (or k s ) are associated with weak (or strong) variations as a function of k. Further, we merely assume that the spectrum can be written using an Ansatz of separation of variables: . 7 l l l s s 5 3 The functions ( ) f k l l and ( ) f k s s are both dimensionless. Their exact forms are to be determined later. We note that similar functional forms have been postulated to discuss the energy spectrum in a steady-state flow [9,50]. When f l =f s =1, equation (7) is nothing but the celebrated k −5/3 spectrum. We chose the boundary conditions for f l and f s to be = With these boundary conditions, E(k, t) behaves as k −5/3 over the wavenumber range 1/L=k=1/η. We note that the total kinetic energy is

Pao's closure
The equation for E(k, t) involves, in an essential way the energy transfer term Π(k, t), which reflects the effect of nonlinearity in the Navier-Stokes equations, equation (1). This term is in general unknown. Among many attempts to model the energy flux (see e.g. [10,43,51]), the one proposed by Pao [43] appears as the simplest, and possibly the best in terms of reproducing the properties of the energy spectra in steady-state flow, for inertial range scales, up to kη≈0.5 [13]. The main assumption by Pao is that the energy transfer is local in k-space, hence it involves only E(k, t) and k. Its expression can then be obtained from dimensional considerations: , , where C is a universal dimensionless constant. The determination of its precise value will be discussed later, see section 5.

Time scales
In a turbulent flow, eddies of size ℓ are evolving with a characteristic time e » ℓ ( ) ℓ t 2 1 3 , resulting in a wide range of time scales. Here, we briefly discuss the time scales associated with the spectral equation, equation (5).
To proceed, we notice that as seen in figures 1(a) and (c), quantities such as the total energy,˜( ) E t , and the instantaneous dissipation rate, ẽ( ) t , both evolve over a time of the order of the large eddy turnover time, is the velocity correlation length defined by equation (4) and e = -Ẽ t d d . This observation is consistent with earlier reasoning [11,12], and implies that the Kolmogorov  . Our statement is merely that the length scale η(t), which is defined by an averaging process, varies with a slow characteristic time scale T(t). A plausible physical picture is that as time goes on, there are fewer and fewer eddies at a given very small scale because the energy transfer rate is slowly decaying, so statistically the length scale η(t) slowly grows.

Decay regimes
In the absence of external forcing, the kinetic energy of the flow decays, typically with a power law [13]. The existence of power laws suggest that T t t are all constants of order unity, where the dot '˙' means time derivative. We will use this as an assumption in the analysis below. Note that the quantity K(t) introduced in equation (7) is µḰ E L, up to numerical constants, which implies that , several regimes have been proposed in the literature, with different values for the exponents n and m. We note that, with our notation, In general, n and m are related by the relation n+2m=2, which results from e = -Ẽ t d d and the well-known relation e μẼ L 3 2 (consistent with our own data, see figure 5 below). As a consequence, the knowledge of the ratio A L /A K =−p is sufficient to fully determine the exponents m and n.
The results in section 2 clearly correspond to m=0 since the large-scale of the flow is bounded by the size of the simulation box. Therefore p=0, which implies that m=0 and n=2. This is consistent with our own data, see the dashed line corresponding to µ -E t 2 in figure 1(d). From the considerations above, when looking for solutions of equation (5) corresponding to the DNS presented in section 2, we therefore assume that  ( ) ( ) TL t L t 1, a condition very well satisfied in our numerical work. In wind tunnels, the exponents obtained experimentally [36,52] show that the decay corresponds to the regime originally analyzed by Saffman [53,54]: n=6/5 and m=2/5. In this case, A L =2/5 and A K =−4/5 (p=1/2). The classic prediction of Kolmogorov [15] for the decay of turbulence, n=10/7 and m=2/7, corresponds to p=1/4. The analysis presented below applies to all these cases.
We stress that in all cases of interest, A K <0. This can be seen by writing and noticing that < A 0 E since energy is decaying and n1, as expected from theoretical considerations and found experimentally [36,52].

General formalism
We now turn to the spectrum equation equation (5). Dividing both sides of equation (5) by E(k, t) leads to: Using the functional form equation (7) leads to: Pao's closure hypothesis, equation (10), can be used to evaluate the expression of ¶ P( ) ( ) k t E k t , ,  Equation (14) can be rearranged in the following form: - and another one involving only k s : These equations are first order in k l and k s , and are completely determined once the functions of time, K(t), L(t), η(t) and Φ(t) are known. We now justify that the function Φ(t) turns out to be 0.
, and we note that in equation (17), f s ≈1 and ¢  f 1 s , due to the boundary condition equation (9). We also use the fact that h ḣ ( ) ( ) t t is of order 1/T. Therefore, the first two terms on the LHS of equation (17) are both of the order . The last term on the LHS is of the order n h , again much smaller than 1/T. Hence Φ(t) must be much smaller than 1/T, while one can easily see that all the terms on the LHS of equation (16) are of order 1/T. Combining equations (16) and (17), we therefore set, as a first approximation, Φ(t)=0, since it appears on the RHS of both equations.
Equations (16) and (17), together with the boundary conditions equations (8) and (9), can be explicitly solved, by taking into account Φ(t)=0. We begin by re-writing equation (16)  which has the solution: We note that ( ) f k l l is an increasing function of k l , as A K <0. Similarly, equation ( whose solution reads:

Solution for A L =0
Motivated by the numerical results on decaying turbulence discussed in section 2, for which the integral length scale, set by the size of the system, does not vary much with t, we begin by presenting the solution in the case A L =0. By imposing the condition that the total energy is , we obtain the following expression for the spectrum: The constant I is obtained from the integral of the expression (25). We find that = In conclusion, the energy spectrum can be written as and when  Re 1, i.e., L/η?1, it reduces to a simpler form: We recall that this expression for the spectrum is obtained using the assumptions listed in section 3, including the fact that the velocity correlation length L(t) hardly varies in a decaying flow. These assumptions and the spectrums equations (26) and (27) will be checked in section 5 against DNS data.

Solution for
We now consider the more general case, with ¹ A 0 L , which is relevant in particular to turbulence produced behind a passive grid in a wind tunnel, where experiments demonstrate that L grows as a power law [36,52], consistent with the picture proposed in [53,54], which implies a value of A L =−A K /2.
The solution for ¹ A 0 L and L/η?1 can be obtained by substituting equations (19) and (24) into equation (7). The result reads: which can be re-expressed as: , whereas turbulence decay in a wind tunnel leads to an exponent −A K /A L =2, consistent with previous work [37].
We note that the limit A L → 0 can be recovered directly from equation (28). Namely: which is exactly equation (27). In the last step of the derivation of equation ( Equations (27) and (29) show that in all cases, the spectrum is an increasing function of k in the very small k limit.

Comparison with other analytic forms
The corrections predicted in [31,46] constant. This simple form can be obtained by purely dimensional considerations, in the spirit of [56,57]. Our work provides an explicit form for these corrections. We note that an interesting form of the energy spectrum has been proposed by Davidson [13] using Pao's conjecture:  [13], where m is flow dependent and its value is in principle determined by the boundary conditions such as~-E k m 2 3 5 3 when k → 0. Typical values are m=17/2 or 11/2 for the Batchelor or Saffman spectra.
The derivation of the solution in [13], sketched above, differs from ours in a number of ways. We use the time-dependence of the total kinetic energy and of the velocity correlation length of the flow as input parameters, and construct the spectrum by hypothesizing a (generalized) self-similar form. In contrast, the solution in [13] does not explicitly take into account the variation of the total kinetic energy of the flow. This explains the different functional forms for the spectrum in the inertial range. Also, the determination of the constants in equation (31) is not straightforward. In comparison, the constants involved in our expressions of the spectrum can be directly determined from the numerical results, as we show below. This enables a quantitative comparison of our predictions with the DNS results.

Verification of the assumptions and the model spectrum
Among the assumptions listed in section 3, Pao's closure equation (10) is certainly the most difficult to justify from first principles. In order to check its validity, we show in figure 4 the ratio e P ( ) ( ) k E k t k t , , 1 3 5 3 obtained from our numerical study, which, according to Pao's conjecture, should be a constant and equal 1/C. Figure 4 shows that the ratio from DNS is not a constant, but varies slightly as a function of kη(t). The variation of the ratio, which is found around 3, remains within ±50%, over the entire range of wavenumbers considered.
As already pointed out, the constant C in equation (10) is not known from first principles. Its value can be estimated by looking at the decay rate of the energy, as we now show. versus wavenumber kη at different t, which is expected to be a constant, and the model prediction is 1/C≈2.9, shown as the horizontal dashed line.
With our definition, the energy dissipation is simply: As in our simulations, we consider the case where the velocity correlation length L, defined by equation (4) which leads to: Moreover, a relation between A K and C can be obtained from the definition of L, equation (4), which gives, in the case of L/η?1: When evaluating L to obtain equation (36), we have set to 1 the dissipative term in the expression of the spectrum E(k, t), equation (27), which amounts to setting h  ¥ L or ¥ Re . A more precise expression for the relation between A K and C can be obtained in terms of special functions when L/η is finite, as explained in appendix A. Now, combining equations (35) and (36) allows us to determine the constant C from the dimensionless ratio e º e˜C L E 3 2 , which, in stationary flows, is independent of the Reynolds number when the Reynolds number is high enough [58,59]. We note that C is the only constant in our theory that needs to be determined from empirical data. In addition, as it is the constant in Pao's conjecture on local energy transfer, it is expected to be universal. We therefore determine it from our DNS data with L fixed and assume it remains the same for other decaying flows even if ¹ A 0 L . As we will see later, the results from our theory agree well with experimental data, suggesting that this underlying assumption is acceptable.
The DNS results shown in figure 5 indicate that after a short transient, for t/T 0 5, the ratio e »L E 0.21 3 2 is almost constant, which, when substituting into equation (35), gives −A K ≈0. 35. We note that the values of A L in DNS are less than 2×10 −3 over the range where A K is approximately constant. We then use the expression in appendix A for R λ ∼O(100), as appropriate in the DNS shown in section 2, to determine the constant C to be C≈0. 35. The corresponding value of 1/C≈2.9 is close to the value obtained directly from checking Pao's conjecture, as shown in figure 4 6 .
Next, we demonstrate that the expressions using our generalized self-similar assumption provide a good approximation of our numerical data. We recall that figures 2(b) and 3(b) show good collapses of the spectra, E (k, t), and of the transfer, Π(k, t), as a function of k×η(t), i.e., in the dissipative range, which is an justification for using the functional form equation (7). Our analysis in section 4 motivates us to superpose E(k, t) and Π(k, t) in the large scales using the variable k×L(t). Figure 6(a) shows the energy spectrum, E(k, t), divided by  Comparison between the predictions of the model and the numerical results for the spectrum E(k, t) (a) and for the energy transfer, Π(k, t) (b). A perfect agreement between the theoretical predictions and the numerics would correspond to a ratio of 1, shown as a thick horizontal line in both panels. At t/T 0 5, the agreement revealed by the figure is within 40% for the spectrum and for the energy transfer spectrum, except at the largest values of k where approximation equation (10) is not justified. 6 We note that there is a subtle relation between the constant C in our model and C K , the so-called 'Kolmogorov constant' for the spectrum in the inertial range. Pao's conjecture, applied to the inertial range only, leads to C=1/C K ≈0.6 or 1/C≈C K ≈1.5, using the generally accepted value of C K in statistically stationary conditions [12]. Here, however, we apply Pao's conjecture throughout the wavenumber range, and determine the constant C using the rate of kinetic energy dissipation, i.e., the 'zeroth-law' constant C ε . Using Pao's conjecture throughout results in a spectrum that differs from the true one in the dissipative range [13], which in turn affects the accuracy of the relation between C and C ε derived from our theory. The value of 1/C shown in figure 5, determined using our approach, could be understood as an average of the normalized flux from Pao's conjecture over the entire wavenumber range. As a result of the shortcomings of the Pao's assumption, the value of 1/C differs from the generally accepted value of C K in the inertial range. On the other hand, it can be expected that at Reynolds numbers larger than those in the present DNS, the effect of the error due to dissipative range would be smaller, and the value of C determined using our approach would become closer to the value C=1/C K , obtained from the spectrum in the inertial range. ( ) E t L t , as a function of kL(t), at different times. The spectra superpose quantitatively well. Note that, although the various spectra at t/T 0 5 seem to be well-represented by a straight line once plotted on a log-log scale, this straight line differs significantly from a k −5/3 spectrum. Similarly, the energy transfer function at different times, Π(k, t), normalized byẼ L 3 2 and plotted as a function of kL, superpose well at small k, see figure 6(b).
An even more stringent test is the direct comparison between the spectrum from DNS and the functional form, equation (26), obtained in section 3. Figure 7(a) shows the ratio between the spectrum predicted by the model, E M , and the spectrum determined numerically, E DNS . After a time of the order of » T 5 0 , the ratio varies between 0.8 and 1.2 in the inertial range. Figure 7(b) plots the ratio between the energy transfer function predicted by the model, Π M , and Π DNS determined numerically, which are also very close to 1 for  t T 5 0 . We note that the predictions of the model significantly differ from the numerical results in the dissipative range. This is not a surprise, as it is known that Pao's approximation does not provide accurate results for the energy spectrum in the dissipative range [13]. The predictions, however, are quantitatively quite good throughout the inertial range [13]. 6. Implications of the model 6.1. Predictions of the second and third order structure functions Our analysis so far has been carried out in wavenumber space. A discussion of the corresponding problem in terms of the structure functions in real space is particularly useful for comparing with experimental measurements. In principle, these functions can be determined by solving directly the Kármán-Howarth equation [3] (see also the recent work of [60] and references therein).
With the spectrum E(k, t) and the transfer function Π(k, t) determined, we can deduce the second and third order structure functions, S 2 (r) and S 3 (r), defined as  3 . In contrast, the structure functions of order 4 and higher cannot be determined from the present analysis.
In  approaches, but has not yet reached 4/5 even at R λ =2000. In model predictions, we have assumed that L is approximately constant. the generalized self-similar decay phase. This is due to the deviations from the steady state regime, and also to the reduced Reynolds number as the turbulence decays. Figure 8(b) shows e -( ) (˜) S r r 3 obtained from equation (39) using the model prediction of the spectrum E(k), equation (26), at different Reynolds numbers, which we incorporated into the model through the scale ratio L/η. We observe that, consistent with the agreement in spectral space reported in the previous section, the structure function e -( ) (˜) S r r 3 from DNS is well reproduced by our theoretical prediction. The predictions from the model at much higher Reynolds number, R λ =1000 and R λ =2000, shown in figure 8, indicate that the plateau at 4/5 is hardly reached even at these very large values of R λ [61]. which give direct access to the scaling exponents and can be used to evaluate how the structure functions approach the expected scaling behavior. Figure 9 shows the exponents e 2 and e 3 as a function of  r at several Reynolds numbers, as obtained from our model. The calculation was carried out by considering A L =−A K /2, corresponding to decaying turbulence in a wind tunnel [36,52]. The deviations from the expected exponents e 2 =2/3 and e 3 =1 when r is comparable to  are the consequence of the deviations of the spectra from the k −5/3 regime, see section 2.3. As the Reynolds number increases, the model prediction for the local slope e 3 moves downwards towards 1, and the curve becomes flatter. However, we still could not observe a true plateau at 1 even at the highest Reynolds number tested. A similar situation occurs in the high-Reynolds number experiments at the VDTT in Göttingen [44]. The results of our analysis are semi-quantitatively consistent with the experimental measurements, at least when r is not too small, i.e., when r remains in the inertial range (compare figure 9(b) with figure 4.3 of [44]). When r is in the dissipative range, for reasons already explained, the predictions of the model are not correct. Figure 9 shows that as R λ increases, the values of e 3 return faster to the predicted constant form, 1, than e 2 returns to its limiting value, 2/3.
To quantify how far are the behaviors of the structure functions from a power law in the inertial range, we evaluate the quantity defined in [44,61]: When the structure function S n behaves as a true power law, δ n =0. Figure 10 shows experimental data points from wind tunnel measurements in [44], together with the model prediction of δ 3 as a function of R λ , which were determined by taking = - , as it was the case for figure 9. Both wind tunnel data and model prediction show that δ 3 approaches 0 as the Reynolds number increases. It is interesting to note numerical calculations from our model prediction show that both −δ 2 and −δ 3 decay with R λ apparently as power-laws:  figure 11, which provides a quantitative measure that S 3 (r) approaches the power-law scaling faster than S 2 (r). dlog dlog 3 , for R λ =300, 610, 860, 1320 from above to below, respectively. These structure functions have been obtained for a decaying turbulence with A L /A K =−1/2. Note that here we normalize the length by the integral scale  as done in [44], which is related to the velocity correlation scale L by  p = ( )L 3 4 .

On the dissipation rate constant
It has been found empirically that the ratio e º is sometimes called 'the zeroth law of turbulence' [62]. However, a compilation of available data showed that C ε seems to be flow dependent and is not universal [58]. Based on experimental observation, Pearson et al [59] suggested that instead of using the integral scale  in the definition of C ε , using L p , which is the length scale corresponding to the peak of the compensated energy spectrum ( ) kE k , gave e º e˜C L E p p , 3 2 that is more universal. Note that Pearson et al [59] actually used the compensated one-dimensional spectrum ( ) k E k 1 11 1 as that is experimentally accessible in wind tunnel measurements using hot-wires. Here we show that our model spectrum for generalized self-similar decay provides theoretical supports to the empirical observations of [59].
We first note that from e = -= -˜Ẽ  Figure 10. d 3 versus R λ , blue line is the model prediction, green circles and red diamonds are the experimental data from the VDTT tunnel using different measurement techniques [44]: P11 are data measured with a conventional hot-wire while NSTAP are data from a 'nano-probe' manufactured using micro-fabrication technology. Figure 11. Variation of −δ 2 (red crosses), and −δ 3 (blue circles) with R λ , shown on a log-linear plot. The inset shows the same data, but on a log-log plot.
which gives It shows that even when R λ ? 1, C ε depends on the ratio -= p A A L K , i.e., on how exactly the turbulence is decaying. This is consistent with the observations reviewed in [58].
, using the spectrum found for the general case ¹ A 0 L , equation (28), a straightforward derivation shows that kE(k) peaks at Thus if we define L p =1/k p and e º e˜C L E p p , which depends only on the universal constant C and is independent of A L /A K . In other words, e C p , is the same for all decaying turbulent flows, at least in the self-similar regime.

Conclusions
In summary, our numerical simulations of decaying turbulence, consistent with recent measurements of the energy spectrum in a wind tunnel, show deviations from the k −5/3 energy spectrum at small k (large scales). These deviations are related to the scale dependence of the energy flux in decaying turbulence. Numerically, we observe that the energy spectra and energy transfer spectra collapse very well both at large wavenumbers when expressed as a function of k×η(t), where η(t) is the Kolmogorov length (h n e = (˜( )) t 3 1 4 ), and at small wavenumbers when expressed as a function of k×L(t), where L(t) is the velocity correlation length, defined by equation (4).
To describe these effects, we have proposed a generalized self-similar formalism for decaying HIT flows. Instead of using only one single length scale, we proposed to write the energy spectrum with two length scales: one for the energy-containing range and the other for to the dissipation range.
Analytical expressions for ( ) E k t , and Π(k, t) can be obtained by further closure assumptions. Here, we adopt Pao's hypothesis. The second and third-order structure functions can also be predicted by transforming E(k, t) and Π(k, t) into real space. The predictions of the model reproduce semi-quantitatively the DNS results of decaying turbulence. Preliminary comparisons of our models also show a very encouraging agreement with measurements in wind tunnel experiments at very high Reynolds numbers [44,45]. We expect that further comparisons will provide illuminating information on the turbulent flows in wind tunnels and in other experimental facilities. Finally, we conclude by pointing out that the approach developed in this work may also have application for a broader class of systems involving many length scales.

Acknowledgments
We are very thankful to E Bodenschatz for his support, and for many fruitful discussions. We have also benefited from useful discussions with CKüchler and CVassilicos. The simulations have been performed on the cluster of the Max-Planck Institute for Dynamics and Self-Organization in Göttingen, Germany. PY and HX acknowledge partial support from the National Science Foundation of China (NSFC) under grant number 11672157. AP gratefully acknowledge the support of the Foundation Alexander von Humboldt.  Appendix B. An alternative approach to the prediction of structure functions The determination of the structure functions shown in the previous subsection was based on the spectrum, obtained by solving the von Karman-Lin equations (5). The structure functions are solutions of the von Karman-Howarth equations [3]; see e.g. [60] and references therein. We follow here another approach, based on the following remark. While the second-order structure function S 2 is traditionally used to describe the distribution of kinetic energy among different size of eddies, Davidson [63] argued that the structure function mixes large-and small-scale information, as well as information about energy and enstrophy. The signature function V(r) is introduced to reduce the influence of large scale. It is a real-space function which plays the role of an energy density, similar to E(k). The evolution equation for V(r) has been derived in [63] (see equation (12)): In equation (47), Π V plays the same role as the spectrum transfer function Π in wavenumber space. It is denoted the 'real space kinetic energy flux'. To solve equation (47), we use a generalized self-similar Ansatz, as it was used to solve equation (7), as well as a closure similar to equation (10). Namely, we look for a solution of the form: Together with equation (47), we can obtain two equations for V 1 and V 2 (which are parallel to equations (16) and (17)): then analytical expressions for V(r, t) and Π V (r, t) can be solved from these two equations.