Simulation of High-Enthalpy Turbulent Shock Wave/Boundary Layer Interaction Using a RANS Approach

In the era of space exploration, the scientific community is strongly focusing on the analysis of hypersonic flows in the presence of shock wave/boundary layer interaction. In these conditions, the flow field presents a complex shock structure due to the interaction of different shock waves with the boundary layer. The strong adverse pressure gradient makes the boundary layer separate, giving rise to a separation bubble. In the reattachment zone, the temperature can reach very high values, inducing thermochemical non-equilibrium effects. This research field is recently achieving more and more relevance in aerospace research, as the analysis of turbulent shock wave/boundary layer interaction so far has been mainly focused on perfect gas flows. In this manuscript, a Reynolds averaged Navier–Stokes (RANS) approach is considered, the shear stress transport (SST) model being coupled with the multitemperature approach proposed by Park to investigate thermochemical non-equilibrium effects in hypersonic turbulent shock wave/boundary layer interaction. The first part of the manuscript is devoted to the validation of the solver, and results for low enthalpy flat plate and compression ramp flows are presented. The numerical results are shown to be in good agreement with numerical solutions and experimental measurements. Afterward, the free stream conditions are modulated to make non-equilibrium relevant and analyze a reacting flow.


Introduction
The scientific community is currently devoting a lot of efforts to the numerical investigation of turbulent boundary layers, especially in the context of hypersonic flight, due to its impact for several breakthrough applications such as military or space exploration [1].At hypersonic speeds, most of the kinetic energy of the flow is converted to internal energy, leading to excitation of internal modes of energy (translational, rotational, vibrational and electronic).Moreover, the high temperature reached ( ≈ 10000 K) induces molecu- lar dissociation and ionization, making the prediction of the flow feature very challenging [2].Most of the times, experiments are very hard to carry out due to the difficulties encountered in the physical reproduction of the flight conditions [3], making numerical simulations a valid alternative for the accurate evaluation of the main flow field characteristics (chemical composition, shear stresses, heat flux).
In this scenario, turbulence introduces more and more difficulties when solving the Navier-Stokes equations for simulating hypersonic re-entry.A direct numerical simulation (DNS) represents the most accurate and high-fidelity approach for the correct evaluation of the flow features, as all the turbulence scales are resolved.However, the huge computational cost and high order schemes required by such simulations often limit a DNS approach to restricted Reynolds numbers and simple configurations such as flat plate or wall-bounded flows.A wide field of numerical studies is available in the literature for both supersonic [4][5][6][7] and hypersonic [8][9][10][11][12] regimes, also accounting for hightemperature effects [13][14][15].These works focused on the validity of Morkovin's hypothesis [16], which states that Francesco Bonelli and Giuseppe Pascazio contributed equally to this work.compressibility effects can be taken into account considering the variation of mean fluid density; the velocity transformation obtained by means of the law of the wall proposed by Van Driest [17][18][19] is coherent with the incompressible version also for moderate wall and boundary layer edge temperature ratios [9].However, the assumption of an isothermal, cooled surface is essential in the context of hypersonic flight as it is more realistic [20,21].Many authors found that Morkovin's hypothesis breaks down in the presence of relevant wall cooling [9,22].
One of the most relevant phenomena in turbulence is the shock wave/boundary layer interaction (SWBLI), whose analysis is fundamental for a proper design of the aerodynamic shape of space vehicles.It is typical of complex configurations such as ramp, double-wedge or double-cone geometries: the boundary layer interacts with the shock wave inducing the laminar-turbulent transition [23][24][25].The SWBLI can be forced by employing a shock generator which makes the wave interact with the flow over a flat plate, provoking the turbulent transition [25][26][27][28].However, natural SWBLI interaction is induced by the wall deflection, typical of compression ramp geometries.Many researchers have analyzed the turbulent shock wave/boundary layer interaction in high-speed compression ramp flows [29][30][31][32].Given the high computational cost required by compression ramp simulations, these works are mainly based on the use of turbulent models, with some exception [33,34].In the framework of hypersonic flows, non-equilibrium conditions have been investigated for different configurations [35][36][37], also by means of modal analysis [38,39].
This work aims at assessing the capabilities of an inhouse GPU based solver to investigate hypersonic flows in the presence of thermochemical non-equilibrium and turbulent shock wave/boundary layer interaction.Turbulence is modeled by a RANS approach in conjunction with the well-known two-equation shear-stress-transport ( k − SST) turbulence model [40].Such an approach allows to switch from the classical k − model in the far field region to the k − model near the wall, combining the characteristics of both models.Non-equilibrium is modeled through the multitemperature model proposed by Park [41], hence considering the high temperature effects making the calorically perfect gas assumption no longer valid.Therefore, the mixture is composed of five species (N 2 , O 2 , NO, N, O) interacting with each other.A continuity equation is solved for each species in the mixture and, also, a transport equation for the vibrational energy of molecules is needed.For this reason, an MPI-CUDA environment is developed to exploit multi-GPU executions and drastically reduce the computational cost [42].
The paper is organized as follows: firstly, the methodology is illustrated, giving details about the governing equations, models and numerical scheme.Afterward, the numerical results are discussed for perfect gas flows, to isolate the turbulence phenomena and assess the effectiveness and accuracy of the model.More in detail, a flat plate and compression ramp flows are investigated.In the final part, thermochemical non-equilibrium conditions are analyzed for the flow over a compression ramp.These must be considered as preliminary results, useful as benchmark for future works.

Governing Equations
The system of Navier-Stokes equations is coupled with the SST turbulence model [40] and with the multitemperature model proposed by Park [41] to deal with thermochemical non-equilibrium.The SST model combines two classical turbulence models, namely, the k − and the k − model, and it is widely used in the literature [37,43,44].In the following, the equations describe the mean flow field obtained by performing the Favre (density-weighted) averaging: for a given variable , = − � indicates the standard time average, with ′ the corresponding fluctuation, whereas ̃ = − �� = ∕ denotes the density-weighted Favre averaging, with ′′ the Favre fluctuation.However, to simplify the notation, overbars and tildes are omitted.To model turbulence high-temperature effects, the model proposed by Jiang et al. [37] has been implemented with few modifications.
Dealing with a multispecies mixture, a continuity equation is solved for each species: where subscript s indicates the s-th species in the mixture and j s is the diffusive flux which accounts for mass diffusion phenomena.This is modeled by means of Fick's law: with T the turbulent viscosity, D s the equivalent diffusion coefficient [48] and Sc T = 0.9 the turbulent Schmidt number (equal for all the species).The turbulent viscosity is defined as follows [45]: with a * = 0.31 .k T is the turbulent kinetic energy, the spe- cific dissipation rate and S the mean strain rate tensor, given by the symmetric part of the gradient of the velocity vector u In the definition of T , the term F 2 is with d being the distance from the wall and * = 0.09 .The source term in Eq. ( 1), ωs , is expressed using the law of mass action.More details are given in [42,46].Momentum balance reads: where L i,j and T i,j denote the viscous and turbulent stress tensors, given by the following expressions: leading to the definition of the "total" stress tensor, given by: where i,j is the Kronecker delta and k T the turbulent kinetic energy.The molecular viscosity L is computed by means of Wilke's mixing rule [47], which combines single species properties evaluated using Gupta's fitting [48].The introduction of the turbulent viscosity arises from the Boussinesq assumption: the stresses deriving from the turbulence behave as the molecular ones, but are linked to the flow property ( T ) instead of depending on a fluid property ( L ).The term − 2 3 k T i,j represents a closure term for a consistent definition of the trace of the Reynolds stress tensor [40].
The total energy conservation accounts for the turbulent stress tensor, the turbulent heat flux and the turbulent transport: with k = 0.85 [49] and E being the total energy content, accounting also for the turbulent kinetic contribution.Before going into the details of the turbulence modeling, one should recall that, dealing with non-equilibrium flows, the total heat flux is given by three contributions: (5) where N m and N s are the number of molecules and species in the mixture, whereas h s is the total specific enthalpy of spe- cies s.Molecular (laminar) thermal conductivity of the mixture ( L ) is calculated through classical mixing rules [50], whereas Eucken formulation [51] is employed for the evaluation of the vibrational conductivity of the molecule m ( vib m,L ).Concerning the turbulent heat flux, q T , a definition simi- lar to the one of the turbulent stress tensor is used : it still depends on the temperature gradient of the mean field, but a "turbulent conductivity" is introduced as a function of the turbulent viscosity.First of all, recall that the thermal conductivity is = c p Pr , c p being the constant pressure specific heat and Pr the Prandtl number (in this definition, the energy modes are neglected for simplicity).In this way, one can define the "turbulent conductivity" as: assuming Pr T = 0.9 .The vibrational conductivity is written for each molecule as a convention.Furthermore, note that c vib m does not refer to constant pressure or constant volume specific heat, since e vib = h vib .Moreover, to account for tur- bulent mass diffusion, one can introduce the "turbulent diffusivity" as in the species continuity equations.In this way, the total turbulent heat flux is given by: Since a multitemperature model is employed to deal with thermochemical non-equilibrium, a transport equation is solved for each vibrational energy: where the subscript m refers to the m-th molecule in the mixture.In this equation, the heat fluxes are expressed as: Details about the source term ωvib m are given in [42,46].Finally, two additional equations are integrated for modeling the turbulence.These read: (10) and represent the evolution of the turbulent kinetic energy ( k T ) and specific dissipation rate ( ).Formally, these equa- tions are very similar to those previously described, but they present source terms which are responsible for promoting the turbulence: term P identifies production terms, whereas term D identifies dissipation terms.More in detail, [49].The limiter on P k avoids the buildup of turbulence in the stagnation regions [45].Sometimes, , which is correct for incompressiblw flows, but is still a good approximation also for compressible flows [49].Concerning Eq. ( 16), source terms are expressed as where T = T ∕ is the turbulent kinematic viscosity.and appearing in the production and dissipation terms are functions blended by F 1 , namely the blending function (varying from 0 to 1) [45]: with d the distance from the wall and the constants 1 = 5∕9 , 2 = 0.44 , 1 = 0.075 , 2 = 0.0828 and ,2 = 0.856 [37,49].Moreover, in Eq. ( 16), an additional source term appears, namely the cross diffusion C .This term is respon- sible for making the model adapt itself to k − model in the far field or to k − model in near-wall regions.This term reads:

Boundary Conditions
Supersonic outflow and inflow boundary conditions are applied on the conservative variables of the Navier-Stokes equations.Concerning turbulence variables, their values must be guessed.A typical way is to compute the turbulent kinetic energy as a function of the turbulent intensity [40] where the double quotes indicate the fluctuations of the velocity.Hence, the free stream inflow condition for k T is: V being the velocity magnitude and I the turbulent intensity.The latter can assume different values, depending on the case [40,49]: as a first attempt, I = 1% has been used in this work.Based on the relation between k T and , the free stream value of the specific dissipation rate is: Concerning the boundary conditions, the turbulent kinetic energy k T must be null at wall which implies that T = 0 at wall.The wall boundary condi- tion for is [52]: L being the kinematic viscosity, d the distance of the first point from the wall, and 1 = 0.075 [40].

Numerical Scheme
The governing equations are integrated using a finite-volume scheme implemented in a body-fitted solver.Advection terms are discretized by means of the well-known Steger-Warming flux vector splitting [53], with the right and left state reconstructed through a second-order MUSCL reconstruction [54,55].Viscous fluxes are evaluated by means of a second-order centered scheme.Lastly, the time integration is performed in two steps: in the first step, an explicit third-order Runge-Kutta scheme is employed for advancing the homogeneous equations [56]; in the second step, source terms are computed thanks to an implicit iterative Gauss-Seidel scheme [57].Such an approach (splitting approach) is suitable for reacting flows and it was shown to provide results comparable to fully coupled implicit schemes [46].The main idea is to employ the most suitable algorithm for treating the different time scales of fluid dynamics and thermochemistry.Solving the coupled equations would drastically reduce the time step size due to the stiffness of the thermochemistry source terms.Instead, the use of a splitting approach allows one to employ an explicit time integration for the fluid dynamics, which involves fewer floating point operations and results in an efficient strategy for GPU implementation.On the other hand, the implicit evaluation of the thermochemistry source terms makes the convergence faster and, involving point operations, maintains the effectiveness of GPU implementation.

Implicit Formulation for the Source Terms
The solution vector Ũ = k T ,  T is taken as an example.
The tilde indicates that this vector is advanced in time starting from the solution obtained by the first step, namely, from the homogeneous equations, to be updated with the source terms.Hence, one needs to solve the system of equations given by with P and D a vector and a diagonal matrix whose com- ponents are nonnegative and represent the production and dissipation terms for a certain quantity.For turbulence equations, these read: ,2 1 k T x j x j , Note that the cross-diffusion term behaves as a production term and has been incorporated into P.For an implicit formulation, the solution vector in the right-hand side of Eq. ( 28) must be evaluated at the 'new time'.So, indicating with k the k-th inner iteration and with n the n-th time step, one has: The same approach is applied to evaluate thermochemical source terms.In this work, eight inner iterations are considered for turbulence and four iterations for thermochemistry.This set of inner iterations turned out to provide an effective convergence.Moreover, since no threshold value is imposed (31) .
for the residuals, GPU threads are ensured to be synchronized during the calculation.

Hypersonic Flow over a Flat Plate (Perfect Gas)
Perfect gas assumption helps in isolating turbulence, neglecting high-temperature effects.Hence, a Mach 7 flow over a flat plate is taken as an example [37,58,59].The sample is 3.2 m long, and two grids were tested to analyze the convergence of the solution.The first computational domain is composed of 100 × 200 nodes, with a wall resolution of 5 × 10 −6 m; the second one is composed of 200 × 400 nodes, with a wall resolution of 2.5 × 10 −6 m.They provide the same results as shown in Fig. 1(a), though a little deviation appears in the transition region.However, this model is not able to capture the laminar-to-turbulent transition, making more relevant the downstream region.The comparison of the wall heat flux with the numerical reference [37] and experimental findings [59] for both laminar and turbulent regimes is shown in Fig. 1(b).The agreement is very satisfactory, even if some discrepancies can be observed.These are attributed to some difference in the turbulence model employed in [37].

Hypersonic Flow over a Compression Ramp
In this section, the results for the flow over a compression ramp is discussed.The first test aims at the validation of the code.For this purpose, the experiment conducted by Coleman and Stollery [23] is simulated.In such a case, the low free stream enthalpy makes the perfect gas assumption suitable.Once the validation is assessed, the free stream enthalpy is increased to investigate the high-temperature effects.The geometry considered is the same for both cases: a sketch reporting the dimensions of the domain and boundary conditions is illustrated in Fig. 2. The wall is isothermal at T w = 295 K.
In accordance with references [37,43,44], the solver employed a mesh with 100 × 200 control volumes, distrib- uted so that the height of the first cell at wall is 5 × 10 −7 m.

Low-Enthalpy Flow (Perfect Gas)
The test gas is pure nitrogen with free stream density of ∞ = 0.133894 kg∕m 3 and temperature T ∞ = 64.5 K .The free stream Mach number is M ∞ = 9.22 , whereas the free stream Reynolds number (per unit length) is Re ∞ = 4.7 × 10 7 1∕m .The aforementioned wall resolution results in a value of y + of 0.08 near the compression corner.The results are illustrated in Fig. 3, where the profiles of the normalized pressure ( p w ∕p ∞ ) and heat flux ( q w ∕q ∞ , with q ∞ = 6.29 W/cm 2 [58]), are reported and compared with numerical references from the literature [37,43,44] and experimental findings.
The results reported by the present simulations are in good agreement with references, with slight underestimation of pressure and heat flux downstream of the corner.Specifically, the pressure profiles are very smooth and present an overall constant value downstream of the shock; contrarily, observing the experimental measurements a peak would be expected.The most critical aspect is the absence of the heat flux jump near the corner ( x ≈ −0.3 m ): the trend found by Zhang [43] is even opposite to the experiments, with the heat flux decreasing toward the wall deflection.Besides these comments, the results are in a good agreement with the expectations.
It is noteworthy that these results were obtained imposing P � k = T S 2 in Eq. ( 17), which provides good agreement with the experiments.In fact, within the context of the closure problem in the RANS approach, one of the issues in highly compressible flows is represented by the modeling of the production term of the turbulent kinetic energy equation [60] that is particularly critical in hypersonic flows with shock waves [61].In Fig. 4, wall quantity profiles are shown for the aforementioned grid and a coarser one (nodes halved in both directions) with P � k = T S 2 (source 1) and (source 2).Important differences are observed.It is fundamental to highlight that such differences do not appear for the flat plate test case: as a matter of fact, deviations between the two approaches take Table 1 Free stream conditions for the high-enthalpy test case place from the separation bubble on, whereas they provide the same results in the upstream field (flat plate zone).Furthermore, the authors have experienced similar behavior on finer meshes, but not on coarser ones where the separation point remains unchanged, as reported in the right-hand side of Fig. 4.This behavior, probably linked to numerical modeling issues, needs further investigations which go beyond the scope of the preliminary assessment of this work.
To better understand the wall quantity profiles, Schlieren images were extrapolated from the flow field.They are shown in Fig. 5(a) together with the zoom of the x-velocity contour in the compression corner (Fig. 5(b)) and the contour of the source term P k , normalized by ∞ u 3 ∞ (Fig. 5(c)) (medium mesh).The black lines in figures (b) and (c) identify the null values.A separation shock is formed near the compression corner due to the interaction of the ramp oblique shock with the boundary layer on the flat plate.Hence, a separation region forms downstream of the so-called separation point, where the flow is subsonic.While moving along the deflected wall, the flow reattaches and is again supersonic.The prediction of separation point location is fundamental in view of a full understanding of the involved phenomena.
The source 2 induces a stronger compressibility of the flow due to a stronger separation shock: this is evident also from the pressure peak in Fig. 4(a).Moreover, the separation extent is much larger with respect to the source 1 case: the zoom of the x-velocity puts in evidence how the boundary layer separation is markedly larger and starts more upstream in case of the source 2, pushing the reattachment toward the trailing edge of the ramp.Finally, a completely different behavior is observed for P k in the source 2 case, which presents large zones with negative values.Also, turbulent kinetic energy production is higher in the case of source 2.
It is fundamental to underline that the production term should be non-negative.De facto, P � k = T S 2 is non-negative, but in the case of source 2 . Performing the algebra yields where the second term is responsible for the change of sign in the global production term.The comparison of these two contributions is reported in Fig. 6.The larger separation encountered in the case of the source 2 induces a larger turbulent viscosity ratio ( T ∕ L ), as depicted in Fig. 7.

High-Enthalpy Flow (Non-equilibrium)
The numerical simulation of the high-enthalpy flow exploits the same computational domain of the low-enthalpy case.However, to induce relevant non-equilibrium, the free stream conditions are changed as reported in table 1. Adiabatic and no-slip boundary conditions are imposed at wall.It is worth highlighting that the adiabatic condition makes the temperature reach a very high value at the wall: specifically, for this test case it is about 6000 K.Such high values are never (32) Schlieren images, x-velocity and normalized P k contours are shown in Fig. 8.The separation bubble is bigger if predicted by the source 2 (see Fig. 8(b)).However, the difference is less relevant with respect to the low enthalpy case: indeed, the contour maps of the turbulent kinetic energy production term are quite similar in the two cases (see Fig. 8(c)).This is potentially attributable to the adiabatic thermal boundary condition: indeed, wall temperature gradients can induce high compressibility effects, which are weaker in case of adiabatic walls.Numerical investigation of isothermal, cooled wall flows is a challenging task due to the uncertainty about turbulence term closure.This goes beyond the scope of the preliminary analysis of this work and is still an ongoing task.
Similarly to the low-enthalpy case, the two contributions of the turbulent kinetic energy production term are reported in Fig. 9 for the source 2, confirming the change of sign is given by P k2 .Also, the turbulent viscosity ratio is shown in Fig. 10.
Wall quantity profiles (pressure and shear stress) are shown in Fig. 11.Again, the source 2 leads to an overestimation of the separation bubble, while underestimating the wall shear stress.
Given this opposite behavior of the two expressions of P ′ k , it is fundamental to understand which are the consequences on the chemical dissociation activity.Contours of N 2 and O 2 mass fractions and vibrational temperatures are shown in Fig. 12.The pictures are zoomed near the compression corner to better observe the differences between the two approaches.First of all, the stronger compressibility effects promoted by the source 2 are evident also on the chemical activity, as the dissociation is more pronounced in such a case.Both molecular nitrogen and oxygen dissociation are more relevant in the case of source 2; moreover, they are anticipated given the bigger size of the separation bubble.
However, no relevant differences are observed in the contours of the vibrational temperatures: besides the different separation bubble extent, the values are mostly comparable in both cases.

Conclusions
In this paper, preliminary numerical results for hypersonic turbulent boundary layer have been presented.Turbulent shock wave/boundary layer interaction has been simulated, underlining the influence of high-enthalpy effects.Turbulence has been modeled by means of the standard SST k − model, which well describes the physics of the boundary layer in the presence of relevant curvature of the geometry.The results showed good agreement with experimental findings and numerical references in the literature.However, they are strongly affected by the expression of the turbulent kinetic energy production term.Depending on how this term is modeled, compressibiltiy effects and chemical activity change, leading to a different separation bubble extent in the case of a compression ramp flow.However, this assessment needs further investigations.
Finally, in the context of hypersonic turbulent boundary layer studies, a tuning of the classical turbulence models has to be devised, in order to improve the quality of the solution obtained by numerical simulations with no need of performing a DNS, whose practical applications are limited due to their huge computational cost.Needless to say DNS studies are essential to provide a sufficiently wide range of data to make turbulence models training feasible in the most common configurations, such as flat plate and shock wave/ boundary layer interaction.

Fig. 1
Fig. 1 Heat flux profiles along the wall of the flat plate obtained on the two grids (a) and comparison of the results with reference profiles (b)

Fig. 2 Fig. 3
Fig. 2 Computational domain with dimensions and boundary conditions

Fig. 4
Fig. 4 Comparison of pressure (a), shear stress (b) and heat flux (c) profiles for different expressions of the turbulent kinetic energy production term: medium mesh (left) and coarse mesh (right)

Fig. 5
Fig. 5 Schlieren images (a), x-velocity contour in the compression corner (b) and contour of the normalized source term P k (c) for the low-enthalpy case.Black lines in (b) and (c) identify the null values

Fig. 6 Fig. 7
Fig. 6 Contour maps of the two contributions of the turbulent kinetic energy production term for the low-enthalpy case

Fig. 8
Fig. 8 Schlieren images (a), x-velocity in the compression corner (b) and contour of the normalized source term P k (c) for the highenthalpy case.Black lines in (b) and (c) identify the null values

Fig. 9 Fig. 10
Fig. 9 Contour maps of the two contributions of the turbulent kinetic energy production term for the high-enthalpy case

Fig. 11
Fig. 11 Comparison of pressure (a) and shear stress (b) profiles for different expression of the turbulent kinetic energy production term

Fig. 12
Fig. 12 Comparison of N 2 mass fraction (a), O 2 mass fraction (b), T v,N 2 (c) and T v,O 2 (d) for different expressions of the turbulent kinetic energy production term