Luca Sciacovelli

Donatella Passiatore

Paola Cinnella

Giuseppe Pascazio

Assessment of a high-order shock-capturing central-difference scheme for hypersonic turbulent flow simulations

High-speed turbulent flows are encountered in most space-related applications (including exploration, tourism and defense fields) and represent a subject of growing interest in the last decades. A major challenge in performing high-fidelity simulations of such flows resides in the stringent requirements for the numerical schemes to be used. These must be robust enough to handle strong, unsteady discontinuities, while ensuring low amounts of intrinsic dissipation in smooth flow regions. Furthermore, the wide range of temporal and spatial active scales leads to concurrent needs for numerical stabilization and accurate representation of the smallest resolved flow scales in cases of under-resolved configurations. In this paper, we present a finite-difference high-order shock-capturing technique based on Jameson's artificial diffusivity methodology. The resulting scheme is ninth-order-accurate far from discontinuities and relies on the addition of artificial dissipation close to large gradient flow regions. The shock detector is slightly revised to enhance its selectivity and avoid spurious activations of the shock-capturing term. A suite of test cases ranging from 1D to 3D configurations (namely, perfect-gas and chemically reacting shock tubes, Shu-Osher problem, isentropic vortex advection, under-expanded jet, compressible Taylor-Green Vortex, supersonic and hypersonic turbulent boundary layers) is analyzed in order to test the capability of the proposed numerical strategy to handle a large variety of problems, ranging from calorically-perfect air to multi-species reactive flows. Results obtained on underresolved grids are also considered to test the applicability of the proposed strategy in the context of implicit Large-Eddy Simulations.

Introduction

Hypersonic flight has gained renewed attention in recent years, due to its importance for multiple breakthrough applications in the defense and military fields, as well as in the areas of spatial tourism and trans-atmospheric flight [START_REF] Leyva | The relentless pursuit of hypersonic flight[END_REF]. The accurate prediction of hypersonic flow fields is a challenging task, due to the massive conversion of kinetic energy from the hypersonic free stream into internal energy as the fluid approaches the body. The shock waves generated in such flight conditions produce a sharp increase in the fluid temperature, possibly causing vibrational excitation and gas dissociation, and resulting in a nonequilibrium thermochemical state. These processes have major effects on aerodynamic performance, heat transfer rates, fluid-surface interaction (e.g., ablation), and hydrodynamic instabilities leading to boundary layer transition and breakdown to turbulence [START_REF] Candler | Rate effects in hypersonic flows[END_REF]. Their accurate prediction is of crucial importance for the design of the thermal protection system and the prediction of the overall force and heat transfer coefficients, and requires advanced numerical solvers and models. In this paper, the focus is put on so-called high-fidelity numerical schemes for Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES), two major enablers for a deeper understanding of out-of-equilibrium flow regions dominated by laminar-to-turbulent transition and turbulent regimes.

A major difficulty in DNS and LES of high-speed flows is the extreme sensitivity of small flow scales to numerical approximation errors. The occurrence of shock waves, with physical thicknesses of the order of a few mean free paths, leads to unfeasible resolution requirements for numerical simulations, at least in the strict DNS sense. On the other hand, velocity fluctuations of the order of the sound speed [START_REF] Lee | Eddy shocklets in decaying compressible turbulence[END_REF] may lead to the formation of eddy shocklets, embedded in the turbulent flow. Usually, both kinds of structures are dealt with by using shock-capturing techniques. These consist in locally injecting controlled amounts of numerical dissipation to spread shocks over a few mesh cells wherever the mesh is not sufficiently fine to resolve the shock thickness. This technique, corresponding to a regularization method, contrasts with the need of using low-dissipation schemes not to alter the fine-scale turbulent motions. These opposite requirements become even more critical as the Mach and Reynolds numbers increase, because of the growing difficulty to distinguish strong gradients due to shocks from those related to turbulent fluctuations. In such a challenging framework, much effort has been done to devise numerical methods able to handle strong shock waves, while ensuring minimal amounts of dissipation elsewhere.

Two great families of discretization method for the non-linear terms in the Euler and Navier-Stokes equations can be distinguished [START_REF] Pirozzoli | Numerical methods for high-speed flows[END_REF]: the first one, originally designed to handle inviscid flows with strong discontinuities, relies on some form of upwinding along with flux or slope limiters ensuring non-linear stability; the second one -in principle more suited for smooth flows -uses central schemes supplemented with some form of filtering or artificial dissipation or, alternatively, ensuring discrete conservation of solution invariants such as the overall kinetic energy. Despite the large number of comparative studies [START_REF] Johnsen | Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves[END_REF][START_REF] Lo | High-order shock capturing schemes for turbulence calculations[END_REF][START_REF] Lusher | Assessment of low-dissipative shock-capturing schemes for transitional and turbulent shock interactions[END_REF], a global consensus on the ''best'' numerical strategy for high-speed turbulent flows has not been reached, since each method proposes a different compromise among concurrent needs: namely, high accuracy, robustness, low computational cost, few tuning parameters and suitability for different configurations.

A well-known drawback of the first class of schemes is the excess of numerical viscosity introduced in the solution, leading to spurious entropy generation and kinetic energy losses in the low Mach number limit [START_REF] Thornber | On entropy generation and dissipation of kinetic energy in high-resolution shock-capturing schemes[END_REF]. Weighted essentially non-oscillatory (WENO) schemes are probably the most popular upwind schemes in the context of LES and DNS of compressible flows. First introduced by Liu et al. [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF] and later improved by Jiang & Shu [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF], they rely on the assembly of high-order numerical fluxes from linear combination of lower-order polynomial reconstructions using suitable weighting coefficients. Many variants have been proposed, e.g., to improve their dispersion and dissipation properties [START_REF] Martín | A bandwidth-optimized WENO scheme for the effective direct numerical simulation of compressible turbulence[END_REF][START_REF] Hu | An adaptive central-upwind weighted essentially non-oscillatory scheme[END_REF] and to reduce the nonlinear dissipation [START_REF] Borges | An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[END_REF][START_REF] Acker | An improved WENO-Z scheme[END_REF]. Coupling with purely central schemes has led to hybrid methods [START_REF] Pirozzoli | Conservative hybrid compact-WENO schemes for shock-turbulence interaction[END_REF][START_REF] Ren | A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws[END_REF] and enhanced weighting strategies and smoothness sensors [START_REF] Van Lith | Embedded WENO: A design strategy to improve existing WENO schemes[END_REF][START_REF] Fu | Targeted ENO schemes with tailored resolution property for hyperbolic conservation laws[END_REF][START_REF] Fu | A very-high-order TENO scheme for all-speed gas dynamics and turbulence[END_REF][START_REF] Zhao | On shock sensors for hybrid compact/WENO schemes[END_REF]. Contrary-wise, the family of central schemes generally introduces very low dissipation, provided that a selective enough numerical filter or artificial viscosity term is used, but is generally limited to compressible flows at moderately supersonic Mach numbers, i.e. with weak shocks. These schemes must be supplemented with selective nonlinear filtering [START_REF] Kim | Adaptive nonlinear artifical dissipation model for Computational Aeroacoustics[END_REF][START_REF] Visbal | Shock capturing using compact-differencing-based methods[END_REF][START_REF] Hixon | Shock-capturing dissipation schemes for high-accuracy Computational Aeroacoustics (CAA) codes[END_REF][START_REF] Bogey | A shock-capturing methodology based on adaptative spatial filtering for high-order non-linear computations[END_REF], artificial diffusive fluxes [START_REF] Kim | Adaptive nonlinear artifical dissipation model for Computational Aeroacoustics[END_REF][START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes[END_REF][START_REF] Tadmor | Convergence of spectral methods for nonlinear conservation laws[END_REF], or localized artificial diffusivity (LAD) under the form of modified transport coefficients [START_REF] Cook | A high-wavenumber viscosity for high-resolution numerical methods[END_REF][START_REF] Cook | Artificial fluid properties for large-eddy simulation of compressible turbulent mixing[END_REF][START_REF] Kawai | Assessment of localized artifical diffusivity scheme for large-eddy simulation of compressible turbulent flows[END_REF] for damping grid-to-grid oscillations in smooth flow regions and to ensure shock capturing. The amount of numerical dissipation introduced at a point of the computational domain is adjusted by means of properly-devised sensors, allowing to switch on shock-capturing capabilities where needed. Shock-capturing highorder central-difference schemes have been successfully applied, for instance, to overexpanded jet flows with shock cells [START_REF] De Cacqueray | Noise of an overexpanded mach 3.3 jet: Nonlinear propagation effects and correlations with flow[END_REF] and highspeed boundary layers of perfect gases up to Mach 6 [START_REF] Franko | Breakdown mechanisms and heat transfer overshoot in hypersonic zero pressure gradient boundary layers[END_REF], as well as to the direct and large-eddy simulations of high-speed flows of singlespecies, molecularly complex gases at thermodynamic conditions close to the liquid/vapor critical point up to Mach 6 [START_REF] Sciacovelli | Dense gas effects in inviscid homogeneous isotropic turbulence[END_REF][START_REF] Sciacovelli | Direct numerical simulations of supersonic turbulent channel flows of dense gases[END_REF][START_REF] Sciacovelli | Small-scale dynamics of dense gas compressible homogeneous isotropic turbulence[END_REF][START_REF] Sciacovelli | Numerical investigation of high-speed turbulent boundary layers of dense gases[END_REF]. However, their suitability for the numerical simulation of severe hypersonic, chemically reacting flows with strong shocks and stiff chemical reactions has not yet been assessed.

For LES, where only the dynamics of the large scales is computed while the effects of sub-grid scales (SGS) are modeled, the choice of the numerical scheme is possibly even more critical than in DNS. Scale separation is indeed difficult to establish since the cut-off between resolved and modeled scales is not sharp and arises from a complex combination of implicit filtering by the grid and the discretization schemes. The intricate interactions between numerical errors and SGS modeling errors has been investigated by numerous authors (a recent discussion can be found in [START_REF] Gloerfelt | Large eddy simulation requirements for the flow over periodic hills[END_REF]), leading once again to two separate modeling strategies: one relying on the explicit introduction of a SGS model, and the other one using the dissipative part of the discretization scheme for ensuring regularization of the unresolved SGS scales [START_REF] Rizzetta | Numerical simulation of supersonic flow over a three dimensional cavity[END_REF][START_REF] Boris | New insights into large eddy simulation[END_REF][START_REF] Mathew | An explicit filtering method for large eddy simulation of compressible flows[END_REF][START_REF] Mathew | A new approach to LES based on explicit filtering[END_REF][START_REF] Bogey | Turbulence and energy budget in a self-preserving round jet: Direct evaluation using large eddy simulation[END_REF]. The latter has become increasingly spread in the scientific community, due to the good tradeoff between computational cost and accuracy offered for a wide range of applications, provided that a high-resolution scheme is used along with sufficiently resolved computational grids [START_REF] Gloerfelt | Large eddy simulation requirements for the flow over periodic hills[END_REF][START_REF] Aubard | Comparison of subgrid-scale viscosity models and selective filtering strategy for large-eddy simulations[END_REF].

The goal of the present study is twofold: (i) to assess the capability of a high-order shock-capturing central scheme, used in our previous works [START_REF] Sciacovelli | Dense gas effects in inviscid homogeneous isotropic turbulence[END_REF][START_REF] Sciacovelli | Numerical investigation of high-speed turbulent boundary layers of dense gases[END_REF], to robustly predict compressible flows with shock waves and chemical nonequilibrium effects while accurately resolving fine-scale turbulent structures; and (ii) to demonstrate the suitability of the non-linear numerical dissipation of the scheme to act as a SGS regularization in under-resolved turbulent flow simulations. The scheme uses tenth-order accurate finite-difference approximations of the non-linear fluxes, supplemented with a higher-order extension of Jameson's adaptive artificial dissipation [START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes[END_REF]. The order of accuracy of the artificial viscosity term is chosen to obtain an overall dissipative-dominant truncation error, which reduces the appearance and amplification of spurious oscillations [START_REF] Thomée | Stability of difference schemes in the maximum-norm[END_REF] and limits the activation of lower-order nonlinear viscosity. The latter is triggered by a highlyselective shock sensor, built on a combination of the original Ducros' extension [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] of Jameson's pressure-based sensor with the Bhagatwala & Lele [START_REF] Bhagatwala | A modified artificial viscosity approach for compressible turbulence simulations[END_REF] modification proposed in the context of LAD methods (more details are given in Section 3). The scheme is applied to a suite of welldocumented test cases of increasing complexity, ranging from 1D and 2D inviscid flow problems to the 3D simulation of a fully turbulent boundary layer at Mach 10 in chemical nonequilibrium conditions. The results are systematically assessed against exact or numerical reference solutions.

The paper is structured as follows. Section 2 reports the governing equations under investigation, along with the models used for thermodynamics, transport properties and chemical source terms. The high-order shock-capturing scheme is presented in Section 3. Sections 4 and 5 contain a selection of preliminary validation cases and multi-dimensional turbulent configurations, respectively. Conclusions are drawn in Section 6.

Governing equations

Our goal is to simulate flows governed by the compressible Navier-Stokes equations for multicomponent chemically-reacting gases, written in differential form:

𝜕𝜌 𝜕𝑡 + (𝜕𝜌𝑢 𝑗 ) 𝜕𝑥 𝑗 = 0 (1) 𝜕(𝜌𝑢 𝑖 ) 𝜕𝑡 + 𝜕 ( 𝜌𝑢 𝑖 𝑢 𝑗 + 𝑝𝛿 𝑖𝑗 ) 𝜕𝑥 𝑗 = 𝜕𝜏 𝑖𝑗 𝜕𝑥 𝑗 (2) 𝜕(𝜌𝐸) 𝜕𝑡 + 𝜕 [ (𝜌𝐸 + 𝑝) 𝑢 𝑗 ] 𝜕𝑥 𝑗 = 𝜕(𝑢 𝑖 𝜏 𝑖𝑗 -𝑞 𝑗 ) 𝜕𝑥 𝑗 - 𝜕 𝜕𝑥 𝑗 ( NS ∑ 𝑛=1 𝜌 𝑛 𝑢 𝐷 𝑛𝑗 ℎ 𝑛 ) ( 3
)
𝜕𝜌 𝑛 𝜕𝑡 + 𝜕 ( 𝜌 𝑛 𝑢 𝑗 ) 𝜕𝑥 𝑗 = - 𝜕𝜌 𝑛 𝑢 𝐷 𝑛𝑗 𝜕𝑥 𝑗 + ω𝑛 (𝑛 = 1, … , NS -1) (4)
In the preceding equations, 𝜌 𝑛 is the density of the 𝑛th species, 𝜌= ∑ NS 𝑛=1 𝜌 𝑛 the mixture density, NS the total number of species, 𝑢 𝑖 the velocity vector components in a Cartesian coordinate system, 𝑝 the pressure, 𝛿 𝑖𝑗 the Kronecker symbol and 𝜏 𝑖𝑗 the viscous stress tensor. In Eq. ( 3), 𝐸 = 𝑒 + 1 2 𝑢 𝑖 𝑢 𝑖 represents the specific total energy (with 𝑒 the mixture specific internal energy) and 𝑞 𝑗 the heat flux; moreover, 𝑢 𝐷 𝑛𝑗 , ℎ 𝑛 , and ω𝑛 denote the 𝑛th species diffusion velocity, specific enthalpy and rate of production, respectively. For temperature values lower than 9000 K, ionization and electronic processes can be usually neglected; air is then modeled as a five-species mixture of N 2 , O 2 , NO, O and N (NS = 5 [START_REF] Park | Nonequilibrium hypersonic aerothermodynamics[END_REF]). To ensure total mass conservation, we solve for the mixture density and NS-1 species conservation equations, the NS-th species being computed as 𝜌 NS = 𝜌 -∑ NS-1 𝑛=1 𝜌 𝑛 . This species is chosen to be Nitrogen since it is the most abundant species (i.e., having the largest mass fraction) throughout the computational domain for each case. The viscous stress tensor is modeled as:

𝜏 𝑖𝑗 = 𝜇 ( 𝜕𝑢 𝑖 𝜕𝑥 𝑗 + 𝜕𝑢 𝑗 𝜕𝑥 𝑖 ) - 2 3 𝜇 𝜕𝑢 𝑘 𝜕𝑥 𝑘 𝛿 𝑖𝑗 , (5)
with 𝜇 the mixture dynamic viscosity. The heat flux is modeled by means of Fourier's law, 𝑞 𝑗 = -𝜆 𝜕𝑇 𝜕𝑥 𝑗

, 𝜆 being the mixture thermal conductivity and 𝑇 the temperature. Each species is assumed to behave as a perfect gas; Dalton's pressure mixing law leads then to the thermal equation of state:

𝑝 = NS ∑ 𝑛=1 𝑝 𝑛 = 𝜌𝑇 NS ∑ 𝑛=1 𝑌 𝑛  𝑛 = 𝑇 NS ∑ 𝑛=1 𝜌 𝑛 𝑅 𝑛 , ( 6
)
𝑌 𝑛 = 𝜌 𝑛 ∕𝜌, 𝑅 𝑛 and  𝑛 being the mass fraction, gas constant and molecular weight of the 𝑛th species, respectively, and  = 8.314 J/mol K the universal gas constant. The thermodynamic properties of high-𝑇 air species are computed considering the contributions of translational, rotational and vibrational modes (named TRV formulation [START_REF] Gnoffo | Conservation equations and physical models for hypersonic air flows in thermal and chemical nonequilibrium[END_REF]); specifically, the internal energy reads:

𝑒 = NS ∑ 𝑛=1 𝑌 𝑛 ℎ 𝑛 - 𝑝 𝜌 , with ℎ 𝑛 = ℎ 0 𝑓 ,𝑛 + ∫ 𝑇 𝑇 ref (𝑐 tr 𝑝,𝑛 + 𝑐 rot 𝑝,𝑛 ) d𝑇 ′ + 𝑒 vib 𝑛 . ( 7
)
Here, ℎ 0 𝑓 ,𝑛 is the 𝑛th species enthalpy of formation at the reference temperature (𝑇 ref = 298.15 K), 𝑐 tr 𝑝,𝑛 and 𝑐 rot 𝑝,𝑛 the translational and rotational contributions to the isobaric heat capacity of the 𝑛th species, computed as

𝑐 tr 𝑝,𝑛 = 5 2 𝑅 𝑛 and 𝑐 rot 𝑝,𝑛 = { 𝑅 𝑛 for diatomic species 0 for monoatomic species, (8)
and 𝑒 vib 𝑛 the vibrational energy of species 𝑛, given by

𝑒 vib 𝑛 = 𝜃 𝑛 𝑅 𝑛 exp (𝜃 𝑛 ∕𝑇 ) -1 , (9)
with 𝜃 𝑛 the characteristic vibrational temperature of each molecule (3393 K, 2273 K and 2739 K for N 2 , O 2 and NO, respectively [START_REF] Park | Nonequilibrium hypersonic aerothermodynamics[END_REF]).

After the numerical integration of the conservation equations, temperature is computed from the specific internal energy by means of Newton-Raphson iterations. Pure species' viscosity and thermal conductivities are computed using curve-fits by Blottner et al. [START_REF] Blottner | Chemically Reacting Viscous Flow Program for Multi-component Gas Mixtures[END_REF] and Eucken's relations [START_REF] Vincenti | Introduction to physical gas dynamics[END_REF]:

𝜇 𝑛 = 0.1 exp[(𝐴 𝑛 ln 𝑇 + 𝐵 𝑛 ) ln 𝑇 + 𝐶 𝑛 ], 𝜆 𝑛 = 𝜇 𝑛 ( 5 2 𝑐 tr 𝑣,𝑛 + 𝑐 rot 𝑣,𝑛 + 𝑐 vib 𝑣,𝑛 ) (10)
where 𝐴 𝑛 , 𝐵 𝑛 and 𝐶 𝑛 are fitted parameters. The corresponding mixture properties are evaluated by means of Wilke's mixing rules [START_REF] Wilke | A viscosity equation for gas mixtures[END_REF]: [START_REF] Martín | A bandwidth-optimized WENO scheme for the effective direct numerical simulation of compressible turbulence[END_REF] where

𝜇 = NS ∑ 𝑛=1 𝑋 𝑛 𝜇 𝑛 ∑ NS 𝑚=1 𝑋 𝑚 𝜙 𝑛𝑚 , 𝜆 = NS ∑ 𝑛=1 𝑋 𝑛 𝜆 𝑛 ∑ NS 𝑚=1 𝑋 𝑚 𝜙 𝑛𝑚
𝑋 𝑛 = 𝑌 𝑛 𝑅 𝑛 ∑ NS 𝑚=1 𝑌 𝑚 𝑅 𝑚
denotes the molar fraction of species 𝑛 and

𝜙 𝑛𝑚 = 1 √ 8 ( 1 +  𝑛  𝑚 ) -1 2 ⎡ ⎢ ⎢ ⎣ 1 + ( 𝜇 𝑛 𝜇 𝑚 ) -1 2 (  𝑚  𝑛 ) 1 4 ⎤ ⎥ ⎥ ⎦ 2 . ( 12
)
Mass diffusion is modeled by means of Fick's law:

𝜌 𝑛 𝑢 𝐷 𝑛𝑗 = -𝜌𝐷 𝑛 𝜕𝑌 𝑛 𝜕𝑥 𝑗 + 𝜌 𝑛 NS ∑ 𝑛=1 𝐷 𝑛 𝜕𝑌 𝑛 𝜕𝑥 𝑗 , ( 13
)
where the first term on the r.h.s. represents the effective diffusion velocity and the second one is a mass corrector term that should be taken into account in order to satisfy the continuity equation when dealing with non-constant species diffusion coefficients [START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. Specifically, 𝐷 𝑛 is an equivalent diffusion coefficient of species 𝑛 into the mixture, computed following Hirschfelder's approximation [START_REF] Hirschfelder | Molecular theory of gases and liquids[END_REF] as

𝐷 𝑛 = 1 -
where 𝐷 𝑚𝑛 is the binary diffusion coefficient of species 𝑚 into species 𝑛, and 𝐴 1,𝑚𝑛 , … , 𝐴 4,𝑚𝑛 are curve-fitted coefficients computed as in Gupta et al. [START_REF] Gupta | A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000 K[END_REF].

The five species interact with each other through a reaction mechanism consisting of five reversible chemical steps, according to Park's model [START_REF] Park | Nonequilibrium hypersonic aerothermodynamics[END_REF][START_REF] Park | A review of reaction rates in high temperature air[END_REF]:

R1 ∶ N 2 + M ⟺ 2N + M R2 ∶ O 2 + M ⟺ 2O + M R3 ∶ NO + M ⟺ N + O + M ( 15
) R4 ∶ N 2 + O ⟺ NO + N R5 ∶ NO + O ⟺ N + O 2
being M the third body (any of the five species considered). Dissociation and recombination processes are described by reactions R1, R2 and R3; whereas the shuffle reactions R4 and R5 represent rearrangement processes. The mass rate of production of the 𝑛th species is governed by the law of mass action:

ω𝑛 =  𝑛 NR ∑ 𝑟=1 ( 𝜈 ′′ 𝑛𝑟 -𝜈 ′ 𝑛𝑟 ) × [ 𝑘 𝑓 ,𝑟 NS ∏ 𝑛=1 ( 𝜌𝑌 𝑛  𝑛 ) 𝜈 ′ 𝑛𝑟 -𝑘 𝑏,𝑟 NS ∏ 𝑛=1 ( 𝜌𝑌 𝑛  𝑛 ) 𝜈 ′′ 𝑛𝑟 ] , (16)
where 𝜈 ′ 𝑛𝑟 and 𝜈 ′′ 𝑛𝑟 are the stoichiometric coefficients for reactants and products in the 𝑟th reaction for the 𝑛th species, respectively, and NR is the total number of reactions. Lastly, 𝑘 𝑓 ,𝑟 and 𝑘 𝑏,𝑟 denote the forward and backward reaction rates of reaction 𝑟, modeled by means of Arrhenius' law.

For configurations in which air is modeled as a single-species calorically-perfect gas, a constant specific heat ratio is considered (𝛾 = 1.4), such that 𝑐 𝑝 = 𝛾𝑅 𝛾 -1

. Since NS = 1 and there is no chemical activity, Eq. ( 4) is not solved and the diffusion velocity is zero. Unless otherwise stated, viscosity is computed by means of Sutherland's Law, 𝜇(𝑇 ) = 𝐶𝑇 3∕2 𝑇 + 𝑆 , with 𝐶 = 1.457933 × 10 -6 kg∕m s K 1∕2 and 𝑆 = 110.4 K, whereas the thermal conductivity follows a constant-Prandtl assumption, for which 𝜆 = Pr∕(𝜇𝑐 𝑝 ).

Numerical method

In this section we describe the high-order shock-capturing centraldifference numerical scheme under investigation. We first present the most important ingredient, i.e. the spatial discretization scheme for the nonlinear terms, for a 1D system of hyperbolic conservation laws:

𝜕𝑤 𝜕𝑡 + 𝜕𝑓 (𝑤) 𝜕𝑥 = 0 ( 17
)
where 𝑤 is the vector of conservative variables and 𝑓 (𝑤) the flux function, such that 𝐴 = 𝜕𝑓 ∕𝜕𝑤 is a diagonalizable matrix with real eigenvalues. Extension to multidimensional cases is straightforwardly carried out by applying the scheme in each direction. Introducing the classical difference and cell-average operators over one cell:

(𝛿•) 𝑗 ∶= (•) 𝑗+ 1 2 -(•) 𝑗-1 2 , (𝜇•) 𝑗+ 1 2 ∶= 1 2
[ (•) 𝑗+1 + (•) 𝑗 ] (18)
and considering a regular Cartesian grid with constant mesh spacing 𝛿𝑥 (so that 𝑥 𝑗 =𝑗 𝛿𝑥), a conservative semi-discrete approximation of the spatial derivative writes:

( 𝜕𝑤 𝜕𝑡 ) 𝑗 + (𝛿 ) 𝑗 𝛿𝑥 = 0 (19)
The numerical flux at cell interface

𝑗+ 1 2 ,  𝑗+ 1 2
, can be calculated using a simple upwind scheme, written hereafter as the sum of a central approximation  and a dissipative term  [START_REF] Lerat | Schémas dissipatifs précis à l'ordre trois pour les systèmes hyperboliques = Dissipative third-order accurate schemes for hyperbolic systems[END_REF]:

 𝑗+ 1 2 = 𝑗+ 1 2 - 𝑗+ 1 2 (20)  𝑗+ 1 2 = 𝑓 𝑗+1 + 𝑓 𝑗 2 = (𝜇𝑓 ) 𝑗+ 1 2 (21)
 𝑗+ 1 2 = 1 2 |𝑄 𝑗+ 1 2 |(𝑤 𝑗+1 -𝑤 𝑗 ) = 1 2 (|𝑄|𝛿𝑤) 𝑗+ 1 2 (22)
with 𝑄 a dissipation matrix. For instance, 𝑄 = 1 2 𝐼 (with 𝐼 the identity matrix) gives Lax-Friedrichs' scheme [START_REF] Leveque | Numerical methods for conservation laws[END_REF], 𝑄 = 𝜌(𝐴)𝐼 Rusanov's scheme [START_REF] Rusanov | Calculation of intersection of non-steady shock waves with obstacles[END_REF] (𝜌(𝐴) being the spectral radius of the flux Jacobian matrix 𝐴), and 𝑄 = 𝐴 𝑅 Roe's scheme [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] (with 𝐴 𝑅 the Roe matrix). All such schemes are monotonicity preserving according to Godunov's theorem [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF], but only first-order-accurate, which makes them unsuitable for turbulence resolving simulations.

A standard way for increasing accuracy consists in using MUSCL extrapolations [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF]. Here we adopt an alternative approach, first introduced in [START_REF] Lerat | Schémas dissipatifs précis à l'ordre trois pour les systèmes hyperboliques = Dissipative third-order accurate schemes for hyperbolic systems[END_REF] to construct a third-order scheme and generalized to any order of accuracy in [START_REF] Lerat | Approximations d'ordre élevé pour les écoulements compressibles[END_REF], which consists in recursively correcting the truncation error of Eqs. ( 21), [START_REF] Visbal | Shock capturing using compact-differencing-based methods[END_REF]. For a scheme of order 2𝑃 + 3 with a stencil of 2(𝑃 + 2) + 1 points, the numerical flux is of the form:

 𝑗+ 1 2 = [( 𝐼 - 𝑃 ∑ 𝑝=0 𝑎 𝑝 𝛿 2+2𝑝 ) 𝜇𝑓 - 𝑎 𝑃 2 |𝑄|𝛿 2𝑃 +3 𝑤 ] 𝑗+ 1 2 (23)
The coefficients 𝑎 𝑝 have alternate negative and positive signs as 𝑃 increases. The interested reader may refer to [START_REF] Lerat | Approximations d'ordre élevé pour les écoulements compressibles[END_REF] for details about their calculation.

The preceding high-order, constant coefficient schemes are not total variation diminishing (TVD) nor monotonicity preserving, and slope or flux limiters should be introduced to avoid the appearance of spurious oscillations. Note that only schemes of odd order accuracy are included in the family, with a leading truncation error term given by:

𝜀 = 𝑎 𝑃 2 𝛿𝑥 2𝑃 +3 𝜕 2(𝑃 +2) 𝑤 𝜕𝑥 2(𝑃 +2) (24)
proportional to an even-order derivative, i.e. of dissipative nature. Explicit odd-order schemes with constant coefficients have been shown to remain stable in the maximum norm 𝐿 ∞ under the CFL constraint for problems with non-smooth initial conditions [START_REF] Thomée | Stability of difference schemes in the maximum-norm[END_REF][START_REF] Desprès | Uniform asymptotic stability of Strang's explicit compact schemes for linear advection[END_REF]. This implies that, although the considered schemes do generate spurious oscillations around discontinuities, such oscillations remain bounded in space and time. For the sake of clarity, we give hereafter the expressions of the schemes of order 3 to 9 (i.e. 𝑃 = 0, 1, 2 and 3) of the preceding family:

 𝑗+ 1 2 = [( 𝐼 - 1 6 𝛿 2 ) 𝜇𝑓 + 1 12 |𝑄|𝛿 3 𝑤 ] 𝑗+ 1 2 (order 3) ( 25
)
 𝑗+ 1 2 = [( 𝐼 - 1 6 𝛿 2 + 1 30 𝛿 4 ) 𝜇𝑓 - 1 60 |𝑄|𝛿 5 𝑤 ] 𝑗+ 1 2 (order 5) (26)
 𝑗+ 1 2 = [( 𝐼 - 1 6 𝛿 2 + 1 30 𝛿 4 - 1 140 𝛿 6 ) 𝜇𝑓 + 1 280 |𝑄|𝛿 7 𝑤 ] 𝑗+ 1 2 (order 7) (27)
 𝑗+ 1 2 = [( 𝐼 - 1 6 𝛿 2 + 1 30 𝛿 4 - 1 140 𝛿 6 + 1 630 𝛿 8 ) 𝜇𝑓 - 1 1260 |𝑄|𝛿 9 𝑤 ] 𝑗+ 1 2 (order 9) (28)
The same schemes can be derived in the finite volume framework by applying MUSCL reconstruction to the physical fluxes [START_REF] Lerat | Approximations d'ordre élevé pour les écoulements compressibles[END_REF]. Specifically, for 𝑄 = 𝐴 𝑅 , one recovers a flux-extrapolation higher-order extension of Roe's scheme. This choice has the advantage of introducing the minimal amount of numerical damping along each characteristic field. However, the extension of Roe's average to real-gas flows is not unique and can introduce significant overcost (see, e.g. [START_REF] Glaister | An approximate linearised Riemann solver for the Euler equations for real gases[END_REF][START_REF] Cinnella | Roe-type schemes for dense gas flow computations[END_REF]).

With the purpose of simplifying the application to real-gas flows and reducing computational cost, a different approach is adopted. We first observe that the preceding schemes can be considered as standard central-difference approximations of order 2(𝑃 + 2) on 2(𝑃 + 2) + 1 points of the flux derivative, plus a high-order artificial viscosity of order 2𝑃 + 3 on the same stencil depending on matrix 𝑄. Afterwards, we choose 𝑄 = 𝜌(𝐴), as in Rusanov's scheme, which is less optimal than Roe's matrix but avoids complexities associated with the extension of the approximate Riemann solver to real gases. This leads to a scalar numerical dissipation term of the form:

 𝑗+ 1 2 = [ 𝑎 𝑃 2 𝜌(𝐴)𝛿 2𝑃 +3 𝑤 ] 𝑗+ 1 2 (29)
Finally, we nonlinearly combine the preceding high-order dissipation with a lower-order term activated in the vicinity of flow discontinuities by means of a highly selective shock sensor. The dissipation  then becomes:

 𝑗+ 1 2 = 𝜌(𝐴) 𝑗+ 1 2 [ 𝜀 2 𝛿𝑤 + (-1) (𝑃 +1) 𝜀 2(𝑃 +2) 𝛿 (2𝑃 +3) 𝑤 ] 𝑗+ 1 2 (30)
with

𝜀 2 𝑗+ 1 2 = 𝑘 2 max(𝜑 𝑗 , 𝜑 𝑗+1 ), 𝜀 2(𝑃 +2) 𝑗+ 1 2 = max(0, 𝑘 2(𝑃 +2) -𝑘 𝜀 𝜀 2 𝑗+ 1 2 ), (31)
where 𝑘 2 and 𝑘 2(𝑃 +2) are adjustable dissipation coefficients and 𝑘 𝜀 is a constant equal to 𝑎 0 ∕𝑎 𝑃 , determining the threshold below which the higher-order dissipation is switched off.

For schemes of order 3 to 9, this gives the following expressions:

 𝑗+ 1 2 = 𝜌(𝐴) 𝑗+ 1 2 [ 𝜀 2 𝛿𝑤 -𝜀 4 𝛿 3 𝑤 ] 𝑗+ 1 2 𝜀 4𝑗+ 1 2 = max(0, 𝑘 4 -𝜀 2𝑗+ 1 2 ) (32)
 𝑗+ 1 2 = 𝜌(𝐴) 𝑗+ 1 2 [ 𝜀 2 𝛿𝑤 + 𝜀 6 𝛿 5 𝑤 ] 𝑗+ 1 2 𝜀 6𝑗+ 1 2 = max(0, 𝑘 6 - 1 5 𝜀 2𝑗+ 1 2 ) (33)
 𝑗+ 1 2 = 𝜌(𝐴) 𝑗+ 1 2 [ 𝜀 2 𝛿𝑤 -𝜀 8 𝛿 7 𝑤 ] 𝑗+ 1 2 𝜀 8𝑗+ 1 2 = max(0, 𝑘 8 - 3 70 𝜀 2𝑗+ 1 2 ) (34)
 𝑗+ 1 2 = 𝜌(𝐴) 𝑗+ 1 2 [ 𝜀 2 𝛿𝑤 + 𝜀 10 𝛿 9 𝑤 ] 𝑗+ 1 2 𝜀 10𝑗+ 1 2 = max(0, 𝑘 10 - 1 105 𝜀 2𝑗+ 1 2 ) (35)
For 𝑘 2 =0 and 𝑘 2(𝑃 +2) = 𝑎 𝑃 2 one recovers the upwind schemes of Eq. ( 23). The activation of the low-order dissipation component rests on the value of the shock-capturing sensor 𝜑 𝑗 , which consists in a combination of different terms. Specifically, one has:

𝜑 𝑗 = 1 2 [ 1 -tanh ( 2.5 + 10 𝛿𝑥 𝑐 ∇ ⋅ 𝐮 )] ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ I × (∇ ⋅ 𝐮) 2 (∇ ⋅ 𝐮) 2 + |∇ × 𝐮| 2 + 𝜖 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ II × | | | | | 𝑝 𝑗+1 -2𝑝 𝑗 + 𝑝 𝑗-1 𝑝 𝑗+1 + 2𝑝 𝑗 + 𝑝 𝑗-1 | | | | | ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ III (36)
The second and third terms denote the classical Ducros' [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] and Jameson's pressure-based [START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes[END_REF] shock sensors, respectively, 𝜖 being a small positive value (𝜖 = 10 -16 ) to avoid division by zero. Their combination palliates to some of the deficiencies related to the stand-alone application of the Ducros' sensor, which takes into account only the relative magnitudes of dilation and vorticity and may result in unwanted activations of the shock-capturing term in regions where both of these two quantities are small (e.g., in the irrotational flow outside boundary layers or mixing layers). To correct this deficiency, the constant 𝜖 can be parametrized by introducing suitable characteristic velocity and length scales depending on the flow under investigation [START_REF] Pirozzoli | Numerical methods for high-speed flows[END_REF]. The main drawback of this method resides in the loss of generality, 𝜖 being transformed in a configuration-dependent parameter. The introduction of the pressure-based sensor allows one to bypass the activation of Ducros' sensor, strongly reducing the amount of low-order dissipation injected and leaving the acoustic perturbations crossing the domain much less affected. The first term of Eq. ( 36) takes into account the Ducros' sensor modification of Bhagatwala & Lele [START_REF] Bhagatwala | A modified artificial viscosity approach for compressible turbulence simulations[END_REF], initially proposed to enhance the selectivity of the artificial bulk viscosity in the Localized Artificial Diffusivity (LAD) technique. In regions of positive dilation 𝜑 𝑗 is switched off, whereas its value increases slowly with the magnitude of the negative dilation. Moreover, the scaling factor 10𝛿𝑥∕𝑐 has the twofold role of (i) normalizing the grid-dependent numerical dilation and (ii) making it invariant with the mesh-size. The sensor is (1) in high-divergence regions and tends to zero in vortex-dominated regions, allowing the capture of flow discontinuities with minimal damping of the vortical structures inside the flow.

In the following, we mostly focus on the ninth-order accurate scheme of the preceding family. Far from flow discontinuities, such scheme has low phase and dissipation errors. Its leading truncation error term is of the form 𝑘 10 𝛿𝑥 9 𝜕 10 𝑤 𝜕𝑥 10 , i.e. it is consistent with a tenthorder viscosity. Such viscosity term acts differently according to the wavenumber, dissipating scales characterized by reduced wavenumbers of about 0.35𝜋 or higher (i.e. wavelengths that are discretized with less than 6 mesh points), while leaving larger scales almost unaffected. In Fig. 1 we report the dissipation and phase errors of the ninthorder scheme for the approximation of a linear advection problem, as a function of the reduced wavenumber 𝑘𝛿𝑥. Lower-order schemes of the same family are also reported to illustrate the effect of increasing accuracy. For all schemes, the dispersion error is exactly the same as for the standard central scheme of order 2(𝑃 +2). Thanks to its selectiveness in the wavenumber space, the ninth-order dissipation constitutes a suitable implicit subgrid regularization term for LES simulations [START_REF] Gloerfelt | Large eddy simulation requirements for the flow over periodic hills[END_REF], with the capability of seamlessly converging to DNS in smooth flow regions as the grid is refined. Unless otherwise stated, we set 𝑘 2 =1 and 𝑘 10 = 1 1260 for all computations. In Navier-Stokes calculations, the viscous flux derivatives are approximated by fourth-order-accurate central formulae, if not specified differently. Finally, in all of the following calculations time advancement is carried out by means of an explicit third-order TVD Runge-Kutta scheme [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]. The non-uniformity of the wall-normal mesh spacing is taken into account by a suitable 1-D coordinate transformation. Near the non-periodic boundaries, the finite-difference stencil for the convective terms is progressively reduced down to the fourth order (and, correspondingly, the numerical dissipation term); then, both the convective and viscous fluxes are evaluated from the interior points by using fourth-order backward differences.

Preliminary validations

The ninth-order shock-capturing central scheme under investigation is first applied to selected inviscid test cases, in order to verify its convergence order in smooth flow regions and to assess its shock-capturing capabilities.

Isentropic vortex advection problem

The accuracy of the discretization scheme in smooth inviscid flow is quantified for the well-known two-dimensional isentropic vortex advection problem [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF][START_REF] Davoudzadeh | Accuracy evaluation of unsteady CFD numerical schemes by vortex preservation[END_REF][START_REF] Yee | Entropy splitting and numerical dissipation[END_REF], in which an inviscid vortex is superimposed to an uniform, perfect-gas (𝛾 = 1.4) air flow. The perturbations in velocity and temperature are given by

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝛿𝑢 = - 𝑦 𝑅 𝛺 𝛿𝑣 = 𝑥 𝑅 𝛺 𝛿𝑇 = - 𝛾 -1 2 𝛺 2 with ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝛺 = 𝛽 exp ( - 1 2𝜎 [ ( 𝑥 𝑅 ) 2 + ( 𝑦 𝑅 ) 2 ]) 𝛽 = 𝑀 ∞ 5 √ 2 4𝜋 e 1 2 (37)
These allow to define, along with the isentropic relations, the initial flow conditions of the primitive variables as [START_REF] Boris | New insights into large eddy simulation[END_REF] where the subscript (•) 0 denotes a quantity given at 𝑡 = 0, and the symbol (•) indicates a nondimensional quantity (𝜌 ∞ , 𝑇 ∞ and 𝑎 ∞ being the characteristic density, temperature and velocity, respectively). In the classical case [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF], one has 𝑅 = 𝜎 = 1, 𝑀 ∞ = √ 2 𝛾 and 𝛼 = 45 • . Periodic conditions are applied at the boundaries. The length of the computational domain has been increased from [𝐿 𝑥 , 𝐿 𝑦 ] = [-5, 5] to [-10, 10] in order to reduce the influence of the small artificial shear layers generated near the boundaries by the non-zero velocity perturbations, which can pollute the results when considering smaller domains [START_REF] Spiegel | A survey of the isentropic Euler vortex problem using high-order methods[END_REF].

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ρ0 = (1 + 𝛿𝑇 ) 1 𝛾-1 ũ0 = 𝑀 ∞ cos 𝛼 + 𝛿𝑢 ṽ0 = 𝑀 ∞ sin 𝛼 + 𝛿𝑣 p0 = 1 𝛾 (1 + 𝛿𝑇 ) 𝛾 𝛾-1
The error with respect to the exact solution (pure advection of the initial vortex) is measured as

𝜀 𝜏,ℎ = 𝐿 2 (𝛹 𝜏,ℎ (𝑖,𝑗) ) = √ 1 ℎ ∑ 𝑖,𝑗 ( 𝛹 𝜏,ℎ (𝑖,𝑗) -𝛹 𝑒𝑥 (𝑖,𝑗) ) 2 , ( 39
)
with 𝛹 𝜏,ℎ (𝑖,𝑗) and 𝛹 𝑒𝑥 (𝑖,𝑗) the computed and exact values of a generic flow variable at the grid point (𝑖, 𝑗), and ℎ = 𝐿 𝑥 ∕(𝑁 -1) the spatial grid size, 𝑁 = 𝑁 𝑥 = 𝑁 𝑦 being the number of grid points. For an unsteady problem, the numerical error may be written as 𝜀 𝜏,ℎ = 𝐶 𝜏 𝜏 𝑝 + 𝐶 ℎ ℎ 𝑞 , where 𝐶 𝜏 and 𝐶 ℎ are some constants, 𝜏 the time-step, ℎ the grid size, and 𝑝 and 𝑞 the order of the temporal integration and spatial discretization schemes, respectively. The ratio of error decay between numerical solutions using time steps 𝜏 and 𝑚𝜏 and grid sizes ℎ and 𝑛ℎ writes:

𝑟 𝑒 = 𝜀 𝑚𝜏,𝑛ℎ 𝜀 𝜏,ℎ = 𝐶 𝜏 (𝑚𝜏) 𝑝 + 𝐶 ℎ (𝑛ℎ) 𝑞 𝐶 𝜏 𝜏 𝑝 + 𝐶 ℎ ℎ 𝑞 = 𝛿𝑚 𝑝 + 𝑛 𝑞 𝛿 + 1 with 𝛿 = 𝐶 𝜏 𝜏 𝑝 𝐶 ℎ ℎ 𝑞 . ( 40
)
Assuming 𝑚, 𝑛 > 1, the overall accuracy ranges between 𝑞 for small 𝛿 (meaning that the temporal error is negligible with respect to spatial error) and 𝑝 for large 𝛿.

Several runs are carried out by varying the number of grid points (between 100 2 and 500 2 ) and the CFL number (from 0.6 to 0.00625); the errors measured after one cycle (i.e. when the vortex returns to the initial position for the first time) are shown in Fig. 2, where each symbol denotes a run. It is worth pointing out that each simulation has been performed twice, with and without the shock-capturing term (𝑘 2 =1 and 𝑘 2 =0, respectively). As expected, differences were found to be negligible, with an influence on the 6th significant digit of the computed error norm. We then report hereafter only the results obtained for 𝑘 2 =1. Fig. 2a displays the value of the numerical error as a function of the CFL number for several grids (color plot). In order to recover the formal order of accuracy of the spatial discretization scheme for a given grid, the CFL number should be small enough not to affect 𝜀 𝜏,ℎ ; i.e., 𝜀 𝜏,ℎ should reach a plateau for sufficiently small values of the CFL number. This is clearly visible in Fig. 2a, which also highlights that even smaller CFL numbers should be considered for grids finer than 500 2 . In Fig. 2b we report the error as a function of 𝑁 for several CFL numbers. The order of accuracy of the numerical solution varies as expected between 3 to 9 (formal orders of the time integration and spatial discretization numerical schemes, respectively) as the CFL number gets smaller. For CFL numbers of 0.2 or smaller, a more or less extended region with ninth-order slope is observed; the slope tends to decrease for finer grids, due to the larger relative importance of the temporal error. The local slopes are shown in Fig. 2c versus the number of grid points. The presence of a maximum can be explained by considering the temporalto-spatial error ratio, 𝛿. In order to keep a constant 𝛿 when refining the grid, the timestep should be decreased as 𝜏 ∝ ℎ 𝑞∕𝑝 = ℎ 9∕3 = ℎ 3 ; whereas a constant CFL imposes 𝜏 ∝ ℎ (i.e., a varying 𝛿). At CFL = 1.25 × 10 -2 Fig. 2. Analysis of the order of convergence of the error norm for the isentropic vortex convection configuration. (a): error of the 𝐿 2 -norm as a function of the CFL number for computational grids ranging from 100 2 to 500 2 ; (b): error of the 𝐿 2 -norm as a function of the number of grid points for CFL numbers ranging from 0.6 to 0.00625; (c): overall spatio-temporal order of accuracy as a function of the number of grid points for CFL numbers ranging from 0.6 to 0.00625. and 6.25 × 10 -3 the formal order of the spatial discretization scheme is recovered over the whole range of considered grid resolutions; this confirms previous studies suggesting the use of CFL ≲ 0.01 [START_REF] Renzo | HTR solver: An open-source exascale-oriented taskbased multi-GPU high-order code for hypersonic aerothermodynamics[END_REF].

A parametric study on the influence of the maximum strength of the perturbation was also carried out. The value of 𝛽 in Eq. ( 37) was increased by a factor in the range [1, 1.9], with the minimum density value in the vortex core passing from 𝜌 min ≈ 0.5 to 𝜌 min ≈ 4 × 10 -3 , as shown in Fig. 3a. Fig. 3b shows that an increase of the perturbation strength leads to larger absolute values of the error norm, but does not alter substantially the formal order of accuracy of the numerical scheme.

Ideal-gas shock tube problems

The next test case is a classical benchmark for the assessment of shock-capturing capabilities. Specifically, we consider the well-known Sod [START_REF] Sod | A survey of several finite difference methods for stystems of nonlinear hyperbolic conservation laws[END_REF] and Lax [START_REF] Lax | Weak solutions of nonlinear hyperbolic equations and their numerical computation[END_REF] 1D shock tube problems, corresponding to the following Riemann problems: Numerical results are compared to the exact solutions at the nondimensional time 𝑡 = 0.2 (𝛥𝑡 = 5 × 10 -4 ) and 𝑡 = 0.13 (𝛥𝑡=10 -3 ), respectively. For these test cases, the artificial viscosity coefficients are chosen equal to 𝑘 2 = 2, 𝑘 10 = 1∕630. The numerical domain is discretized with 200 evenly-spaced grid points; air is considered either as a singlespecies, calorically-perfect gas (𝛾 = 1.4) or as a two-species non-reacting mixture of Nitrogen (𝑌 N 2 = 0.79) and Oxygen (𝑌 O 2 = 0.21). Fig. 4 shows the results for density, velocity and pressure, whose profiles are compared to the corresponding exact solutions. A good agreement is shown for both cases and the single-species and multi-species numerical solutions are perfectly superposed. The slight smearing of the numerical solution across the contact discontinuity can be attributed to the highorder dissipation term, the shock-capturing component being switched off in that region due to the constant value of the pressure. Note that a similar behavior is observed for other high-order schemes [START_REF] Fu | Targeted ENO schemes with tailored resolution property for hyperbolic conservation laws[END_REF].

(𝜌, 𝑢, 𝑝) SOD = { (1,

Shu-Osher problems

The Shu-Osher problem [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes, II[END_REF] consists of a 𝑀=3 shock propagating in a perturbed density field and allows to evaluate the behavior of the numerical scheme for a simplified shock-turbulence interaction configuration. The extent of the computational domain is [-5, 5] and the initial conditions are (𝜌, 𝑢, 𝑝) SHU = { (3.857143, 2.629369, 10.33333) 𝑥 < -4

(1 + 0.2 sin(5𝑥), 0, 1)

𝑥 ≥ -4 (42)
The 1D Euler equations are solved on a uniform mesh with 𝑁 = 200 and a reference solution is computed with the same numerical scheme on a mesh with 𝑁 = 2000. Fig. 5 shows results for density, pressure, velocity and entropy at the final time 𝑡 = 1.8. The profiles of pressure and velocity are in good agreement with the reference solution, with some oscillations registered near the shock (which is captured reasonably well) that remain bounded to small values throughout the simulation. On the contrary, the density and entropy profiles exhibit a stronger damping after the shock train passage. This is partly due to the use of a scalar artificial viscosity, which introduces the same amount of dissipation for all characteristic fields, whether they be of acoustic or entropic nature. Present results are in agreement with those of other classical numerical schemes [START_REF] Johnsen | Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves[END_REF]. Of note, the introduction of the Bhagatwala & Lele correction to the shock-capturing term strongly enhances the entropy waves resolution (contrary to the shock tube cases, where no appreciable differences were found), mitigating the spurious activation of Jameson's pressure-based sensor. A possible fix for this overlydissipative behavior would be a selective reduction of the numerical dissipation by using different 𝑘 2 values for each equation. However, this approach would increase the number of tuning parameters, and was thus discarded in order to keep the numerical strategy as simple and general as possible. Lastly, it should be pointed out that the core regions of turbulent flows are intrinsically rotational: the Ducros sensor, which always degenerates to unity in 1D flows, plays therefore a fundamental role in minimizing the numerical dissipation outside of shocked-flow turbulent regions. For the purpose of comparison with other numerical strategies we also considered a two-dimensional generalization of the Shu-Osher problem, which consists in a vorticity/entropy wave interacting with a normal shock [START_REF] Johnsen | Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves[END_REF]. The extents of the computational domain are 𝑥 ∈ [0, 4𝜋), 𝑦 ∈ [-𝜋, 𝜋), with 𝛥𝑥 = 𝜋∕50 and 𝛥𝑦 = 𝜋∕16; The vorticity content at 𝑡 = 25 is reported in Fig. 6 and compared to results from the codes ADPDIS3D, Hybrid and Stan, digitalized from the Ref. [START_REF] Johnsen | Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves[END_REF] (the reader may refer to the original paper for details about the different shock-capturing approaches). Each test has been repeated twice, with 𝑘 2 =0 and 𝑘 2 =1, to assess the influence of the shock-capturing term.

Both the instantaneous and mean quantities exhibit profiles of vorticity (Fig. 6) and kinetic energy (not shown) comparable to the reference. As expected, the shock-capturing activation reduces the postshock oscillations and slightly damps the mean kinetic and vortical contents, which remain close to the reference linear analysis (LIA) results [START_REF] Mahesh | The interaction of a shock wave with a turbulent shear flow[END_REF]. Interestingly, the post-shock mean vorticity profiles shows the same oscillating trend of the Stan code (which uses a LAD-based shock capturing), keeping however a better response in each case. We also note that for such a mild configuration, a lower value of 𝑘 2 would be more appropriate, resulting in a larger response.

Real-gas shock tube problem

The last inviscid test case consists of a multi-species, hightemperature shock tube designed by Grossman & Cinnella to study thermochemical effects [START_REF] Grossman | Flux-split algorithms for flows with non-equilibrium chemistry and vibrational relaxation[END_REF]. For such a configuration, however, thermal nonequilibrium effects are rather mild and can be neglected, taking into account only chemical nonequilibrium as shown in previous validation studies [START_REF] Renzo | HTR solver: An open-source exascale-oriented taskbased multi-GPU high-order code for hypersonic aerothermodynamics[END_REF]. At 𝑡 = 0, the initial conditions are the following:

(𝑝, 𝑢, 𝑇 ) GROSSMAN = { 1.95256 × 10 5 Pa, 𝑢 𝐿 = 0 m∕s, 9000 K 𝑥 < 0 10 4 Pa, 𝑢 𝑅 = 0 m∕s, 300 K 𝑥 ≥ 0 (45)
Chemically-reacting air is modeled by means of Park's 5-species model (N 2 , O 2 , N, O, NO); the initial values for the species mass fractions correspond to the mixture equilibrium composition at the given pressure and temperature for the right and left states, respectively. Due to the stiffness of the chemical source terms, the solution is advanced in time using a CFL number of 0.02, as suggested in [START_REF] Grossman | Flux-split algorithms for flows with non-equilibrium chemistry and vibrational relaxation[END_REF] for explicit Runge-Kutta time integrations. The simulation is stopped when the shock-wave reaches the location 𝑥 = 0.11 m. Given the severe conditions developed by flow in the present problem, the Bhagatwala & Lele's modification of the shock sensor was turned off for better numerical robustness. Results at the final time are shown in Fig. 7. The numerical scheme is able to capture correctly the rarefaction wave, the contact discontinuity and the shock; moreover, the distributions of the species mass fractions agree very well with data from [START_REF] Grossman | Flux-split algorithms for flows with non-equilibrium chemistry and vibrational relaxation[END_REF].

Applications to multiscale turbulent flows

In this section we analyze the performance of the shock-capturing scheme for viscous compressible flows with shocks and fine-detail vortical structures. The applications range from a two-dimensional under-expanded jet flow to a hypersonic turbulent boundary layer in chemical nonequilibrium.

Two-dimensional underexpanded jet

A N 2 -O 2 inert, highly underexpanded jet has been first considered to test the suitability of the numerical scheme to deal with strong discontinuities and vortical layers in multi-species flows. This configuration represents the starting point for the study of reactive jets and has been widely analyzed in the past years, both experimentally and numerically. In addition, it has been shown that 2D simulations are able to capture reasonably well some detailed features of the problem, such as the location and dimension of the Riemann wave, the scale of the jet and the timewise averages of the thermodynamic variables.

The present setup is similar to the one reported in Martínez Ferrer et al. [START_REF] Ferrer | A detailed verification procedure for compressible reactive multicomponent Navier-Stokes solvers[END_REF]. Specifically, the pressure 𝑃 0 and temperature 𝑇 0 of the mixture in the injection plane are set to 15 atm and 1000 K, whereas the ambient values are 1 atm and 300 K, respectively, resulting in a nozzle pressure ratio (NPR) of 15. The inflow jet Mach number is equal to 1 and a slow coflow at 𝑀 = 0.05 is imposed consistently with previous studies [START_REF] Ferrer | A detailed verification procedure for compressible reactive multicomponent Navier-Stokes solvers[END_REF][START_REF] Su | Numerical simulations of inert and reactive highly underexpanded jets[END_REF]. The height of the injector exit is equal to 𝐷 = 3 cm and the extent of the computational domain is 𝐿 𝑥 × 𝐿 𝑦 = 50𝐷 × 25𝐷. At the inflow, the jet velocity is prescribed by means of an hyperbolic-tangent profile (see ''Profile 2'' of Michalke [START_REF] Michalke | Survey on jet instability theory[END_REF]):

𝑢(𝑟) 𝑈 = 1 2 { 1 + tanh [ 0.25 𝑅 𝜃 ( 𝑅 𝑟 - 𝑟 𝑅 )]} ( 46
)
where 𝑈 is the inflow centerline velocity, 𝑅 the inflow half-height, 𝜃 the initial momentum thickness of the shear layer and 𝑟 the local distance from the jet axis. The ratio 𝑅∕𝜃 is an important parameter influencing the jet stability, which is mainly related to the introduction of vorticity in the jet shear layer. In our study, we consider 𝑅∕𝜃 = 12.5. At the remaining outflow boundaries, characteristic conditions are imposed, along with sponge regions and grid stretching to avoid spurious wave reflections.

Five computational grids were considered, with dimensional mesh sizes equal to 𝛥𝑥 = 𝛥𝑦 = 4, 2, 1, 0.5 and 0.25mm. Fig. 8 shows an instantaneous snapshot of several quantities for the most refined grid used in the current study, revealing the challenging physics of the problem (for a thorough description of the flow physics of free underexpanded jets, the reader may refer to Franquet et al. [START_REF] Franquet | Free underexpanded jets in a quiescent medium: A review[END_REF]). The influence of the grid resolution is shown in Fig. 9, where the axial profiles of the Mach number (panel a), temperature and density (panel b) are shown. The location and width of the Mach disk is in agreement with experimental observations [START_REF] Sheeran | Observations on jet flows from a two-dimensional, underexpanded, sonic nozzle[END_REF], as well as recent numerical results obtained for similar NPR values [START_REF] Su | Numerical simulations of inert and reactive highly underexpanded jets[END_REF]. Averaged profiles of the thermodynamic quantities are in accordance with results shown by Su et al. [START_REF] Su | Numerical simulations of inert and reactive highly underexpanded jets[END_REF]; the small post-shock discrepancies in density may be attributed to the slight difference in the inflow N 2 and O 2 mass fractions. The r.m.s. temperature values extracted at the vertical line 𝑥∕𝐷 = 20 reveal nearly-converged profiles for the 𝛥𝑥 = 0.5 mm grid, whereas 𝛥𝑥 = 2 mm grid is already sufficiently fine to capture correctly the mean profiles of first-order quantities. In any case, the numerical scheme is shown to handle satisfactorily under-resolved grids without the insurgence of numerical instabilities nor non-physical results.

In Section 3, we stressed the importance of using a well-constructed shock-capturing sensor for the numerical simulation of high-speed turbulent flows. The current 2D jet represents a suitable configuration for benchmarking such a sensor, due to the concurrent presence of strong steady shocks, turbulent shear layers and propagation of acoustic waves. The influence of each one of the three components of the sensor presented in Eq. ( 36) is shown in Fig. 10, displaying their isocontours on an instantaneous snapshot. Values below 10 -2 have been cropped in order to highlight the flow regions marked by each term as potential shocks. As expected, the Ducros' sensor (panel a) is not able to properly identify strong gradients if applied alone; its values are shown to be close to unity even in smooth-flow regions, essentially because of the absence of vorticity (shown in Fig. 8). The Bhagatwala & Lele correction and the Jameson's pressure-based sensors (panel b and c, respectively) allow one to obtain a more satisfactory large-gradients tracing, the latter being slightly more selective. One may object the utility of using both of them; however, their combined use palliates each other's deficiency. Specifically, the Jameson's sensor could hardly tell the difference between a shock and a high-vorticity, purely-solenoidal region; on the other hand, the Bhagatwala & Lele correction can result in excessive damping of all dilatational motions when dealing with slightly under-resolved simulations. The combination of the three sensors (panel d) results in an excellent localization of the regions in which low-order numerical dissipation should be injected; that is, along the lateral intercepting shocks, the Mach disk and the reflected shock waves. A few wave fronts are also marked downstream of the Riemann wave; however, the sensor magnitude is almost negligible there (i.e., lower than 0.02 everywhere), and such is the amount of injected dissipation.

Two-dimensional viscous shock tubes

A two-dimensional chemically-reacting viscous shock tube is now investigated. In such a configuration, a complex unsteady interaction is generated between the incident boundary layer and the shock wave reflected at the end wall. This case was first studied by Daru & Tenaud [START_REF] Daru | Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme[END_REF], who performed numerical simulations for a perfect gas at different Reynolds numbers by means of a high-resolution monotonicity-preserving scheme. The computational domain is a square shock tube with a unit side length and adiabatic walls at each boundary. The initial conditions for the perfect gas case read: with 𝛾 = 1.4. Three different resolutions are analyzed, using respectively 2000, 4000 and 6000 grid points per direction. Reference data are given for 𝑁 = 4000. Fig. 11 shows the density distribution along the bottom wall and at 𝑦 = 0.05 for the case 𝑅𝑒 = 1000 at the final time (𝑡 = 1). Results are in perfect agreement with reference data and grid convergence is already obtained at 𝑁 = 4000. Such a configuration was recently extended to chemical nonequilibrium conditions by Chen et al. [START_REF] Chen | Numerical study of chemically reacting flow in a shock tube using a high-order point-implicit scheme[END_REF], who considered a relatively low initial temperature (𝑇 ref = 800 K) in order to trigger only oxygen dissociation-recombination reactions and therefore simulate a simplified O 2 /O/N 2 mixture. Here, we select 𝑇 ref = 2000 K such that NO production is enabled by means of the Zeldovich mechanism. The left and right pressure values are 12 kPa and 0.12 kPa, respectively, and the initial mass fractions are 𝑌 N 2 = 0.767, 𝑌 O 2 = 0.233. The computational domain lengths are 𝐿 𝑥 = 1 m and 𝐿 𝑦 = 0.25 m; non-catalytic adiabatic walls are considered at the left, bottom and right boundaries, whereas purely non-reflective conditions are applied at the top end. Three computational grids are considered (4000 × 1000, 6000 × 1250 and 8000 × 2000, respectively), with a constant CFL number of 0.04. Results at the final time 𝑡 = 1 ms are shown in Fig. 12. After the shock reflection on the end wall, a lambda-shock pattern is generated and extends into the bottom boundary layer, where a separation bubble entraining hot gases develops. Kelvin-Helmholtz instabilities are also observed along the slip line originating at the triple point; for a detailed description of the influence of high-temperature effects one might refer to Ref. [START_REF] Chen | Numerical study of chemically reacting flow in a shock tube using a high-order point-implicit scheme[END_REF]. The large temperature values (𝑇 > 4000 K) and small flow velocity near the end wall promote a strong chemical activity; almost all molecular oxygen dissociates and the nitric oxide amount is close to the equilibrium composition at the same 𝑝 and 𝑇 values. Figs. 12e-f show density and temperature profiles along the bottom wall, where the footprints of the vortical structures contained inside the separation bubble are observed; grid-converged results are obtained starting from the 6000 × 1500 grid.

(𝜌, 𝑢, 𝑣, 𝑝) DARU = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ( 120, 0, 0, 120 𝛾 ) 𝑥 < 0.5 ( 1.2, 0, 0, 1.2 𝛾 ) 𝑥 ≥ 0.5 (47)

Compressible Taylor-Green vortex

We consider a compressible extension of the canonical Taylor-Green Vortex (TGV) problem in order to assess the behavior of the numerical scheme when applied to highly under-resolved cases and the suitability of the non-linear artificial dissipation to act as a regularization term in the context of low-and high-𝑀 ILES. The first authors to extend the TGV problem to strongly compressible configurations were Peng & Yang [START_REF] Peng | Effects of the Mach number on the evolution of vortex-surface fields in compressible Taylor-Green flows[END_REF], with the aim of investigating the evolution of vortexsurface fields at large Mach numbers. They considered Mach numbers ranging from 0.5 to 2, albeit at a lower Reynolds number with respect to the one considered in most of the incompressible and low-𝑀 studies (i.e., 400 instead of 1600). Recently, Lusher & Sandham [START_REF] Lusher | Assessment of low-dissipative shock-capturing schemes for transitional and turbulent shock interactions[END_REF] carried out TGV simulations up to 𝑀=1. [START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes[END_REF], in order to evaluate the behavior of different shock-capturing schemes for high-𝑀 turbulent flows. The initial conditions for velocity and pressure fields are:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑢(𝑥, 𝑦, 𝑧) = sin ( 𝑥 𝐿 ) cos ( 𝑦 𝐿 ) cos ( 𝑧 𝐿 ) 𝑣(𝑥, 𝑦, 𝑧) = -cos ( 𝑥 𝐿 ) sin ( 𝑦 𝐿 ) cos ( 𝑧 𝐿 ) 𝑤(𝑥, 𝑦, 𝑧) = 0 𝑝(𝑥, 𝑦, 𝑧) = 𝑝 ∞ + 𝜌 ∞ 𝑈 2 ∞ 16 [ cos ( 2𝑥 𝐿 ) + cos ( 2𝑦 𝐿 )] [ 2 + cos ( 2𝑧 𝐿 )] (48)
When considering the compressible form of the TGV problem, different choices are possible to set the initial conditions for the thermodynamic quantities [START_REF] Peng | Effects of the Mach number on the evolution of vortex-surface fields in compressible Taylor-Green flows[END_REF]. We selected a constant density initial condition, for which the initial density field is equal to 𝜌 ∞ = 1.291834 kg∕m 3 everywhere, whereas the initial temperature is computed from density and pressure (𝑝 ∞ = 1 atm) according to the perfect-gas equation of state. The Reynolds number, Prandtl number and the specific heat ratio are set to 1600, 0.71 and 1.4, respectively. In our test, computational grids ranging from 64 3 to 1024 3 and Mach numbers equal to 0.1, 0.5 and 1 have been considered, in conjunction with a fixed CFL number of 0.5. Periodic conditions are applied in all directions. Comparisons are made by considering the total kinetic energy 𝐾 and the resolved enstrophy 𝛺 integrated over the computational domain:

𝐸 𝑘 = 1 𝜌 ref  ∫  1 2 𝜌𝑢 𝑖 𝑢 𝑖 d, 𝛺 = 1 𝜌 ref 𝑅𝑒 ∫  𝜇 ( 𝜀 𝑖𝑗𝑘 𝜕𝑢 𝑘 𝜕𝑥 𝑗 ) 2 d. ( 49
)
Note that the enstrophy correspond to the solenoidal part of the total viscous dissipation rate; the dilatational component will not be discussed since its contribution has been shown to be negligible also at Mach numbers higher than unity [START_REF] Lusher | Assessment of low-dissipative shock-capturing schemes for transitional and turbulent shock interactions[END_REF]. Fig. 13 shows the evolution of the volume-integrated kinetic energy, solenoidal enstrophy and turbulent kinetic energy spectra at selected times (arbitrarily chosen after the time at which enstrophy peaks), as a function of the initial Mach number and of the grid resolution. Several considerations are in order. First, we observe that the low-𝑀 case is perfectly superposed to classical spectral results, which confirms the good behavior of the numerical scheme for low-speed configurations. The low-order artificial nonlinear dissipation is correctly turned off, whereas the higher-order one guarantees an amount of dissipation such that numerical stability is ensured and virtually all active scales are almost untouched. Increasing the Mach number, the flow dynamics are heavily altered: in the early stages of the decay, internal energy is converted into kinetic energy, such that the average pressure-work can exceed the viscous dissipation rate and the total kinetic energy does not exhibit any more the classical monotonous decline of the incompressible case. The enstrophy peak slightly moves towards later times, whereas its magnitude is almost unchanged. The enstrophy is insufficiently resolved on the 64 3 and 128 3 grids, albeit the largest scales are correctly captured as displayed by the kinetic energy spectra and the 𝑄-criterion visualizations (Fig. 14). In any case, no energy pileup is observed at the smaller scales. The numerical scheme is therefore able to handle efficiently highly under-resolved configurations, while ensuring a correct representation of the largest scales. Grid-converged results are observed for the 512 3 grids; whereas the 256 3 grids allow one to observe a clear trend; that is, a faster convergence of results for larger initial Mach numbers. Enstrophy values are superposed to the reference solution at 𝑀=1, whereas they are slightly underestimated at 𝑀=0.1. This is coherent with the small-scales representation given by the turbulent kinetic energy spectra: although the cutoff wavenumber 𝜅 𝑐 (i.e., the wavenumber at which the dissipation term of the numerical scheme starts to have a large influence on the flow field) does not change for a given grid, the spectral energy contents for 𝜅 > 𝜅 𝑐 tend to be closer to the DNS distributions for larger 𝑀. This behavior can be attributed to the time step constraints imposed by acoustic waves in low-𝑀 configurations: achieving the same final nondimensional time (𝑡 = 16) requires a larger number of time steps, which also implies a larger number of applications of the numerical dissipation operator. This effect (visible only at the smallest, underresolved scales beyond 𝑐 ) underlines the satisfactory selectivity properties of the numerical method.

To illustrate the importance of using a high-order in Fig. 15 we report enstrophy profiles and kinetic energy spectra for a 256 3 grid obtained by changing the order of the non-linear dissipation term of the scheme, from 3rd up to 9th one. Results indicate that (1) the 3rd-order dissipation on the 256 3 grid has an accuracy similar to the 9th-order one on a 64 3 grid, (2) the 9th-order dissipation allows to obtain almost grid-converged results, and (3) the gain of further increasing the dissipation order is almost negligible, even for coarser grids. Of note, results obtained with 3rd-order dissipation in conjunction with higher-order discretizations for the convective fluxes (up to 10th-order) are almost superposed, which confirms the fundamental role of the nonlinear dissipation strategy in obtaining high-quality data. Finally, the effect of the order of the numerical scheme for the viscous fluxes discretization was also assessed; results are reported in Fig. 16 for the 256 3 and 128 3 grids. Consistently with the analysis in the incompressible limit of DeBonis [START_REF] Debonis | Solutions of the Taylor-Green vortex problem using high-resolution explicit finite difference methods[END_REF], a negligible influence is observed when considering discretization orders higher than 4. On the contrary, the use of a 2nd-order discretization produces excessive amounts of viscous dissipation, which for the 256 3 grid is shown to be larger than the one predicted by the DNS. Furthermore, the 128 3 grid with 2nd-order viscous fluxes seems to closely follow the DNS profile. This is rather fortuitous, since the larger enstrophy values come from an incorrect (overestimated) spectral energy repartition at the smallest scales. These results should warn about possible misleading interpretations of the flow physics in case of unsuitable numeric ingredients and/or strongly under-resolved simulations, as well as the need to resorting to spectral analyses for the study and characterization of turbulent flows, the classical volume-integrated quantities being not necessarily accurate enough.

Turbulent boundary layer flows

A crucial aspect for the numerical simulation of wall-bounded turbulent flows is the ability of the numerical scheme to retain good dissipation and dispersion properties near the boundaries. In this section we consider a compressible, calorically-perfect turbulent boundary layer and we carry out a grid-sensitivity analysis to compare results obtained by means of wall-resolved ILES-like grids with respect to reference DNS data. Next, the same numerical strategy is applied to the simulation of an hypersonic turbulent boundary layer of a multi-species, chemically out-of-equilibrium mixture.

Supersonic turbulent boundary layer

The flow conditions for such a configuration are similar to those investigated by several authors [START_REF] Guarini | Direct numerical simulation of a supersonic turbulent boundary layer at Mach 2.5[END_REF][START_REF] Pirozzoli | Direct numerical simulation and analysis of a spatially evolving supersonic turbulent boundary layer at M = 2.25[END_REF][START_REF] Pirozzoli | Characterization of coherent vortical structures in a supersonic turbulent boundary layer[END_REF][START_REF] Pirozzoli | Turbulence in supersonic boundary layers at moderate Reynolds number[END_REF][START_REF] Wenzel | DNS of compressible turbulent boundary layers and assessment of data/scaling-law quality[END_REF]; specifically, a nominal Mach number equal to 2.25, a free-stream temperature of 65 K and a freestream density of 0.13 kg∕m 3 . The fluid considered is calorically-perfect air; Sutherland's law is used to compute dynamic viscosity, along with a constant Prandtl number equal to 0.72. The rectangular computational domain is discretized with even spacing in the streamwise (𝑥) and spanwise (𝑧) direction, whereas a grid stretching is used in the wall-normal (𝑦) direction, following the profile:

𝑦(𝑗) 𝐿 𝑦 = (1 -𝛼) ( 𝑗 -1 𝑁 𝑦 -1 ) 3 + 𝛼 𝑗 -1 𝑁 𝑦 -1 (50)
where 𝐿 𝑦 and 𝑁 𝑦 denote the domain length and the number of grid points, respectively; 𝑗 ∈ [1, 𝑁 𝑦 ] and 𝛼=0.08. The total extent of the domain is 𝐿 𝑥 × 𝐿 𝑦 × 𝐿 𝑧 = 1600𝛿 * in × 100𝛿 * in × 20𝜋𝛿 * in , the initial boundary layer displacement thickness 𝛿 * in being used as length scale. No-slip and isothermal boundary conditions at a temperature close to the laminar adiabatic value (𝑇 𝑤 = 120.18 K) are applied at the wall, whereas characteristic-based conditions are used for the top and right boundaries and periodic conditions in the spanwise direction. A similarity profile is imposed as inlet condition at a distance 𝑥 0 from the leading edge, corresponding to Re 𝛿 * = 1700. Transition to turbulence is triggered by means of a suction-and-blowing forcing method; this excitation technique consists in applying a time-and-space-dependent wall normal velocity disturbance of the form:

𝑣 𝑤 𝑢 ∞ = 𝐴𝑓 (𝑥)𝑔(𝑧) [cos(𝜔𝑡 + 𝛽𝑧) + cos(𝜔𝑡 -𝛽𝑧)] ( 51
)
where 𝐴, 𝜔= ω𝛿 * in ∕𝑐 ∞ and 𝛽= β𝛿 * in represent the amplitude, nondimensional pulsation and spanwise wave number, respectively. Here, the symbol (

•) denotes dimensional values and 𝑐 ∞ is the free-stream speed of sound; 𝑓 (𝑥) and 𝑔(𝑧) are two perturbation-modulation functions defined as in Franko & Lele [START_REF] Franko | Breakdown mechanisms and heat transfer overshoot in hypersonic zero pressure gradient boundary layers[END_REF]. The forcing strip is located near the inlet, at Re 𝛿 * =2000; the parameters prescribed for the current set of simulations are 𝐴=0.025, 𝜔=0.12 and 𝛽=0.2. Of note, the spanwise extent of the domain has been selected in order to contain exactly two oblique waves. Statistics are computed by averaging both in time and in the periodic direction; their collection spans approximately two turnover times and starts after that the initial transient has been evacuated and a statistically-steady state is reached. The sampling time interval is constant and equal to 𝛥𝑡 = 𝛥 t𝑐 ∞ ∕𝛿 * in ≈ 3.93 × 10 -3 .

Table 1

List of the computational grids (and associated legend) selected for the supersonic boundary layer, along with the total number of grid points and resolutions in inner units. 𝛥𝑦 + 𝑤 and 𝛥𝑦 + 𝛿 refer to the wall-normal resolutions at the wall and at the boundary layer edge, respectively. Five different computational grids have been considered and are listed in Table 1. A reference solution is generated by means of a wellresolved DNS and is compared to data obtained on coarser, ILES-like grids, built by progressively halving the total number of points in the streamwise and spanwise directions. On the contrary, the wall-normal grid point distribution is kept unaltered such that discrepancies with respect to DNS data can directly be attributed to excessively coarse discretizations in the other two directions. The least-refined grid has less than 20 million grid points; that is, it is 64 times coarser than the grid used for DNS, which counts more than 1.2 billion grid points. Note that the (•) + quantities listed in Table 1 have been evaluated at the point where the skin friction coefficient peaks, corresponding to the Fig. 16. Influence of the order of the numerical scheme for the discretization of the viscous fluxes for the TGV case at 𝑀=1. Left: enstrophy evolution; right: turbulent kinetic energy spectra at 𝑡=11.1 (zoom). Red color: 256 3 grids; blue color: 128 3 grids. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 17. Resolutions of the computational grids along the streamwise direction in wall units. Line legend as in Table 1. most stringent requirements in terms of resolution. The overall streamwise evolutions of 𝛥𝑥 + and 𝛥𝑧 + for the five computational meshes, reported in Fig. 17, clearly show that their values rapidly grow during transition and become approximately constant in the fully turbulent region. The influence of resolution can be appreciated from the lateral and frontal instantaneous views of the 𝑄-criterion shown in Fig. 18 for all cases but ILES-C, which bears a strong resemblance with case ILES-B. DNS (panels a and e) and ILES-A (panels b and f) are almost indistinguishable; fine-scale structures are properly resolved and the streamwise development of the boundary layers is alike in terms of integral thicknesses. Some visual discrepancies start to appear for case ILES-B, in which both the transition region and the smallest features of the flow are shown to be more grainy. Lastly, the coarser grid of case ILES-D (panels d and h) results in excessive damping of the turbulent motions: hairpin-like structures are smeared out and the occurrence of coherent structures becomes much more sporadic.

Legend

The streamwise profile of the skin friction coefficient 𝐶 𝑓 is displayed in Fig. 19, along with a close-up in the fully-turbulent region of the same quantity as a function of the momentum-thickness-based Reynolds number 𝑅𝑒 𝜃 . DNS and ILES-A results are almost superposed both in the laminar-to-turbulent transition and fully-developed turbulent regions, confirming that grid resolutions of 𝛥𝑥 + < 15 and 𝛥𝑧 + < 10 are fine enough to capture the main flow features. On the contrary, ILES-B and ILES-C cases present some discrepancies in the breakdown to turbulence, the 𝐶 𝑓 peak being slightly smaller and moved towards the domain inlet. A comparison with ILES-A underlines the importance of keeping 𝛥𝑧 + < 10 for the spanwise resolution, whereas the influence of doubling the streamwise mesh size (from 𝛥𝑥 + ≈ 11.5 to 𝛥𝑥 + ≈ 22.5) is shown to be much less pronounced. Lastly, the ILES-D grid results in wrong predictions of both transitional and fully turbulent regions. Despite the much coarser resolutions with respect to the DNS, the fully-turbulent regions tend to match DNS predictions as 𝑅𝑒 𝜃 increases, mainly due to boundary-layer thickening. The 𝐶 𝑓 evolution is indeed in good agreement with values extracted from numerical simulations of similar configurations available in literature [START_REF] Pirozzoli | Characterization of coherent vortical structures in a supersonic turbulent boundary layer[END_REF][START_REF] Pirozzoli | Turbulence in supersonic boundary layers at moderate Reynolds number[END_REF][START_REF] Wenzel | DNS of compressible turbulent boundary layers and assessment of data/scaling-law quality[END_REF], the slight discrepancies being related to different free-stream Mach numbers and forcing locations. In Fig. 20a, the Van-Driest-transformed velocity profiles collapse reasonably well for all the computational meshes except for ILES-D, in which the log-region values are shown to be largely overestimated (in accordance with observations of previous studies about the resolution limits for turbulent boundary layers [START_REF] Sayadi | Large eddy simulation of controlled transition to turbulence[END_REF][START_REF] Poggie | Resolution effects in compressible, turbulent boundary layer simulations[END_REF]). Root-main-square values for each component of velocity, temperature and density are displayed in Figs. 20b, c andd, respectively; here, the same conclusions hold: all the computational grids predict wallnormal profiles very close (or superposed) to the DNS solution, the only exception being ILES-D case which registers large mismatches in proximity of the inner peak. The current grid resolution assessment confirms then the reliability of the numerical strategy in a wall-resolved ILES framework, pointing out that (i) 𝛥𝑥 + < 15 and 𝛥𝑧 + < 10 should be used to obtain DNS-like results, and (ii) slightly coarser resolutions may be considered when focusing the analysis only on fully-turbulent regions.

Hypersonic turbulent boundary layer with finite-rate chemistry

The last test case here considered highlights the capabilities of the numerical code in reproducing turbulent configurations of hightemperature flows, by considering a hypersonic, chemically out-ofequilibrium boundary layer undergoing laminar-to-turbulent transition. It is worth noting that a similar configuration has been recently analyzed from the physical standpoint [START_REF] Passiatore | Finite-rate chemistry effects in turbulent hypersonic boundary layers: A direct numerical simulation study[END_REF]. The thermodynamic conditions are similar to those adopted in several stability studies (see Refs. [START_REF] Malik | Real gas effects on hypersonic boundary-layer stability[END_REF][START_REF] Hudson | Linear stability of hypersonic flow in thermochemical nonequilibrium[END_REF][START_REF] Perraud | Studies of laminar-turbulent transition in hypersonic boundary layers at ONERA[END_REF][START_REF] Franko | Effects of chemistry modeling on hypersonic boundary layer linear stability prediction[END_REF][START_REF] Marxen | Direct numerical simulations of hypersonic boundary-layer transition with finite-rate chemistry[END_REF][START_REF] Miró | Diffusion and chemical non-equilibrium effects on hypersonic boundary-layer stability[END_REF]); specifically, the imposed free-stream values are 50) is used to generate the wall-normal grid distribution, with 𝛼 = 0.13. Outside the boundary layer, air in equilibrium at its free-stream conditions is prescribed (namely, 𝑌 N 2 = 0.767082 and 𝑌 O 2 = 0.232918, similar to Marxen et al. [START_REF] Marxen | Direct numerical simulations of hypersonic boundary-layer transition with finite-rate chemistry[END_REF]), whereas a non-catalytic boundary condition is imposed at the wall. The locally self-similar profile, computed under finite-rate chemistry assumption, is prescribed as inlet condition at a distance corresponding to Re 𝛿 * = 6375. Further details about the computation of the local self-similar solution are reported in Appendix A. Near the inlet, a sponge layer is applied to prevent abrupt distortions caused by the high Mach number. Differently from the supersonic case presented in Section 5.4.1, the forcing is located at Re 𝛿 * = 13880 and the suction-and-blowing forcing function reads:

𝑣 𝑤 𝑢 ∞ = 𝑓 (𝑥)𝑔(𝑧) 𝑁 mode ∑ 𝑚=1 𝐴 𝑚 sin(𝜔 𝑚 𝑡 + 𝛽 𝑚 𝑧), ( 52
)
19. Skin friction coefficient 𝐶 𝑓 as a function of (𝑥 -𝑥 0 )∕𝛿 * in . Sub-figure: 𝐶 𝑓 as a function of 𝑅𝑒 𝜃 in the fully turbulent zone. Line legend as in Table 1. 1.

the collection of 𝐴 𝑚 , 𝜔 𝑚 and 𝛽 𝑚 for each mode being listed in Table 2. The forcing strategy adopted, similar to the one described by Franko & Lele [START_REF] Franko | Breakdown mechanisms and heat transfer overshoot in hypersonic zero pressure gradient boundary layers[END_REF], consists in exciting the two dimensional Mack mode with a strong amplitude. Additional oblique and stationary waves are also imposed at a lower amplitude to promote secondary instabilities. This leads to a late breakdown to turbulence with respect to the supersonic boundary layer case presented in Section 5.4.1, requiring an adequately large computational domain in the streamwise direction. Fig. 21 shows the streamwise evolution of molecular oxygen mass fraction Y O 2 ; it can be clearly seen that the regime is laminar and transitional in almost one half of the entire domain. Statistics have been collected for approximately one turnover time, with a sampling time interval equal to 𝛥𝑡 s 𝑐 ∞ ∕𝛿 * in = 3.45 × 10 -2 . Fig. 22 shows the streamwise evolution of the skin friction coefficient 𝐶 𝑓 along with the compressible extension of Blasius' laminar breakdown to turbulence occurs. On the other hand, the similarity solution, which does not take into account the external pressure gradient created by the viscous/inviscid interaction (see Appendix A), tends to overestimate the skin friction coefficient. Averaged distributions of normalized temperature and species mass fractions along the wall are displayed in Fig. 23a-b, respectively. The extremely high temperature triggers the O 2 dissociation and exchange reactions, and dissociation of N 2 in a smaller amount. As a result, the amount of O 2 decreases as the flows evolves in the laminar region, whereas O and NO increase.

After the breakdown to turbulence, 𝑌 n does not vary significantly due the asymptotic behavior of wall temperature. It is worth noting the inversion of the trend in the fully turbulent region with respect to the laminar evolution, due to the enhanced mixing with the colder, external layers. Although molecular oxygen dissociation continues downstream, the overall mixture composition remains far from the equilibrium one.

Wall-normal distributions of selected first-and second-order statistics are also investigated. The profiles are extracted at a streamwise station close to the domain outlet; the corresponding parameters are reported in Table 3. In Fig. 24a, the Van-Driest transformation for the streamwise velocity is displayed; the logarithmic law is well described by 𝑢 + = 1 0.41 ln(𝑦 + ) + 6.0 (see Ref. [START_REF] Duan | Direct numerical simulation of hypersonic turbulent boundary layers. Part 4. Effect of high enthalpy[END_REF]) and a perfect match with the linear law is found. The wall-normal evolution of averaged mass fractions, shown in Fig. 24b, reveals that most of the atomic oxygen is produced within the viscous sublayer, and tends towards zero at 𝑦 + ≈ 250; nitric oxide (mainly produced by shuffle reactions) follows a similar trend. On the other hand, the mass fraction of atomic nitrogen is almost two order of magnitude smaller than the other reaction products, since the temperatures reached do not trigger significant N 2 dissociation. The registered amount of N is therefore mainly due to the recombination reaction between N 2 and O. Second-order statistics are displayed in Fig. 25; panel (a) shows that root-main-square value of temperature doubles the free-stream value at the peak of production, whereas density fluctuations stay relatively small. Accordingly, the fluctuating distributions of species mass fractions (Fig. 25b) exhibit the largest values where 𝑇 rms peaks, although the strongest chemical activity is concentrated at the wall where the temperature reaches the largest values.

Conclusions

In this study, an efficient methodology for the numerical simulation of turbulent flows at high speeds has been presented. The main difficulty associated to such flows is related to the amount of numerical dissipation to be added, which must satisfy the concurrent needs of obtaining a reliable picture of small-scale turbulent motions, while being able to handle strong discontinuities.

The proposed numerical scheme represents an high-order extension of the original artificial-viscosity-based method of Jameson [START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes[END_REF]. It relies on a blending of second-and tenth-order derivatives equipped with a suitable shock detector, which results in a ninth-order-accurate scheme in smooth flow regions. Several test cases have been considered to assess the capability of such technique to deal with shocks, shock/entropy perturbations and turbulence, both for calorically- perfect gases and high-temperature multicomponent flows. Different grid refinements (representative of typical resolutions used to carry out implicit Large-Eddy simulations) were analyzed for multi-scale configurations; results were systematically compared to DNS data to quantify the error introduced by the adoption of under-resolved meshes.

A parametric study on the isentropic vortex advection was carried out and the formal orders of accuracy of the spatial discretization and temporal integration were correctly retrieved. The four shock tube configurations considered (Sod, Lax, Grossman-Cinnella and Daru-Tenaud) confirm that the scheme handles properly strong shocks and contact discontinuities, even in presence of chemical nonequilibrium processes. Furthermore, the two-dimensional N 2 -O 2 under-expanded jet, the compressible Tayor-Green Vortex and the supersonic boundary layer were used to prove the reliability and robustness of the numerical method when applied on coarse grids for ILES. It was shown that results tend seamlessly to DNS predictions as the resolution increases, and computations were stable also for severely under-resolved simulations. The last configuration investigated, a turbulent boundary layer at hypersonic speeds, demonstrates the good behavior of the numerical with 𝑓 ′ = 𝑢∕𝑈 𝑒 , one can manipulate the boundary-layer equations and obtain their formulation in the self-similar coordinate system: Here, 𝑔 = ℎ∕ℎ 0,𝑒 is the self-similar parameter for the enthalpy, ℎ 0,𝑒 being the edge stagnation enthalpy, 𝜃 = 𝑇 ∕𝑇 𝑒 and 𝐶 = 𝜌𝜇∕𝜌 𝑒 𝜇 𝑒 . For the sake of clarity, in case of frozen chemistry or chemical equilibrium, the only 𝜉-dependent term in Eq. (A.9) vanishes and the resulting equations become globally self-similar. The previous system of equations can be integrated numerically, subjected to the following boundary conditions:

( 𝐶𝑓 ′′ ) ′ + 𝑓 𝑓 ′′ = 0 (A.
𝑓 ′ = 𝑓 = 0, 𝑌 ′ 𝑛 =

Fig. 1 .

 1 Fig. 1. Dispersion and dissipation errors as a function of the reduced wavenumber 𝑘𝛿𝑥 up to the ninth-order scheme. (a) dispersion error 𝛷, (b) dissipation error 𝐷.

0 ( 41 )Fig. 3 .

 0413 Fig. 3. Analysis of the order of convergence of the error norm for the isentropic vortex convection configuration, for different values of the perturbation amplitude (from 1 to 1.9) at CFL = 0.025 and 𝑘 2 = 1. (a): initial density profile at 𝑦=0; (b): error of the 𝐿 2 -norm as a function of the number of grid points.

Fig. 4 .

 4 Fig. 4. Reference and numerical solutions for Sod's (left) and Lax's (right) shock tubes. Triangles: single species, circles: N 2 -O 2 mixture. The gray solid line denotes the reference solution.

Fig. 5 .

 5 Fig. 5. Numerical solution for the Shu-Osher problem at 𝑡 = 1.8. Top: profiles of density, pressure, velocity and entropy; bottom: zoom on the post-shock regions. In each subfigure, the black solid lines denote the reference solution, and the colored solid and dashed lines the present solution, respectively with and without Bhagatwala & Lele's correction.

Fig. 6 .

 6 Fig. 6. Results for the shock-vorticity/entropy wave interaction at 𝑡 = 25. (a) Instantaneous vorticity profiles for 𝑦 = 0 (𝑘 𝑦 = 1); averaged vorticity content, normalized with respect to the upstream value, for (b) 𝑘 𝑦 = 1 and (c) 𝑘 𝑦 = 2. The shock location is 𝑥 𝑠 = 3𝜋∕2; the dashed horizontal lines in panels (b) and (c) denote the linear analysis solution [75].

Fig. 7 .

 7 Fig. 7.Profiles of density, pressure, velocity and mass fractions for the reacting, real gas shock tube case. In panels (a), (b) and (c), quantities are made non-dimensional with respect to the corresponding initial values on the right side of the tube. Lines: current simulation; symbols: reference solution[START_REF] Grossman | Flux-split algorithms for flows with non-equilibrium chemistry and vibrational relaxation[END_REF].

Fig. 8 .

 8 Fig. 8. Instantaneous snapshot of an N 2 -O 2 inert underexpanded jet. From left to right, top to bottom: isocontours of density 𝜌, temperature 𝑇 , numerical Schlieren log(∇𝜌∕𝜌 𝑗 ), vorticity 𝜔, Mach number 𝑀 and pressure fluctuations 𝑝∕𝑝 𝑗 -1.

Fig. 9 .Fig. 10 . 2 [

 9102 Fig. 9. Influence of the grid resolution on the underexpanded jet. Axial profiles of Mach number (a); axial profiles of normalized temperature and density (b); rms profiles of temperature along the vertical line 𝑥∕𝐷 = 20 (c). In panel (a), the vertical lines denote the streamwise locations of the Riemann wave observed by Martínez Ferrer et al. [77] ( ), Su et al. [78] ( ) and Sheeran & Dosanjh [80] ( ). In panel (b), black circles represent data extracted from Su et al. [78].

Fig.

  Fig. Results for the two-dimensional perfect-gas viscous shock tube, for 𝑅𝑒 = 1000 at 𝑡 = 1. (a): density distribution along the bottom adiabatic wall (𝑦 = 0); (b): density distribution at 𝑦 = 0.05.

Fig. 12 .

 12 Fig. 12. Results for the two-dimensional chemically-reacting viscous shock tube at 𝑡 = 1 ms. Instantaneous visualizations of (a) normalized density, (b) normalized temperature, (c) O 2 mass fraction and (d) NO mass fraction. Panels (e) and (f) show the normalized density and temperature distributions, respectively, along the bottom adiabatic wall for the different computational grids considered.

Fig. 13 .

 13 Fig. 13. Influence of the grid resolution for the Taylor-Green Vortex case. Top: 𝑀=0.1; center: 𝑀=0.5; bottom: 𝑀=1. Left and center columns: timewise evolutions of volumeintegrated kinetic energy 𝐾 and resolved solenoidal enstrophy 𝛺; right column: turbulent kinetic energy spectra 𝐸(𝜅) at selected times (𝑡 = 8.9 for 𝑀=0.1, 9.4 for 𝑀=0.5 and 11.1 for 𝑀=1).

Fig. 14 .

 14 Fig. 14. Isocontours of the 𝑄-criterion colored with local Mach number for computational grids ranging from 64 3 to 512 3 (𝑀 = 1 case).

Fig. 15 .

 15 Fig. 15. Influence of the order of the nonlinear dissipation term for the TGV case at 𝑀=1 for 256 3 grid. Left: enstrophy evolution; right: turbulent kinetic energy spectra at 𝑡=11.1 (zoom).

Fig. 18 .

 18 Fig. 18. Lateral and frontal views of instantaneous snapshots for the 𝑀 = 2.25 supersonic boundary layer. Isosurfaces of the 𝑄-criterion colored by the distance from the wall. (a, e): DNS; (b, f): ILES-A; (c, g): ILES-B; (d, h): ILES-D.

Fig. 20 .

 20 Fig. 20. Wall-normal profiles of the Van-Driest-transformed longitudinal velocity (a), r.m.s. of velocity (b), temperature (c) and density (d), extracted at 𝑅𝑒 𝜃 = 3500. Line legend as in Table1.

Fig. 21 .Fig. 22 .

 2122 Fig. 21. Instantaneous visualization of a 𝑥𝑧-plane extracted at 𝑦∕𝛿 * in = 12, colored by molecular oxygen mass fraction Y O 2 .

Fig. 23 .

 23 Fig. 23. Streamwise evolutions of (a) normalized wall temperature and (b) mean mass fractions 𝑌 𝑛 along the wall. In panel b, 𝑌 N 2 is not shown being outside of the prescribed range.

Fig. 24 .

 24 Fig. 24. Wall-normal profiles extracted at 𝑅𝑒 𝜃 = 10000 of (a) Van-Driest transformed longitudinal velocity, and (b) mass fractions of selected species.

Fig. 25 .

 25 Fig. 25. Wall-normal profiles extracted at 𝑅𝑒 𝜃 = 10000 of (a) normalized r.m.s. temperature and density, and (b) r.m.s. species mass fractions.

Fig. A. 26 .

 26 Fig. A.26. Evolution of normalized and normalized (left) and species mass fractions (right), at the inflow. 𝑌 N 2 and 𝑌 O 2 are not shown being outside the range.Independent variable transformations are introduced as follows (see, for instance, Ref.[START_REF] Lees | Laminar heat transfer over blunt-nosed bodies at hypersonic flight speeds[END_REF]):𝜉 = 𝜌 𝑒 𝜇 𝑒 𝑈 𝑒 𝑥 = 𝜉(𝑥), 𝜂 = 𝑢 𝑒 √ 2𝜉 ∫ 𝑦

  configuration reported in Section 5.4.2, the prescribed edge variables are 𝑈 𝑒 = 3754 m∕s, 𝑇 𝑒 = 350 K and ℎ 0,𝑒 = 7.1 MJ∕kg. Considering a distance from the leading edge equal to 𝑥 = 0.01 m, we obtain the inflow profiles of normalized temperature, normalized velocity and species mass fractions, shown in Fig. A.26 as a function of the incompressible similitude variable 𝜂 𝑖 = 𝑦 √ 𝑅𝑒 𝑥 ∕𝑥.

Table 2

 2 Parameters of the modes excited by the forcing function in Eq. (52): nondimensional amplitude, frequency and spanwise wave number.

	Mode	1	2	3	4	5	6 7	8 9	10 11 12 13
	𝐴 𝑚 (×10 -3 ) 50	2.5	2.5 2.5 2.5 2.5	2.5	2.5	2.5
	𝜔 𝑚	1.71 0.855 0	0	0	1.71	0.855	1.71	1.71
	𝛽 𝑚	0	0	0.2 0.4 0.6 ±0.2	±0.2	±0.4	±0.6
	correlation, 𝐶 𝑓 ,𝑙𝑎𝑚 =	0.664 √ 𝑅𝑒 𝑥	√	𝜌 𝑤 𝜇 𝑤 𝜌 ∞ 𝜇

∞ , as well as the results predicted by the locally self-similar theory. After the forcing strip (not visible in the figure) and despite the presence of chemical reactions, the evolution of 𝐶 𝑓 stays close to the laminar correlation up to (𝑥-𝑥 0 )∕𝛿 * in ≈ 3300, where

Table 3

 3 Boundary layer properties at a selected station in the fully turbulent region. In the table, 𝑅𝑒 𝜏 = 𝜌 𝑤 𝑢 𝜏 𝛿 99 ∕𝜇 𝑤 denotes the friction Reynolds number, 𝑀 𝜏 = 𝑢 𝜏 ∕𝑐 𝑤 stands for the friction Mach number and 𝐻 = 𝛿 * ∕𝜃 is the boundary layer shape factor.

	(𝑥 -𝑥 0 )∕𝛿 * in	𝑅𝑒 𝜃	𝑅𝑒 𝛿 * × 10 -4 𝑅𝑒 𝜏	𝑀 𝜏	𝐻	𝛥𝑥 +	𝛥𝑦 + 𝑤	𝛥𝑦 + 𝛿	𝛥𝑧 +
	7377	10000 29.83	236 0.16 31 4.51 0.51 3.82 4.07

Acknowledgments

This work was granted access to the HPC resources of IDRIS and TGCC under the allocation 2019-2B10947 made by GENCI (Grand Equipement National de Calcul Intensif). D.P. and G.P. were partially supported by the Italian Ministero dell'Istruzione, dell'Università e della Ricerca under the Programme Department of Excellence Legge 232/2016 (Grant No. CUP-D94I18000260001).

scheme also with the concomitant occurrence of shocklets, broadband turbulence and chemical nonequilibrium processes.

Future works have been planned concerning the extension of the current scheme to curvilinear meshes, the coupling with skewsymmetric formulations, and the numerical investigations of high-speed configurations encompassing thermal relaxation phenomena by means of multi-temperature models.

CRediT authorship contribution statement

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Locally self-similar solution for laminar boundary layers

The boundary-layer equations for a steady, compressible, multicomponent, reacting, two-dimensional flow without pressure gradient write: