L Sciacovelli

P Cinnella
email: paola.cinnella@ensam.eu

C Content

F Grasso

Dense gas effects in inviscid homogeneous isotropic turbulence

Keywords: compressible turbulence, isotropic turbulence, gas dynamics

A detailed numerical study of the influence of dense gas effects on the large-scale dynamics of decaying homogeneous isotropic turbulence is carried out by using the van der Waals gas model. More specifically, we focus on dense gases of the Bethe-Zel'dovich-Thompson type, which may exhibit non-classical nonlinearities in the transonic and supersonic flow regimes, under suitable thermodynamic conditions. The simulations are based on the inviscid conservation equations, solved by means of a ninth-order numerical scheme. The simulations rely on the numerical viscosity of the scheme to dissipate energy at the finest scales, while leaving the larger scales mostly unaffected. The results are systematically compared with those obtained for a perfect gas. Dense gas effects are found to have a significant influence on the time evolution of the average and root mean square (r.m.s.) of the thermodynamic properties for flows characterized by sufficiently high initial turbulent Mach numbers (above 0.5), whereas the influence on kinematic properties, such as the kinetic energy and the vorticity, are smaller. However, the flow dilatational behaviour is very different, due to the non-classical variation of the speed of sound in flow regions where the dense gas is characterized by a value of the fundamental derivative of the gas dynamics (a measure of the variation of the speed of sound in isentropic compressions) smaller than one or even negative. The most significant differences between the perfect and the dense gas case are found for the repartition of dilatation levels in the flow field. For the perfect gas, strong compressions occupy a much larger volume fraction than expansion regions, leading to probability distributions of the velocity divergence highly skewed toward negative values. For the dense gas, the volume fractions occupied by strong expansion and compression regions are much more balanced; moreover, strong expansion regions are characterized by sheet-like structures, unlike the perfect gas which exhibits tubular structures. In strong compression regions, where compression shocklets may occur, both the dense and the perfect gas exhibit sheet-like structures. This suggests the possibility that expansion eddy shocklets may appear in the dense gas. This hypothesis is also supported by the fact that, in dense gas, vorticity is created with equal probability in strong compression and expansion regions, whereas for a perfect gas, vorticity is more likely to be created in the strong compression ones.

Introduction

The present work aims at investigating the large-scale dynamics of compressible homogeneous isotropic turbulence in the presence of strong dense gas effects. Dense gases (DG) are usually defined as single-phase fluids, characterized by complex molecules, working in pressure and temperature conditions of the same order of magnitude as their thermodynamic critical point (see, e.g. [START_REF] Cinnella | Inviscid and viscous aerodynamics of dense gases[END_REF] and references cited therein). Dense gas fluid dynamics is governed by a thermodynamic parameter known as the fundamental derivative of gas dynamics [START_REF] Thompson | A fundamental derivative in gasdynamics[END_REF], defined as

Γ := v 3 2c 2 ∂ 2 p ∂v 2 s = 1 + ρ c ∂c ∂ρ s , (1.1)
where ρ is the density, v the specific volume, p the pressure, s the entropy and c = √ ∂p/∂ρ| s the sound speed. The first equality in the preceding definition shows that Γ is related to the curvature of isentropic lines in the p-v plane; according to the second definition, Γ is a measure of the rate of change of the sound speed in isentropic transformations. From (1.1), Γ < 1 implies ∂c/∂ρ| s < 0, meaning that the sound speed grows in isentropic expansions and drops in isentropic compressions, unlike the case of perfect gases, where:

Γ = γ + 1 2 , (1.2)
γ = c p /c v being the ratio of the isobaric specific heat c p to the isochoric specific heat c v . For thermodynamic stability, γ is always greater than unity, hence Γ > 1 for perfect gases. On the contrary, heavy gases characterized by high c v /R ratios (R being the gas constant), exhibit extended ranges of density and pressures in which Γ is less than unity, or even negative, the perfect gas behaviour being recovered in the low-density limit.

Dense gas flows exhibit interesting physical behaviours that are not yet fully understood. The most interesting phenomena are expected for a family of heavy polyatomic compounds, named Bethe-Zel'dovich-Thompson (BZT) fluids [START_REF] Bethe | The theory of shock waves for an arbitrary equation of state[END_REF][START_REF] Zel'dovich | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF][START_REF] Thompson | A fundamental derivative in gasdynamics[END_REF], which exhibit a region of negative Γ values (called the 'inversion zone', see [START_REF] Cramer | On the propagation of waves exhibiting both positive and negative nonlinearity[END_REF] in the vapour phase close to the liquid/vapour coexistence curve. This is theoretically predicted to result in non-classical compressibility effects in the transonic and supersonic flow regimes, such as rarefaction shock waves, mixed shock/fan waves and shock splitting (see e.g. [START_REF] Cramer | On the propagation of waves exhibiting both positive and negative nonlinearity[END_REF]Cramer 1989b, 1991[START_REF] Rusak | Transonic flow of dense gases around an airfoil with a parabolic nose[END_REF]. Furthermore, the entropy change across a weak shock, given by [START_REF] Bethe | The theory of shock waves for an arbitrary equation of state[END_REF]:

s = - c 2 Γ v 3 ( v) 3 6T + O(( v) 4 ) (1.3)
with T the absolute temperature, is much weaker than usual for dense gases with Γ 1, leading to reduced shock losses. Note that, if Γ < 0, then s > 0 only if v > 0, meaning that expansion shocks are admissible according to the second principle of thermodynamics.

For dense gases, the perfect gas (PFG) model is no longer valid, and more complex equations of state (EoS) are used to account for their peculiar thermodynamic behaviour. One of the simplest models for dense gases is the Van der Waals (VDW) EoS (Van der Waals 1873), which takes into account the co-volume occupied by the molecules and the two-body collision interactions. The VDW EoS allows us to predict the inflexion of the critical isotherm at the critical point (p/p c = 1, v/v c = 1; subscript (•) c indicates the critical conditions) and the change of curvature of the isotherms in the Clapeyron plane pv for Γ = 0. As a consequence, VDW is also the simplest possible EoS allowing us to account for the BZT effects [START_REF] Bethe | The theory of shock waves for an arbitrary equation of state[END_REF].

Turbulent flows of dense gases are of interest for a wide range of engineering applications, such as energy conversion cycles [START_REF] Monaco | Supersonic flows of dense gases in cascade configurations[END_REF][START_REF] Brown | Application of Bethe-Zel'dovich-Thompson fluids in organic Rankine cycle engines[END_REF][START_REF] Horen | Numerical simulation of real-gas flow in a supersonic turbine nozzle ring[END_REF], high Reynolds number wind tunnels [START_REF] Anderson | A numerical study on the use of sulfur hexafluoride as a test gas for wind tunnels[END_REF], chemical transport and processing [START_REF] Kirillov | Analysis of modern natural gas liquefaction technologies[END_REF]) and refrigeration [START_REF] Zamfirescu | Performance investigation of high-temperature heat pumps with various BZT working fluids[END_REF].

Considerable progress has been made in the past 30 years about the study of inviscid dense gas flows. Inviscid flows of dense gases have been extensively studied analytically, with focus on the generation of non-classical compressibility effects such as expansion shocks, sonic and double-sonic shocks and shock splitting [START_REF] Thompson | Negative shock waves[END_REF][START_REF] Cramer | On the propagation of waves exhibiting both positive and negative nonlinearity[END_REF][START_REF] Cramer | Exact solutions for sonic shocks in van der Waals gases[END_REF]Cramer 1989a,b;[START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF][START_REF] Cramer | Nonclassical dynamics of classical gases[END_REF][START_REF] Cramer | Transonic flows of Bethe-Zel'dovich-Thompson fluids[END_REF][START_REF] Cramer | Nozzle flows of dense gases[END_REF][START_REF] Rusak | Transonic flow of dense gases around an airfoil with a parabolic nose[END_REF][START_REF] Wang | Numerical studies of transonic BZT gas flows around thin airfoils[END_REF][START_REF] Rusak | Low-drag airfoils for transonic flow of dense gases[END_REF], both for one-dimensional unsteady and two-dimensional steady configurations. For the case of viscous dense gases, analytical studies were carried out to investigate the dissipative structure of shock waves and the impact of dense gas effects on laminar boundary layers and laminar shock/boundary layer interactions [START_REF] Cramer | The dissipative structure of shock waves in dense gases[END_REF][START_REF] Cramer | On the suppression of shock-induced separation in Bethe-Zel'dovich-Thompson fluids[END_REF][START_REF] Kluwick | Internal flows of dense gases[END_REF][START_REF] Kluwick | Shock boundary layer interactions in dense gases[END_REF][START_REF] Kluwick | Shock regularization in dense gases by viscous-inviscid interactions[END_REF], 2011). These studies pointed out that it is possible to greatly reduce, or even suppress, the separation induced by shock/boundary layer interactions for some optimal choice of the operating thermodynamic state. Moreover, dense gases are generally characterized by heat capacities much larger than those of gases of less molecular complexity, leading to small characteristic Eckert numbers (Ec = U ref /(c p T ref ), with U ref and T ref a reference velocity and temperature, respectively). Hence, the coupling between the dynamic and thermal boundary layers is much weaker than in flows of standard gases like air.

Theoretical findings have been confirmed by numerical simulations. Computations of inviscid dense gas flows have been presented, among others, by [START_REF] Argrow | Computational analysis of dense shock tube flow[END_REF], [START_REF] Brown | Nonclassical dense gas flows for simple geometries[END_REF], [START_REF] Cinnella | Aerodynamic performance of transonic Bethe-Zel'dovich-Thompson flows past an airfoil[END_REF]. [START_REF] Cinnella | Inviscid and viscous aerodynamics of dense gases[END_REF] performed the first numerical simulation of turbulent dense gas flows. Specifically, they investigated compressible turbulent boundary layers and transonic turbulent flows around airfoils by using the Reynolds-averaged Navier-Stokes (RANS) equations closed by standard turbulence models. The simulations were carried out under a set of working assumptions, namely: (i) the flow conditions being sufficiently far from the thermodynamic critical point, dramatic variations of the fluid specific heats and compressibility coefficients, as well as density fluctuations were neglected; (ii) the compressible RANS equations supplemented by an eddy viscosity turbulence model were supposed to provide an adequate description of the mean flow for the configurations under study; (iii) the turbulent heat transfer was modelled though a 'turbulent Fourier law', by introducing a turbulent Prandtl number, assumed to be roughly constant and of O(1) throughout the flow, as commonly done in PFG flows.

More recently, several other authors have carried out RANS simulations of turbulent dense gas flows, in particular for turbomachinery applications (e.g. [START_REF] Harinck | Computational study of a high-expansion ratio radial organic Rankine cycle turbine stator[END_REF][START_REF] Wheeler | The role of dense gas dynamics on organic rankine cycle turbine performance[END_REF], 2014;[START_REF] Sciacovelli | Numerical study of multistage transcritical organic rankine cycle axial turbines[END_REF]. However, the accuracy of RANS models for these flows has not been properly assessed up to now, due to the lack of reference data either experimental or numerical. In particular, reliable experimental results are extremely difficult to obtain -most attempts of constructing dense gas shock tube and nozzle experiments are at the preliminary design stage as yet -and all of the past and present attempts aim at getting experimental proofs of the existence of non-classical shock waves, rather than providing detailed turbulence data. For single-phase gases close to saturation conditions, the most impressive non-classical phenomenon, i.e. the formation of a rarefaction shock wave, has been observed experimentally by [START_REF] Borisov | Rarefaction shock waves near the critic liquid-vapour point[END_REF] and [START_REF] Kutateladze | Rarefaction waves in liquid and gas-liquid media[END_REF]; however, the interpretation of such results has been challenged by several authors [START_REF] Cramer | Shock formation in fluids having embedded regions of negative nonlinearity[END_REF][START_REF] Fergason | Theory for producing a single-phase rarefaction shock-wave in a shock tube[END_REF]. More recently, dense gas shock tubes and test rigs have been devised by [START_REF] Fergason | Theory for producing a single-phase rarefaction shock-wave in a shock tube[END_REF], [START_REF] Colonna | Design of the dense gas flexible asymmetric shock tube[END_REF], [START_REF] Spinelli | Design of a test rig for organic vapours[END_REF], but no experimental results have been published to the present date.

High-fidelity numerical datasets such as those provided by direct numerical simulations (DNS) represent an effective tool for gaining insight into the physics of turbulent dense gas flows. In the literature, DNS have been carried out for low Mach number turbulent flows of light supercritical fluids such as carbon dioxide and nitrogen, with a focus on heat-transfer properties, the so-called piston effect and dissipation effects in mixing layers [START_REF] Zappoli | Anomalous heat transport by the piston effect in supercritical fluids under zero gravity[END_REF][START_REF] Okong'o | Direct numerical simulation of a transitional supercritical binary mixing layer: heptane and nitrogen[END_REF][START_REF] Bae | Direct numerical simulation of turbulent supercritical flows with heat transfer[END_REF][START_REF] Okong'o | Small-scale dissipation in binary-species, thermodynamically supercritical, transitional mixing layers[END_REF][START_REF] Tanahashi | DNS of supercritical carbon dioxide turbulent channel flow[END_REF]. [START_REF] Selle | Large-eddy simulation of single-species flows under supercritical thermodynamic conditions[END_REF] have studied the influence of supercritical initial conditions on the decay of homogeneous turbulence for nitrogen, assumed to obey a cubic EoS, even though their analysis is limited to the study of the kinetic energy decay evolution at low turbulent Mach number (M t = 0.1). At these conditions, compressibility effects are negligible and the authors conclude that real gas thermodynamic effects have a negligible influence on the turbulence dynamics. [START_REF] Battista | Turbulent mixing of a slightly supercritical Van der Waals fluid at Low-Mach number[END_REF] extended the low Mach number Navier-Stokes equations for reactive perfect gas flows to a generic real gas equation of state. They analysed the turbulent mixing of slightly supercritical fluids at low Mach number by means of the VDW EoS, and found differences in the ligament structures. Much less is known about turbulent dense gas flows in the transonic and supersonic regime. The impact of dense gas effects on compressible homogeneous isotropic turbulence (CHIT) and, more generally, highly compressible dense gas turbulent flows, represents an open research field.

In this work, we investigate for the first time the impact of dense gas effects on the structure of compressible turbulence, with focus on high Reynolds number flows that are of interest for typical dense gas industrial applications. As a natural starting point, we select a configuration that has been well documented for perfect gases, namely, the decay of high turbulent Mach number CHIT in a periodic box, which may exhibit eddy shocklets. For this kind of flow, the peculiar behaviour of the speed of sound when Γ is close to or less than zero is expected to modify the turbulence structure significantly. Specifically, for turbulent flows of BZT fluids, non-classical eddy shocklets may appear, leading to significant changes in the probability distribution functions of several flow properties. In decaying CHIT, the turbulent Mach number decreases with time. As a consequence, dense gas effects are expected to have an impact most of all on the initial stages of the decay, which are dominated by large to medium scales of turbulence.

Much of the numerical work on incompressible or perfect gas CHIT is restricted to low/moderate Taylor-based Reynolds numbers (generally less than 500, see [START_REF] Wang | Effect of compressibility on the small-scale structures in isotropic turbulence[END_REF], [START_REF] Donzis | Fluctuations of thermodynamic variables in stationary compressible turbulence[END_REF]), restricting severely the range of scales in the inertial range undergoing a cascade mechanism unimpeded by dissipative processes [START_REF] Porter | Inertial range structures in decaying compressible turbulent flows[END_REF]. Addressing high Reynolds numbers is more difficult as the turbulent Mach number of the simulation increases. Dissipative processes are, of course, fundamental at the smallest scales, where the turbulent kinetic energy is converted into heat. For high Reynolds number incompressible isotropic decaying turbulence, [START_REF] Lesieur | Turbulence in Fluids[END_REF] identified two distinct stages: in the first stage, viscous effects are negligible, worm-like structures due to sheet roll-up development, and vortex stretching causes a large increase of enstrophy; in the second stage of the decay, the energy transferred to the smallest scales by means of the energy cascade mechanism start to be dissipated by viscosity. The transition from the first to the second stage corresponds to a peak in the time history of the total enstrophy, since no smaller structures can be created in the second phase and the existing ones are destroyed due to viscous effects. Lesieur then concludes that the initial stages of the decay are well represented by an inviscid model.

Euler-based turbulent simulations have been largely employed in the literature to investigate the behaviour of high Reynolds number supersonic turbulent flows, such as wakes of objects moving at supersonic speed or those encountered in the interstellar medium. High-resolution numerical simulations of both forced and decaying inviscid supersonic turbulence have been performed by Porter, Pouquet & Woodward (1992a,b, 1994), [START_REF] Porter | Inertial range structures in decaying compressible turbulent flows[END_REF][START_REF] Porter | Turbulence in compressible flows[END_REF], [START_REF] Porter | Measures of intermittency in driven supersonic flows[END_REF], [START_REF] Sytine | Convergence tests for the piecewise parabolic method and Navier-Stokes solutions for homogeneous compressible turbulence[END_REF], who provided information about the spectra, intermittency, energy transfers and scaling exponents for structure functions of velocity, density and entropy. For decaying turbulence, the impact of dilatational modes on the different temporal phases has been addressed by [START_REF] Kritsuk | The statistics of supersonic isothermal turbulence[END_REF], who performed large-scale three-dimensional (3-D) simulations of supersonic inviscid turbulence. They measured the velocity scaling exponents and the slope of the velocity power spectrum, finding substantial differences with respect to incompressible turbulence, and proposed an extension of Kolmogorov's theory taking into account compressibility effects. [START_REF] Benzi | Intermittency and universality in fully developed inviscid and weakly compressible turbulent flows[END_REF] performed high-resolution numerical simulations of homogeneous isotropic inviscid weakly compressible turbulence. They showed that the scaling properties of the velocity field in the inertial range were in excellent agreement with those observed in DNS of the Navier-Stokes equations. One important result of the studies cited above is that the dynamics of the inertial range is independent not only of viscosity but also of the detailed dynamics of the dissipative mechanism. As a consequence, the latter does not affect the statistical properties in the inertial range. Similar considerations have led several investigators to model the small-scale dissipative processes using hyperviscosity algorithms, e.g. [START_REF] Cook | Hyperviscosity for shock-turbulence interactions[END_REF], [START_REF] Lamorgese | Direct numerical simulation of homogeneous turbulence with hyperviscosity[END_REF].

Based on the above-mentioned results, in the present paper we adopt the Euler equations, written for gases governed by a general equation of state. The governing equations are discretized by using a ninth-order accurate scheme, which introduces a numerical dissipation based on tenth-order derivatives in smooth flow regions. This dissipates wavelengths close to the grid cutoff, while leaving larger scales almost unaffected.

The paper is organised as follows. The governing equations and the numerical method are described in § 2. Details of the numerical set-up (initial and boundary conditions) adopted for the CHIT simulations, as well as a validation against literature results for a PFG, are provided in § 3. In addition, in § 3 we also discuss the validity of the inviscid assumption by comparison with viscous results at various Reynolds numbers. Dense gases are characterized by values of the specific heat ratio closer to unity than diatomic gases like air or triatomic gases like carbon dioxide (see, e.g. [START_REF] Harinck | The influence of molecular complexity on expanding flows of ideal and dense gases[END_REF]. As a consequence, deviations with respect to classical PFG simulations for a diatomic gas can be due both to the different specific ratio and to the EoS. In order to distinguish the contributions of these two effects, in § 4.1 we present a preliminary study of the sensitivity of PFG solutions to the specific heat ratio of the gas. In § 4.2, dense gas compressible homogeneous isotropic turbulence at various initial turbulent Mach numbers is investigated by using the VDW model, and the results are compared with those of a PFG characterized by the same heat ratio. The main findings of the study are summarized in § 5.

Governing equations and numerical scheme

The flow model is based on the compressible Euler equations, whose validity for large-scale dynamics is discussed in § 3.2. The equations, written in integral form for a control volume V with boundary ∂V , are:

d dt V w dV + ∂V f • n dS = 0, (2.1)
with w = [ρ, ρv, ρE] T the conservative variable vector, f = [ρv, ρvv + p Ī, ρvH] T , the inviscid flux, v the velocity vector, E = e + |v| 2 /2 the specific total energy, H = E + p/ρ the specific total enthalpy, e the specific internal energy, Ī the unit tensor and n the outer normal to ∂V . The system is supplemented by thermal and caloric equations of state, respectively: p = p(ρ, T) and e = e(ρ, T).

(2.2a,b)

These satisfy the compatibility relation:

e = e ref + T T r c v,∞ (T ) dT - ρ ρ r T ∂p ∂T ρ -p dρ ρ 2 , (2.3)
where the subscript (•) ref indicates a reference state, c v,∞ (T) is the specific heat at constant volume in the ideal gas limit and the superscript (•) denotes auxiliary integration variables. For a thermally and calorically perfect gas, (2.2) become:

p = ρRT and e = c v T, (2.4a,b)
with c v = R/(γ -1) = constant and R = R/M (being R the universal gas constant and M the gas molecular weight).

In this paper we use the thermal Van der Waals equation of state, which models the gas molecules as rigid spheres of given radius and subject to inter-particle attractive forces. Such EoS captures the qualitative aspects of dense gas dynamics without the complexities and computational overhead associated with more complex, and accurate, thermodynamic models [START_REF] Thompson | Negative shock waves[END_REF]. This is written: (2.5) with: and Z c = p c v c /RT c the critical compressibility factor. Assuming a calorically perfect (polytropic) gas and combining (2.5) and ( 2.3), one obtains an explicit relationship for the pressure in terms of density ρ and internal energy e: p = (γ -1)ρe + aρ 2 1 -bρ -aρ 2 .

p = RT v -b - a v 2 ,
a = 9p c v 2 c 8Z c b = v c 3 . ( 2
(2.7)

For a polytropic gas obeying the VDW EoS, the fundamental derivative of gas dynamics becomes:

Γ = 1 2 γ (γ + 1)RT (1 -bρ) 3 -6aρ γ RT (1 -bρ) 2 -2aρ .
(2.8)

A region of negative values of the fundamental derivative appears if the specific heat ratio is in the range 1 < γ 1.06 [START_REF] Thompson | Negative shock waves[END_REF]. Figure 1 shows some iso-Γ lines and the inversion zone for a VDW gas with γ = 1.0125, which is roughly representative of heavy fluorocarbons.

In the following simulations, we consider the perfluoro-perhydrophenanthrene, (chemical formula C 14 F 24 ), called hereafter PP11 (commercial name); its thermodynamic properties are provided in table 1. This fluid has been often used in the literature since it exhibits a wide inversion zone and, as a consequence, significant dense gas effects. For such a complex gas, the specific heat ratio γ varies with the temperature. However, since the flow conditions investigated in this work are characterized by small temperature variations, we assume hereafter that γ is constant and equal to 1.0125, which is representative of heavy fluorocarbons [START_REF] Brown | Nonclassical dense gas flows for simple geometries[END_REF].

Time integration is carried out by means of a low-storage six-stage Runge-Kutta scheme [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF]. The flux derivatives in the governing equations are approximated in space by means of centred tenth-order accurate differences. In order ), critical density (ρ c ), critical pressure (p c ), critical compressibility factor (Z c ), acentric factor (ω), dipole moment of the gas phase (d.m.), boiling temperature (T b ) and ratio of ideal gas specific heat at constant volume over the gas constant (c v (T c )/R) at the critical point.

M T c ρ c p c Z c ω d.m. T b c v (T c )/R g mole -1 K kg m -3 MPa - - D K - PP11 (C
to introduce a minimal amount of numerical dissipation while ensuring computational robustness for compressible flow simulations, the centred scheme is supplemented by a high-order artificial viscosity term, similar to the approach of [START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF] and [START_REF] Kim | Adaptive nonlinear artificial dissipation model for computational aeroacoustics[END_REF]. In our case, we use a blend of second-and tenth-order derivatives approximated by standard central differences.

For a 1-D problem and a regular Cartesian grid with constant mesh spacing δx (so that x j = j δx), the semi-discrete tenth-order scheme in space writes:

(w t ) j + (δF ) j δx = 0, (2.9)
where δ is the classical difference operator over one cell, δ(•

) j := (•) j+1/2 -(•) j-1/2
and F j+1/2 is the numerical flux at cell interface j + 1/2: (2.10) where f = f (w) is the physical flux and µ is the cell average operator, µ(•) j+1/2 := ((•) j+1 + (•) j )/2. A high-order dissipation term is added to prevent high-frequency oscillations. Furthermore, the scheme is modified by introducing a nonlinear artificial viscosity term to capture flow discontinuities. The numerical flux of the scheme supplemented by the dissipation term finally becomes:

F j+1/2 = I -1 6 δ 2 + 1 30 δ 4 -1 140 δ 6 + 1 630 δ 8 µf j+1/2 ,
F j+1/2 = I -1 6 δ 2 + 1 30 δ 4 -1 140 δ 6 + 1 630 δ 8 µf -D j+1/2 (2.11)
with the numerical dissipation term D:

D j+1/2 = ρ(A) j+1/2 [ε 2 δw + ε 10 δ 9 w] j+1/2
(2.12) with ρ(A) the spectral radius of the flux Jacobian matrix A = ∂f /∂w, and

ε 2j+1/2 = κ 2 max(ν j Φ j , ν j+1 Φ j+1 ), ε 10j+1/2 = max(0, k 10 -ε 2j+1/2 ), (2.13a,b)
with k 2 and k 10 adjustable dissipation coefficients. In the above, ν j is the well-known pressure-based shock sensor of [START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF] and Φ is the discrete form of Ducros's sensor [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF],

Φ = (∇ • u) 2 (∇ • u) 2 + ω 2
(2.14) 1) in high-divergence regions and tends to zero in vortex dominated regions, thus allowing capturing of flow discontinuities sharply while minimizing the effect of numerical dissipation on vortical structures. Note that when k 2 = 0 and k 10 = 1/1260, the preceding method becomes equivalent to a ninth-order accurate upwind scheme.

If not specified, the present compressible computations were carried out using k 2 = 2 and k 10 = 1/420. The preceding scheme has very low phase and dissipation errors. Far from flow discontinuities, its leading truncation error term is of the form k 10 δx 9 ∂ 10 w/∂x 10 , i.e. it is consistent with a tenth-order viscosity. Such a viscosity term acts differently according to the wavenumber. It dissipates scales characterized by reduced wavenumbers of approximately 0.35π or higher, i.e. wavelengths that are discretized with less than 6 mesh points, leaving larger scales essentially unaffected. By comparing the Fourier-space representation of a given mode of the numerical solution to the exact transfer function of an advection-diffusion equation, the effective viscosity ν eff introduced by the present scheme is found to be:

ν eff = CFL δx 2 2 10 k 10 sin 10 ξ /2 ξ 2 , (2.15)
where the Courant-Friederichs-Lewy (CFL) number is of O(1) in the present simulations and ξ is the reduced wavenumber. Plots of the dissipation function, of the normalized phase error and of the effective viscosity introduced by the present scheme are provided in figure 2. The scheme being of odd order, dispersion errors are smaller than dissipation ones. For the smallest wavenumbers resolved on a given grid, the effective viscosity is of the order of 10 -5 .

Set-up of CHIT simulations

Initial and boundary conditions

The CHIT decay problem is solved on a cubic computational domain with extension [0, 2π] 3 . Periodic boundary conditions are imposed in the three Cartesian directions.

The issue of the initial conditions for compressible isotropic turbulence has been addressed by various authors [START_REF] Blaisdell | Compressibility effects on the growth ano structure of homogeneous turbulent shear flow[END_REF][START_REF] Ristorcelli | Consistent initial conditions for the DNS of compressible turbulence[END_REF][START_REF] Samtaney | Direct numerical simulation of decaying compressible turbulence and shocklet statistics[END_REF]. In general, the shape of the initial three-dimensional spectrum and the r.m.s. level of each flow variable must be provided. Furthermore, different spectra can be assigned to the solenoidal and dilatational components. The main difference between an incompressible and compressible initialization is the presence of an intrinsic velocity scale for the latter, which is related to the speed of sound. The initial r.m.s. velocity u rms is defined by prescribing the turbulent Mach number (u rms = M t c , where c is the prescribed average speed of sound). Temperature and pressure fluctuations are specified in accordance with the velocity fluctuations. Several simplifying hypotheses are usually made in the compressible case [START_REF] Blaisdell | Compressibility effects on the growth ano structure of homogeneous turbulent shear flow[END_REF]): usually, the same spectrum shape is imposed for all the fluctuating fields; furthermore, fluctuations of the thermodynamic quantities are neglected. A review of initialization methods for CHIT can be found in [START_REF] Samtaney | Direct numerical simulation of decaying compressible turbulence and shocklet statistics[END_REF], as well as an analysis of their influence on the evolution of turbulent quantities. The latter authors show that for turbulent Mach numbers up to 0.5 (the maximum value considered in their analysis), initial conditions have moderate influence on the turbulence decay. Incompressible-like initialization leads to a fast transient in which the dilatational velocity component grows and becomes coherent with the solenoidal one, while fluctuations of thermodynamic quantities develop. This transient leads to higher values of the r.m.s. of the velocity divergence. Note that, for PFG, neglecting fluctuations of the thermodynamic quantities in the initial conditions has been found to lead to weaker compressibility effects on the turbulence decay [START_REF] Blaisdell | Compressibility effects on the growth ano structure of homogeneous turbulent shear flow[END_REF].

For PFG, the initialization is similar to the one described in [START_REF] Sarkar | The analysis and modelling of dilatational terms in compressible turbulence[END_REF] and [START_REF] Pirozzoli | Direct numerical simulations of isotropic compressible turbulence: influence of compressibility on dynamics and structures[END_REF]. In the case of dense gases that use a complex EoS, additional assumptions are required. Since the thermodynamic variables are nonlinearly related, contrary to the perfect gas case, it is not possible to assume a direct proportionality between density (or temperature) and dilatational velocity fluctuations.

In all of the following simulations (PFG and DG), we have imposed an initial velocity spectrum of the Passot-Pouquet type [START_REF] Passot | Numerical simulation of compressible homogeneous flows in the turbulent regime[END_REF]):

E(k) = Ak 4 exp -2 k k 0 2 , (3.1)
where k 0 is the peak wavenumber and A is a constant that depends on the initial amount of kinetic energy. For small values of k 0 , this distribution associates most of the energy to the largest scales and practically none to the smallest ones. As a consequence, the large-scale dynamics (which is of interest here) is enhanced, while the impact of the small scales (which are not meaningful in an inviscid simulation) is reduced.

To ensure fair comparisons between perfect and dense gas simulations, fluctuations of the thermodynamic quantities are set to zero, and the turbulent velocity field is assumed to be purely solenoidal in all of the following simulations. With this choice, turbulent scales set at the beginning of the simulation can be controlled accurately. The expression for the initial turbulent kinetic energy spectrum is used to compute analytically the initial turbulent kinetic energy K 0 , the initial enstrophy Ω 0 , the initial integral length scale L I and the initial large-eddy turnover time τ LE :

K 0 = 3A 64 √ 2πk 5 0 , Ω 0 = 15A 256 √ 2πk 7 0 , L I = √ 2π k 0 , τ LE = 32 A (2π) 1/4 k -7/2 0 . (3.2a-d)
A sensitivity study to the initial conditions, and specifically to the initial compressibility ratio (χ 0 , the ratio of dilatational to solenoidal energy) and initial peak wavenumber has been conducted in the case of a perfect gas. The results, omitted for brevity, show in particular that the main effect of setting compressibility ratio greater than zero is to enhance compressibility effects, especially at the initial stages of the decay, as already pointed out by [START_REF] Blaisdell | Compressibility effects on the growth ano structure of homogeneous turbulent shear flow[END_REF]. For the investigation of dense gas effects in CHIT, we have made the conservative choice of setting the initial compressibility ratio to zero.

Concerning the choice of the initial peak wavenumber, we have chosen k 0 = 2, thus concentrating approximately 99 % of the initial energy at large scales (the ones that are well resolved on the grids used in the present study). As a consequence, the initial stages of the decay are practically unaffected by numerical dissipation. Then, the kinetic energy spectrum broadens due to vortex stretching mechanisms and, for non-dimensional times greater than approximately 2.5, a well-identified inertial range is observed.

Assessment of the numerical strategy

The numerical method used in this work is first assessed versus results available in the literature for inviscid CHIT of a perfect gas. Results provided by the present ninth-order scheme, implemented within an in-house code [START_REF] Outtier | The high-order dynamic computational laboratory for CFD research and applications[END_REF]) are compared to those obtained by [START_REF] Garnier | On the use of shockcapturing schemes for large-eddy simulation[END_REF], who studied the capability of various shock-capturing schemes, including Jameson's [START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF], TVD-MUSCL (Van Leer 1979) and essentially non-oscillatory (ENO) [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF][START_REF] Shu | Numerical experiments on the accuracy of ENO and modified ENO schemes[END_REF][START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF], to behave as implicit subgrid-scale models for several schemes. For the sake of comparison, we also run computations with the classical lower-order scheme by Jameson available in our code. Computational grids made of 64 3 and 128 3 cells have been used, and all of the results are expressed in terms of the non-dimensional time of [START_REF] Garnier | On the use of shockcapturing schemes for large-eddy simulation[END_REF] based on a unit length and the initial r.m.s. velocity. Note that for the sake of comparison, we use the same non-dimensional time throughout our study.

Figure 3 shows results obtained for case 1 of [START_REF] Garnier | On the use of shockcapturing schemes for large-eddy simulation[END_REF], corresponding to an initial turbulent Mach number equal to 0.2 and initial compressibility ratio equal to 0, but similar results were obtained for the other cases (not reported for the sake of brevity). For comparison, we show in figure 3 the time histories of kinetic energy, enstrophy and pseudo-Taylor microscale λ [START_REF] Porter | Kolmogorov-like spectra in decaying three-dimensional supersonic flows[END_REF], as well as the kinetic energy spectrum at the non-dimensional time t = 10. Note that λ, defined in the appendix A, should tend to zero in a purely inviscid simulation. Here, due to the numerical dissipation introduced by the scheme, it tends to a small but non-zero value. λ can be seen as a measure of the smallest scales resolved by a numerical method, and its limit value becomes smaller as the grid is refined or dissipation errors are reduced. As one would expect, the ninth-order scheme improves the results dramatically not only with respect to the low-order scheme (which exhibits a very dissipative behaviour), but also with respect to a fifth-order accurate MENO scheme. We observe that: (i) the kinetic energy is preserved for longer times (figure 3a);

(ii) the enstrophy peak is shifted to the right and the maximum value is more than five time greater than the one obtained with Jameson's scheme and more than 2.5 times the one provided by the MENO scheme (figure 3b); (iii) the pseudo-Taylor microscale is reduced by approximately a factor 6 for times greater than 4 (where the solution is dominated by the smaller scales) with respect to the low-order scheme, and more than a factor 2 with respect to the fifth-order scheme (figure 3c); (iv) the cutoff in the kinetic energy spectrum is also moved toward wavenumbers that are closer to the grid aliasing limit (figure 3d). Finally, figure 4 reports the time evolution of the r.m.s. density values at various initial turbulent Mach numbers. Comparisons with the results provided by [START_REF] Garnier | On the use of shockcapturing schemes for large-eddy simulation[END_REF] for a ENO scheme match to within plotting accuracy. The sensitivity of the present numerical simulations to mesh resolution has been assessed for both a perfect gas (γ = 1.4) and a VDW gas, at M t,0 = 0.8. For that purpose, we have considered four grids with the number of cells equal to 64 3 , 128 3 , 256 3 and 512 3 . The time histories of λ, ρ rms and of Ω normalized with the initial enstrophy, as well as the kinetic energy spectrum and the probability distribution function of the velocity divergence at t = 2.5, are reported in figure 5 for the PFG.

Figure 5 shows that ρ rms is little affected by grid resolution and is well captured even on a 64 3 grid (a). The same is true for the average and r.m.s. of the other thermodynamic quantities. Differences measured between the two finest grids are found to be of the order of 1 % or lower. Resolution has a stronger influence on the pseudo-Taylor microscale (b), which is reduced by a factor of approximately 1.6 when halving the mesh size [START_REF] Garnier | On the use of shockcapturing schemes for large-eddy simulation[END_REF]. We observe that the enstrophy should become unbounded within a finite time in a purely inviscid simulation and should exhibit a peak increasing with the Reynolds number in viscous simulations. In the present calculations, due to numerical dissipation, the enstrophy peaks at a non-dimensional time of the order of 5. Due to the creation of smaller scales, its amplitude tends to double when doubling the mesh (c). Note that the solutions on the two finest grids are superposed up to t = 2 and tend to depart from each other at later times. At time t = 2, a significant inertial range has not developed yet, and most of the energy is still contained in scales much larger than the grid size. At time t = 2.5 (considered for the analyses discussed later), a well-identified inertial range is observed already on the 128 3 grid (see d), and enstrophy values computed on the two finest grids differ by approximately 5 %. The inertial range extends approximately over more than a decade for a grid resolution of 256 3 , and slightly less than two decades for a resolution of 512 3 (d). The exponential roll off of the spectra at the end of the inertial range is due to the numerical viscosity cutoff, and moves toward higher wavenumbers as grid resolution increases. The observed cutoff wavenumber is rather close to the grid resolvability limit (four points per wavelength), in agreement with the spectral properties of the current numerical method (see § 1). With the aim of assessing the influence of the numerical cutoff on local turbulence properties, the probability density function (p.d.f.) of the velocity divergence θ /θ rms normalized with respect to M 2 t,0 at t = 2.5 is also reported (figure 5e). Grid resolution modifies mainly the tails of the distribution, which are however characterized by rather low probability values (less than 10 -3 ). This leads to some differences between the various grids in terms of integrated probabilities, like the volume fractions occupied by flow structures with given divergence levels. Nevertheless, the orders of magnitude are well conserved, and the associated percentages are reliable at least to one significant digit or more, according to the considered divergence level. Similar trends are observed for different initial turbulent Mach numbers (not reported) and for VDW simulations (only results for the more sensitive quantities on the two finest grids are shown in figure 6). The p.d.f. exhibits a shape different than the one observed for a perfect gas, independently of the grid. The physical mechanism leading to this different behaviour is discussed in the next section.

The 256 3 grid leads to reasonably small numerical errors both on the mean and r.m.s. properties, as well as an adequate description of the inertial range over more than a decade at the non-dimensional time of interest (t = 2.5), and it has been selected for the parametric study discussed in later sections.

As pointed out in the introduction, due to the high-order dissipation of the scheme, the present Euler-based simulations can be considered as a form of implicit large-eddy simulation (LES) of high Reynolds number CHIT. To further assess the numerical model, we have also carried out a series of Navier-Stokes-based simulations (both for PFG and DG), whereby the viscous fluxes are discretized by standard fourth-order central differences. Comparisons were first made for PFG CHIT at M t,0 = 0.5 and Re λ = 175, 1750 and 17 500 (Re λ being based on the initial Taylor microscale; all simulations were carried out on the 256 3 grid (and setting k 0 = 2 and χ 0 = 0). For the Euler-based simulation, the effective Reynolds number (Re eff λ ) is estimated by using the numerical viscosity of the scheme. In particular, the numerical viscosity associated with the smallest wavenumbers is of O(10 -5 ) and Re eff λ = O(10 4 ), which is indeed of the same order as the highest Re λ considered in the Navier-Stokes-based simulations. We observe that at t = 0, the smallest wavelength that can be resolved on the given grid, k max η (η being the Kolmogorov length scale), is equal to 4.8, 1.5 and 0.48, respectively for the three Re λ . Hence, for the lowest Reynolds number case, all scales are well resolved up to the Kolmogorov scale. The other viscous simulations do not resolve all of the scales, and represent instead implicit LES, with the numerical dissipation acting as a subgrid model. Figure 7 compares the time histories of the Euler-based and Navier-Stokes-based results in terms of ρ rms , c rms , kinetic energy and θ rms . The figure shows that, for the highest Reynolds number cases, the initial stages of the decay are well represented by inviscid simulations. At t = 2.5, the selected grid resolves accurately large and medium scales up to a cutoff wavenumber well within the inertial range for the high-Re λ cases. This is clearly shown in figure 8(a), where we compare the kinetic energy spectra. As the Reynolds number increases, the inertial range becomes larger and it extends up to a wavenumber k v of approximately 10 for Re λ = 1750 and approximately 20 for Re λ = 17 500 (for Re λ = 175 there is no scale separation and the inertial range is virtually non-existent). We observe that the inviscid model captures essentially the same energy spectrum as the viscous one at Re λ = 17 500, the effective Reynolds number being of the same order.

The comparison of the p.d.f. of the velocity divergence (figure 8b) shows that the Euler-based and the Navier-Stokes-based results are almost superposed for all cases (except for the lowest Re λ , for which the probability content of the tails is smaller; however, even in this unfavourable case the overall shape of the p.d.f. remains similar).

Similar comparisons were also carried out for dense gas VDW simulations. The kinetic energy spectrum and the p.d.f. of normalized dilatation at t = 2.5 are reported in figure 9, which shows that the Euler-based model gives results very close to the high Reynolds simulations. Note again the peculiar p.d.f. that exhibits a more symmetric shape than that of a PFG (this point will be further discussed in § 4).

Finally, the influence of the artificial dissipation coefficients k 2 and k 10 on the turbulence decay has been investigated considering a 256 3 grid and M t,0 = 0.8 and χ 0 = 0. Four different parameter combinations were considered (see table 2). Sample results are shown in figure 10. The evolutions of kinetic energy (a) and r.m.s. of the thermodynamic properties (not reported for brevity) are essentially insensitive to k 2 and k 10 . Small changes in the enstrophy time history (b) are observed for t greater than ≈3. At t = 2.5, the turbulent kinetic energy spectrum (c) is slightly dependent on k 10 at high wavenumbers, whereas it is found to be largely insensitive to k 2 . Finally, the p.d.f. of the normalized velocity divergence at the same time (d) shows that the coefficients have only a small influence on the tails of the distribution.

Compressible homogeneous isotropic turbulence

This section presents a parametric study of the influence of dense gas effects on CHIT decay. Due to their greater molecular complexity, these gases are generally characterized by values of the specific heat ratio closer to unity than diatomic gases like air or triatomic gases like carbon dioxide (see, e.g. [START_REF] Harinck | The influence of molecular complexity on expanding flows of ideal and dense gases[END_REF]. It is useful to recall that, for gases with γ close to one, the specific heat coefficients tend to infinity (this is easily seen, e.g. for calorically perfect gases); hence, isentropic transformations are also approximately isothermal.

To the authors' knowledge, all of the results regarding CHIT available in the literature up to now have been obtained for air, modelled as a perfect gas with γ = 1.4. Since changing the specific heat ratio has an impact on the coupling of the kinematic and thermodynamic fields, we first perform a parametric investigation of the effect of γ on perfect gas CHIT simulations, and generate reference results for γ = 1.0125 that is a value representative of the PP11 gas of interest here. We then run dense gas simulations using the VDW EoS with γ = 1.0125, and setting k 0 = 2 and χ 0 = 0.

Perfect gas simulations: influence of the specific heat ratio

In table 3 we report the different values of γ considered, as well as the values of the specific heat coefficients and fundamental derivative. Simulations were performed on the 256 3 grid and assuming an initial turbulent Mach number M t,0 = 0.8.

Using dimensional analysis and momentum conservation, the pressure fluctuations can be related to the velocity fluctuations:

p rms ∼ ρ 0 u 2 rms 0 ∼ ρ 0 M 2 t,0 c 2 0 ; (4.1) hence: p rms p 0 ∼ γ M 2 t,0 . (4.2)
On the other hand, the normalized density fluctuation is related to the normalized pressure fluctuations through the sound speed definition ρ rms /ρ 0 ∼ p rms /c 2 0 ; hence using (4.2), one obtains

ρ rms ρ 0 ∼ M 2 t,0 . (4.3)
Recalling the definition of the fundamental derivative of gas dynamics (1.1) and its expression for PFG (1.2)), one obtains:

c rms c 0 ∼ ρ rms ρ 0 (Γ -1) ∼ ρ rms ρ 0 γ -1 2 ∼ γ -1 2 M 2 t,0 . (4.4)
The time evolution of the normalized r.m.s. of the speed of sound, pressure and density at various γ are reported in figure 11, where p rms and c rms are rescaled according to (4.2) and (4.4). The results show that changing the specific heat ratio has a strong influence on c rms , whose maximum value varies between 0.002 (for the low-γ case) and 0.84 (for the high-γ case). This is due to both direct and indirect effects. On the one hand, the speed of sound depends on √ γ for a perfect gas; on the other hand, lower values of γ correspond to higher values of the specific heat, and hence to smaller heating effects and lower values of the mean temperature. For the case of a polyatomic heavy gas (case G4 of table 3), characterized by high values of the specific heats, the mean sound speed remains roughly constant with time. The opposite behaviour is observed for monoatomic gases that are characterized by specific heats lower than air. Intermediate behaviours are observed for cases G2 and G3. The p rms exhibits trends similar to c rms , with higher fluctuations for case G4, in which ρ rms shows an opposite behaviour even though it is much less dependent on the value of γ .

The time evolution of M t , K, ω rms and θ rms , as well as the turbulent kinetic energy spectra and the p.d.f. of the local Mach number at 2.5 are reported in figure 12, at various γ . The figure shows that the vorticity fluctuations, the turbulent kinetic energy and its spectra are not affected by γ . The figure also shows that the probability of finding local Mach numbers greater than unity is greater for lighter gases.

In table 4 we report the volume fractions occupied by weak (W), moderate (M) and strong (S) compressions (C) and expansions (E), which are defined in terms of normalized dilatation rescaled by M 2 t,0 [START_REF] Wang | Effect of compressibility on the small-scale structures in isotropic turbulence[END_REF]. We observe that, independently of γ , weak expansions are more probable than compression ones (approximately 52 %-38 %). On the contrary, moderate and strong compression regions are more likely to occur than expansion ones. However, overall weak compressions and expansions occupy more than 90 % of the flow volume, and the probability of the occurrence of eddy shocklets (defined as regions for which θ /θ rms < -3) exists, but it is rather small.

Dense gas simulations

In order to assess dense gas effects on CHIT, we have considered the perfluoroperhydrophenanthrene that is assumed to obey the VDW EoS with γ = 1.0125. The results are compared to those of a perfect gas having the same specific heat ratio. Since the PFG and VDW models lead to different intrinsic scales (the sound speed being very different), it is not possible to set the same initial kinetic energy without changing the thermodynamic conditions. Consequently, comparisons were performed by renormalizing the calculated quantities with respect to their initial values. The initial thermodynamic state is reported in the p-v plane (figure 13) for the VDW model; the initial thermodynamic conditions for both perfect and dense gas simulations are summarized in table 5, where all variables are normalized with respect to their critical value. Figure 13 shows the location of the (uniform) initial thermodynamic state in the p-v plane for the VDW model. The initial state is located close to the transition line and corresponds to a negative initial value of the fundamental derivative of gas dynamics. As a consequence, non-classical gas dynamic behaviours are likely to appear in supersonic flow regions. We have then carried out a parametric study by considering various initial turbulent Mach numbers (M t,0 = 0.2, 0.5, 0.8 and 1).

The computed kinetic energy spectra at t = 2.5, reported in figure 14, exhibit an inertial range for all values of the initial turbulent Mach number and the perfect and dense gas solutions do not show any difference. However, the thermodynamic variables as well as the turbulent Mach number are strongly affected by the dense gas effects as shown in figure 15, where the time histories of ρ rms , M t , c and c rms (normalized by the respective initial values) are reported at various M t,0 . In particular, we observe that for the lower M t,0 cases (M t,0 = 0.2, 0.5) there are no significant differences between PFG and VDW solutions. For the PFG, the average and the r.m.s. FIGURE 14. PFG and DG simulations. Kinetic energy spectra at t = 2.5 at various initial turbulent Mach numbers. ----, M t,0 = 0.2; -----, M t,0 = 0.5; --, M t,0 = 0.8; ---, M t,0 = 1. E , VDW; A , PFG. speed of sound remains almost constant in time independently of M t,0 . In the VDW dense gas case these quantities increase strongly with M t,0 . As a consequence the mean and r.m.s. fundamental derivative of gas dynamics (which is directly related to sound speed variations) also exhibits a strong dependence on M t,0 (figure 16). This behaviour is due to the fact that for VDW gases both c and Γ depend on temperature and density, which are changing. For the lowest values of M t 0 , fluctuations of the thermodynamic quantities are small enough, so that the thermodynamic states remain close to the initial one throughout the flow. This is confirmed by inspection of figure 17 that shows the representation of the thermodynamic states in the p-v diagram at t = 2.5 for various M t,0 . The thermodynamic states are distributed along the isentropic line corresponding to s = s 0 (s 0 being the initial entropy). At M t,0 = 0.2 (a), the thermodynamic variables exhibit weak variations and Γ is rather constant in time. As the initial turbulent Mach number increases, thermodynamic states attained by the flow lay along the isentropic line and the thermodynamic variables exhibit strong variations, thus explaining the strong dependence of Γ and Γ rms observed in figure 16. The rapid growth of both quantities (attaining a maximum at about 1/1.5) is likely due to the appearance of strong local compression phenomena (namely, eddy shocklets). A close up of the instantaneous isocontours of Γ , taken along a mesh plane at t = 2.5, for a dense gas simulation with M t,0 = 1 is shown in figure 18. The figure shows that there exist regions of the plane where Γ varies abruptly from negative to positive values (from approximately -0.5 to values greater than 4). The representation in the p-v plane of the thermodynamic states at the same time (figure 17) shows that, indeed, for the high turbulent Mach numbers, strong compressions drive the thermodynamic states far away from the initial one, up to supercritical pressures and temperatures. In the high-pressure region, the second nonlinearity parameter (the rate of change of Γ with respect to density at constant entropy) is positive, and Γ increases rapidly with pressure, leading to the formation of ----, M t,0 = 0.2; -----, M t,0 = 0.5; --, M t,0 = 0.8; ---, M t,0 = 1. E , VDW; A , PFG.

spots characterised by very high values. As a consequence, the sound speed decreases abruptly, leading to large variations of the local Mach number.

To further analyse compressibility effects, figure 19 displays the time evolution of the r.m.s. velocity divergence normalized by the initial r.m.s. vorticity at various M t,0 . During the initial transient phase, the development of the dilatational modes exhibits similar trends for both perfect and dense gas. Afterwards, compressibility effects are enhanced for dense gas cases.

Regardless of the gas model in use, the observed dynamics is similar. The almost-incompressible initial distribution of large-scale eddies is modified at the beginning of the decay with the appearance of compressible modes, leading to the generation of eddy shocklets that start interacting at a later stage. As observed by [START_REF] Passot | Numerical simulation of compressible homogeneous flows in the turbulent regime[END_REF], when shocklets collide they may either coalesce or 'ignore' each other, pursuing their path with a slightly different velocity due to momentum exchange. Because of turbulence kinetic energy decay, the number and intensity of shocklets tend to decrease with time. Nevertheless, the natures of the eddy shocklets for the perfect and dense gas are different. Figure 20 ----, M t,0 = 0.2; -----, M t,0 = 0.5; --, M t,0 = 0.8; ---, M t,0 = 1. E , VDW. divergence greater than 3, which identifies strong expansions. We observe that strong compression zones are somewhat more frequent and wider in PFG than in dense gas simulations. Conversely, strong expansion regions, which occupy a significant portion of the domain in the dense gas solution, are nearly absent in PFG.

θ/θ rms • M 2 t,0 (%) M t0 [-∞, -2] [-2, -1] [-1, 0] [0, 1] [1, 2] [2,
The isosurfaces θ/θ rms = 3 and θ/θ rms = -3 at t = 2.5 are reported in figures 21 and 22, respectively, both for PFG and DG. The figures show that in the PFG case, strong compression regions correspond to sheet-like structures, while strong expansions have tubular shapes. On the contrary, the dense gas exhibits sheet-like and tubular structures, both in compression and in expansion regions.

Table 6 reports the volume fractions occupied by compression and expansion regions at t = 2.5, as defined in § 4.1. At the lowest M t,0 , only weak expansion and compression regions are present both for PFG and dense gas cases, with a balanced distribution. The percentage of moderate and strong compression and expansion regions increases as the initial turbulent Mach number increases. The volume distribution for the PFG case with M t,0 = 1 is similar to the one reported in [START_REF] Wang | Effect of compressibility on the small-scale structures in isotropic turbulence[END_REF]. Strong compression regions become significant starting from M t,0 = 0.8. At M t,0 = 1, strong dilatation regions (including compressions and expansions) occupy approximately 5 % of the total volume both for perfect and dense gases. However, in the latter case, expansions and compressions are equally probable and the volume fraction shows a more symmetric statistical behaviour (which is preserved for all M t,0 due to the non-classical sound speed behaviour in zones with Γ < 0). For PFG the volume fraction distribution is symmetric only at the lowest M t,0 . At M t,0 = 1 strong expansion zones occupy a volume fraction which is approximately 400 % less than strong compressions. In dense gases strong expansions occupy a volume fraction that is approximately 10 % less than strong compressions at the same M t,0 . In the PFG case, moderate and weak expansion occupy 57 % of the flow volume, while the weak and moderate compressions occupy 38 %. However, in VDW case the volume occupied by weak and moderate compressions and expansions is approximately equal (approximately 47 % each). To further elucidate the influence of dense gas and compressibility, in figure 23 we report the p.d.f. of the normalized dilatation at various M t,0 (t = 2.5) both for PFG and VDW. The dashed vertical lines identify the regions of strong compression bounded by the area of negative (respectively, positive) values of θ /θ rms to the left (respectively, to the right) of θ /θ rms = -3 (respectively, θ /θ rms = 3). In the case of the perfect gas, the p.d.f. of the dilatation is highly skewed, as already observed in weakly and moderately compressible turbulent flows at high M t,0 (and γ = 1.4) by [START_REF] Porter | Measures of intermittency in driven supersonic flows[END_REF], [START_REF] Pirozzoli | Direct numerical simulations of isotropic compressible turbulence: influence of compressibility on dynamics and structures[END_REF]. In the dense gas case, strong expansions are enhanced, resulting in a more symmetric p.d.f. (i.e. strong compressions and strong expansions are equally probable). Table 7 shows the volume fractions occupied by flow regions having different Γ values conditioned on the dilatation for the dense gas case at M t,0 = 1. We observe that 40 % of the flow volume is characterized by thermodynamic states falling into the inversion zone (exhibiting Γ < 0), whereas only 4 % is characterized by Γ > 1. Strong expansion regions characterized by a negative Γ occupy a volume fraction of approximately 1 %, which represents 50 % of all Γ values, whereas strong expansion regions having Γ > 1 represent a negligible fraction. ----, M t,0 = 0.2; -----, M t,0 = 0.5; --, M t,0 = 0.8; ---, M t,0 = 1. E , VDW; A , PFG. Strong compression regions having positive Γ occupy approximately 70 % of the total volume. This suggests that in the dense gases the occurrence of both eddy shocklets and eddy 'expansion shocklets' is admissible.

In figure 24 we report the conditioned distributions of Γ at various M t,0 (a) and the p.d.f. of the normalized local speed of sound at M t,0 = 1 (b). In compression regions, Γ is positive and varies considerably with θ for the higher M t,0 values; accordingly, the local thermodynamic state moves towards the supercritical conditions (see also figure 17). The changes in the concavity are due to the fact that, for the chosen initial thermodynamic conditions, weak isentropic compressions drive the flow toward thermodynamic regions with lower Γ . Conversely, weak expansions drive the flow outside the inversion region. For highly positive levels of the dilatation, the flow is statistically far away from the inversion zone. However, variations associated with high positive dilatation are smaller since Γ varies more slowly when moving to the right part of the pv plane (see again figure 17).

In figure 25 we report the joint p.d.f. of Γ with θ /θ rms at various initial turbulent Mach numbers at time t = 2.5. For M t,0 = 0.2, the weak dilatation levels cause the p.d.f.s to be concentrated in the neighbourhood of the initial condition (Γ ≈ -0.045). At M t,0 = 0.5, the p.d.f. becomes broader, but most of the states are characterized by Γ < 1 and θ /θ rms ∈ [-1, 1]. At higher M t,0 , due to strong dilatation effects, the thermodynamic states spread over a large range of thermodynamic conditions, and a secondary peak is obtained around Γ ≈ 4 and small positive values of θ /θ rms . When strong compressions occur, flow particles becomes increasingly difficult to compress

Γ |θ /θ rms [-∞, -2] [-2, -1] [-1, 0] [0, 1] [1, 2] [2, +∞] [-∞, +∞] SC (%) MC (%) WC (%) WE (%) ME (%) SE (%) Total [-∞, 0] 1 4 16 14 4 1 40 [0, 1] 2 5 21 22 5 1 56 [1, +∞] 0 1 1 2 1 0 4 TABLE 7.
Volume fraction occupied by flow regions characterized by different Γ values, conditioned on dilatation levels (dense gas simulations at M t,0 = 1, t = 2.5). First row corresponds to regions characterized by Γ < 0 (thermodynamic states falling in the inversion zone); second row to regions with 0 < Γ < 1 (dense gas region); third row to regions having Γ > 1. The last column provides the total volume fraction for each thermodynamic region. as density grows. This results in a clustering of the thermodynamic states in the highpressure region, close to the critical isotherm.

Figure 26 shows that, in the case of dense gas, the p.d.f. of the normalized local speed of sound c/ c for the dense gas spans over a wider range than PFG for a Dirac-like distribution is recovered (i.e. the speed of sound is nearly constant). This is true for the total p.d.f. (conditioned on all possible values of θ /θ rms ) and for p.d.f.s conditioned on either strong expansions either strong compressions.

The heavy right tail of the dense gas p.d.f.s, for c/ c > 1.5, is due to high compression regions characterized by high Γ values, for which the sound speed increases considerably. In fact, the p.d.f. conditioned on strong compressions exhibits a heavy right tail as well, whereas the p.d.f. conditioned on strong expansions attributes lower probability to high values of the speed of sound. On the other hand, in compression regions, low values of c are also very probable when Γ < 0, since in this case the sound speed initially decreases. This explains the higher probability associated with values of c/ c below 0.7 in the case of strong compressions, with respect to strong expansions.

The peculiar behaviour of the speed of sound deeply modifies the p.d.f. of the local Mach number (figure 27). However, the probability of having high local Mach number is higher in dense gases. In figure 28 we report the p.d.f. of the local Mach number (a) and of the normalized vorticity ω/ω rms (b) conditioned on the strong compression and expansion regions at M t,0 = 1. The dense gas and the perfect gas case exhibit a bellshaped form with an exponential right tail, as also found by [START_REF] Moisy | Geometry and clustering of intense structures in isotropic turbulence[END_REF] and [START_REF] Wang | Effect of compressibility on the small-scale structures in isotropic turbulence[END_REF] for a perfect gas with γ = 1.4. Restricting our attention to the conditional p.d.f. we note that for dense gas, unlike the PFG, the p.d.f.s conditioned on strong dilatations and on strong compressions do not exhibit significant differences. In other terms, the probability of producing a given vorticity is similar in both regions. A possible explanation is that, on the one hand, extra vorticity may be generated in strong expansion regions through expansion eddy shocklets; on the other hand, shocklets of any kind are expected to be weaker in dense gases with Γ close to zero and specifically compression shocklets are weaker than in a PFG. For the PFG case, vorticity production in high compression regions is more likely than in high dilatation ones, as observed also by [START_REF] Wang | Effect of compressibility on the small-scale structures in isotropic turbulence[END_REF]. We then infer that, statistically, even if the overall vorticity behaviour is similar for dense and perfect gases, the local production mechanisms may be different.

Conclusion

In the present work, the influence of dense gas effects on the large-scale dynamics of decaying homogeneous isotropic turbulence (CHIT) has been investigated by means of Euler-based numerical simulations assuming the van der Waals thermodynamic model. Specifically, our study has addressed turbulent flows of a BZT gas, for which non-classical nonlinearities, such as expansion shocks, may occur at particular thermodynamic conditions in the transonic and supersonic flow regions. The simulations have been carried out by means of a ninth-order accurate numerical method, validated beforehand against inviscid CHIT results available in the literature for perfect diatomic gases. BZT gases are characterized by specific heat ratios γ very close to one, unlike standard gases such as air. Hence, we have first carried out a preliminary parametric study of the influence of γ for perfect gas CHIT simulations. At high turbulent Mach numbers (M t,0 0.8), the flow thermodynamic properties are found to be highly dependent on γ . Specifically, we find that the r.m.s. of the pressure and sound speed scale respectively with γ and (γ -1)/2, while the r.m.s. of the density is weakly affected by γ . These findings have been shown on theoretical grounds dimensional analysis and momentum conservation. On the contrary, kinematic properties like the kinetic energy and the vorticity are almost insensitive to γ . The influence of the specific heat ratio on local velocity divergence levels is also found to be weak.

The influence of the gas model has been investigated at different initial turbulent Mach numbers considering a polytropic gas characterized by a given constant heat ratio, typical of heavy fluorocarbons. Once again, at high M t,0 , the thermodynamic model has a significant influence on the time evolutions of the average and r.m.s. of the thermodynamic properties, whereas the influence on kinematic properties is smaller. However, the flow dilatational behaviour is deeply different. For a dense gas, the r.m.s. values of the thermodynamic properties show higher amplitudes than for a PFG. The simulations show that the most affected quantity is the speed of sound. In the PFG case, sound speed is nearly constant with time, independently of the initial turbulent Mach number, whereas for the dense gas, it varies over a range of values that increases with M t,0 . Accordingly, the fundamental derivative of gas dynamics varies abruptly from negative to positive values. The peculiar behaviour of the speed of sound for a dense gas has a strong influence on the local Mach number and flow dilatation levels. The most significant differences between the perfect and the dense gas case are found for the repartition of dilatation levels. For PFG, strong compression regions occupy a much larger volume fraction than strong expansion regions. As a consequence, the probability distribution of the velocity divergence is highly skewed toward negative values, even for values of γ typical of dense gases. For the dense gas, the volume fraction occupied by strong compression and expansion regions is much more balanced: strong compression regions are reduced by 30 % and strong expansion regions are increased by 80 % with respect to the perfect gas, and strong expansions and compressions are found to be equally probable. More than a half of the strong expansion regions are characterized by negative values of Γ , suggesting the possibility that expansion shocklets may occur in the dense gas. The presence of expansion shocklets is also suggested by the fact that, for the dense gas, strong expansion regions exhibit a sheet-like structure, rather than the tubular structure observed in the PFG case. In strong compression regions, both the dense and the perfect gas exhibit sheet-like structures, characteristic of eddy shocklets. Furthermore, for PFG the probability of producing a given vorticity is greater in strong compression regions, where compression shocklets possibly occur.

For the dense gas, the production of vorticity in strong compression and expansions is equally probable. This is due to the vorticity generated across expansion eddy shocklets forming in strong expansion regions, and to the fact that shocklets of any kind are expected to be weaker in dense gases with Γ close to zero.

The findings of the present study show that, for flows at sufficiently high turbulent Mach number, compressible homogeneous isotropic turbulence is significantly modified by dense gas effects. Further work is being carried out to elucidate the role of dense gas effects in the small-scale dynamics of viscous homogeneous isotropic turbulence and in wall-bounded flows.

FIGURE 1 .

 1 FIGURE 1. Clapeyron diagram for a BZT Van der Waals gas with γ = 1.0125. (a) The iso-Γ (grey contours and solid lines) are reported as a function of v/v c . (b) The iso-Γ (grey contours) and the iso-c (solid lines) are reported as a function of v/v c . The dashdotted line represents the critical isotherm T/T c = 1. In the figure the dark grey zone represents the inversion zone, corresponding to the region between the transition line and the saturation curve where Γ < 0.

FIGURE 2 .

 2 FIGURE 2. Numerical errors (a) and effective viscosity (b) of the ninth-order scheme as a function of the reduced wavenumber. In (a): --, dissipation function; -• -, normalized phase error.

FIGURE 3 .

 3 FIGURE 3. Comparison between the present results for CHIT decay at M t 0 = 0.2 and χ 0 = 0 using the ninth-order scheme and Jameson's second-order scheme and those of[START_REF] Garnier | On the use of shockcapturing schemes for large-eddy simulation[END_REF] based on a fifth-order modified ENO scheme (MENO) and Jameson's scheme. The grid in use is composed by 128 3 cells. Time histories of the turbulent kinetic energy (a), enstrophy (b), pseudo-Taylor microscale (c) and spectrum of the turbulent kinetic energy at t = 10 (d).

FIGURE 4 .

 4 FIGURE 4. Comparison of the present results for a ninth-order accurate scheme with those provided by Garnier et al. (1999) using a fourth-order ENO scheme, for different choices of the initial conditions. Lines denote the present results, symbols results from the reference.

FIGURE 5 .

 5 FIGURE 5. Influence of grid resolution for a PFG CHIT simulation. Here, the working fluid is a perfect gas with γ = 1.4, and the initial turbulent Mach number is M t,0 = 0.8. (a-c) Show the time histories of the normalized r.m.s. density (a), the pseudo-Taylor microscale (b) and the normalized enstrophy (c); (d,e) show the turbulent kinetic energy spectrum (d) and the p.d.f. of velocity divergence (e) at t = 2.5.

FIGURE 6 .

 6 FIGURE 6. Influence of grid resolution for a DG CHIT simulation. Here, the working fluid is PP1, modelled as a polytropic VDW gas, and the initial turbulent Mach number is M t,0 = 0.8. (a) Shows the time history of the normalized enstrophy; (b,c) show the turbulent kinetic energy spectrum (b) and the p.d.f. of the velocity divergence (c) at t = 2.5.

FIGURE 7 .

 7 FIGURE 7. Comparison of Euler-based and Navier-Stokes-based results for PFG. (a-d) Show the time histories of the normalized r.m.s. density (a), the normalized r.m.s. sound speed (b), the turbulent kinetic energy (c) and the normalized r.m.s. velocity divergence (d) at various Reynolds numbers.

FIGURE 8 .FIGURE 9 .

 89 FIGURE 8. Comparison of Euler-based and Navier-Stokes-based results for PFG simulations. (a) Turbulent kinetic energy spectrum and (b) p.d.f. of the normalized velocity divergence evaluated at t = 2.5 at various Reynolds numbers.

FIGURE 10 .

 10 FIGURE 10. Influence of the artificial dissipation coefficients k 2 and k 10 (256 3 grid, M t,0 = 0.8, perfect gas with γ = 1.4). (a,b) Show the time histories of the normalized turbulent kinetic energy (a), the normalized enstrophy (b). (c,d) Show the turbulent kinetic energy spectrum (c) and the p.d.f. of velocity divergence (d) at t = 2.5.

  FIGURE 11. PFG simulations: influence of the specific heat ratio. Time evolution of the r.m.s. of the sound speed (a), of the pressure (b) and of the density (c). p rms and c rms are rescaled according to (4.2) and (4.4), respectively.

FIGURE 12 .

 12 FIGURE 12. PFG simulations: influence of the specific heat ratio. Time evolution of the turbulent Mach number (a), turbulent kinetic energy (b), normalized r.m.s. of the vorticity ω rms /ω 0 (c), normalized r.m.s. of the dilatation θ rms /ω 0 (d), kinetic energy spectrum (e) and p.d.f. of the local Mach number at t = 2.5 ( f ).

FIGURE 13 .

 13 FIGURE 13. Representation in the p-v plane of the initial thermodynamic conditions used for the dense gas simulations using the VDW model. I.C. represents the initial thermodynamic state; the dashed line s = s i represents the isentropic curve containing the initial thermodynamic state; T = T c the critical isotherm, the dark grey region bounded by Γ = 0 (the transition line) is the inversion zone; the grey region to the right of the Γ = 1 is the dense gas region.

  FIGURE 15. PFG and DG simulations. Time evolution of the r.m.s. density (a), turbulent Mach number (b), average speed of sound (c) and r.m.s. speed of sound (d).----, M t,0 = 0.2; -----, M t,0 = 0.5; --, M t,0 = 0.8; ---, M t,0 = 1. E , VDW; A , PFG.

  FIGURE 16. PFG and DG simulations. Time histories of the average (a) and r.m.s. (b) fundamental derivative of gas dynamics for the dense gas at various initial turbulent Mach numbers.----, M t,0 = 0.2; -----, M t,0 = 0.5; --, M t,0 = 0.8; ---, M t,0 = 1. E , VDW.

FIGURE 17 .

 17 FIGURE 17. Representation in the p-v diagram of the thermodynamic states at t = 2.5. Each symbol represents the thermodynamic conditions computed at a given mesh points at various initial Mach numbers. The spread of thermodynamic states (approximately distributed along an isentropic curve) increases with the Mach number, due to compressibility effects.

FIGURE 18

 18 FIGURE 18. (Colour online) Instantaneous isocontours (t = 2.5) of the fundamental derivative of gas dynamics (a), the normalized speed of sound c/c 0 (b), and normalized density ρ/ρ 0 (c) along a mesh plane for a dense gas computation at M t,0 = 1. At high Mach number dense gas computations, these thermodynamic properties exhibit significant large-scale fluctuations. Specifically, Γ can take both positive and negative values. In negative regions, non-classical nonlinearities may appear.

  FIGURE 19. PFG and DG simulations. Time evolution of the r.m.s. divergence of the velocity field normalized with the initial vorticity at various initial turbulent Mach numbers.----, M t,0 = 0.2; -----, M t,0 = 0.5; --, M t,0 = 0.8; ---, M t,0 = 1. E , VDW; A , PFG.

  FIGURE 20. (Colour online) Instantaneous iso-contours of the normalized velocity divergence θ/θ rms at various times, for M t 0 = 1. Results for the PFG case are represented in (a), dense gas in (b). Red regions correspond to strong compressions (θ/θ rms < -3), green regions to strong expansions (θ /θ rms > 3).

  FIGURE 21. (Colour online) Snapshots at t = 2.5 of the isosurfaces θ /θ rms = 3, coloured by the local Mach number, for perfect and dense gas simulations at M t,0 = 1.

  FIGURE 23. P.d.f. of the normalized velocity divergence scaled by M 2 t,0 at various initial turbulent Mach numbers (and t = 2.5). The vertical dashed lines indicate the values θ /θ rms = ±3. E , VDW; A , PFG.

FIGURE 24 .

 24 FIGURE 24. Distribution of the conditioned mean of the fundamental derivative of gas dynamics Γ on θ /θ rms at various M t,0 .

  FIGURE 25. (Colour online) Contours of the joint probability density function (log scale) of the fundamental derivative of gas dynamics Γ with the scaled normalized local velocity divergence θ/θ rms at various initial turbulent Mach numbers.

FIGURE 26 .

 26 FIGURE 26. Probability density function of the normalized sound speed conditioned on θ/θ rms at M t,0 = 1. ----, [-∞, -2]; ---, [2, ∞]; --, [-∞, ∞]. E , VDW; A , PFG.

FIGURE 27 .

 27 FIGURE 27. Probability density functions of the local Mach number at various initial turbulent Mach numbers at t = 2.5. VDW, ( ); PFG, ( ).

FIGURE 28 .

 28 FIGURE 28. Total probability density function of the local Mach number (a) and the normalized local vorticity (b) and p.d.f. conditioned on strong compressions (θ /θ rms < -2) and strong expansions (θ /θ rms > 2) for M t,0 = 1 at time t = 2.5. ----, [-∞, -2]; ---, [2, ∞]; --, [-∞, ∞]. E , VDW; A , PFG.

TABLE 1 .

 1 Thermodynamic properties of PP11: molecular weight (M ), critical temperature (T c

	14 F 24 )	624.11	650.2 627.14 1.46 0.2688 0.4963 0.0 488.15	97.3

TABLE 2 .

 2 Cases considered to study the influence of the artificial dissipation coefficients.

TABLE 3 .

 3 PFG simulations: influence of the specific heat ratio. Values of the specific heat ratio (γ ), specific heat coefficient at constant volume (c v ) and at constant pressure (c p ) and fundamental derivative of gas dynamics Γ .

TABLE 5 .

 5 Initial thermodynamic conditions used for the dense gas simulations using the PFG and VDW equations of state.

TABLE 6 .

 6 Volume fraction occupied by compression and expansion regions at time t = 2.5.

	+∞]

Acknowledgement

High performance computing resources were allocated by the GENCI French supercomputing center under project GENCI-ldf7332.

Appendix A. Characteristic turbulent quantities

In this appendix we recall some useful relations for the analysis of compressible homogeneous isotropic turbulence. Since the mean velocity field is zero, we define the turbulent fluctuating velocity as u = u = u 2 i /3 1/2 , where the operator • refers to a volume average over the computational domain at a fixed time instant and the index i = 1, 2, 3 represents the three Cartesian coordinates. The r.m.s. values are computed as η rms = η 2η 2 , where η is a generic fluctuating quantity. The integral length scale L I is defined as:

where k = k is the wavenumber and E(k) is the spectrum of the turbulent velocity u i , integrated over shells of radius k. The Taylor microscale λ is defined as in Hinze (1975)

According to [START_REF] Jiménez | The structure of intense vorticity in isotropic turbulence[END_REF], the Taylor microscale can be computed directly from integral quantities, by using the relation:

where E = ∞ 0 E(k) dk is the integrated energy and Ω = ∞ 0 k 2 E(k) dk is the enstrophy. The turbulent Mach number M t and the r.m.s. Mach number M rms are defined as:

In the definition of M rms in (A 4) it is assumed M loc = 0 since the mean flow is zero and only the fluctuating components of the kinetic energy are computed.

In the spectral space, the kth harmonic of the turbulent velocity u(k) can be properly decomposed in a divergence-free component, u s (k), and a dilatational component, u c (k). The latter is computed as

Considering a discrete representation of the spectrum over N uniformly spaced grid points, the total and compressible kinetic energies K and K c are computed as follows:

where n k is the number of Fourier modes in the kth bin, satisfying

and the spectra of the total kinetic energy and its dilatational component are computed as, respectively:

(A 7a,b)

Finally, the compressibility ratio χ is defined as χ = K c /K.