About Universality and Thermodynamics of Turbulence

This paper investigates the universality of the Eulerian velocity structure functions using velocity fields obtained from the stereoscopic particle image velocimetry (SPIV) technique in experiments and direct numerical simulations (DNS) of the Navier-Stokes equations. It shows that the numerical and experimental velocity structure functions up to order 9 follow a log-universality (Castaing et al. Phys. D Nonlinear Phenom. 1993); this leads to a collapse on a universal curve, when units including a logarithmic dependence on the Reynolds number are used. This paper then investigates the meaning and consequences of such log-universality, and shows that it is connected with the properties of a “multifractal free energy”, based on an analogy between multifractal and thermodynamics. It shows that in such a framework, the existence of a fluctuating dissipation scale is associated with a phase transition describing the relaminarisation of rough velocity fields with different Hölder exponents. Such a phase transition has been already observed using the Lagrangian velocity structure functions, but was so far believed to be out of reach for the Eulerian data.


Introduction
A well-known feature of any turbulent flow is the Kolmogorov-Richardson cascade by which energy is transferred from large to small length scales until the Kolmogorov length scale below which it is removed by viscous dissipation. This energy cascade is a non-linear and an out-of-equilibrium universal process. Moreover, the corresponding non-dimensional energy spectrum E(k)/ 2/3 η 5/3 is an universal function of kη, where η = (ν 3 / ) 1/4 is the Kolmogorov length scale, the mean energy dissipation rate per unit mass, and ν the kinematic viscosity. Every used quantity is identified with its definition in a nomenclature available in Table 1. However, there seems to be little dependences on the Reynolds number, boundary, isotropy or homogeneity conditions [1]. In facts, the energy spectrum is based upon a quantity, the velocity correlation that is quadratic in velocity. Nevertheless, it is now well admitted that the universality does not carry over for statistical quantities that involve higher order moments. For example, the velocity structure functions of order p, given by S p ( ) = u(x + r) − u(x) p x, r = are not universal, at least when expressed in units of the Komogorov scale η and velocity u K = (ν ) 1/4 (see below, Section 3.2 for an illustration).

Symbol Mathematical Definition Interpretation
u(x, t) ∈ R 3 × R → R 3 Velocity field k ∈ R + Wavenumber E(k) FT u i (x + r, t)u i (x, t) x, r = ,t Energy spectrum k f ∈ R * + Forcing wavenumber N x ∈ N Grid size in direction x ν ∈ R * + Kinematic viscosity ∈ R * + Mean dissipation power per unit mass Taylor Reynolds number ∆x ∈ R * + SPIV spatial resolution p ∈ [1,9] Power ∈ R * + Scale L ∈ R * + Inertial large scale δ u(x, t) u(x + r, t) − u(x, t) r = Velocity increment at scale Φ(x) exp(− x 2 /2)/(2π) 3 2 Wavelet filter In theory (δW(x, , t)) p x,t For data analysis Velocity structure functioñ Rescaled length General intermittency correction τ(p, ) τ(p, θ( )) General intermittency correction γ(Re),β(Re) R + → R Fitting functions G R 2 → R General function from Castaing [2] Universal parameters Parabolic fit log( /L) for in Inertial range Intermittency correction from general rescaling Free energy The mechanism behind this universality breakage is identified in [3], where a generalization of the Kolmogorov theory is introduced, based on the hypothesis that a turbulent flow is multifractal. In this model, the velocity field is locally characterized by a Hölder exponent h, such that δ u(x) ≡ u(x + r) − u(x) r = ∼ h(x) ; here h is a stochastic function that follows a large deviation property [4] P (log(|δ u|/u 0 ) = h log ( /L 0 )) ∼ ( /L 0 ) C(h) , where u 0 (resp. L 0 ) is the characteristic integral velocity (resp. length), and C(h) is the multifractal spectrum. Velocity fields with h < 1 are rough in the limit → 0. Indeed they are at least not differentiable. In real flows, any rough field with h > −1 can be regularized at sufficiently small scale (the "viscous scale") by viscosity. The first computation of such dissipative scale was performed by Paladin and Vulpiani [5], who showed that it scales with viscosity like η h ∝ ν 1/(1+h) , thereby generalizing the Kolmogorov scale, which corresponds to h = 1/3. Such a dissipative scale fluctuates in space and time (along with h), resulting in non-universality for high order moments, at least when expressed in units of η and u K .
A few years later, Frisch and Vergassola [6] claimed that the universality of the energy spectrum can be recovered, if the fluctuations of the dissipative length scale are taken into account by introducing a new non-dimensionalisation procedure. The new prediction was that log E(k) − 2 should be a universal function of log(kη)/ log(Re), where Re is the Reynolds number. This claim was examined by Gagne et al., later using data from the Modane wind tunnel experiments [7]. They further suggested that the prediction can be extended to the velocity structure functions S p , so that log(S p ( )/u p K )/ log(Re) should be a universal function of log( /η)/ log(Re), at any given p. They found good agreement for p up to 6. The velocity measurements, in the above experiments, were performed using hot wire anemometry, which provide access to only one component of velocity. To our knowledge, no further attempts have been made to check the claim with more detailed measurements.
The purpose of the present paper is to reexamine this claim. However, now using the velocity fields obtained from the Stereoscopic Particle Image Velocimetry (SPIV) in experiments and the direct numerical simulations (DNS) of the Navier-Stokes equations (NSE). We show that the numerical and experimental velocity structure functions up to order 9 follow a log-universality [7]; they indeed collapse on a universal curve, if we use units that include log(Re) dependence. We then investigate the meaning and consequences of such a log-universality, and show that it is connected with the properties of a "multifractal free energy", based on an analogy between multifractal and thermodynamics (see [8] for summary). This framework uses co-existing velocity fields with different Hölder exponents which are regularized at variable scales. We show that in such a framework, this fluctuating dissipation length scale is associated with a phase transition describing the relaminarisation of velocity fields.

Experimental Facilities and Parameters
We use experimental velocity field described in [9]. The radial, axial and azimuthal velocity are measured in a Von Kármán flow, using Stereoscopic Particle Image Velocimetry technique at different resolutions ∆x. The Von Kármán flow is generated in a cylindrical tank of radius R = 10 cm through counter-rotation of two independent impellers with curved blades. The flow was maintained in a turbulent state at high Reynolds number by two independent impellers, rotating at various frequencies. Figure 1 shows the sketch of the experimental setup. The five experiments are performed in conditions so that the non-dimensional mean energy dissipation per unit mass is constant. The viscosity is monitored using mixture of water and glycerol, so as to vary the Kolmogorov length η. Table 2 summarizes the different parameters; R λ = λu rms /ν is the Reynolds number based on the Taylor length scale λ = u 2 ∇u 2 , the root mean squared velocity u rms and the kinematic viscosity ν. All velocity measurements are performed in a vertical plane that contains the rotation axis. The case (A) corresponds to measurements over the whole plane contained in between the two impellers, and extending from one side to the other side of the cylinder. Its resolution is 5 to 10 times coarser than similar measurements performed by zooming on a region centered around the symmetry point of the experiment (on the rotation axis, half way in between the two impellers), over a square window of size 4 cm × 3 cm. Since the flow is not homogeneous, statistics in this central region may differ from statistics computed over the whole tank. This explains the strong difference of R λ between (A) and (B,C). The little differences between (B) and (C) are explained by the different experimental resolutions used.

Direct Numerical Simulation
The direct numerical simulations (DNS), based on pseudo-spectral methods, are performed in order to compare with our experimental data. The DNS runs with R λ = 25, R λ = 80, R λ = 90 and R λ = 138 are performed using the NSE solver VIKSHOBHA [10], whereas the run with R λ = 53 is carried out using another independent pseudo-spectral NSE solver. The velocity field u is computed on a 2π triply-periodic box.
Turbulent flow in a statistically steady state is obtained by using the Taylor-Green type external forcing in the NSE at wavenumber k f = 1 and amplitude f 0 = 0.12, the value of viscosity is varied in order to obtain different values of R λ (see Ref. [10] for more details).

Velocity Increments vs. Wavelet Transform (WT) of Velocity Gradients
The classical theories of Kolmogorov [11,12] are based on the scaling properties of the velocity increment, defined as δ u(x, t) = u(x + r, t) − u(x, t) r = where = r is the distance over which the increment is taken. As pointed out by [8], a more natural tool to characterize the local scaling properties of the velocity field is the wavelet transform of the tensor ∂ j u i , defined as: In what follows, we choose a Gaussian function Φ(x) = exp(− x 2 /2)/(2π) 3 2 such that Φ(r)dr = 1. We then compute the wavelet velocity increments as This formulation is especially well suited for the analysis of the experimental velocity field, as it naturally allows to average out the noise. It has been verified that the wavelet-based approach yields the same values for the scaling exponents as those computed from the velocity increments [10].

K41 and K62 Universality
In the first theory of Kolmogorov [11], the turbulence properties depend only on two parameters: the mean energy dissipation per unit mass and the viscosity ν . The only velocity and length unit that one can build using these quantities are the Kolmogorov length η = (ν 3 / ) 1/4 and velocity u K = ( ν) 1/4 . The structure functions are then self-similar in the inertial range η L 0 , where L 0 is the integral scale, and follow the universal scaling: which can also be recast into:S where C p is a (non universal) constant. This scaling is typical of a global scale symmetry solution, and was criticized by Landau, who considered it incompatible with observed large fluctuations of the local energy dissipation. Kolmogorov then built a second theory (K62), in which fluctuations of energy dissipation were assumed to follow a log-normal statistics, and taken into account via an intermittency exponent κ and a new length scale L, thereby breaking the global scale invariance. The resulting velocity structure functions then follow the new scaling: which implies a new kind of universality involving the relative structure functionsS p as: where τ(p) = κ p(3 − p) and A p is a constant. Such a formulation already predicts an interesting universality, if L = L 0 , as we should have: Therefore, we should be able to collapse all structure functions, at different Reynolds number by plotting ( L 0 η ) τ(p)S p as a function of η , given that L 0 /η ∼ Re 3/4 . There is however no clear prediction about the value of L and we show in the data analysis (Section 4) that L differs from L 0 .
The relation (7) shows that log is a linear function of log( η ). In principle, such universal scaling is not valid outside the inertial range, i.e., for example when < η. To be more general than previously thought, it can however be shown using the multifractal formalism as first shown by [6].

Multifractal and Fluctuating Dissipation Length
For the multifractal (MFR) model, it is assumed that the turbulence is locally self-similar, so that there exists a scalar field h(x, , t), such that for a range of scales in a suitable "inertial range" η h L, where L is a large inertial scale, η h a cut-off length scale, and u 0 a characteristic large-scale velocity. This scale η h is a generalization of the Kolmogorov scale, and is defined as the scale where the local Reynolds number |δ u|/ν is equal to 1. Writing δ u = u 0 ( /L) h leads to the expression of η h as a function of the global Reynolds number Re = u 0 L/ν as η h ∼ LRe −1/(1+h) . This scale thus appears as a fluctuating cut-off which depends on the scaling exponent and therefore on x. This is the generalization of the Kolmogorov scale η ∼ LRe −3/4 ≡ η 1 3 , and was first proposed in [5]. Below η h , the velocity field becomes laminar, and δ u ∝ . When the velocity field is turbulent, h ≡ log(δ u/u 0 )/ log( /L) varies stochastically as a function of space and time. Also, if the turbulence is statistically homogeneous, stationary and isotropic, h only depends on , the scale magnitude. Therefore, formally, h can be regarded as a continuous stochastic process labeled by log( /L). By Kramer's theorem [13], one sees that as in the limit → 0, log(L/ ) → ∞, we have: where C(h) is the rate function of h, also called multifractal spectrum. Formally, C(h) can be interpreted as the co-dimension of the set where the local Hölder exponent at scale is equal to h. Using Gärtner-Elis theorem [13], one can connect C and the velocity structure functions as: To proceed further and make connection with previous section, we set = u 3 0 /L so that S p ( ) can now be written: This shows that τ(p) is the Legendre transform of the rate function C(h + 1/3), i.e., τ(p) = min h (p(h − 1/3) + C(h)), and equivalently, that C(h) is the Legendre transform of τ(p). Because of this, it is necessarily convex. The set of points where C(h) ≤ 3, represents the set of admissible or observable h, is therefore necessarily an interval, bounded by −1 ≤ h min and h max ≤ 1.
As noted by [6], the scaling exponent ζ(p) = p/3 + τ(p) defined via Equation (11) is only constant in a range of scale where > η h for any h ∈ [h min , h max ]. For small enough , this condition is not met anymore, since as soon as < η h , all velocity fields corresponding to h are "regularized", and do not contribute anymore to intermittency since they scale like . This results in a slow dependence of ζ(p) with respect to the scale, which is obtained via the corrected formula: To understand the nature of the correction, we can compute the value of h such that = η h . This gives: h( ) = −1 + log(Re)/ log(L/η h ). With θ = log(L/ )/ log(Re), Equation (12) can be rewritten as:S where . As discussed by [6], this implies a new form of universality that extends beyond the inertial range, into the so-called extended dissipative range, as; If the scale L is constant and equal to L 0 , the integral scale, then we have Re = (L 0 /η) 4/3 and the multifractal universality implies that log(S p )/ log(L 0 /η) is a function of log( /η)/ log(L 0 /η). When the function is linear, we thus recover the K62 universality. The multifractal universality is thus a generalization of the K62 universality.
This form of universality is however not easy to test, as the scale L is not known a priori, and may still depend on Re. In what follows, we demonstrate a new form of universality that allows more freedom upon L and encompass both K62 and multifractal universality.

General Universality
Using the hypothesis that turbulence maximizes some energy transfer in the scale space, Castaing [2] suggested a new form of universality for the structure functions, that reads: where A p and K 0 are universal constants and β and G are general functions, G being linear in the inertial range, G(p, x) ∼ τ(p)x. The validity of this universal scaling was checked by Gagne and Castaing [7] on data obtained from the velocity fields measured in a jet using hot wire anemometry. They found good collapse of the structure functions at different Taylor Reynolds R λ , provided γ(Re) is constant at low Reynolds numbers and follows a law of the type: γ(Re) ∼ γ 0 / log(R λ /R * ), where R * is a constant, whenever R λ > 400. Since we have R λ ∼ Re 1/2 and (L 0 /η) ∼ Re 3/4 , we can rewrite Equation (15) as: where S 0p are some constants and β and H are general functions. Compared to the K62 or MFR universality Formulas (7) or (14), we see that Formula (16) is a generalization of these two universality with L = L 0 . It allows however more flexibility than K62 or MFR universality through the function β(Re), which is a new fitting function. We test these predictions in Section 4 and provide a physical interpretation of (16) in Section 5.

Check of Universality Using Data Analysis
The various universality are tested using the velocity structure functions based on the wavelet velocity increments Equation (2), in order to minimize the noise in the experimental data. We define: We then apply this formula to both experimental data ( Table 2) and numerical data (Table 3), to get wavelet velocity structure functions at various scales and Reynolds numbers. Table 3. Parameters for the DNS. R λ is the Reynolds based on the Taylor micro-scale. η is the Kolmogorov length. The third column gives resolution of the simulation through k max η, where k max = N x /3 is the maximum wavenumber. The fourth column gives the grid size; notice that the dimensionless length of the box is 2π. Here, min is the smallest scale available for the calculations of the wavelets. k f is the forcing wavenumber. The Sample columns gives the number of points (frames × grid size) over which the statistics are computed.
This is shown in Figure 2 for both experimental and numerical data. Obviously, the data do not collapse on a universal curve, meaning that K41 universality does not hold. This is well known, and is connected to intermittency effects [14].

Check of K62 Universality
The K62 universality (7) can be checked by plotting: The collapse depends directly on τ(p), the intermittency exponents. Obtaining the best collapse of all curves is in fact a way to fit the best scaling exponents τ(p). We thus implement a minimization algorithm that provides the values of τ(p) that minimized the distance between the curve and the line of slope τ(p). The values of τ(p) are reported in Table 4. The best collapse is shown on Figure 3a for the DNS, and Figure 3b for the experiment. The collapse is better for experiments than for the DNS. However, in both cases, there are significant differences in between points at different R λ , at larger scales, showing that universality is not yet reached.   [9]. The exponents τ EXP (p)(red square) and τ DNS (blue circle) have been computed through a least square algorithm.

Check of General Universality
We can now check the most general universality, by plotting: In this case, best collapse is obtained by fitting two families of parameters: S 0p , β(Re) that are obtained through a procedure of minimization. We take the DNS at R λ = 138 as the reference case, and find for both DNS and experiments, the values of β(Re) and S 0p that best collapse the curves. The corresponding collapses are provided in Figure 5. The collapses are good for any value of Re, except for the DNS at the lowest Reynolds number, which does not collapse in the far dissipative range.   Table 4), while the black squares correspond to ζ EXP (p)/ζ EXP (3).

Function β(Re)
Motivated by earlier findings by [7], we plot in Figure 6 the value 1/β as a function of R λ . Our results are compatible with 1/β ∼ β 0 / log(R λ ), with β 0 ∼ 4/3 over the whole range of Reynolds number. For comparison, we provide also on Figure 6 the values found by Gagne and Castaing [7] in jet of liquid Helium, shifted by an arbitrary factor to make our values coincide with them at large Reynolds number. This shift is motivated by the fact that β(Re) is determined up to a constant, depending upon the amplitude of the structure functions used as reference. At large Reynolds, our values are compatible with theirs. At low Reynolds, however, we do not observe the saturation of 1/β that is observed in the jet experiment of [7]. An interpretation of the meaning of β(Re) is provided in Section 5.

Scaling Exponents
Our Collapse method enables us to obtain the scaling exponents of the structure functions ζ(p) by the following two methods: (i) Using the K62 universality, we get τ(p), and then ζ(p) = ζ(3)p/3 + τ(p). These estimates still depend on the value of ζ(3), which is not provided by the K62 universality plot. To obtain them, we use a minimization procedure on both experimental log(S 3 /u 3 K ) from the one hand, and the numerical log(S 3 /u 3 K ) on the other hand (see Figure 4a), to compute ζ(3) as the value that minimizes the distance between the curve and a straight line of slope ζ (3). The values so obtained are reported in Table 4, and are used to compute ζ(p) from τ(p). In Table 4, two different methods are used to process the experimental data. The subscript SAW refers to the values obtained by [9] on the same set of experimental data, using velocity increments and Extended Self-Similarity technique [15]. The quantities with subscript EXP are computed through a least square algorithm upon τ(p), minimizing the scatter of the rescaled structure functions log L 0 η τ(p)S p with respect to the line ( /η) τ(p) . DNS data have been processed the same way as EXP.
(ii) Using the general universality, we may also get τ p,univ by a linear regression on the collapse curve. Please note that since the data are collapsed, this provides a very good estimates of this quantity, with the lowest possible noise. In practice, we observe no significant differences with the two estimates; therefore, we only report the values obtained by following the first method.
The corresponding values are plotted in Figure 4 and summarized in Table 4. Please note that for both DNS and experiments, the value of ζ(3) is different from 1, which is apparently incompatible with the famous Kolmogorov 4/5th law that predicts ζ(3) = 1. This is because we use absolute values of wavelet increments, while the Kolmogorov 4/5th law uses signed values. We have checked that using unsigned values, we obtain a scaling that is closer to 1, but with larger noise. Note also that when we consider the relative value ζ(p)/ζ (3), we obtain values that are close to the values obtained [9] on the same set of experimental data, using velocity increments and Extended Self-Similarity technique [15].

Multifractal Spectrum
From the values of ζ(p), one can get the multifractal spectrum C(h) by performing the inverse Legendre transform: Practically, this allows to use the following formula: To estimate C, we thus first perform a polynomial interpolation of order 4 on ζ(p), then derivate the polynomial to estimate d ζ(p) d p , thus get C through Equation (22). The result is provided in Figure 6b for both the DNS and the experiment.
The curve looks like a portion of parabola, corresponding to a log-normal statistics for the wavelet velocity increments. Specifically, fitting by the shape: we get a = 0.35 and b = 0.045. This parabola also provides a good fit of the scaling exponents, as shown in Figure 4 by performing Legendre transform of C(h) given by Equation (23).

Thermodynamical Analogy
Multifractal obeys a well-known thermodynamical analogy [8,16,17] that will be useful to interpret and extend the general universality unraveled in the previous section. Indeed, considering the quantity: By definition µ is positive definite and µ = 1 for any . It therefore can be interprated as a scale dependent measure. It then also follows a large-deviation property as: where S(E) is the large deviation function of log(µ ) and has the meaning of an energy while log( /η) has the meaning of a volume, and log(µ )/ log( /η) is an energy density. With the definition of µ , it is easy to see that S is connected to C, the large deviation function of |δW |. In fact, since in the inertial range where |δW | 3 ∼ ζ(3) , we have S(E) = C(3h − ζ (3)). By definition, we also have: so thatS 3p is the partition function associated with the variable log(µ ), at the pseudo-inverse temperature p = 1/k B T. Taking the logarithm of the partition functionS 3p , we then get the free energy F as: By the Gärtner-Elis theorem, F is the Legendre transform of S: F = min E (pE − S(E)). The free energy a priori depends on the temperature T = 1/k B p, on the volume V = log( /η) and on the number of degrees of freedom of the system N. If we identify N = (1/β(Re)) log(L 0 /η), we see that the general universality means: i.e., can be interpreted as extensivity of the free energy. The thermodynamic analogy is thus meaningful and is summarized in Table 5. It can be used to derive interesting prospects. Table 5. Summary of the analogy between the multifractal formalism of turbulence and thermodynamics.

Thermodynamics Turbulence
Temperature Free energy F log(S 3p )

Multifractal Pressure and Phase Transition
Given our free energy, F = log(S 3p ), we can also compute the quantity conjugate to the volume, i.e., the multifractal pressure as: P = ∂F/∂V. In the inertial range, whereS p ∼ τ(p) , we thus get P = τ(p), which only depends on the temperature. Outside the inertial range, P has the meaning of a local scaling exponent that also depends upon the scale, i.e., on the volume V and on N (Reynolds number). Using our universal functions derived in Figure 5, we can then compute empirically the multifractal pressure P and see how it varies as a function of T, V and N. It is provided in Figure 7 for R λ = 25 and R λ = 53, and in Figure 8 for R λ = 90 and R λ = 138. We see that at low Reynolds number, the pressure decreases monotonically from the dissipative range, reaches a lowest points and then increases towards the largest scale. There is no clear flat plateau that would correspond to an "inertial" range. In contrast, at higher Reynolds number, a plateau appears for p = 1 to p = 4 when going towards the largest scale, the value of the plateau corresponding to τ DNS . The plateau transforms into an inflection point for p ≥ 5 making the derivative ∂P/∂V change sign. This is reminiscent of a phase transition occurring in the inertial range, with coexistence of two phases: one "laminar" and one "turbulent". We interpret such a phase transition as the result of the coexistence of region of flows with different Hölder exponents, with areas where the flow has been regularized due to the action of viscosity, because of the random character of the dissipative scale (see below).  Table 4.

Conclusions
We show that a deep analogy exists between multifractal and classical thermodynamics. In this framework, one can derive from the usual velocity structure function an effective free energy that respects the classical extensivity properties, provided one uses several degrees of freedom (given by N = 1/β(Re)) that scales like log(R λ ). This number is much smaller than the classical N ∼ Re 9/4 that is associated with the number of nodes needed to discretize the Navier-Stokes equation down to the Kolmogorov scale. It would be interesting to see whether this number is also associated with the dimension of a suitable "attractor of turbulence". Using the analogy, we also find the "multifractal" equation of state of turbulence, by computing the multifractal "free energy" F and "pressure" P = ∂F/∂V. We find that for large enough R λ and p (the temperature), the system obeys a phase transition, with coexistence of phase like in the vapor-liquid transition. We interpret this phase transition as the result of the coexistence of region of flows with different Hölder exponents, with areas where the flow is relaminarized due to the action of viscosity, because of the random character of the dissipative scale. We note that this kind of phenomenon has already been observed in the context of Lagrangian velocity increments, using the local scaling exponent ζ(p, ∆t) = d(log(S p (∆t)))/d(log(∆t)) [18]. The phase transition is then associated with the existence of a fluctuating dissipative time scale. It is further shown that in a multifractal without fluctuating dissipative time scale, the local exponent decreases monotonically from dissipative scale to large scale, implying a disappearance of the phase transition [19]. Funding: H.F. is supported by a CFR. F.N. is supported by a fellowship from the ENS. This work has been supported by the Labex PALM (project Interdist) and by the ANR EXPLOIT, grant agreement no. ANR-16-CE06-0006-01. Part of this work was granted access to the resources of IDRIS under the allocation 2A310096 made by GENCI (Grand Equipement National de Calcul Intensif).

Conflicts of Interest:
The authors declare no conflict of interest.