Finite Reynolds number corrections of the 4/5 law for decaying turbulence

Finite Reynolds number contributions from the unsteady and viscous terms to the inertial range solution of the third-order structure functions are examined.


I. INTRODUCTION
From the Kármán-Howarth equation [1], Kolmogorov [2] derived an exact equation for the second order structure function D 2,0 = ( u 1 ) 2 , u 1 = u 1 (x 1 + r,x 2 ,x 3 ) − u 1 (x 1 ,x 2 ,x 3 ) with the separation vector r i = (r,0,0). This procedure was generalized by Hill [3], who systematically derived structure function equations for arbitrary orders directly from the Navier-Stokes equations. These equations written for the second order structure functions are the basis for the following analysis and are given by where ε is the mean dissipation, ν the kinematic viscosity, and D 0,2 = ( u 2 ) 2 the transverse second order structure function with u 2 = u 2 (x 1 + r,x 2 ,x 3 ) − u 2 (x 1 ,x 2 ,x 3 ). In the derivation of Eqs. (1) and (2), isotropy has been used which results in vanishing pressure correlations (cf. [4]). D 3,0 = ( u 1 ) 3 and D 1,2 = u 1 ( u 2 ) 2 are the third order longitudinal and transverse structure functions, respectively. These two equations are analogous to the equation derived by Kolmogorov. Kolmogorov then introduced an inertial range situated in between the large and small scales, for which he postulated a vanishing influence of the viscosity ν and neglected the unsteady term. Integrating Eq. (2) under these assumptions then yields and using this in Eq. (1) one obtains Kolmogorov's famous result, the so-called 4/5 law: Both Eqs. (3) and (4) are exact results under the assumption of isotropy and very large (infinite) Reynolds numbers and have been subject to many experimental and numerical studies.
However, the range for which Eqs. (3) and (4) are found to hold is rather small for experiments at finite Reynolds numbers (see, e.g., [5][6][7][8]). For that reason, modifications of the asymptotic results which include the finite Reynolds number effects were proposed for different kinds of flows.
Lindborg [9] considered the isotropic second order equations of Kolmogorov and kept the unsteady term ∂D 2,0 /∂t. He then proceeded to express it using Kolmogorov's second similarity hypothesis, i.e., he assumed that Employing the k-ε model to solve for ∂ε/∂t and assuming a decay of the kinetic energy k ∼ t −n enabled him to obtain solutions for different kinds of flows and Reynolds numbers with good qualitative agreement of measurements and his model. Lundgren [10,11] used asymptotic expansions to derive the longitudinal third order structure function in the inertial range. He found the same Reynolds number dependence as Lindborg.
Qian [12,13] examined the approach of the 4/5 law in the inertial range using the energy spectrum equation. He found that the asymptotic results are approached rather slowly. Danaila et al. [14] examined the inhomogeneous second order structure function equations of Hill [4] adapted for grid turbulence and looked at the balance of the respective terms. They found that the inhomogeneities contribute significantly for larger r. Similar conclusions were drawn by Zhou et al. [15]. Danaila et al. measured the second order structure function balance for channel flows [16] as well as homogeneous shear turbulence [17] and obtained similar results as for grid turbulence, in the sense that the inhomogeneities are important at large r and also in the inertial range.
Here, we use the isotropic equations, Eqs. (1) and (2), and examine the influence of the unsteady and viscous terms, i.e., their contribution to the inertial range solutions for the third order structure functions. While we keep the unsteady term as did Lindborg, our approach differs inasmuch as we transform Eqs. (1) and (2) into a self-preserving form depending only on a normalized length scale r similarly to the work of Schaefer et al. [18] on the velocity correlation and the unsteady term is reformulated assuming a decay of the kinetic energy k ∼ t −n . Lundgren used the same coordinate transform but neglected the unsteady term. Rather, he matched the leading order terms of the asymptotic expansions of the structure functions for an outer and an inner region. Corrections to Eq. (4) then follow from the second order terms of the expansion.
The paper is organized as follows: After presenting the data used for the analysis in Sec. II, the (normalized) second order structure function equations used here are derived in Sec. III. We then close the resulting system of equations using two different approaches as outlined in Sec. IV. This allows us both to examine the Reynolds number scaling as well as to make predictions about large Reynolds number behavior. This is discussed in more detail in Sec. V, where the balances of the (normalized) system of equations are examined using direct numerical simulation (DNS) data in order to check for the validity of the assumptions made in Sec. IV.

II. DATA SET DESCRIPTION
Direct numerical simulation of decaying homogeneous isotropic turbulence has been performed on the supercomputer JUQUEEN at the research center in Juelich, Germany. The incompressible Navier-Stokes equations are solved in a triply periodic cubic box with size 2π by a pseudospectral 064403-2 method. For numerical stability, the nonlinear term of the momentum equation is rewritten in rotational form. In a pseudospectral method the nonlinear term is computed in real space and transformed to spectral space for temporal integration. Temporal integration is carried out by a low-storage stability preserving third order Runge-Kutta method. The viscous term is treated exactly by using an integrating factor technique. A standard isotropic truncation procedure in combination with a random phase-shift technique is used to eliminate aliasing effects, allowing us to keep all wave numbers with κ < √ 2N/3. The grid resolution is N 3 = 2048 3 , which adequately resolves the smallest scales during the simulation. For pseudospectral methods, the resolution requirement can be written in terms of the nondimensional number κ max η, where κ max is the largest wave number appearing in the truncated Fourier series, and η is the Kolmogorov length. The resolution condition κ max η for the four time steps under consideration is indicated in Table I and has been shown to be sufficient to compute second order velocity gradient statistics [19]. The flow is initialized by a prescribed isotropic energy spectrum of the form where κ is the wave number and κ p is the location at which the initial energy spectrum peaks. In this work we aim at reaching high Reynolds numbers to obtain a well established inertial range. For this reason we set κ p to a comparable small value of 3.5 and tolerate small confinement effects due to the finite size of the computational domain [20]. Following Ishida et al. [20], the initial state of freely decaying turbulence can be characterized by a Reynolds number defined as where k(t = 0) = ∞ 0 E(κ,t = 0)dκ denotes the initial turbulent kinetic energy, and ν is the kinematic viscosity. From this definition and with k(t = 0) = 10 and ν = 0.00021 we obtain an initial Reynolds number of Re = 4302.
We show the temporal evolution of k ∼ t −n and ε ∼ t −n−1 in Figs. 1(a) and 1(b), where the red lines correspond to a decay exponent n = 1.45. This value is slightly larger than the theoretical value n = 10/7 obtained for a κ 4 spectrum as in Eq. (7) (cf., e.g., the discussion in Refs. [21] or [22]). The times used for the present analysis are indicated by the dotted vertical black lines in the decaying regime. We use all four times to examine Reynolds number dependencies of the closures presented below and the highest (leftmost dotted black line, Re λ = 254.75) and lowest (rightmost dotted black line, Re λ = 121.39) Reynolds number for more detailed analysis.

III. UNSTEADY TERMS
To examine the influence of the unsteady terms ∂D 2,0 /∂t and ∂D 0,2 /∂t, we use the identity where f (r,t) is the longitudinal correlation function and similarly for the transverse structure function where g(r,t) is the transverse correlation function Under the assumption of isotropy, u 2 1 = u 2 2 = u 2 . Taking the derivative of Eq. (9) gives and a similar expression for the transverse structure function, Eq. (11). In the following, r is normalized with a large scale in the spirit of [1,11,18], defined as where u = √ u 2 1 and the dissipation ε = 2ν S ij S ij . Then, the normalized length scale and therefore 064403-4 Consequently, and where ∂f ( r, t)/∂ t and ∂g( r, t)/∂ t vanish if f ( r, t) and g( r, t) are self-similar. We see in Sec. IV below that neglecting ∂f/∂ t and ∂g/∂ t has very little impact on the balance equations, i.e., that the assumption is well justified. From Eqs. (9) and (11), we then have as u 2 1 and u 2 2 do not depend on r and where D 2,0 = D 2,0 /u 2 and D 0,2 = D 0,2 /u 2 are the normalized second order structure functions. Similarly for the third order structure functions, For decaying homogeneous isotropic turbulence, the energy balance reduces to where ε = 2ν S ij S ij is the mean energy dissipation. In the self-similar decay state, u 2 = u 2 0 (t/t 0 ) −n , where n is a decay exponent (cf. Fig. 1) and consequently ε = ε 0 (t/t 0 ) −n−1 in agreement with Eq. (20).
Substituting Eqs. (20) and (10) into Eq. (13) and normalizing with ε yields 1 ε and similarly, Finally, from Eq. (14) 1 u Dividing Eqs. (1) and (2) by ε then gives the following equations: where Re L = uL/ν is a large scale Reynolds number. That is, the derivative with respect to t of the second order structure functions has been reformulated in terms of spatial derivatives and the decay of kinetic energy expressed by the decay exponent n. Therefore, the partial differential equations are reduced to ordinary differential equations. This allows the integration of Eqs. (24) and (25) in r, if D 2,0 and D 0,2 are known as functions of r.
Equations (24) and (25) are also valid for grid turbulence, when invoking Taylor's hypothesis. Specifically, one obtains with U 1 as mean velocity and X 1 = [(x 1 + r) + x 1 ]/2 and where the x 1 coordinate corresponds to the streamwise direction. This leads to the equations also considered by Danaila et al. [14], where the large scales are now determined by inhomogeneities in the x 1 direction.

IV. CLOSURES
As Eqs. (24) and (25) are unclosed, anything besides computing and comparing the individual terms from DNS requires additional closure assumptions. For this, we introduce two different approaches to close the system of equations in the following. First, we assume that the second order structure functions follow a power law, which allows us to directly integrate the two equations, finding explicit expressions for the third order structure functions in the inertial range. The drawback of this approach is of course that the scaling range of the second order structure functions is small at low Reynolds numbers and therefore the range for which the resulting third order expressions hold is also quite limited. However, interestingly enough, the same results derived in different ways by Lindborg [9] and Lundgren [11] are then recovered. Second, we close the equations by employing an eddy-viscosity ansatz as presented by Oberlack and Peters [23], which relates the second and third order structure functions. This allows us to also compare overall agreement and extrapolate the results to higher Reynolds numbers. We cannot rule out that the decay exponent has some influence on the parameters of the closures discussed in the following. However, we have varied the decay exponent while keeping all other parameters constant and have found no significant influence on the numerical solutions, as long as the decay exponent is not unrealistically large or small. Specifically, we compared n = 1.4, n = 1.45, and n = 1.5. For that reason, we are confident that the deviation from n = 10/7 is negligible.

A. Power law closure
Here, we assume that the normalized second order structure functions D 0,2 and D 2,0 follow a power law of the form in the inertial range, where both the prefactors C 2,0 and C 0,2 as well as the exponents ζ 2,0 and ζ 0,2 are assumed to be independent of the separation distance r. Power laws for the second order structure functions were introduced by Kolmogorov [2,24]. In his theory (K41 theory), he assumed 064403-6 that in the inertial range the structure functions are only dependent on the mean dissipation ε and the scale r: Since the inertial range is situated in between the small scales and the large scales, the solution should not depend on the viscosity ν, the dissipative length η, or the integral length L. From dimensional arguments, one then obtains power laws with ζ 2,0 = 2/3 and ζ 0,2 = 2/3 as exponents. However, the use of the mean dissipation as a scaling parameter has been questioned in an argument rooted in a remark of Landau (cf., e.g., the discussion in Ref. [25]), namely, that fluctuations of the dissipation may influence the scaling of the structure functions and consequently the scaling exponents may differ from 2/3. From experiments (see, e.g., [26,27]) and DNS (e.g., [28,29]), ζ 2,0 > 2/3 and ζ 0,2 > 2/3 have been found.
The numerical values of C 2,0 , C 0,2 , and ζ 2 for the data we use here are shown in Table II. Noticeably, there is a small trend for the prefactors to increase with the Reynolds number, although we cannot say whether they would approach a constant for very large Re λ . As the Kolmogorov constant has been found to vary only slightly (if at all) with increasing Reynolds number (cf. [30]), this might be due to the fact that μ = 0.

B. Eddy-viscosity closure
Another approach to close the coupled system is to directly relate the second and third order structure functions, for which there are different approaches in the literature. One way to close the equations is to employ an eddy-viscosity ansatz, as, e.g., discussed recently by Thiesset et al. [31]. Here, the formula of Oberlack and Peters [23] is used. They proposed an eddy-viscosity closure of the form with the eddy viscosities ν t,30 = κ 1 r D 2,0 , ν t,12 = κ 2 r D 0,2 .
Since both κ 1 and κ 2 as well as D 2,0 and D 0,2 are positive (and consequently also ν t,30 0 and ν t,12 0), the closure implies that ( u 1 ) 2 and ( u 2 ) 2 are transported towards smaller r, in agreement with the notion of the energy cascade towards smaller scales. Together with Eqs. (24) and (25) we then have a closed set of equations for D 2,0 and D 0,2 , with the Reynolds number Re L , the decay exponent n, and the coefficients κ 1 and κ 2 as parameters. The solution of D 3,0 and D 1,2 is then obtained by inserting the computed D 2,0 and D 0,2 into Eqs. (36) and (37). It is readily checked that this closure gives ζ 2 = 2/3 in the inertial range if the unsteady terms are neglected.
We need boundary conditions for r → 0 to solve the system of equations at hand. Specifically, four boundary conditions are needed as we have two second order equations. For homogeneous isotropic turbulence, Kolmogorov [2] showed that for r → 0 and therefore also for r → 0, which provides the four required boundary conditions. Consequently, D 3,0 ∼ r 3 , D 1,2 ∼ r 3 in the viscous range; i.e., the model reproduces the correct r scaling for r → 0. Finally, the model parameters κ 1 and κ 2 have to be specified. 064403-8 The values of κ 1 and κ 2 used here are shown in Table III. Noticeably, κ 1 and κ 2 do not vary much with the Reynolds number and no trend is observable, where it needs to be kept in mind that the considered range of Reynolds numbers is not that large. This observation is in line with the original formulation of Oberlack and Peters [23], who related κ 1 to the Kolmogorov constant C 2,0 . Here, κ 1 and κ 2 are rather directly computed from DNS via Eqs. (36) and (37).
Combining both equations, one obtains Again assuming power laws [Eq. (27)] in the inertial range and writing D 3,0 = −(4/5) r + 3,0 and D 1,2 = −(4/15) r + 1,2 , where 3,0 and 1,2 are corrections due to the unsteady and viscous terms, then if the corrections are small. Using the values of Table II with Eq. (41), one finds that κ 1 /κ 2 ≈ 4.6, which is close to the value κ 1 /κ 2 ≈ 4.5 from Table III. Consequently, κ 1 and κ 2 (and their ratio) weight the prefactors of the longitudinal and transverse structure functions. Noticeably, the continuity equation [Eq. (28)] constrains the ratio C 0,2 / C 2,0 . Then, as a function of the exponent ζ 2 only, which is thought to be independent of the Reynolds number (but not necessarily the flow configuration).

A. DNS
To quantify the influence of the unsteady terms, the balance of Eqs. (24) and (25) 24) and (25) are contributions by the unsteady terms, while the remaining term(s) are transport terms in r space. On the right-hand side (rhs), the first term (the −4/3) is the dissipative term, as we normalized the equation by ε, while the terms in square brackets are viscous terms.
For small r in the dissipative range, the dissipative term is balanced by the viscous terms, while the transport and the unsteady terms are negligible. Solving Eqs. (24) and (25) and neglecting these terms then leads to Eq. (38). The viscous terms are negligible for r outside the dissipative range. At intermediate r in the inertial range, the largest terms are the transport terms, which give the leading order solutions, Eqs. (4) and (3), after integration. However, the contribution of the unsteady terms is not negligible, as the transport terms alone are not sufficient to balance the dissipative term (the 4/3). For larger r, the transport terms become smaller and the unsteady terms larger. For very large scales r > 1 the unsteady terms are dominant and balance the dissipation. We also plot the sum 064403-9 of the unsteady, transport, and viscous terms as indicated by the dashed black line and find that it balances the 4/3 very well (i.e., the dashed lines nearly coincide with the 4/3). In other words, the assumption that the temporal changes ∂f/∂ t and ∂g/∂ t in Eq. (19) are negligible is well justified.
With increasing Reynolds number, the range of r for which the transport terms are larger than the viscous and the unsteady terms (i.e., the inertial range) increases (note that because of the normalization with large scale quantities, the inertial range is shifted to smaller r). Consequently, the scaling range of the transport terms for Re λ = 254.75 is larger than at Re λ = 121.39, but still very limited. Thus, the unsteady terms may not be neglected, as they contribute significantly at intermediate to large r.

Power law
Let us first look at the power law closure. The first term on the rhs of Eqs. (29) and (32) corresponds to Kolmogorov's asymptotic results in the inertial range for very large Reynolds numbers, while the second term on the rhs is the contribution due to the unsteady term and the third term on the  Table IV. rhs stems from the viscous terms. Indeed, the unsteady corrections (the second terms on the rhs) become smaller with increasing Reynolds number. However, they increase with increasing r/η. As A 3,0 and A 1,2 are positive, | D 3,0 | and | D 1,2 | then become smaller than 4/5 and 4/15 at a fixed Reynolds number; the deviations are not negligible for r/η larger than a certain threshold and the higher the Reynolds number, the higher the threshold value of r/η. This behavior is exactly the same as observed from our DNS data (cf. Figs. 2 and 3). Noticeably, the influence of μ, i.e., deviations from the K41 value ζ 2 = 2/3, play only a marginal role, as μ is small. As μ is found to be positive (see Table II for our DNS and, e.g., [5,29,32] in the literature) and small, 0 < μ 4/3, the viscous terms decrease much faster with increasing r and are therefore negligible as expected (r/η 1 in the inertial range). For K41 scaling, μ = 0 and there is no Reynolds number dependence of the viscous terms. Physically, μ = 0 corresponds to the statement that the second order structure functions are determined in the inertial range solely by the scale r (with dimension [m]) as well as another quantity with dimensions [m 2 /s 3 ] which is usually taken to equal the mean dissipation ε (cf. K41 theory). In this spirit, μ = 0 implies that there are more (albeit a priori unknown) quantities which influence the second order structure functions in the inertial range.
We compare Eqs. (32) and (29) Table IV. 064403-11 n determined from the DNS are given in Table IV. Because C 2,0 and C 0,2 vary with the Reynolds number, the coefficients A and B do so as well. The closure agrees better with the transverse D 1,2 than the longitudinal D 3,0 . This is probably due to the fact that D 1,2 feeds into D 3,0 , so that any errors of Eq. (29) are carried over to Eq. (32). Nevertheless, we find good qualitative agreement, also for the lower Reynolds number. As expected, the closure improves with increasing Reynolds number, because the scaling range of the second order structure functions increases. However, the deviations from Kolmogorov's results, Eqs. (3) and (4)

Eddy-viscosity closure
The results of the eddy-viscosity closure are shown in Figs. 6 and 7, where we compare the numerical solutions of D 3,0 and D 1,2 for the system of equations both with and without the unsteady terms to the DNS data. The model solution has been computed using an explicit Runge-Kutta solver. We use constant values of κ 1 and κ 2 throughout the numerical integration, which have been determined from DNS in the (presumed) inertial range. As the transport terms are small in the dissipative range, the deviations of the constants κ 1 and κ 2 from their true dissipative range values do not play a crucial role for the model performance, although they lead to small deviations in the dissipative range. The model parameters used here are listed in Table III,  between DNS and model increases with r. This is not that surprising, because the contributions of the unsteady terms to the balances as seen by Figs. 2 and 3 increase with increasing r. As discussed above, the model gives ζ 2 = 2/3 and ζ 3 = 1 [i.e., Eqs. (3) and (4)] if the unsteady terms are neglected. Consequently, their absence at the intermediate and large scales then results in an infinitely long inertial range for the model, even at finite Reynolds numbers. This is in agreement with the observation that the unsteady term is the only term in Eqs. (24) and (25) which explicitly contains the large scales.
After having established that the eddy-viscosity closure agrees very well with the DNS data when the unsteady corrections are included, we proceed to extrapolate towards higher Reynolds numbers. As κ 1 and κ 2 evaluated from our DNS (cf . Table III) Fig. 9. We find that the range for which the transport terms dominate increases with increasing Reynolds number, in agreement with Kolmogorov's classical notion of the inertial range and Fig. 8. The transport terms peak close to the intersection of the unsteady and viscous terms, which we briefly discuss in the following. We find that the intersection point ( r C , y C ), i.e., the crossover after which the unsteady terms are larger than the viscous terms, scales with the Reynolds number as indicated with the dashed black line in Figs. 9(a) and 9(b). We may thus write r C, = A r, Re where indicates the crossover of the terms of Eq. (24) and ⊥ the crossover of the terms of Eq. (25). Using a least squares fit, the model then gives the parameters shown in Table V and the longitudinal and transverse exponents are found to be approximately equal for both r C and y C . Noticeably, only the prefactors of the crossover length but not the corresponding values of the ordinate differ. That is, the inertial range of the longitudinal and transverse structure function is approached equally fast, but its location in r differs. In particular, the inertial range of the transverse structure function D 1,2 is shifted towards smaller r more than the corresponding inertial range of D 3,0 and we have r C, / r C,⊥ ∼ 1.55 and y C, / y C,⊥ ∼ 1.15 independent of the Reynolds number. We also find good agreement of the scaling as given by Table V with the lower Reynolds numbers of our DNS. As r = r/L, for both the longitudinal and transverse non-normalized crossover length where the prefactor is O (1). In other words, the transport terms of Eqs. (24) and (25) peak at the Taylor scale λ. That is, Kolmogorov's inertial range assumption that both the viscous and unsteady terms are small is best fulfilled at r on the order of the Taylor scale λ. This result is in agreement with the observation that λ is an intermediate length scale, smaller than the integral length L and larger than the dissipative scale η.

VI. CONCLUSION
After normalizing the unsteady terms in the second order equations of Hill with the large scale L in the spirit of [1], they can be written as a function of r = r/L only. Evaluating the balance of the second order using DNS, the unsteady terms contribute significantly to the inertial range solution of the third order structure functions D 3,0 and D 1,2 in agreement with previous results in the literature. Using a power law closure, we find that the contribution of the unsteady terms increases with increasing r/η, but decreases with increasing Reynolds number in agreement with the notion of an inertial range. If the second order structure functions follow K41 scaling, the same result as previously reported by Lundgren and Lindborg is recovered. Closing the system of equations by directly coupling the second and third order structure functions using an eddy-viscosity ansatz as proposed by Oberlack and Peters gives very good agreement with our DNS when the unsteady terms are included. This model also allows solving the equations for higher Reynolds number for which no DNS is available while retaining the influence of the unsteady and viscous terms. From the model, we find that the intersection of these two terms scales with the Taylor scale λ; i.e., λ is situated in the inertial range.