Shock-induced energy transfers in dense gases

Dense gases are characterised by molecules featuring large numbers of active degrees of freedom (quantified by the cv/R ratio). The isentropes in such gases have the distinct property of following rather closely the isotherms (the two become identical in the limit of cv/R going to infinity). Near the liquid-vapour critical point, this makes the isentropes very shallow and possibly concave (in the pressure-specific volume diagram). Whilst shallow isentropes are desirable when designing expanders (i.e. a large specific-volume increase may be achieved for virtually no pressure drop), could such extreme compressibility effects modify turbulence in a profound manner? This paper discusses two particularly interesting aspects: (i) shock-refraction properties (i.e. the way a shock can redistribute the energy of incoming perturbations), (ii) enstrophy production in homogeneous turbulence. A linear interaction analysis (LIA) is conducted on the shock configuration for which the incoming perturbation is decomposed into linear modes of the compressible Euler equations. The transmission coefficients relative to each eigen modes are solved analytically and results are compared against fully non-linear compressible direct numerical simulation reproducing the weak perturbation of an isolated two-dimensional compression shock wave. The linear analysis is found to be capable of predicting the shock-induced redistribution of the energy of the incoming perturbation between the different eigen modes. Non-ideal gas effects are observed both analytically and numerically with especially an unusual selective response for some particular choice of incoming Mach number. A two-dimensional isotropic turbulence configuration is then numerically investigated for the case of an inviscid compressible dense-gas flow close to the liquid-vapour critical point. Strong non-ideal-gas effects on enstrophy production are observed with the formation of eddy shocklets. In both cases non-convex isentropes close to the liquid-vapour critical point are extremely influential in letting both the shock and the turbulence redistribute any supply of turbulence kinetic energy in ways which are simply not observable in ideal gases. This will hopefully spark enthusiasm amongst turbulence modellers (and their end users?).


Introduction
Turbulence is the process by which any supplied amount of energy is redistributed across a broad range of length and velocity scales, constituting the turbulence kinetic energy (TKE). It is well accepted that the TKE cascades from large to small scales, down to scales for which friction forces dominate, ultimately dissipating the TKE into heat. For incompressible flows, this is the complete picture. If in addition fluid parcels can compress (i.e. change volume), the work done by the pressure forces converts some amount of TKE into turbulence internal energy (TIE). In the absence of viscosity, the conversion is isentropic and the TIE can be converted back to TKE. When the velocity fluctuations of the turbulence are of the order of the sound speed, localised shockwaves will appear (called eddy shocklets). Eddy shocklets immediately and irreversibly redistribute the energy across its various forms. Both vorticity and entropy are produced as the result of shocklets and their interactions. The entropy generated represents a leakage of TKE which is no longer recoverable -this is a direct loss. The vorticity production at small scales induced by the shocklets precipitates the cascade and speeds up the TKE dissipation (through friction forces at small scales) -this is an indirect loss.
The ability of eddy shocklets to significantly alter the cascade is conditioned to the actual thermo-physical properties of the fluid. For fluids comprised of molecules solely undergoing elastic collisions, i.e. ideal gases, pressure is linearly related to density and temperature. Thus, density fluctuations necessarily scale with that of pressure. Some fluids of practical interests (in the area of energy production, aerodynamics and propulsion) are made of molecules undergoing long-range interactions and therefore exhibit inelastic collisions. A particular occurrence is near the gas-liquid thermodynamic critical point (TCP), which is characterised by the presence of a saddle point on the critical isotherm (in the pressure/specific-volume phase diagram), a feature which is impossible to achieve in ideal gases. The presence of a saddle point on the isotherms is foot-printed on the isentropes (the two coincide for molecules with an infinite number of active degrees of freedom). Thus the isentropes can be made arbitrarily shallow near the TCP. Shallow isentropes translate in low sound speeds and mechanically promote the formation of eddy shocklets. Thus, for a given amount of energy, TKE transfers near the TCP will exhibit exacerbated compressibility effects relative to an ideal gas.
This paper aims at exploring some of the mechanisms by which shocks and shocklets affect the evolution of the turbulence kinetic energy when the fluid is governed by locally-concave isentropes. Only elementary interactions are considered: (i) the refraction properties of a compression shock following the impingement by an entropy (density) perturbation, (ii) the production of vorticity (and enstrophy) in two-dimensional turbulence.

Compression-shock refraction properties 2.1. Admissible compression shocks
Discontinuous solutions to the Euler equations are found by solving the well-known Rankine-Hugoniot (RH) algebraic system. The solutions are conveniently seen in the pressure/specificvolume diagram as the intersection points between the Rayleigh line (R-line, along which mass is conserved and momentum balanced by the pressure forces) and the Hugoniot line (H-line, along which changes in internal energy are balanced by the work done by the averaged pressure), therein joining all reachable thermodynamic states across the discontinuity: Notations (.) 1 and (.) 2 refer to pre-and post-shock quantities (in a reference frame attached to the shock), respectively, with p the pressure, ϑ the specific volume, and j = ρu the mass flow rate (u being the fluid velocity normal to the shock). In classical gas dynamics (e.g. when an ideal-gas equation of state is assumed) the H-lines are convex and only two intersection points are at most solutions to equation (1). Admissible discontinuities are selected based on the classical Lax condition [1]. Dense gases modelled with a non-ideal equation of state can have non-convex H-lines close to the thermodynamic critical point and equation (1) admits more

Linear perturbations in one dimension
Linear perturbations of a uniform solution to the Euler equations can be projected on three different modes: (1) an entropy mode (pure temperature and density disturbance fields); (2) two irrotational acoustic modes and; (3) a solenoidal vorticity mode [4]. In the linear regime all three modes propagate independently and may interact only at boundaries such as shockwaves.
The ability of compression shocks to transfer energy from one mode to another is investigated close to the TCP.

Analytical results
An entropy mode is injected on the supersonic side of a compression shock close to the TCP. Perturbations normal to the shock front are considered. Thus, no vortical mode may contribute and our focus is on energy transfers between the entropy and acoustic modes at the shock.
Reference thermodynamic values are taken at the TCP (p c , ρ c , T c ), providing a velocity scale (u = p c /ρ c ) and allowing the Euler equations to be written in non-dimensional form. The reference length scale (λ ) is prescribed by the perturbation. Dimensionless quantities are respectively for the entropy mode, noted with (.) s , and the upstream/downstream propagating acoustic modes, noted with (.) a + /(.) a − . Harmonic disturbances of pulsation ω are thus expressed as: with ξ j and k j the corresponding eigen-modes amplitudes and wave numbers defined as σ j = ω/k j with j ∈ {s, a + , a − }. In equation (3) the variable x is the space coordinate in the direction aligned with the mean flow.
The particular choice of equation of state (EOS) appears in equation (2) through the sound speed (noted c) which depends on the thermodynamic variables. Figure 2 (left) illustrates the convergence of the critical isentrope towards the critical isotherm in the (p, ϑ)-diagram for fluids composed of molecules featuring a large number of degrees of freedom [5]. Thus, a significant decrease of the sound speed is expected near the isotherm saddle point. This may come with a degeneracy of Kovasznay's eigen-mode basis in which case a strong bias towards density fluctuation is expected, see figure 2 (right). This is a unique feature of non-ideal gas dynamics.
Compression-shock refraction properties are investigated using the aforementioned linear framework. The incoming supersonic pure-entropy disturbance is projected on Kovasznay's eigen-mode basis, giving: where η = x − x s is the space coordinate centred around the mean shock position. Refracted entropy and acoustic disturbances are then produced and projected onto Kovasznay's eigen-mode basis: with similar notations as in equation (3). The refraction properties are discussed through the entropy-mode amplification (χ s = ξ s 2 /ξ s 1 ) and acoustic-mode generation (χ a + = ξ a + 2 /ξ s 1 ) coefficients. Analytical expressions are determined by solving a boundary value problem at the shock front (using the linearised RH conditions). These expression are obtained for a general equation of state: This result will be illustrated here using the van der Waals equation of state close to the TCP, and compared with fully nonlinear numerical simulations.

Numerical method
All numerical results presented in this paper are based on shock-resolved numerical simulations of the Euler equations. The equations are integrated in time explicitly using a 3 rd -order TVD Runge-Kutta scheme. Spatial derivatives are evaluated with a 4 th -order centred finite difference dispersion-relation-preserving scheme [6] using a 13-point stencil (a 4 points per wavelength cut-off is used for the optimisation). Numerical stability is ensured using a centred 13-point 8 thorder explicit discrete filter optimised in spectral space following the methodology developed by Bogey et al. [7] (a 5 points per wavelength cut-off is used). This allows for accurate propagation of large-wavelength content (both from a dispersion and dissipation point of view) whilst minimising small-scale dissipation. Stretching of the grid is performed using a coordinate transform.
Shockwaves are captured using a localised artificial diffusivity method [8]. An additional viscous term is added to the Euler equations where a discontinuity is detected. This method captures shocks on four grid points with pre-and post-shock spurious oscillations reduced to about 1% of the shock strength [9]. Whilst being robust and easy to implement, this formulation does not explicitly set the shock thickness and a reduction of oscillations below a level of one percent requires tedious case-dependant adjustments of the model parameters. For simulations of isolated-shock refraction properties, an alternative and more flexible strategy has been developed to directly control shock-induced spurious oscillations down to a level of 0.01% of the shock strength. This was achieved by explicitly fixing the artificial viscosity magnitude (evaluated from a weak-shock approximation theory [5]) allowing a direct control of the shock resolution (set up to 30 points in the present work).   [10,11], specific to non-ideal gas dynamics, are also shown.

Results
For both pre-shock conditions χ s has a large local maximum with respect to M 1 (corresponding to conditions for which Kovasznay's acoustic eigen-vectors in the post-shock states are mostly biased towards the density component). Furthermore, χ s can be several orders of magnitude higher than that of the equivalent ideal gas. This makes such shockwaves highly selective with respect to small changes in upstream conditions (e.g. the Mach number M 1 ), particularly when the H-line is discontinuous. The analytical results are in very good agreement with numerical simulations of fully nonlinear one dimensional entropy-mode shock perturbation (for wich ξ s 1 = 0.001ρ 1 ).
These peculiar one dimensional shock-refraction properties match the M 1 range for which the two dimensional DK spontaneous emission is expected. The multidimensional shock-perturbation problem is then expected to be far from the classical ideal-gas case.

Refraction of a density pulse
Two-dimensional refraction properties are numerically investigated through the interaction of a density pulse with shocks of varying strengths. Special attention is paid to the shock-selectivity properties found in the one dimensional analysis. Grid stretching at the shock location is applied so as to enforce a large separation of scales between the viscous-shock structure (shock capturing) and the initial pulse width (a ratio of 10 2 is used).
Results for base flows with pre-shock conditions at p 1 = 0.55, T 1 = 1.00 for γ = 1.013 and for which the H-line is discontinuous at M 1 = 1.29 are shown in figure 4. Three different values of upstream Mach numbers (M 1 ) are considered around the discontinuity. The corresponding ideal-gas case (same γ) at M 1 = 1.30 is shown for comparison, and found to be in accordance with the linear interaction analysis (LIA) performed by Fabre et al. [12].
The linear interaction generates a cylindrical acoustic wave, a vorticity (shear) wave through baroclinic effects and a stretched refracted density pulse. The critical LIA angle for which the refracted shear wave reaches a maximum amplitude can be seen in all vorticity fields for the cases presented. A stronger acoustic intensity (shown by the dilatation field) is expected near the shock position due to the additional contribution of the non-propagative acoustic [12]. The generated acoustic wave for the M 1 = 1.29 case has an unusual circular shape due to the nearly sonic post-shock Mach number (M 2 ) close to the H-line discontinuity. The vorticity production can be significantly larger in the non-ideal gas case, relative to the corresponding ideal-gas case (same pre-shock Mach number and thermodynamic conditions).
The strong-selectivity effect of the discontinuous H-lines is made clear when comparing the refracted fields between the three van der Waals cases. Indeed, unlike for ideal gases, it is now possible to channel the incoming-perturbation energy preferably into the acoustic or the vorticity modes for virtually no change in the upstream Mach number. This new degree of freedom and associated selectivity effect should be considered in turbulence models and/or for practical application involving non-ideal gas compressible dynamics of dense vapours close to critical conditions.

Decaying isotropic turbulence
This section reports early results illustrating the importance of shocks in dense-gas turbulence in close vicinity to the critical point. It rests on the assumption that the slowing down of the sound speed in this region of the phase diagram allows for the formation of eddy shocklets, which are effective at redistributing the turbulence kinetic energy. The objective here is not to produce a realistic flow but rather to explore how some of the mechanisms associated with shock-refraction properties fit within a turbulent-like flow. In particular, we are interested in the production of enstrophy. For that purpose, simulations of decaying two-dimensional isotropic turbulence are performed in the absence of viscosity and conductivity apart from those associated with the explicit filter (thereby ensuring the stability of the simulation) and shock capturing.
The flow is initialised with a solenoidal velocity field with energy spectrum E(k) = Ck exp[−(k/k 0 ) 2 /2], where C is a constant to obtain a maximum velocity of U 0 = 1. The pressure fluctuations are computed by solving the Poisson equation for pressure, and added to the chosen pressure p 0 . The temperature is set uniform to T 0 and the density field computed so as to satisfy the thermal equation of state (van der Waals or ideal). The computational domain is a double periodic box of size L × L with L = 2π, discretised with 1920 × 1920 points. This choice of initial condition sets the dilatational velocity to zero and is of interest for studies at constant initial turbulence intensities but varying thermal equations of state. Figure 5 gives the vorticity fields in the early stage of transition to "turbulence" for both the ideal-gas and van der Waals equations of state for a gas with γ = 1.001. The same solenoidal field was used in both cases and explains why large-scale structures in both cases are similar early on in the simulations. This approach highlights the production of many more vortical structures in the van der Waals gas. In the ideal gas, small-scale vortical structures are found to emerge primarily as a consequence of shock-shock interactions. Indeed, shock-shock interactions produce slip lines which develop vortical structures through Kelvin-Helmholtz instabilities. In the van der Waals gas, the exacerbated compressibility of the gas near the critical point gives birth to more shocks (diverse in nature, e.g. expansion shocklets). When curved, these shocks produce significant amounts of vorticity through the baroclinic production term. When impinged by disturbances (e.g. vortex filaments), these shocks are also effective at producing vorticity as demonstrated with the shock/pulse interactions in the previous section. Figure 6 gives the time-evolution of both the space-averaged enstrophy and baroclinic-production term. Given the incompressible initial solution, an overall expansion is expected at the early stage of the simulation. In the absence of the baroclinic-production term and the limit of no dissipation, the space-averaged enstrophy is then expected to decay in time. This is indeed observed for tU 0 /L < 0.1. From tU 0 /L ≈ 0.1, eddy shocklets have formed in the flow. In the ideal-gas case, the shocks fail to deform and produce significant amounts of vorticity. In the van der Waals case however, the emergence of the shocks coincides with a sharp rise of the enstrophy-production term, two orders of magnitude larger than in the ideal gas. Consequently, the space-averaged enstrophy grows in time and exceeds its initial value.
Although limited in scope, this section confirms that dense gases near the critical point can experience an unusual production of vorticity. We will provide more details about the precise mechanisms responsible for the baroclinic-production term in future communications but hope that the present results demonstrate the need to derive new and suitable turbulence models, either for LES or RANS purposes, so as to predict TKE levels in component designs such as ORC expanders.

Summary
This paper explores some of the non-isentropic energy-transfer mechanisms associated with shockwaves in fluids modelled with a non-ideal equation of state. The focus is on dense gases formed by molecules featuring a large number of degrees of freedom when brought near their thermodynamic critical point.
An analysis of increasing complexity is proposed. First, the refraction properties of compression shockwaves are investigated for supersonic (in the shock reference frame) entropy perturbations normal to the shock front. Analytical expressions for the amplitudes of both the refracted entropy and the generated acoustic contents are derived for a general equation of state. A selective response is demonstrated for which significant changes in the refracted content can be achieved through minute changes of the incoming Mach number. This property is exacerbated for discontinuous Hugoniot lines featuring a non admissible part. The aforementioned onedimensional selective response is then confirmed numerically for the two-dimensional linear interaction of a density pulse with an isolated shockwave near a discontinuous adiabat. The ability to redistribute a given supply of internal energy into either acoustic or vorticity modes in a way not possible in ideal gas is elucidated in multiple dimensions. Eventually, these energytransfer mechanisms are investigated in the context of decaying two dimensional "turbulence" in a dense gas modelled with a van der Waals equation of state. The production of enstrophy through baroclinic effects associated with eddy shocklets is found to be significantly larger than for the equivalent ideal gas with same initial turbulence intensity.
This study is a first step towards the understanding of the complex energy transfers in turbulence in non-ideal gases. It sheds light on new non-isentropic energy pathways associated with shockwaves leading to a richer dynamic than in standard ideal-gas turbulence, eventually requiring a revamp of existing turbulence models.