Dissipation range of the energy spectrum in high Reynolds number turbulence

We seek to understand the kinetic energy spectrum in the dissipation range of fully developed turbulence. The data are obtained by direct numerical simulations (DNS) of forced Navier-Stokes equations in a periodic domain, for Taylor-scale Reynolds numbers up to $R_\lambda=650$, with excellent small-scale resolution of $k_{max}\eta \approx 6$ for all cases (and additionally at $R_\lambda=1300$ with $k_{max}\eta\approx3$), where $k_{max}$ is the maximum resolved wavenumber and $\eta$ is the Kolmogorov length scale. We find that, for a limited range of wavenumbers $k$ past the bottleneck in the range $0.1 \lesssim k\eta \lesssim0.5$, the spectra for all $R_\lambda$ display a universal stretched exponential behavior of the form $\exp(-k^{2/3})$, in rough accordance with recent theoretical predictions. The stretched exponential fit in the near dissipation range $1 \lesssim k\eta \lesssim 4$ does not possess a unique exponent, which decreases with increasing $R_\lambda$. This region can be regarded as a crossover between the stretched exponential behavior and the far dissipation range $k\eta>6$, in which analytical arguments as well as DNS data with superfine resolution (S. Khurshid et al., Phys.~Rev.~Fluids 3, 082601, 2018) suggest a $\exp(-k)$ dependence. We remark on the connection to the multifractal model which hypothesizes a pseudo-algebraic law.

We seek to understand the kinetic energy spectrum in the dissipation range of fully developed turbulence. The data are obtained by direct numerical simulations (DNS) of forced Navier-Stokes equations in a periodic domain, for Taylor-scale Reynolds numbers up to R λ = 650, with excellent small-scale resolution of kmaxη ≈ 6 for all cases (and additionally at R λ = 1300 with kmaxη ≈ 3), where kmax is the maximum resolved wavenumber and η is the Kolmogorov length scale. We find that, for a limited range of wavenumbers k past the bottleneck in the range 0.1 kη 0.5, the spectra for all R λ display a universal stretched exponential behavior of the form exp(−k 2/3 ), in rough accordance with recent theoretical predictions. The stretched exponential fit in the near dissipation range 1 kη 4 does not possess a unique exponent, which decreases with increasing R λ . This region can be regarded as a crossover between the stretched exponential behavior and the far dissipation range kη > 6, in which analytical arguments as well as DNS data with superfine resolution (S. Khurshid et al., Phys. Rev. Fluids 3, 082601, 2018) suggest a exp(−k) dependence. We remark on the connection to the multifractal model which hypothesizes a pseudo-algebraic law.
Introduction: Turbulent fluctuations in fluid flows span a wide range of scales and are often characterized by the energy spectrum E(k), where k is the wavenumber whose inverse measures the scale size in real space [1,2]. The integral of E(k) over all k gives the average kinetic energy of turbulence. The pioneering work of Kolmogorov [3], K41 henceforth, theorized that the small scales are universal at sufficiently high Reynolds numbers, depending solely on the viscosity ν and the mean dissipation rate ǫ . In addition, at an intermediate range of scales, the so-called inertial range, the dependence on ν vanishes as well. These considerations imply that in the range of scales much smaller than the energy injection scale, the energy spectrum can be written as E(k) ∼ ǫ 2/3 k −5/3 f (kη), where η = (ν 3 / ǫ ) 1/4 is the Kolmogorov length scale and f is some universal function of kη, tending to a constant in the inertial range. The energy spectrum has been extensively studied by numerous researchers, and the k −5/3 prediction (with some small intermittency correction) seems to have received substantial validation [4][5][6][7]. However, the functional form of f and its universality in the dissipation range, are still not properly understood.
Many attempts have been made over the years to characterize f using both experiments and direct numerical simulations (DNS) [8][9][10][11][12][13][14][15][16], all of which suggest the following general form (1) However, experiments have seldom resolved the range beyond kη ≈ 1 [4,8], and DNS have been either restricted to low R λ [9][10][11]13] or achieved high R λ by sacrificing small-scale resolution [12]. Consequently, there has been no clarity regarding the values of the coefficients in Eq. (1), especially the exponent γ. The direct interaction approximation [17] and other ideas [8,[18][19][20], predict a pure exponential, i.e., γ = 1 for large wavenumbers or very small scales regularized by viscosity. While this prediction was found to hold at low R λ [10,11,13], it could not adequately describe data at higher R λ and often led to conflicting and ad hoc fits [4,8,21,22]. The above issues were addressed in a recent study [15] by means of a DNS with superfine resolution. This study showed that there are two distinct regimes in the dissipation range: a far-dissipation range (FDR) for kη > 6 consistent with a pure exponential; and a near-dissipation range (NDR) in the vicinity of kη 1, where the spectrum is a pure exponential at very low R λ (γ = 1), but evolves into a stretched exponential with decreasing γ < 1 as R λ increases. This analysis in [15] was restricted to R λ 100, which invites the question as to whether some asymptotic high-R λ limit for the NDR (and hence γ) exists.
Our goal is to assess the full picture by means of a well-resolved DNS of isotropic turbulence based on highly accurate Fourier pseudo-spectral methods, going up to grids of 12288 3 and Taylor-scale Reynolds number R λ ranging from 140 to 1300. The largest R λ here is more than an order of magnitude larger than in [15]. We shall also interpret the findings in terms of two recent theoretical predictions; the first, resulting from ideas based on distributed chaos, predicts γ = 3/4 or 2/3 depending on a particular choice of parameters [23]; and the second, emerging from a nonperturbative renormalization group (NPRG) approach, predicts that γ = 2/3 [14]. Both references have claimed an agreement in their respective inspections with experimental or DNS data [14,16,23,24] but, as mentioned earlier, the data were restricted to either low R λ or limited resolution. We assess these claims and show that there exists an intermediate bridging region between the stretched exponential and the FDR, on which we shall remark only briefly.
DNS data: The data, summarized in Table I, are an extension of those utilized in a recent work [25]; we have also extended the runs at R λ = 390 and 650 for longer computational times. In addition, we have performed a new run at R λ = 1300, with a small-scale resolution kη = 3. The totality of the data allows us to demonstrate that the behavior of the spectrum in the dissipation range, while being consistent with [15], is more complex at higher R λ than was anticipated there. The stretched exponential region of NDR: In Fig. 1a we show the compensated energy spectra for various R λ , as a function of kη. Consistent with earlier results at lower R λ [15], a systematic enhancement in the high-wavenumber spectral density is observed with respect to R λ . The curves clearly show that f (kη) is not a universal function of its argument. Earlier studies such as [4,12,26], which inferred a spectral collapse consistent with K41 phenomenology, were limited by technical reasons: the scatter in [26] was sufficiently large that possible trends could have been easily obscured; and the spectral resolution in yet others was limited to kη < 1. In order to explore the behavior further, we consider its log-derivative of Eq. (1), given by This form allows us to isolate the stretched exponential behavior in a meaningful way. Figure 1b shows φ(kη) for various R λ . The curves clearly suggest that the f (kη) is non-universal and exhibits concave curvatures, confirming that γ < 1. In contrast to the results of [16,24], Fig. 1b at higher R λ shows that the energy spectra in NDR cannot be described by one single value of γ. We now undertake a more detailed analysis to extract γ and its dependence on R λ . We also make a preliminary note that the multifractal formalism should yield a nearly constant form for φ(kη) [27], quite unlike the data (details are discussed later). As noted in [15] and other similar contexts [25], extracting γ through a direct curve fit of Eq. (2) results in a complex nonlinear regression, which is strongly dependent on initial seeds and does not guarantee proper convergence. Hence, alternative strategies must be utilized. We adopt a modified version of the strategy utilized in [15]. In order to evaluate γ, the authors of [15] compensated φ(k) by (kη) γ for different γ values, until a reasonable plateau was observed in the chosen fitting range. Furthermore, they noted that the precise value of α was inconsequential for the fit, and one can set it to zero without any loss of fidelity. Consequently, Eq. (2) reduces to −φ(kη) ∼ βγ(kη) γ , and one can obtain γ by simply fitting a power law for −φ(kη) in the desired range. This procedure is similar to that of [15], but as we will see, it has the added benefit of also identifying the appropriate ranges of power law behaviors. Additional methods of extracting γ are also possible but are not stable as described in the Appendix, and can lead to incorrect conclusions. Figure 2a shows log-log plots of −φ(kη) versus kη and confirms that the log-derivative exhibits two regions of distinct power laws. (An expanded version is provided in Fig. 2b.) In the first, corresponding to the region immediately past the bottleneck (known to occur around kη ≈ 0.1 [28,29]) to kη 0.5, data for all R λ exhibit a spectral collapse, with the exponent ranging from 0.68 ± 0.03 for R λ = 140 to 0.67 ± 0.01 for R λ = 1300, effectively 2/3. This value of γ is in agreement with the theoretical prediction from NPRG [14] (though a precise wavenumber range is not obtainable from the theory). However, the analysis also predicts a strong R λ −3 dependence of the coefficient β in the this range. Given the collapse obtained in Fig. 2, it follows that β is independent of R λ in this range-which invites a possible refinement of the underlying theoretical arguments in [14].
Another prediction in [23] on the basis of distributed chaos yields γ = 3/4, which seems to be ruled out. However, the same author provided an alternative argument that yields γ = 2/3, which is consistent with the present results. Incidentally, some support for γ = 2/3 was also provided in a recent work [24] in the range R λ = 60 − 240, though the fitting range included part of the wavenumber range (0.2 kη 4) that lies outside this range of universal fit-and thus produced the considerable error bar. Our results show that the prediction from NPRG is valid only in a small region of NDR and the behavior in the remainder of NDR (kη > 1) is quite different, as shown next.
The remainder of NDR: From Fig. 2, the second region where power laws can be fitted is the range 1 < kη < 3 − 4, say, which is similar to that utilized in [15,24]. It is clear that no single value of γ is inadequate consistent with ref. [15] for up to R λ = 100. We have plotted the current data in Fig. 3a together with the data for R λ ≤ 100 from [15]. Evidently γ continues to decrease for the R λ range considered here, with a plausible fit that is logarithmic. Alternatively, Fig. 3b shows an equally plausible weak power law with γ ∼ R λ −0.16 , for R λ > 20, say. Both these dependencies are similar to how the bottleneck flattens [28]. In fact, it seems reasonable that the increase in spectral density (with R λ ) in the dissipation range is connected to a decrease in the bottleneck region. Physically, the bottleneck is thought to develop due to inadequate 'thermalization' of the energy transferred from inertial to dissipation scales, leading to a pileup at their crossover [30]. However, with increasing R λ and the scale-range, the energy transfer across the scales is better facilitated, leading to the diminution of the bottleneck and simultaneous rise in spectral density in NDR. Recent experimental results [29] have confirmed the decay of the bottleneck with R λ but suggest that a very mild bottleneck persists even at R λ ≈ 4000. Assuming this behavior to be correct, and if there is a connection between the diminished bottleneck and decreased exponent, we may infer that the exponent in this part of NDR will likely continue to decrease at least up to R λ = 4000; however, if the trend in Fig. 3b persists for higher R λ , it is clear that the asymptotic value will be zero, albeit achieved at extremely high Reynolds numbers.
A log-dependence of the exponent would be nominally consistent with the multifractal (MF) prediction of [27]. However, this formalism explicitly predicts a power law dependence of the spectra in the NDR, which is not consistent with the current data (see Figure 1b). Also, as noted in [15], the MF prediction does not appear to work for the spectra in the FDR (kη ≫ 1), to which we will return momentarily. In addition, recent work of [25] has suggested that the assumptions built in to the MF model for describing this range of scales is not entirely justified, or may require extremely large R λ (that are practically impossible to simulate, even without the fine resolution used here). Nevertheless, one has to leave open the question of whether the trend observed here for γ holds up for much higher R λ .
We note that β in Eq. (2) is also a parameter of the stretched exponential fit. However, given how the NDR beyond kη > 1 departs systematically from the universal regime with γ = 2/3, it follows that the product βγ, the coefficient that appears in Eq. (2), will emerge as independent of R λ (since these power law fits can be thought to have a common origin but different slopes). This implies that β ∼ 1/γ, which is, in fact, consistent with the observation in [15]. However, the precise value of the product βγ is strongly dependent on the exact fitting range and hence not very important.
The far dissipation region: So far as we know, only the authors of [15] were able to adequately resolve the range kη > 6 (FDR). They could do it because of the comparatively low R λ of their simulations. Their conclusion is that the spectral shape in FDR is exponential, consistent with analytical arguments [8,[17][18][19][20] that require viscosity to regularize the velocity field at very small scales. Thus, it appears reasonable to expect that the spectrum in the FDR would be a pure exponential. It has not been possible for us to have the same resolution and also extend to the R λ values attained in this paper. Thus, we leave open the possibility that a pure exponential occurs for wavenumbers higher than kη > 6 even at very high R λ . It is not lost on us that the increasing demands on resolution could be hinting something important at the analytic structure of the Navier-Stokes equations.
Concluding remarks: We have analyzed the dissipation range behavior of the energy spectra obtained from very well resolved DNS of isotropic turbulence at Taylor-scale Reynolds numbers that are an order of magnitude higher than in earlier studies. In the process, we have extended the work of [15] and also undertaken the verification of various theoretical predictions. Our results indicate that the behavior of the spectra in NDR is more complex than previously realized. In a limited range of NDR, 0.15 ≤ kη ≤ 0.5, our results show a universal stretched exponential fit to the spectra, of the form exp(−k 2/3 ). This result matches the theoretical prediction from NPRG, but the anticipated range of validity is much smaller than that asserted in recent works [14,24]. It is also consistent with one version of the distributed chaos [23]. In the FDR, one can anticipate a pure exponential predicted from analyticity arguments [17]. The range of wavenumbers 1 < kη < 6, in which the spectra possess a stretched exponential behavior with the exponent persistently decreasing with increasing R λ , may be regarded as a crossover between the stretched exponential and exponential behaviors.

ACKNOWLEDGMENTS
We thank Alain Pumir, Diego Donzis, P. K. Yeung and Sualeh Khurshid for their comments on the draft and sustained collaboration over the years. We gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gausscentre.eu) for providing computing time on the GCS supercomputers JUQUEEN and JUWELS at Jülich Supercomputing Centre (JSC), where the simulations were performed. This work was also partly supported by supercomputing resources under the Blue Waters sustained petascale computing project, which was supported by the NSF (awards OCI-5725070 and ACI-1238993) and the State of Illinois.
Appendix A: Robust determination of exponent for a stretched exponential curve fit In determining the exponent γ in Eq. (2) from experimental or numerical data, a few different methods can be employed (other than determining the power-law exponent as done in the current work). One method, also utilized in recent works [14,16,23,24], is to plot the log-derivative φ(kη) vs. (kη) γ for a choice of γ and thereafter compare the curve with a straight line. However, we note that this method relies heavily on a visual comparison, rather than an explicit curve fit, and is inherently error prone [31]. For instance, in Fig. 4 we simply plot the log-derivative of various exponential functions f (x) as a function of x 2/3 . As is evident, all curves can be erroneously matched with a straight line on the basis of a simple visual inspection, leading to the incorrect conclusion such as the exponent being 2/3 in a larger range. Another method is to directly determine the log-derivative of −φ(kη) with respect to kη, which in principle can allow for a direct evaluation of γ from the resulting local slopes plot, without any explicit curve fit [24]. Unfortunately, this method is not very reliable in practice, since calculating the local slope of −φ(kη) amounts to taking two successive log-derivatives of the spectrum E(kη). This leads to substantial numerical noise, making it very difficult to meaningfully extract the exponent. Finally, we note that in the method employed in [15] φ(kη) is compensated by (kη) γ until a reasonable plateau is obtained. While this indeed results in a reasonable fit, it requires an advance knowledge of the fitting range (which perhaps is the reason why the 2/3 region was overlooked in their work).