DNS of turbulent flows of dense gases

The influence of dense gas effects on compressible turbulence is investigated by means of numerical simulations of the decay of compressible homogeneous isotropic turbulence (CHIT) and of supersonic turbulent flows through a plane channel (TCF). For both configurations, a parametric study on the Mach and Reynolds numbers is carried out. The dense gas considered in these parametric studies is PP11, a heavy fluorocarbon. The results are systematically compared to those obtained for a diatomic perfect gas (air). In our computations, the thermodynamic behaviour of the dense gases is modelled by means of the Martin-Hou equation of state. For CHIT cases, initial turbulent Mach numbers up to 1 are analyzed using mesh resolutions up to 5123. For TCF, bulk Mach numbers up to 3 and bulk Reynolds numbers up to 12000 are investigated. Average profiles of the thermodynamic quantities exhibit significant differences with respect to perfect-gas solutions for both of the configurations. For high-Mach CHIT, compressible structures are modified with respect to air, with weaker eddy shocklets and stronger expansions. In TCF, the velocity profiles of dense gas flows are much less sensitive to the Mach number and collapse reasonably well in the logarithmic region without any special need for compressible scalings, unlike the case of air, and the overall flow behaviour is midway between that of a variable-property liquid and that of a gas.


Introduction
The physical understanding of the turbulence dynamics of dense gas flows is relevant for several engineering systems, like high-Reynolds wind tunnels [1], chemical transport and processing [2], energy conversion cycles [3] and refrigeration [4]. Dense gases (DG) are single-phase fluids with a molecular complexity such that the fundamental derivative of gas dynamics [5] Γ := 1 + ρ c ∂c ∂ρ s (where ρ is the density, p the pressure, s the entropy, and c the sound speed), which measures the rate of change of the sound speed in isentropic transformations, is less than one or even negative in a range of thermodynamic conditions close to the saturation curve. In such conditions, the speed of sound increases in isentropic expansions and decreases in isentropic compressions, unlike the case of perfect gases. Dense gas effects are expected to give the most interesting phenomena for a family of heavy polyatomic compounds, named Bethe-Zel'dovich-Thompson (BZT) fluids [5], which exhibit a region of negative Γ values (called the "inversion zone" [6]) in the vapour phase close to the liquid/vapour coexistence curve. This is theoretically predicted to result in nonclassical compressibility effects in the transonic and supersonic flow regimes, like rarefaction shock waves, mixed shock/fan waves, and shock splitting [6,7]. Furthermore, the entropy change across a weak shock is much weaker than usual for dense gases with Γ << 1, leading to reduced shock losses. Numerical simulations of turbulent dense gas flows of engineering interest are based on the well-known RANS equations, which need to be supplemented by a model for the Reynolds stress  Table 1: Thermophysical properties of PP11: molecular weight (M), critical temperature (T c ), density (ρ c ) and pressure (p c ), critical compressibility factor (Z c ), acentric factor (ω), dipole moment (d.m.) of the gas phase, boiling temperature (T b ), ratio of ideal-gas specific heat at constant volume over the gas constant (c v (T c )/R) at the critical point and power-law exponent for the low-density specific heat (n).
tensor. The accuracy of RANS models for dense-gas flows has not been properly assessed up to date, due to the lack of both experimental and numerical reference data. The creation of reliable and thorough DNS databases is then needed to quantify the deficiencies of existing turbulence models and to develop and calibrate improved models. Additionally, DNS may improve our understanding of modifications to compressible turbulence due to the fluid in use. Indeed, for dense gases, the perfect gas (PFG) model is no longer valid, and more complex equations of state (EoS) must be used to account for their peculiar thermodynamic behaviour. Moreover, in the dense gas regime, the dynamic viscosity µ and the thermal conductivity κ depend on temperature and pressure through complex relationships. Similarly, the approximation of nearly constant Prandtl number (P r = µc p /κ ≈ const) is no longer valid. In the following, we present recent DNS results for the decay of compressible homogeneous isotropic turbulence (CHIT), and of turbulent plane channel flow (TCF), representative of free and wall-bounded turbulence, respectively, and we highlight the most striking differences with respect to perfect-gas turbulence.

Governing Equations and Numerical Method
In the present study we consider flows of gases in the single-phase regime, governed by the compressible Navier-Stokes equations. The gas behaviour is modelled through the Martin-Hou (MAH) thermal equation of state [8], which provides a realistic description of the gas behaviour and of the extent of the inversion zone. This equation, involving five virial terms and satisfying ten thermodynamic constraints, is reasonably accurate for the fluid of interest and requires a minimum amount of experimental information. A power law of the form c v∞ = c v∞ (T c )(T /T c ) n is used to model variations of the low-density specific heat with temperature, where n is a material-dependent parameter. In addition to the EoS, thermodynamic models relating the dynamic viscosity µ and thermal conductivity κ to the gas temperature and pressure have to be specified. In the following calculations, the transport properties follow the Chung-Lee-Starling model [9], which incorporates a correction term in the dense-gas region. For air cases, the classical polytropic perfect gas model is used along with a temperature-dependent power law for the viscosity, µ = µ ref (T /T ref ) 0.7 , and a constant Prandtl assumption, P r = 0.7. The working fluid is perfluoro-perhydrophenanthrene (chemical formula C 14 F 24 ), called hereafter with its commercial name PP11; its thermodynamic properties are provided in table 1. This fluid has been often used in the dense-gas literature since it exhibits a relatively wide inversion zone and, as a consequence, may lead to the appearance of non-classical effects. The space derivatives in the governing equations are approximated by means of optimized finite difference schemes. Specifically, the convective terms are approximated by fourth-order optimized finite differences with an eleven-point stencil in each mesh direction [10], and the viscous terms by standard fourth-order finite differences. As a part of the algorithm, an optimized selective sixth-order filter also using an eleven-point stencil is applied in each direction to damp grid-to-grid oscillations in smooth flow regions. In order to ensure computational robustness for compressible flows with shocks, an adaptive nonlinear selective filtering strategy is employed [10]. Finally, time integration is carried out by means of a low-storage six-step Runge-Kutta scheme optimized in the wave number space [11].

Decaying homogeneous isotropic turbulence
The CHIT decay problem is solved on a cubic domain of [0, 2π] 3 . Periodic boundary conditions are applied in the three Cartesian directions. The initial turbulent kinetic energy spectrum is of Passot-Pouquet type, i.e., E(k) = Ak 4 exp −2 (k/k 0 ) 2 , where k 0 is the peak wavenumber and A is a constant that depends on the initial amount of kinetic energy. The initial velocity field is purely solenoidal and the r.m.s. of the thermodynamic variables are set to zero. Compressible turbulence decay is governed by the turbulent Mach number M t = u 2 i /c and the Reynolds number Re λ = u λρ/µ, based on the Taylor microscale λ. We performed direct simulations of dense-gas CHIT at various initial turbulent Mach numbers, namely M t 0 = 0.2, 0.5, 0.8 and 1. The initial Reynolds number is set to Re λ = 200 and the initial peak wavenumber k 0 is equal to 2 for all simulations. The decay is simulated up to t = 4 (with t a non-dimensional time scale based on the initial large eddy turnover time). All simulations were carried out on a regular Cartesian grid of 512 3 . In the most severe case (i.e. M t 0 = 1), the smallest resolved wavelength k max η (being k max the largest wavenumber and η = (µ 3 /ρ 2 ) 1/4 the Kolmogorov lengthscale [12]) is always greater than ≈ 6.5 − 7, and the Kolmogorov lengthscale is discretized by at least 2 mesh points, which corresponds to a very good resolution. Dense gas simulations are reported for two choices of the initial thermodynamic state, corresponding to a point in the dense gas region (ρ 0 /ρ c = 1.618, p 0 /p c = 1.02, T 0 /T c = 1.01, Γ 0 = 0.1252, referred-to as IC1) and one inside the inversion zone (ρ 0 /ρ c = 1.618, p 0 /p c = 0.98, T 0 /T c = 1.001, Γ 0 = −0.093, IC2), for which the flow may exhibit BZT effects. Since entropy can only increase during the flow decay, no BZT effects can appear for the first choice of conditions. Results were systematically compared with those of a diatomic perfect gas (air). Compressibility effects are first investigated by comparing the time evolutions of relevant flow statistics for PFG and dense-gas simulations at various M t 0 (figure 1a). The PFG and DG solutions are practically superposed at M t 0 = 0.2 and progressively deviate as M t 0 increases, due to the different compressibility and transport properties of the two fluids. For the sake of brevity, in the following we restrict our analysis to the most representative case, i.e., M t 0 = 1. The average temperature (figure 1b) increases significantly in time for PFG, and viscosity (figure 1c) follows the same behavior. On the contrary, for DG cases the temperature remains approximately constant due to the much higher specific heat, and the time evolution of viscosity (which depends on density variations) exhibits a peak at the beginning of the decay and decays afterwards, i.e. behaves at the opposite of the PFG. The normalized turbulent kinetic energy K/K 0 (not reported for brevity) and enstrophy Ω/Ω 0  (figure 1d) are similar for PFG and DG cases. An enstrophy peak is observed for t ≈ 1.4 − 1.5, which is slightly higher for the dense gas, due to the lower average viscosity compared to PFG. Similar behaviours are observed for IC1 and IC2. Time histories of the fluctuating thermodynamic quantities are presented in figure 2 (top line). The maximum to minimum instantaneous density ratio increases during the initial supersonic phase, driven by vortex-stretching mechanisms with the production of strong velocity gradients, reaching a maximum for a nondimensional time comprised between 1 and 2. In the dense gas, the maximum value is however about 2.5 times lower than in the perfect gas due to the formation of weaker eddy shocklets. Then the ratio decays and both fluids exhibit similar behaviours at later times. Similarly, pressure fluctuations (figure 2b) in the dense gas are about five times lower for PP11, independently on the choice of the initial thermodynamic state. Conversely, the r.m.s. values of the speed of sound are higher for the dense gas. For this fluid, c = c(ρ, T ) follows essentially the density variations (the temperature being almost constant throughout the decay). This mechanism is deeply different from the one observed for PFG, where the speed of sound increases due to friction heating. Time evolutions of the average and r.m.s fundamental derivative of gas dynamics are also reported for PP11 at various initial conditions. The average value Γ (figure 2d) follows the average density variations. For both IC1 and IC2, the initial development of compressible modes increases the average Γ, which however remains below 1 throughout the evolution (implying that the speed of sound tends to decrease in compressed regions) and then decays. The r.m.s. of Γ is as high as 1 at t ≈ 1 meaning that, for high Mach number flows, the pointwise values of Γ are widely spread around the average and may be locally greater than 1 or, for IC2, lower than 0. Despite the different local behaviour of Γ, the macroscopic flow statistics are little affected and tend to behave similarly for both initial conditions. Figure 3 shows the distribution of the flow thermodynamic states in the p − v diagram at time t = 2, for IC1 at various M t 0 and for IC2 at M t 0 = 1. Viscosity and shock losses being small, all states tend to distribute along the initial isentrope. The thermodynamic region spanned by the p − v states throughout the decay increases with M t 0 . For M t 0 = 1, the slight change in the initial thermodynamic conditions drastically changes the maximum pressure achieved. The specific volume, indeed, is bounded by a lower limit corresponding to the liquid/gas transition (respectively, an upper limit for density), as the flow undergoes strong compressions. On the other hand, strong expansions are favoured by the dense gas behaviour, since the speed of sound is increased and compressibility decreases. The dilatation regions are classified according to the local level of normalized dilatation θ/ √ θ 2 M 2 t 0 [13]. Regions with 0 < |θ/ √ θ 2 M 2 t 0 | < 2 are referred to as weak expansions/compressions, whereas |θ/ √ θ 2 M 2 t 0 | > 2 denotes strong dilatation/compression motions. Additionally, it is commonly acknowledged that compression shocks may be formed for θ/ √ θ 2 M 2 t 0 < −3 [14]. The differences in the Clapeyron diagram are confirmed by the probability distribution functions (p.d.f.s) of the normalized dilatation shown in figure 2f. The dense gas exhibits a much more symmetric p.d.f. than PFG, which is characterized by a long left tail, according to similar results in the literature [15]. The greater balance between the left and right tails of the p.d.f. in the dense gas are due to the significant weakening or suppression of compression shocks in flow regions where Γ is close or less than zero; on the other hand, expansions are favoured, leading to a heavier right tail. This is more developed for IC2, due to the possibility of BZT effects leading to the appearance of expansion shocklets. At t = 2, approximately one half of the flow volume is characterized by thermodynamic states enclosed under the transition line. Among these states, approximately 36% consist of weak dilatation regions and ≈ 6% of strong expansion regions. Globally, strong compression regions decrease from 3% to 2% whereas strong expansions increase from 1.7% to 2.6%. More details about dense gas CHIT and its small-scale behaviour are given in [16].

Supersonic Turbulent channel flow
The TCF configuration is now considered. Periodic conditions are applied in the streamwise (x) and spanwise (z) directions. Isothermal no-slip wall conditions are applied on the lower and upper walls. In order to counteract viscous friction and maintain the target bulk mass flow, a spatially constant body force is applied in the streamwise direction [17]. In the following, the subscripts (·) B , (·) w and (·) CL denote time and space averaged values over the channel crosssection, at the wall and at the centerline, respectively; (·) indicates Reynolds averaging and (·) Reynolds fluctuations; similarly, (·) and (·) denote Favre averages and fluctuations. Flow conditions are defined by imposing the bulk Reynolds number Re Bw := ρ B u B h µ w , and bulk Mach number M Bw := u B cw . Since isothermal conditions are applied at the walls, T w = T w =const. For the PFG model, fixing the reference wall temperature allows to fix the sound speed and transport properties accordingly. For dense gas cases, the latter quantities depend on both temperature and density, and thus their values change during the simulation because ρ w cannot be fixed a priori. For this reason we used an iterative procedure on the wall temperature in order to match the prescribed Reynolds and Mach numbers. A parametric study is carried out at various bulk Reynolds numbers (Re Bw = 3000, 7000 and 12000) and two bulk Mach numbers (M Bw = 1.5 and 3.0). The computational domain has dimensions L x × L y × L z = 8πh × 2h × 2πh (being h the channel half-height); the selected computational grids ensure a spatial resolution comparable to other DNS of TCF [18,19,20,21,17], i.e., ∆x * ∈ [5,16], ∆y * w ∈ [0.5, 0.8], ∆y * CL ∈ [2, 6], ∆z * ∈ [2,6]. The overall number of grid points varies between 33 · 10 7 ÷ 1.2 · 10 9 for the various cases. The main numerical parameters and global DNS results are summarized in table 2. Starred quantities refer to the empirical semi-local scaling initially proposed by Huang et al [19] for compressible flows, which corrects the usual wall scaling with centerline quantities.
Specifically, y * = ρ(y)u * τ y µ (y) and Re τ * = Re τ ρ CL ρ w µ w µ CL , with u * τ := τ w ρ(y) the semi-local friction velocity. This mixed scaling (recently discussed in [22]) has proven to give quite satisfactory results in collapsing first-and second-order moments [21] obtained from a wide range of M Bw . We recall that in case of an incompressible flow, Re τ * = Re τ .  Despite considerable improvement is achieved over the classical incompressible scaling (with friction velocity u τ ), the profiles depart from the incompressible log law (ln y + )/0.4 + 5.5 when M Bw increases. On the contrary, quite satisfactory results are obtained for PP11 at any Mach number. In some sense, the dense gas behaves similarly to a liquid with variable properties. For air, the average temperature T + grows rapidly with M Bw , due to the significant increase of friction heating, especially for the lowest Re Bw (panel b). For PP11, due to the large specific heat of the fluid T + is almost constant across the channel for any choice of the Mach and Reynolds numbers (panel f), and the centerline temperature differs less than 1% from T w . Decoupling of dynamic and thermal effects in the dense gas leads to smaller density variations across the channel. For air, the centerline density is up to 60% lower than ρ w (at M Bw = 3) whereas variations below 20% are observed for PP11. Panels c and d show profiles of the normalized viscosity µ/µ w . This quantity follows temperature variations for air, so that µ CL /µ w ≈ 2 for case AM3R3. For PP11, the viscosity exhibits a liquid-like behaviour at the considered thermodynamic conditions and takes lower values at the centerline. Unlike liquids, however, this is not due to a temperature rise, but to the reduced density. The preceding differences explain the opposite behaviour of the friction Reynolds number for the two fluids (panels d and h). For air, Re τ * decreases rapidly up to the buffer layer and more slowly up to centerline. For a given Re Bw , higher Mach numbers enhance temperature gradients leading to the drop of Re τ * in the inner layer. For PP11, Re τ * increases monotonically up to the buffer layer; afterwards, its value is approximately constant. Dense gas effects become more evident at higher M Bw , since the local thermodynamic states spread over a wider range. In general, Re τ * < Re τ for air (gas-like behaviour) and Re τ * > Re τ for PP11 (liquid-like behaviour). In the wall region, Re τ * is ≈ 20% lower in PP11 (leading to a smaller friction coefficient), whereas it is up to two times higher at centerline, that is then characterized by smaller and more erratic turbulent structures. In spite of the liquid-like behaviour of some properties in the dense gas, compressibility effects are not suppressed. On the contrary, the centerline Mach number is 20 to 30% greater in the dense gas than in air, where friction heating reduces the local Mach number with respect to the bulk one.
The distribution of the p − v thermodynamic states for case PM3R12 is reported in figure 5 as a function of the distance to the wall in semi-local units. The diagram is coloured with the local Prandtl number. Moving from the wall to the centerline, a smaller thermodynamic region is spanned, while the average entropy tends to increase. We observe that, since the flow never crosses the transition line, no BZT effects are involved. Additionally, the maximum turbulent Mach number, approximately 0.5, is low enough to neglect the influence of compressibility on turbulence. On the other hand, the Prandtl number undergoes significant variations, especially in the near-wall region, which makes the use of a constant-Prandtl assumption quite inaccurate to predict dense gas flows characterized by a near-critical wall temperature. Figure 6 displays the the r.m.s. values of density, the Reynolds stresses and shear stresses ρu i u j + = τ −1 w ρu i u j in semi-local scaling and the ratio of the turbulent kinetic energy production to dissipation. Profiles of the r.m.s. density for the dense and the perfect gas (panels a and e) are deeply different. Previous reference results [17] have shown that for air cases,  fluctuations grow more quickly with M B,w in PP11 than in air. Nevertheless, Morkovin's hypotesis remains quite well satisfied, even at the highest Mach number. Despite the striking differences in the thermodynamic behaviour, the Reynolds stress profiles are similar for both fluids. Specifically, the dense gas profiles agree with those observed in other studies [22] for low-Mach TCF with temperature-dependent transport properties. The liquid-like behaviour of the viscosity leads to an increase of the spanwise, wall-normal and shear Reynolds stresses with respect to the corresponding incompressible evolution, whereas the streamwise one decreases, and this effect tends to be reinforced when increasing the Mach number (panels f and g). An opposite behaviour is observed for air (panels b and c). For instance, case AM1R7 (Re τ * ≈ 315) is characterized by a semi-local friction Reynolds number higher than case AM3R7 (Re τ * ≈ 200). Nevertheless, being the temperature-dependent viscosity much higher throughout the channel in the latter case, ρv v , ρw w and ρu v are lower with respect to AM1R7. Finally, panels d and h show the profile of the ratio of production to dissipation term P k / k of the turbulent kinetic energy budget. The semi-local scaling, built with τ 2 w /µ(y), behaves better than both the bulk scaling [20] τ w u B /h and the wall scaling ρ w u 4 τ /ν w . Previous works [24] showed that the main effect of compressibility is a reduction of both production and dissipation; this is confirmed for both air and PP11. All air and PP11 cases exhibit the usual production peak in the inner region, located approximately at y * ≈ 12. For the higher Reynolds PP11 case, a second production peak is observed in the outer region, like in high-Re incompressible flow, which is not the case for air at the same Re Bw . This effect is due to the liquid-like behaviour of viscosity in PP11, leading to reduced dissipation when approaching the centerline.

Conclusions
Direct numerical simulations of the decay of compressible homogeneous isotropic turbulence (CHIT) and plane turbulent channel flow (TCF) of dense gases have been performed. Parametric studies with respect to the characteristic Mach and Reynolds numbers have been carried out for both configurations. The results were systematically compared to those obtained for air. In CHIT simulations, important differences are observed for flows with a high initial turbulent Mach number. Different trend for the r.m.s. values of the thermodynamic quantities and transport properties greatly modify the turbulence evolution. In the dense gas, compressibility effects, and specifically eddy shocklets, are weaker. A tendency to symmetrization is observed for the p.d.f. of the normalized velocity divergence, due to the weakening of compression shocklets and the sharpening of rarefaction waves, leading possibly to the formation of rarefaction eddy shocklets for BZT fluids at suitable operating conditions. For TCF, the classical y + scaling based on the friction velocity fails to collapse thermodynamic profiles and Reynolds stresses when high M Bw are considered and corrected scalings accounting for variations of the fluid properties have to be considered. The semi-local scaling, based on a mix of wall and local thermodynamic quantities, is shown to be the most suitable choice for collapsing Reynolds stresses and energy budgets. Nevertheless, the maximum turbulent Mach number is approximately 0.5 despite the high bulk Mach number, and turbulence structure is indeed little affected by compressibility effects. Additionally, the coupling between dynamic and thermal effects is very small for dense fluid characterized by high specific heats, contrary to air which undergoes significant friction heating (typical of highly compressible flows). For these reasons, for PP11 many of the quantities of interest follow an evolution similar to that observed in incompressible turbulence. For the reference thermodynamic conditions adopted, the transport properties have a liquid-like behavior and the actual Reynolds numbers are found to be much higher with respect to air. Of course, compressibility may become important if higher bulk Reynolds or Mach number are considered.
The numerical simulations reported in the paper have been performed for PP11, but preliminary simulations carried out for other perfluorocarbons (i.e. PP9 and PP10) have shown a similar behavior. Analyses are also underway for siloxanes (D5 and D6) and refrigerants (R245fa and R134a). Also note that simulations for a BZT van der Waals gas (not reported) exhibit a similar qualitative behavior. This suggests that, although the magnitude of the observed dense gas effects depends on the extent of the dense gas region for a given fluid and on the presence (or not) of an inversion zone, the qualitative behavior is expected to be similar for other dense gases. Results with other working fluids will be presented in a forthcoming paper.