Antonia

M Meldi

L Djenidi

R A Antonia

Sensitivity analysis of the second and third-order velocity structure functions to the Reynolds number in decaying and forced isotropic turbulence using the EDQNM model

Keywords: HIT, EDQNM, turbulence modelling

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The study of homogeneous isotropic turbulence (HIT) using two-point statistics has proven to be quite fruitful, particularly as a result of the major contribution by Von Karman & Howarth [START_REF] Von Karman | On the statistical theory of isotropic turbulence[END_REF] who developed a transport equation for the two-point longitudinal velocity correlation f (r, t) defined as f (r, t) = u(x, t)u(x + r, t)/U 2 (x, t), where x is a spatial position, r is the magnitude of the longitudinal separation r and U 2 = 2 3 K(t) (K(t) is the turbulent kinetic energy); the over-bar represents an ensemble average. Here u(x, t) is the longitudinal velocity component of the velocity u i.e. u(x, t) = u(x, t) • r/r The transport equation for f (r, t), now referred to as the Karman-Howarth (hereafter denoted KH) equation and written below ∂f (r, t)U (1) where h(r, t) = u(x, t)u(x, t)u(x + r, t)/(U 2 (x, t)) 3/2 is the third order velocity correlation, marked a milestone for the theory of HIT. The term W f represents the correlation between a volume source term (forcing) and the velocity fluctuation in the Navier-Stokes equations. This term is zero for decaying HIT. Equally useful two-point statistical quantities for the theory of HIT are the so-called longitudinal velocity structure functions defined as S n (r, t) = [u(x + r, t) -u(x, t)] n = (δu) n .

(

) 2
where δu = u(x + r) -u(x) is the velocity increment. Of particular interest are the 2nd and 3rd order structure functions. The former is related to f (r) as follows

S 2 (r, t) = 2 3 K (1 -f (r, t)) , (3)
Its transport equation, which can be derived from the Navier-Stokes equations in similar fashion to equation [START_REF] Von Karman | On the statistical theory of isotropic turbulence[END_REF], can be expressed as

2 3 ∂K ∂t (t) = - 2 3 ε = 1 2 ∂S 2 ∂t + 1 6r 4 ∂r 4 S 3 ∂r - ν r 4 ∂ ∂r r 4 ∂S 2 ∂r + W S ( 4
)
where ε is the ensemble average of the instantaneous turbulent energy dissipation rate and W S is the spatial counterpart of W f . Assuming an infinitely large Reynolds number and considering that turbulence at small-scales is in a statistically steady state, allowed Kolmogorov [START_REF] Kolmogorov | Dissipation of energy in the locally isotropic turbulence[END_REF] to drop the first term on the right side of (4) in his K41 theory. After integration with respect to r, one obtains: S 3 (r, t) = -4 5 εr + 6ν ∂ ∂r S 2 (r, t) + Z.

After neglecting the second term on the right of (5) and Z, which represents the effects of the energy injected at large scales and, for convenience, will be denoted loosely as "turbulence production" hereafter, eq. ( 5) reduces to the 4/5 law [START_REF] Kolmogorov | Dissipation of energy in the locally isotropic turbulence[END_REF] S 3 (r, t) = -4 5 εr.

The dropping of ∂S 2 ∂t , i.e. the first term on the right side of (4), is critical for the establishment of [START_REF] Thiesset | Consequences of selfpreservation on the axis of a turbulent round jet[END_REF]. This term reflects the contribution of the large scales to all scales of motion and if reinstated in (5) yields [START_REF] Pg Saffman | Lectures on homogeneous turbulence[END_REF][START_REF] Danaila | A generalization of Yaglom's equation which accounts for the large-scale forcing in heated decaying turbulence[END_REF]

In decaying HIT, Z = 0 and ∂S 2 /∂t = 0, which reflects the statistical time evolution of HIT features such as the energy spectrum E(k, t) (k is the wave number space of the spectral transform). However, ∂S 2 /∂t = 0 can also be associated with the non-homogeneity either in time (e.g, 3D periodic box turbulence) or in space (e.g. grid turbulence). In direct numerical simulations (DNS) of a forced steady state 3D periodic turbulence, ∂S 2 /∂t = 0 and the term Z represents the imposed forcing at given (mainly large) scales.

Notice that Z and the third term on the right side of [START_REF] Sl Tang | Finite Reynolds number effect on the scaling range behaviour of turbulent longitudinal velocity structure functions[END_REF], which may also be considered as a some sort of forcing caused by the large scales, represent a cumulative effect suggesting that the impact of different forms of forcing on the transport of S 2 (r) at a given scale may be felt differently. This would be consistent with the results of Thiesset et al. [START_REF] Thiesset | Consequences of selfpreservation on the axis of a turbulent round jet[END_REF] and Antonia et al. [START_REF] Antonia | Boundedness of the velocity derivative skewness in various turbulent flows[END_REF], who

showed that, at finite Reynolds numbers, the balance between the longitudinal velocity derivative and the destruction coefficient of enstrophy in the transport equation for ε depends on the type of large-scale forcing and thus differs from flow to flow. Further, Tang et al. [START_REF] Sl Tang | Finite Reynolds number effect on the scaling range behaviour of turbulent longitudinal velocity structure functions[END_REF] showed that the impact of the forcing in various turbulent flows at small and moderate Re λ , the Taylor microscale Reynolds number, felt at scales of the order of λ, the Taylor microscale, depends on the types of flow under consideration. For example, they found that for a given Re λ this impact is largest on the centreline of a fully developed channel flow but smallest for stationary forced periodic box turbulence. For decaying-type flows, the strength of this impact lies between the previous two cases. These observations corroborate the results of Qian [START_REF] Qian | Slow decay of the finite Reynolds number effect of turbulence[END_REF], Antonia & Burattini [START_REF] Antonia | Approach to the 4/5 law in homogeneous isotropic turbulence[END_REF] and Tchoufag et al. [START_REF] Tchoufag | Spectral approach to finite Reynolds number effects on Kolmogorov's 4/5 law in isotropic turbulence[END_REF] who showed that the magnitude of the maximum of C 3 = -S 3 /(εr) is affected differently in forced and decaying turbulent flows; for a fixed Re λ , the maximum is larger for forced turbulence than decaying turbulence. Tchoufag et al. [START_REF] Tchoufag | Spectral approach to finite Reynolds number effects on Kolmogorov's 4/5 law in isotropic turbulence[END_REF] observed that C 3 ≈ 4/5 for Re λ ≥ 10 3 in forced turbulence and argued that Re λ should be larger than about 10 6 before the same maximum is observed in both types of turbulence. These observations and results illustrate the impact of the last term on the right side of [START_REF] Sl Tang | Finite Reynolds number effect on the scaling range behaviour of turbulent longitudinal velocity structure functions[END_REF]. They further indicate that this impact is not only flow dependent but also Reynolds number dependent for a given flow when the Reynolds number is finite.

As noted in [START_REF] Antonia | Finite Reynolds number effect and the 4/5 law[END_REF], Eq. 6 has yet to be verified convincingly on the basis of either experimental or DNS data. The expectation is that Eq. 6 is more likely to be verified when forcing is applied. Indeed, the maximum value of In the present work, the evolution with Re λ of several second and thirdorder statistical moments of HIT are investigated. The physical quantities include the energy spectrum E, the non-linear energy transfer T and the structure functions S 2 (r), S 3 (r). The evolution of the velocity increment skewness

S(r) = S 3 (r)/S 3/2 2
is also considered. The analysis is performed with calculations based on the eddy-damped quasi-normal model (EDQNM). In particular, the EDQNM-LMFA model is employed [START_REF] Bos | Lagrangian Markovianized Field Approximation for turbulence[END_REF]. The EDQNM model has been extensively used in the open literature to perform investigations of energy spectra, non-linear energy transfer and second and third order structure functions (see, among the others, references [START_REF] Cambon | Energy transfer in rotating turbulence[END_REF][START_REF] Bos | Reynolds number dependency of the scalar flux spectrum in isotropic turbulence with a uniform scalar gradient[END_REF][START_REF] Tchoufag | Spectral approach to finite Reynolds number effects on Kolmogorov Õs 4/5 law in isotropic turbulence[END_REF][START_REF] Meldi | On non-self-similar regimes in homogeneous isotropic turbulence decay[END_REF][START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically averaged descriptors[END_REF][START_REF] Meldi | The signature of initial production mechanisms in isotropic turbulence decay[END_REF][START_REF] Briard | Spectral modelling for passive scalar dynamics in homogeneous anisotropic turbulence[END_REF]). The paper is structured as follows. The EDQNM model is described in general terms in Section 2; the distribution of the forcing term F , representing production mechanisms in the spectral space, is also discussed in this Section.

Section 3 contains a comprehensive discussion of the results. Conclusions are given in Section 4. Specific features of the adopted EDQNM model are given in more detail in the Appendix A.

The EDQNM model

The EDQNM model [START_REF] Orszag | Analytical theories of turbulence[END_REF][START_REF] Lesieur | Turbulence in Fluids[END_REF][START_REF] Sagaut | Homogenous Turbulence Dynamics[END_REF] is briefly described in this section, while an extended discussion is provided in the Appendix A. The EDQNM is a turbulence closure in spectral space and it relies on the numerical discretization of the Lin equation:

∂E(k, t) ∂t + 2νk 2 E(k, t) = T (k, t) + F (k, t) (8)
where T (k, t) is the non-linear energy transfer due to triadic interactions and Several proposals for the forcing distribution were investigated for the case of statistically steady forced turbulence. However, forced HIT results are here restricted to two forcing schemes. These two schemes were chosen according to their different distribution (local or global) in the spectral space. In addition, both forcing terms comply well with the very large scale requirements imposed by the adaptive spectral mesh strategy employed in the calculations [START_REF] Meldi | An adaptive numerical method for solving EDQNM equations for the analysis of long-time decay of isotropic turbulence[END_REF] which will be described in the following. The first scheme is a single wavenumber forcing:

F
F (k) = γ F • ε(0) if k = k L (0) 0 if k = k L (0) (9)
where γ F is a positive scalar.

The second forcing scheme employed is derived in order to obtain conservation of the large-scale features i.e. ∂E(k, t)/∂t = 0 for k < k L , where k L (t) = L -1 (t) is the wave number associated with the integral length scale L. Artificial damping is applied at large wavenumbers viz.

F (k) = 2νk 2 E(k) -T (k) if k ≤ k L (0) (2νk 2 E(k) -T (k)) • r -k/k L (0) F if k > k L (0) (10)
The parameter r F = 1.1 has been set to obtain a bump between the integral length scale region and the scaling range. This feature has been observed in several experimental and numerical analyses reported in the literature [START_REF] Meldi | On the emergence of nonclassical decay regimes in multiscale/fractal generated isotropic turbulence[END_REF][START_REF] Laizet | The spatial origin of -5/3 spectra in grid-generated turbulence[END_REF][START_REF] Melina | Vortex shedding effects in grid-generated turbulence[END_REF].

Initial conditions are imposed using an energy spectrum functional form suggested by Pope [START_REF] Pope | Turbulent Flows[END_REF] and Meyers & Meneveau [START_REF] Meyers | A functional form for the energy spectrum parametrizing bottleneck and intermittency effects[END_REF]:

E I (k) = C K ε 2/3 k -5/3 f L (kL)f η (kη) (11)
with This parameter controls the shape of the energy spectrum at large scales, and the value chosen corresponds to the well known case of Saffman turbulence.

f L (kL) = kL [(kL) 1.5 + c L ] 1/1.5 5/3+σ , f η (kη) = exp(-β([(kη) 4 + c 4 η ] 1/4 -c η )) (12
Early tests choosing a value of σ = 4 (Batchelor turbulence) showed that the EDQNM prediction of S 2 and S 3 is not sensitive to the parameter σ. The initial Reynolds number has been set to Re λ = 10 6 for the free decay case and to Re λ ∈ [50, 10 6 ] for the forced cases. The initial transient regime in free HIT decay is governed by the features of the functional form prescribed for t = 0. An increase of the Reynolds number is observed up to Re λ ≈ 2 × 10 6 for t ≈ t 0 , where t 0 = K(0)/ε(0) is the initial turnover time. After this first phase, the statistics progressively lose memory of the initial condition and a classical power law decay is observed [START_REF] Comte-Bellot | The use of a contraction to improve the isotropy of grid-generated turbulence[END_REF][START_REF] Meldi | On non-self-similar regimes in homogeneous isotropic turbulence decay[END_REF]. For the present analysis, data are sampled in the range 10 6 ≥ Re λ ≥ 50. For the forced simulations, E(k) will converge towards a forced solution from the initially prescribed functional form E I (k). The statistically steady solution is clearly governed by the shape of the forcing term F . Finally, calculations are performed using an adaptive spectral mesh strategy [START_REF] Meldi | An adaptive numerical method for solving EDQNM equations for the analysis of long-time decay of isotropic turbulence[END_REF] which conserves the large-scale resolution for the free decay case. The condition k L (t)/k min (t) = 10 3 has been imposed, where k min is the smallest resolved mode. This choice allows confinement effects to be excluded.

The present study focuses on the prediction of the energy spectrum E, the non-linear energy transfer T and the structure functions S 2 and S 3 . These quantities are connected via the following integral relations [START_REF] Bos | Reynolds number effect on the velocity increment skewness in isotropic turbulence[END_REF][START_REF] Sagaut | Homogenous Turbulence Dynamics[END_REF]:

S 2 (r, t) = +∞ 0 4 E(k, t) 1 3 - sin(kr) -(kr) cos(kr) (kr) 3 dk (13) S 3 (r, t) = +∞ 0 12 T (k, t) 3(sin(kr) -(kr) cos(kr)) -(kr) 2 sin(kr) (kr) 5 dk (14)
using the EDQNM prediction for E(k, t) and T (k, t).

The spectra calculated using the EDQNM may exhibit remarkable differ- that, for similar Re λ , the free decay case exhibits a shorter scaling range.

The results for the forced case show an almost perfect plateau just outside the energy production region. The analysis of the compensated third-order velocity structure function -S 3 /(εr) in figure 1 (f) provides similar information. This behaviour is intimately connected with the triadic interactions associated with very energetic modes, as discussed by Meldi & Sagaut [START_REF] Meldi | Investigation of anomalous very fast decay regimes in homogeneous isotropic turbulence[END_REF].

Because of the connection between T and S 3 in equation ( 14), this mechanism is also responsible for the emergence of a clearer plateau at moderate

Re λ for the third-order structure function in the case of forced HIT. Also, results from every EDQNM calculation in figure 1 (e)-(f)collapse at the small scales, indicating that the large scale behaviour does not influence the dynamics of the small scales. In addition, all the computed spectra exhibit a power-law range proportional to k n E , n E ≈ -5/3 in the inertial region, once the Reynolds number is high enough. This is shown in figure 2 with the local calculation via polynomial fitting of the relation E ∝ k n E for the free decay case. For the highest Reynolds number considered, Re λ = 10 6 , the averaged deviation of the power law exponent within the scaling range from the K41 behaviour is n E + 5/3 = 5 × 10 -3 (the overbar here represents the averaged value of n E calculated in the range 10 -6 ≤ kη ≤ 10 -2 ), which is around 0.3%.

This value is sensibly smaller than the intermittency corrections usually reported in the literature for the energy spectrum, for which

E(k) ∝ k -5/3-µ/9 , µ ∈ [0.1, 0.3].
Once the distribution of F is prescribed, the production term Z(F ) = Z(r, t) in ( 7) can be calculated exactly via the integral relation [START_REF] Gotoh | Velocity field statistics in homogeneous steady turbulence obtained using a high-resolution direct numerical simulation[END_REF]:

Z(r, t) = 12r ∞ 0 1 15 + sin(kr) (kr) 3 + 3 cos(kr) (kr) 4 -3 sin(kr) (kr) 5 F (k) dk (15)
Using equations ( 13), ( 14) and ( 15) the budget terms in ( 7) can be analysed using the EDQNM data. This is shown in figure 3. Recall that Z = 0 in decaying HIT while ∂S 2 /∂t = 0 for the statistically steady forced HIT.

The values of S 3 calculated using [START_REF] Sl Tang | Finite Reynolds number effect on the scaling range behaviour of turbulent longitudinal velocity structure functions[END_REF], where each term on the right side of the equation are evaluated, and those using ( 14) are in good agreement; very small differences are observed for r/η ≈ 1 which are due to the discretization procedures used to calculate the derivatives of S 2 . As expected, the viscous term is dominant for r/η ≤ 5 while the large-scale contributions forced HIT (Ishihara et al. [START_REF] Ishihara | Study of high Reynolds number isotropic turbulence by Direct Numerical Simulation[END_REF]) and decaying free turbulence (Lamballais et al. [START_REF] Lamballais | Implicit/explicit spectral viscosity and large-scale SGS effects[END_REF]) are shown in figure 4. Overall, a good agreement is observed. The main differences are associated with the prediction of the bottleneck region of three-dimensional energy spectra. A larger peak in the bottleneck region is observed in the DNS distributions. This difference could be due to the EDQNM modelling as well as to the DNS discretization error due to a lack of small scale resolution. The difference between EDQNM and DNS results appears to be more marked for the data of Ishihara et al, where a mesh resolution k max η ≈ 1 is employed. On the other hand, the data of Lamballais et al. are obtained for k max η ≈ 7.5, which is much closer to the EDQNM resolution (k max η > 10). For the forced cases, a satisfactory agreement is observed for the energy flux in figure 4(b) for kη ≥ 6 × 10 -3 and figure 4(d) for kη ≥ 10 -2 . Similarly, the agreement for -S 3 /(εr) shown in figure A. 16(c) seems adequate in the intermediate scale range. It is perhaps surprising that EDQNM and DNS results show a good match for the quantities determined by the non-linear energy transfer T (k, t), which is the only modelled term in the EDQNM. Once T is obtained, the calculation of E in the EDQNM framework is virtually exact, as the energy spectrum is directly derived from the Lin equation. Thus, the differences observed for the three-dimensional energy spectrum in Figure 4(a)-(c) are most probably due to the forcing schemes employed, which are different for DNS and EDQNM. For the forced cases, one can also see that the variation of the interscale energy flux Π is less abrupt for the multiscale forcing proposed in Equation 10. The large scale behaviour in this case is more similar to the DNS data by Ishihara et al. [START_REF] Ishihara | Study of high Reynolds number isotropic turbulence by Direct Numerical Simulation[END_REF], which were also obtained using a large-scale forcing.

A quantitative measure of the discrepancy between a given spectrum E and the reference spectrum E with which E is compared is provided by the function

N (E , E ) = +∞ 0 (E -E ) 2 dk +∞ 0 (E ) 2 dk ( 16
)
For freely decaying HIT, the comparison of the EDQNM spectrum with results from Lamballais et al. [START_REF] Lamballais | Implicit/explicit spectral viscosity and large-scale SGS effects[END_REF] The EDQNM data are now compared with one-dimensional energy spectra measured in grid turbulence. The three-dimensional EDQNM spectra are manipulated to obtain the 1D spectra via the isotropic relation:

E 11 (k 1 ) = +∞ k 1 E(k) k 1 - k 2 1 k 2 dk (17)
The results for freely decaying HIT calculated via the EDQNM model are compared with the grid turbulence measurements of Mydlarski & Warhaft [START_REF] Mydlarski | On the onset of high-Reynolds-number grid-generated wind tunnel turbulence[END_REF] in figure 5 for Re λ = 448, 199. There is a reasonable agreement between the shapes of the spectra. In particular, both the EDQNM and experimental compensated spectrum

E 11 (k 1 )/(ε 1 2/3 k -5/3 1
) exhibit a maximum at about k 1 η 1 0.05. In this case, the measured discrepancy is to provide an estimation of the complete three-dimensional spectrum from measurements of E 11 in HIT:

N (E EDQN M
E(k) = 1 2 k 3 ∂ ∂k 1 k ∂ ∂k E 11 (18)
The results reported by Comte-Bellot & Corrsin [START_REF] Comte-Bellot | Simple eulerian time correlation of full and narrow-band velocity signals in grid-generated, isotropic turbulence[END_REF] are compared with the three-dimensional EDQNM spectrum in Figure 7 (b). One can see that in this case the comparison is not as good as for the one-dimensional spectra in Figure 7 (a). For 10 -2 ≤ k 1 η 1 ≤ 1, comparisons for compensated E 11 and E provide similar indications. However, significant differences arise outside of this range. Equation 18provides an over prediction of the density of turbulent kinetic energy at the large scales and an under prediction in the small scale region, which were not observed for E 11 . These results are mainly due to the approximation in Equation 18.

The previous comparison at moderate Reynolds numbers highlights the capability of the EDQNM model to provide satisfactory predictions for E(k, t)

and T (k, t), which are the essential elements for the calculation of second and third order velocity structure functions. The principal advantage of the EDQNM model is that it can be performed at very high Reynolds regimes, currently unachievable with either experiments or DNS. In both cases, one can see the emergence of three distinct regions with an increasing Re λ ; they are, in order of increasing r/η, the dissipative range, the scaling range and the large-scale range. The shape and evolution of S 2 /(4K/3) appear identical for both cases. To test if S 2 has a scaling range i.e. (εr) ζ 2 with ζ 2 = 2/3 we plot in Figures 8(c) and (d) S 2 /(εr) 2/3 versus r/η using a linear scale for the vertical axis. While both cases exhibit good collapse in the dissipative range (r/η < 100), the two cases differ substantially in the scaling range. In neither case can we identify a region where S 2 ∼ (εr) 2/3 , which should have been marked by a plateau. The forced HIT, which shows a stronger departure from a plateau than the decaying HIT, exhibits a well defined bump in S 2 /(εr) 2/3 as the large scale region is approached.

Interestingly, the magnitude of this bump for the two largest values of Re λ is the same, suggesting that it may have become independent of Re λ . In the scaling range, the local slope n S 2 = d log(S 2 )/d log(εr) continues to evolve with Re λ although the possibility that it may eventually reach a constant for larger values of Re λ than considered here cannot be excluded. The magnitude of n S 2 is shown in Figure 9. If a power law S 2 ∼ (εr) n S 2 exists, then n S 2 should be constant. The figure reveals that this is not strictly verified, so that the concept of a scaling range is at best only approximate. For Re λ = 10 One can observe a difference in the way the constancy of n S 2 is approached for forced HIT (Figure 9 (b)). In decaying HIT, n S 2 decreases monotonically from 2 to 0 with increasing r. In forced HIT, n S 2 first decreases, then, at least for Re λ larger than 471, increases slightly before decreasing to zero. This behaviour, at least up to the largest Reynolds numbers investigated, prevents n S 2 to reach the same values as in decaying HIT in the scaling range. Similar observations have been reported by Lohse & Muller-Groeling [START_REF] Lohse | Bottleneck effects in turbulence: scaling phenomena in r versus p space[END_REF]. Such a behaviour of n S 2 is associated with modified energy mode interactions which are responsible for an independent interscale energy flux from k in the scaling range, which leads to a steeper energy spectrum for k > k L , discernible in Figure 1 (a) and which results in an energy pile-up. The intensity of this energy mode interactions is governed by the energy production mechanisms, represented here by F (k). Further, the EDQNM results show that, for very large Reynolds number, while the energy pile up is restricted to the borders of the scaling range, it spreads over at least one decade in the spectral space.

This can explain the apparent deviation from the 5/3 law observed at moderate Re λ . It is worth noting that the distribution of n S 2 is the same for both decaying and forced HITs in the range 0 ≤ r/η ≤ 10, for all Re λ , reflecting the universality of the dissipative range.

Third-order structure function S 3

The evolution of S 3 with the Reynolds number is reported in figure 10.

The results are in agreement with previous results reported in the literature.

In particular, the present results confirm the approach towards an establishment of the 4/5 law as Re λ continues to increase. Figure 11, which reports the variation of C 3 = max(-S 3 /(εr)) with Re λ , confirms the results of Qian [8], Antonia & Burattini [START_REF] Antonia | Approach to the 4/5 law in homogeneous isotropic turbulence[END_REF] and Tchoufag et al. [START_REF] Tchoufag | Spectral approach to finite Reynolds number effects on Kolmogorov's 4/5 law in isotropic turbulence[END_REF]. While both cases exhibit a similar evolution, the value of Re λ at which C 3 ≈ 0. and Cambon et al. [START_REF] Cambon | Spectral modelling of homogeneous non-isotropic turbulence[END_REF]. As previously discussed, these spectra are affected by a significant noise and may exhibit as well a non-negligible departure from a freely decaying dynamic behaviour. Thus, one can expect that most flow configurations studied in realistic applications, where production and anisotropy effects cannot be completely eliminated, should fall in the area between the black line corresponding to the forced HIT prediction and the grey line representing the HIT free decay behaviour.

Finally, the variation of the slope n S 3 = d log(S 3 )/d log(εr), is shown in figure 12. For both decaying and forced HIT, a linear dependence of S 3 on r emerges in the scaling range as Re λ increases. There is a remarkable difference between the behaviour of S 3 and S 2 in the scaling range when forcing is applied. Indeed, S 2 exhibits a noticeable bump at the upper end of the scaling range; such a bump is absent in S 3 . While this difference has yet to be explained, one can speculate that the bump observed in S 2 is linked to the energy spectrum which also exhibits a peak at the lower wavenumbers. The product between the geometric function and E in the integrand of Eq. ( 13) exhibits very high values when r/η is in the wavenumber range corresponding to the forcing. This is consistent with the analysis by Lohse & Muller-Groeling [START_REF] Lohse | Bottleneck effects in turbulence: scaling phenomena in r versus p space[END_REF]. On the other hand, the product between the geometric function and T in the integrand of Eq. ( 14) shows an oscillatory tail

Skewness factor of the velocity increment, S(r)

We conclude this EDQNM analysis of S 2 and S 3 by presenting (Figure 13) the distributions of S(r) = S 3 /S . One expects that, with increasing Re λ , the EDQNM distributions of S 2 and S 3 will asymptotically approach a state which is consistent with the first similarity hypothesis of K41.

3. as Re λ increases, S 3 reaches a reasonably convincing power-law behaviour over the scaling range, especially when forcing is applied. This behaviour is fully consistent with the 4/5 law, Eq. 6, but it should be emphasized that, with forcing, this behaviour extends to values of r/L (see Figure A.15(b)) which approach 1 and is therefore difficult to reconcile with the classical concept of an inertial range. Over the same range, S 2 continues to evolve with Re λ . For forced HIT, S 2 exhibits a bump at the upper end of the scaling range where the forcing is introduced. As Re λ increases, and whilst the scaling range extends to larger separations, the bump does not vanish, but its impact on smaller separations decreases, paving the way towards a power-law behaviour.

The overriding conclusion that emerges from the present EDQNM results is that, although S 3 exhibits a behaviour consistent with Eq. 6 when Re λ is sufficiently large, S 2 continues to evolve with Re λ over the same range of scales. While the EDQNM is not able to capture intermittency, the present results suggest that observations of anomalous behaviors for n S 2 may not only be due to intermittency effects, but to the interactions (possibly nonlinear) between several different aspects which include intermittency and finite Reynolds number effects. If these observations are, in the future, confirmed by DNS for progressively higher Re λ , this would raise some doubt as to whether S 2 exhibits an anomalous scaling behaviour, as is usually reported in the literature. This topic has also been recently discussed under the perspective of Holder/Cauchy-Schwarz mathematical constraints [START_REF] Djenidi | Mathematical constraints on the scaling exponents in the inertial range of fluid turbulence[END_REF]. The present Re λ trend in the scaling range does not exclude the possibility that, irrespectively of whether the turbulence is decaying or forced, S 2 → (εr) 2/3 at even larger Reynolds numbers than considered here.

Also, it is difficult to discount the difference in the relative behavior of S 2 and S 3 between decaying and forced HIT. This difference is evident in experimental and DNS data and it is captured by the present EDQNM model.

The model shows that it persists when Re λ is several orders of magnitude larger than are currently achievable via experiment or DNS. An obvious inference is that conclusions, e.g. with regard to the magnitudes of n S 2

and n S 3 , drawn from forced turbulence studies do not necessarily apply to decaying turbulence and vice versa.

S 3

 3 (r, t) = -

  ) where C K ∈ [1.4, 1.6] is the Kolmogorov constant, L is the integral length scale and η is the Kolmogorov length scale. The free coefficients have been set to c η = 0.4, β = 5.3; c L has been chosen in order to obtain L(0) = 1.Calculations are performed by fixing an initial value (= 2) of the parameter σ.

Figure 1 :

 1 Figure 1: (a) Energy spectrum E and (b) forcing F for the case of statistically steady forced HIT using the forcing in Equation 10 for Re λ ≈ 10 5 . (c) Normalized energy spectrum E and (d) non-linear energy transfer T for the cases of decaying HIT and statistically steady forced HIT (equation 9) for Re λ ≈ 10 5 . (e) Second-order velocity structure function and (f) third-order velocity structure function for free decaying HIT and forced HIT.

Figure 2 : 4 r 0 s 4

 244 Figure 2: Local slope, n E (k), of the spectrum E(k) for decaying HIT.

3 . Results 3 . 1 .Figure 3 : 4 r 0 s 4 3 r 4 r 0 s 4

 331344344 Figure 3: Budget terms (normalized by ( r)) in Eq. 7 obtained from the EDQNM results for S 3 . Results are shown at Re λ = 10 6 for (a) decaying HIT and (b) statistically steady forced HIT using the single wavenumber model in equation 9. Continuous black lines represent the terms (a) -3 r 4

  provides a discrepancy of N (E EDQN M , E DN S ) = 0.1044 i.e. a global difference of ≈ 10%. The measure of the discrepancy for the EDQNM forced cases with the DNS data from Ishihara et al. is globally much higher because of the differences observed at the large scales due to the different forcing employed as well as for the tails observed in the DNS at the small scales. However, if these two regions are excluded, the measured discrepancy is N (E EDQN M , E DN S ) = 0.116 for Re λ = 1131 and N (E EDQN M , E DN S ) = 0.1246 for Re λ = 471. The discrepancy is slightly higher for the lower Re λ . This result is somehow expected considering that the EDQNM hypotheses should be better fulfilled for very high Re λ , when a clear scale separation is obtained.

11 , E Exp 11 )

 1111 = 0.092 for Re λ = 448 and N (E EDQN M 11 , E Exp 11 ) = 0.1338for Re λ = 199, which is comparable to the discrepancy measured between DNS and EDQNM data. An adequate comparison is also obtained with the one dimensional energy spectra by Bodenschatz et al.[START_REF] Bodenschatz | Variable density turbulence tunnel facility[END_REF] for Re λ =

Figure 4 :

 4 Figure 4: Comparison of EDQNM and DNS data for moderately high Re λ . The compensated (left column) energy spectrum and (right column) energy flux are shown. DNS data are taken from Ishihara et al. [37] ((a) to (d)) and Lamballais et al. [38] (e).

Figure 5 :

 5 Figure 5: Comparison of EDQNM and grid turbulence results for low Re λ . The compensated one-dimensional energy spectrum is shown. Grid turbulence data are taken from Mydlarski & Warhaft[START_REF] Mydlarski | On the onset of high-Reynolds-number grid-generated wind tunnel turbulence[END_REF]

Figure 6 :

 6 Figure 6: Comparison of EDQNM and grid turbulence results for moderate Re λ . The compensated one-dimensional energy spectrum is shown. Grid turbulence data are taken from Bodenschatz et al.[START_REF] Bodenschatz | Variable density turbulence tunnel facility[END_REF]

Figure 7 : 2 S 2

 722 Figure 7: Comparison of EDQNM and grid turbulence results for moderate Re λ . (a) The compensated one-dimensional energy spectrum and (b) the three-dimensional energy spectrum obtained via Eq. 18 are shown. Grid turbulence data are taken from Comte-Bellot & Corrsin[START_REF] Comte-Bellot | Simple eulerian time correlation of full and narrow-band velocity signals in grid-generated, isotropic turbulence[END_REF]

Figure 8 :

 8 Figure 8: Structure function S 2 . (left column) free HIT decay and (right column) forced HIT obtained with EDQNM.

Figure 9 :

 9 Figure 9: Local value of the slope, n S2 (r) ∼ d log(S 2 )/d log(εr) for (a) decaying HIT and (b) forced HIT.

  behaviour at larger values (almost by a factor of 3) of r/η, which are represented in Figure3(b). This can explain the difference in the scaling range behaviour between S 2 and S 3 .

Figure 10 :

 10 Figure 10: Structure function S 3 for (a) decaying and (b) forced HIT obtained with EDQNM.

Figure 11 :

 11 Figure 11: Maximum value C 3 for the normalized longitudinal structure function S 3 . Lines: present EDQNM; white symbols: reference data for HIT free decay; black symbols: reference data for forced HIT.

Figure 12 :

 12 Figure 12: Local value of the slope, n S3 (r) ∼ d log(S 3 )/d log(εr). (a): decaying HIT ;(b) forced HIT.

3 / 2 2

 32 Figure A.14 in the appendix. Also both exhibit the same sensitivity to Re λ . Secondly, S(r) shows a tendency to form a plateau in the scaling range when Re λ increases. One can expect that S(r) may become constant in the scaling range as Re λ increases to infinitely large values. This behaviour seems to conform with K41.

Figure 13 :

 13 Figure 13: Skewness factor S(r) of δu. EDQNM results for (left column) decaying HIT and (right column) forced HIT.

Figure A. 14 :Figure A. 15 :Figure A. 16 :

 141516 Figure A.14: Comparison of results obtained using the classical version of the EDQNM and the recent EDQNM-LMFA proposal for freely decaying HIT. (a) Velocity derivative flatness S behaviour in the range Re λ ∈ [200, 5×10 5 ]. (b) Compensated structure function S 2 /(ε r) 2/3 for Re λ = 10 5 .

  6 in both cases, n S 2 approaches a value of about 0.674 in the scaling range; 1% larger than 2/3. This value is very close to that (0.679) obtained by McComb et al.[START_REF] Mccomb | Spectral analysis of structure functions and their scaling exponents in forced isotropic turbulence[END_REF] in forced HIT but smaller than that (0.72) found by Iyer et al.[START_REF] Iyer | Reynolds number scaling of velocity increments in isotropic turbulence[END_REF][START_REF] Biferale | Self-similar subgrid-scale models for inertial range turbulence and accurate measurements of intermittency[END_REF][START_REF] Iyer | Scaling exponents saturate in three-dimensional isotropic turbulence[END_REF] also in forced HIT. This discrepancy reflects the different approaches employed to estimate the power law exponent. In Iyer et al.[START_REF] Iyer | Reynolds number scaling of velocity increments in isotropic turbulence[END_REF][START_REF] Biferale | Self-similar subgrid-scale models for inertial range turbulence and accurate measurements of intermittency[END_REF][START_REF] Iyer | Scaling exponents saturate in three-dimensional isotropic turbulence[END_REF] n S 2 is obtained via a global measure over almost two decades for r/L. The present values and those of McComb et al.[START_REF] Mccomb | Spectral analysis of structure functions and their scaling exponents in forced isotropic turbulence[END_REF] are instead estimated locally.

  , in particular for lower Reynolds numbers (Re λ < 10 2 ). This difference is due to the different strategy used to obtain C 3 . Briard et al. performed a large number of simulations with initial spectra prescribed for several Re λ .For each calculation, they sampled C 3 after 20 turn-over times, which may not be long enough to eliminate any lingering effects associated with the initial conditions. This would support the idea that a departure from a pure free decay HIT would result in a larger C 3 yielding a larger rate at which the 4/5 law is reached with increasing Re λ . It is also interesting to observe that the data by Briard et al. match pretty well C 3 values obtained from the reconstructed three-dimensional spectra from Comte-Bellot & Corrsin[START_REF] Comte-Bellot | Simple eulerian time correlation of full and narrow-band velocity signals in grid-generated, isotropic turbulence[END_REF]

	8 is first attained
	is, as expected, much lower for forced than for decaying HIT. The present
	results, which have been obtained by a single model run, match almost per-
	fectly the fit proposed by Antonia & Burattini [9]. On the other hand, the
	EDQNM results of Briard et al. (as reported in [23]) indicate higher values
	of C 3

Acknowledgements: Dr. C. Cambon and Dr. A. Briard are acknowledged for fruitful discussion.

Competing interests: We declare we have no competing interests.

Appendix A. Numerical details about the EDQNM model

The EDQNM model [START_REF] Orszag | Analytical theories of turbulence[END_REF][START_REF] Lesieur | Turbulence in Fluids[END_REF][START_REF] Sagaut | Homogenous Turbulence Dynamics[END_REF] is a turbulence closure in spectral space.

The non-linear energy transfer T in equation 8 is obtained via two important approximations:

• The fourth order cumulants, which represent the deviation of the velocity derivative pdf from a Gaussian distribution, are approximated using an eddy-damping term (eddy damping hypothesis)

• Time integration is simplified via a Markovianization procedure. Markovian processes have been extensively used for the analysis of a large number of applications for turbulent flows (see the review article by Meyer and Saggini [START_REF] Meyer | Testing the markov hypothesis in fluid flows[END_REF] and references therein). Its application within the EDQNM formalism implies that local non-linear times are much shorter than the large-eddy turnover time, which is strictly tied with scale separation observed for very high Reynolds numbers A closed expression for the non-linear energy transfer T (k, t) is obtained via shell integration of radius k:

where Θ kpq is a characteristic time resulting from the Markovian approximation and g is a geometric function. In the present article, a recent proposal by Bos & Bertoglio [START_REF] Bos | Lagrangian Markovianized Field Approximation for turbulence[END_REF] is used to calculate the term Θ kpq in equation A.1.

In the classical EDQNM formulation [START_REF] Orszag | Analytical theories of turbulence[END_REF][START_REF] Lesieur | Turbulence in Fluids[END_REF] the term is calculated as:

where η E is the EDQNM damping factor. This last term is usually calculated following the proposal by Pouquet et al. [START_REF] Pouquet | Evolution of high reynolds number two-dimensional turbulence[END_REF]:

The free coefficient A ∈ [0.3, 0.5] is chosen to optimize the value for

. One of the major drawbacks of such a procedure is that it is not possible to derive an optimal value for which every HIT statistical quantity exhibits a behaviour in agreement with experiments and DNS in the literature. For this reason, the damping factor η E is here calculated resolving an evolution equation for the velocity-displacement cross-correlation spectrum F CC [START_REF] Bos | A single-time two-point closure based on fluid particle displacements[END_REF]:

The terms T F i are calculated via integration in the spectral space, similarly to T EDQN M in equation A.1:

where the terms f i are geometrical functions and the spectral time scale Θ F kpq is defined as:

The algebraic relation between η E and F CC is:

The term νk 2 in (A.14) was included in the earlier work by Bos & Bertoglio [START_REF] Bos | A single-time two-point closure based on fluid particle displacements[END_REF] but was then excluded in their most recent formulation [START_REF] Bos | Lagrangian Markovianized Field Approximation for turbulence[END_REF].

In the present work, it is employed after initial tests provided a more accurate prediction of numerous physical quantities. The term η X E in equation A.13 is instead set to zero. Comparison between the classical EDQNM and the present proposal (EDQNM-LMFA) [START_REF] Bos | Lagrangian Markovianized Field Approximation for turbulence[END_REF] show that results for free decay regimes are in qualitative agreement, as shown in figure A.14 (a) for the velocity skewness S. The high-Reynolds number asymptotic limit for the velocity derivative skewness S is in the second case S ≈ -0.57, which is approximately 7% smaller than the result (S ≈ -0.53) obtained with the classical EDQNM version [START_REF] Meldi | Further insights into self-similarity and selfpreservation in freely decaying isotropic turbulence[END_REF] or with that observed in the majority of experimental data [START_REF] Antonia | Boundedness of the velocity derivative skewness in various turbulent flows[END_REF]. The results obtained with the EDQNM-LMFA model appear to be slightly more sensitive to FRN effects. This larger sensitivity is associated with a slow convergence towards the expected behaviour [START_REF] Gagne | Reynolds dependence of third-order velocity structure functions[END_REF][START_REF] Ishihara | Study of high Reynolds number isotropic turbulence by Direct Numerical Simulation[END_REF][START_REF] Mccomb | Spectral analysis of structure functions and their scaling exponents in forced isotropic turbulence[END_REF][START_REF] Lamballais | Implicit/explicit spectral viscosity and large-scale SGS effects[END_REF]. The agreement is generally adequate. 700