NOMENCLATURE
- RANS
Reynolds-averaged Navier–Stokes
- SWBLI
shock wave boundary layer interaction
- DNS
direct numerical simulation
- EVM
eddy viscosity model
- LSY
two-equation k–
$\varepsilon $ model of Launder and Sharma with the Yap correction- SST
two-equation k–
$\omega $ model of Menter- CLS
two-equation k–
$\varepsilon $ model of Craft et al.- ${Pr_ t}$
turbulent Prandtl number, 0.9
Greek symbols
- $\bar{ \rho }$
mean density
- $\gamma$
ratio of specific heat capacities
- $\tilde{ \sigma }$
molecular viscous stress
- $\tilde{ \tau }$
Reynolds stress
- $ \kappa $
thermal conductivity
- $\mu$
dynamic viscosity
- $\mu_t$
turbulent eddy viscosity
1.0 INTRODUCTION
The need for a reliable and cost-effective approach to travel through the atmosphere at high speed has prompted a renewed interest in supersonic and even hypersonic flight vehicles. The design and optimisation of such vehicles, including accurate predictions of aerodynamic force and aerothermal environment, requires a good understanding of complex turbulent flows, which is one of the most important unsolved problems of classical physics. The development of modern computer technology provides an opportunity for solving turbulent flow problems using Computational Fluid Dynamics (CFD). The most practical approach for the simulations of many complex flows is via Reynolds-Averaged Navier–Stokes (RANS) modelling for the consideration of cost and accuracy.
The classical RANS equations are based on incompressible flow assumptions, where the density is considered as constant. When the density variation is so significant that it cannot be ignored, however, the density should no longer be treated as constant. So, when dealing with compressible turbulent flows, not only must the velocity and pressure fluctuations be considered, but the effects of the density and temperature fluctuations on the mean motion must be accounted for as well. In such an approach, Favre (or mass) averaging is widely used to avoid the additional terms that would otherwise prevent the mean conservation equations from having close analogues to the laminar equations.
When dealing with incompressible turbulent flow, using an eddy viscosity formulation to represent the turbulent stresses, a term of 
$ \frac{2}{3}\bar{ \rho }k $ that arises from the eddy viscosity model (EVM) can be combined with the pressure in the mean momentum equation, and is thus often not included explicitly. Such a treatment is not strictly correct unless the density is uniform; however, together with the neglect of other similar terms involving the turbulent kinetic energy in the mean energy transport equations, it is often employed in low speed, or even supersonic, compressible turbulent flows. Whether this approach is suitable under hypersonic conditions is even more questionable, but nevertheless such approximations are often made, even for very high-speed turbulent flow simulations(Reference Coakley1–Reference Georgiadis, Rumsey and Huang3). Rumsey(Reference Rumsey4) validated compressibility corrections of the k–
$ \omega $ models by examining hypersonic zero-pressure gradient boundary layer flows using the CFL3D(Reference Krist5) code, in which the k terms in both the mean momentum and energy equations are not included. The resulting comparisons revealed that little difference was seen when comparing the computations of the CFL3D and VULCAN(Reference White and Morrison6) codes, the latter of which includes all k terms. However, Wilcox(Reference Wilcox7) pointed out that “at hypersonic speeds, it is entirely possible to achieve conditions under which 
$ \bar{ \rho }k $ is a significant fraction of P”, which suggests that the turbulent kinetic energy should not be ignored in the momentum equation. Since different authors treated the k terms in different ways, sometimes it is difficult to obtain a uniform evaluation from different papers for the same turbulence models on the performance of supersonic and hypersonic turbulent flow predictions.
Good-quality, detailed, experimental data of supersonic and hypersonic turbulent flows, especially those involving SWBLIs, for turbulence modelling are rather scarce. Settles and Dodson(Reference Settles and Dodson8) published a database for turbulence model validation after reviewing hundreds of experiments. Another useful summary for 2D and axisymmetric hypersonic SWBLIs is the report by Marvin et al.(Reference Marvin, Brown and Gnoffo9), in which some baseline CFD solutions were also provided. With the advances in computer technology and numerical methods in the past decades, more and more studies of Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) on high-speed turbulent flows have been published(Reference Martin10
$^-$Reference Ritos, Drikakis and Kokkinakis18), among which some high-quality data of turbulent boundary layers and shock wave boundary layer interactions are very good candidates for future, further validation and calibration of turbulence models for supersonic and hypersonic flows involving SWLBI.
In this paper, to investigate the impact of ignoring the turbulent kinetic energy terms in the RANS equations, a new density-based solver was developed within OpenFOAM (Open source Field Operation and Manipulation), based on the AUSMPW+ convection scheme(Reference Kim, Kim and Rho19), that contains different ways of including, or not including, the relevant k terms in the mean momentum and energy equations. Two linear eddy viscosity models, namely the two-equation k–
$ \varepsilon $ model of Launder and Sharma(Reference Launder and Sharma20) with the Yap correction(Reference Yap21) (LSY model) and the two-equation k–
$ \omega $ model of Menter(Reference Menter, Kuntz and Langtry22) (SST model), and one non-linear eddy viscosity model, the two-equation k-
$ \varepsilon $ model of Craft et al.(Reference Craft, Iacovides and Yoon23) (CLS model), are applied to the computations of five benchmark cases. The first case considered is a 2D zero-pressure gradient hypersonic flat plate boundary layer. Two impinging shock cases are then considered, where the externally fixed shock wave position dictates the location of the shock wave boundary layer interactions. Finally, two compression corner cases are examined, in which the resulting flow separation from the SWBLI itself affects the positioning of the shock wave and interaction. The open-source CFD package OpenFOAM V4.1 has been used to implement the new density-based compressible solver and the non-linear turbulence model.
2.0 PHYSICAL MODEL
Under the continuum mechanics assumption, the Favre-averaged Navier–Stokes equations are often used to describe the unsteady turbulent flows of compressible viscous fluids. Besides the transport equations for the mean variables, additional equations, such as those for the Reynolds stresses, the turbulent variables and other state and heat transfer correlations, are required to close the whole partial differential equation (PDE) system.
2.1 Governing equations
For an unsteady compressible flow, the Favre-averaged transport equations for mean mass, momentum and energy may be written as
\begin{equation}\hspace*{-69pt}\frac{ \partial \bar{ \rho }}{ \partial t}+\frac{ \partial }{ \partial x_{i}} \left( \bar{ \rho }\tilde{u}_{i} \right) =0\hspace*{69pt}\end{equation}
\begin{equation}\hspace*{-63pt}\frac{ \partial }{ \partial t} \left( \bar{ \rho }\tilde{u}_{i} \right) +\frac{ \partial }{ \partial x_{j}} \left( \bar{ \rho }\tilde{u}_{j}\tilde{u}_{i} \right) =-\frac{ \partial \bar{p}}{ \partial x_{i}}+\frac{ \partial }{ \partial x_{j}} \left( \tilde{ \sigma }_{ij}-\overline{ \rho u^{\prime\prime}_{j}u^{\prime\prime}_{i}} \right)\end{equation}\hspace*{70pt}
\begin{align}\frac{ \partial }{ \partial t} \left[ \bar{ \rho }\widetilde{E} \right] +\frac{ \partial }{ \partial x_{j}} \left[ \bar{ \rho }\tilde{u}_{j} \left( \tilde{\!E}+\bar{p} \right) \right] & = \frac{ \partial }{ \partial x_{j}} \left[ \tilde{u}_{i} \left( \tilde{ \sigma }_{ij}-\overline{ \rho u^{\prime\prime}_{j}u^{\prime\prime}_{i}} \right) \right] -\frac{ \partial }{ \partial x_{j}}\tilde{q}_{j} \nonumber\\
&\quad +\frac{ \partial }{ \partial x_{j}} \left[ -c_{p}\overline{ \rho u^{\prime\prime}_{j}T^{\prime\prime}}+\overline{u^{\prime\prime}_{i} \sigma _{ij}}-\overline{ \rho u^{\prime\prime}_{j}\frac{1}{2}u^{\prime\prime}_{i}u^{\prime\prime}_{i}} \right],\end{align}
where 
$\bar{ \rho } $, 
$\tilde{u}_{i} $, 
$\bar{p} $ and 
$\,\tilde{\!E} $ are the mean density, velocity, pressure and total energy with 
$\bar{} $ representing time-averaged values and 
$\tilde{} $ representing mass-averaged values. Obeying the perfect gas assumption, the equation of state and the total energy are
\begin{equation}\bar{p}= \left( \gamma -1 \right) \bar{ \rho }\tilde{e}{ ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \widetilde{E}=\tilde{e}+\frac{1}{2}\tilde{u}_{i}\tilde{u}_{i}+k}\end{equation}
where 
$\tilde{e}=c_{v}T $ is the specific internal energy and k is the turbulent kinetic energy. 
$ \gamma =\frac{c_{p}}{c_{v}} $ is the ratio of specific heat capacities, and 
$c_{p} $ and 
$c_{v} $ are the specific heat capacities at constant pressure and volume, respectively. The term 
$\tilde{ \sigma }_{ij}-\overline{ \rho u^{\prime\prime}_{j}u^{\prime\prime}_{i}} $ is the mean total stress tensor, which includes both the molecular and turbulent parts. For a Newtonian fluid, the mean molecular viscous stress is 
$\tilde{ \sigma }_{ij}= \mu\! \left( \frac{ \partial \tilde{u}_{i}}{ \partial x_{j}}+\frac{ \partial \tilde{u}_{j}}{ \partial x_{i}}-\frac{2}{3}\frac{ \partial \tilde{u}_{k}}{ \partial x_{k}} \delta _{ij} \right) $, in which the dynamic viscosity 
$ \mu $ is determined by Sutherland’s law. The heat flux vector, 
$\tilde{q}_{j} $, is obtained from Fourier’s law so that 
$\tilde{q}_{j} = -\kappa \frac{ \partial \tilde{T}}{ \partial x_{j}} $, where the thermal conductivity of fluid, 
$\kappa$, is calculated by the modified Eucken correlation (see Poling et al.(Reference Poling, Prausnitz and O’Connell24)).
The turbulent Reynolds stress 
$\tilde{ \tau}_{ij}=-\overline{ \rho u^{\prime\prime}_{j}u^{\prime\prime}_{i}} $ and turbulent heat flux 
$ \tilde{q}_{T,j}=c_{p}\overline{ \rho u^{\prime\prime}_{j}T^{\prime\prime}} $ and the so-called molecular diffusion term of 
$\overline{u^{\prime\prime}_{i} \sigma _{ij}} $ and turbulent transport term of 
$\overline{ \rho u^{\prime\prime}_{j}\frac{1}{2}u^{\prime\prime}_{i}u^{\prime\prime}_{i}} $ need to be modelled in order to close this whole PDE system.
2.2 Turbulence modelling
The turbulence models used in this paper are all EVMs, so the Reynolds stress can be modelled with the general form (up to cubic terms) of the following expression:
\begin{eqnarray}\overline{ \rho u^{\prime\prime}_{i}u^{\prime\prime}_{j}} = \frac{2}{3}\bar{ \rho }k \delta _{ij} - \mu _{t}S_{ij} + C_{1}\frac{ \mu _{t}k}{\tilde{ \varepsilon }} \left( S_{ik}S_{jk}-\frac{1}{3}S_{mk}S_{mk} \delta _{ij} \right) + \ C_{2}\frac{ \mu _{t}k}{\tilde{ \varepsilon }} \left( \Omega _{ik}S_{kj}+ \Omega _{jk}S_{ki} \right) \nonumber\\ + \ C_{3}\frac{ \mu _{t}k}{\tilde{ \varepsilon }} \left( \Omega _{ik} \Omega _{jk}-\frac{1}{3} \Omega _{lk} \Omega _{lk} \delta _{ij} \right) + \ C_{4}\frac{ \mu _{t}k^2}{\tilde{ \varepsilon }^{2}} \left( \Omega _{ki}S_{lj}+S_{kj} \Omega _{li} \right) S_{kl} \nonumber\\ + \ C_{5}\frac{ \mu _{t}k^2}{\tilde{ \varepsilon }^{2}} \left( \Omega _{il} \Omega _{lm}S_{mj}+S_{il} \Omega _{lm} \Omega _{mj}-\frac{2}{3}S_{lm} \Omega _{mn} \Omega _{nl} \delta _{ij} \right) \nonumber\\ \hspace*{-25pt}+ \ C_{6}\frac{ \mu _{t}k^2}{\tilde{ \varepsilon }^{2}}S_{ij}S_{kl}S_{kl} + \ C_{7}\frac{ \mu _{t}k^2}{\tilde{ \varepsilon }^{2}}S_{ij} \Omega _{kl} \Omega _{kl},\qquad\end{eqnarray}
where 
$S_{ij}= \left( \frac{ \partial \tilde{u}_{i}}{ \partial x_{j}}+\frac{ \partial \tilde{u}_{j}}{ \partial x_{i}} \right) -\frac{2}{3}\frac{ \partial \tilde{u}_{k}}{ \partial x_{k}} \delta _{ij} $ and 
$ \Omega _{ij}= \left( \frac{ \partial \tilde{u}_{i}}{ \partial x_{j}}-\frac{ \partial \tilde{u}_{j}}{ \partial x_{i}} \right) $ are the strain rate and vorticity, respectively and 
$\mu _{t} $ is the turbulent eddy viscosity. If 
$C_{1} \sim C_{7} $ all equal zero, this form is the so-called Boussinesq approximation and the model under this approximation is often also called a linear EVM. Otherwise, models making use of some, or all, of the non-linear terms in the above stress–strain relation can be referred to as non-linear EVMs. In this paper, all the turbulence models used are in their original forms, so the details of each model can be found in Refs (Reference Launder and Sharma20)–(Reference Craft, Iacovides and Yoon23).
The turbulent heat fluxes 
$\tilde{q}_{T,j}=c_{p}\bar{ \rho u^{\prime\prime}_{j}T^{\prime\prime}} $ are modelled using the simple eddy-diffusivity approximation, which can be written as 
$\tilde{q}_{T,j}=c_{p}\bar{ \rho u^{\prime\prime}_{j}T^{\prime\prime}}=-\frac{ \mu _{t}c_{p}}{\Pr_{t}}\frac{ \partial \tilde{T}}{ \partial x_{j}} $, in which the turbulent Prandtl number, 
$Pr_{t} $, is set to be 0.9. An approximation of a gradient expression for the molecular diffusion and turbulent transport appearing in the mean energy equation is used here as follows:
\begin{equation}\overline{ u^{\prime\prime}_{i} \sigma _{ij}}-\overline{ \rho u^{\prime\prime}_{j}\frac{1}{2}u^{\prime\prime}_{i}u^{\prime\prime}_{i}}= \left( \mu +\frac{ \mu _{t}}{ \sigma _{k}} \right) \frac{ \partial k}{ \partial x_{j}},\end{equation}
where 
$ \sigma _{k} $ is a model coefficient.
2.3 Presence of the turbulent kinetic energy
Following the eddy viscosity assumption, the transport Equations (2) and (3), under the Boussinesq approximation for the turbulent stresses, can be rewritten as
\begin{equation}\frac{ \partial }{ \partial t} \left( \bar{ \rho }\tilde{u}_{i} \right) + \frac{ \partial }{ \partial x_{j}} \left( \bar{ \rho }\tilde{u}_{j}\tilde{u}_{i} \right) = -\frac{ \partial \bar{p}}{ \partial x_{i}}+\frac{ \partial }{ \partial x_{j}} \left[ \underbrace{ \left( \mu + \mu_{t} \right) \left( \frac{ \partial \tilde{u}_{i}}{ \partial x_{j}}+\frac{ \partial \tilde{u}_{j}}{ \partial x_{i}}-\frac{2}{3}\frac{ \partial \tilde{u}_{k}}{ \partial x_{k}} \delta_{ij} \right) }_{ \left( d \right) }-\underbrace{\mathop{\frac{2}{3}\bar{ \rho }k\delta_{ij}}}_{ \left( a \right) } \right],\end{equation}
\begin{eqnarray} \frac{ \partial }{ \partial t} \left[ \bar{ \rho } \left( \tilde{e}+\frac{1}{2}\tilde{u}_{i}\tilde{u}_{i} + \underbrace{k}_{ \left( b \right) } \right) \right] +\frac{ \partial }{ \partial x_{j}} \left[ \bar{ \rho }\tilde{u}_{j} \left( \left( \tilde{e}+\frac{1}{2}\tilde{u}_{i}\tilde{u}_{i}+\underbrace{k}_{ \left( b \right) } \right) +\bar{p} \right) \right] \nonumber \\= \frac{ \partial }{ \partial x_{j}} \left[ \tilde{u}_{i} \left( \underbrace{ \left( \mu + \mu_{t} \right) \left( \frac{ \partial \tilde{u}_{i}}{ \partial x_{j}}+\frac{ \partial \tilde{u}_{j}}{ \partial x_{i}}-\frac{2}{3}\frac{ \partial \tilde{u}_{k}}{ \partial \mathop{x}_{k}} \delta_{ij} \right) }_{ \left( d \right) }-\underbrace{\frac{2}{3}\bar{ \rho }k \delta_{ij}}_{ \left( a \right) } \right) - \tilde{q}_{j} \right. \nonumber \\ \left. -\frac{ \mu _{t}c_{p}}{\Pr_{t}}\frac{ \partial \widetilde{T}}{ \partial x_{j}} +\underbrace{ \left( \mu +\frac{\mu_{t}}{\sigma_{k}} \right) \frac{ \partial k}{ \partial x_{j}}}_{ \left( c \right) } \right].\end{eqnarray}
For the non-linear EVM, the two equations above retain the same form, except for the term (d), which should then contain contributions from the non-linear terms as well. In Equations (7) and (8), there are three different terms related to the turbulent kinetic energy, k, which are often neglected. The term (a) in Equation (7) is often ignored or absorbed into the pressure term when dealing with incompressible turbulent flows. The terms (b) and (c) in Equation (8) are also often ignored for high-speed compressible turbulent flows. Table 1 lists three possible combinations for the presence or absence of these terms which have been found in software or reported in research papers. For combination C1, all three terms are ignored, which is the commonly made assumption in incompressible flow simulations. In OpenFOAM V4.1, the integrated density-based compressible solver, rhoCentralFoam, also takes this approach. The second combination C2, includes term (a) in the mean momentum and energy transport Equations (7) and (8), but neglects terms (b) and (c), which is often a reasonable approximation since 
$ k \ll \,\tilde{\!E} $ in many flows of engineering interest. The last combination C3 includes all the terms related to the turbulent kinetic energy in Equations (7) and (8).
Three combinations of the inclusion of the turbulent kinetic energy terms in the RANS equations
As the turbulent kinetic energy, k, has already been calculated when solving the transport equations of turbulent variables, the inclusion, or exclusion, of the terms (a) and (b), which simply involve the value of k, would be expected to have little impact on computational cost, or solution efficiency. The term (c), in practice, it is calculated as a diffusion term of the turbulent kinetic energy, with diffusivity of 
$\mu +\frac{\mu_{t}}{\sigma_{k} }$, whose cost of calculation, compared with that of the whole mean energy equation, is still relatively low. Thus, the use of these three combinations of the k-related terms would be expect to introduce only very minor differences in the cost and efficiency of the overall process of numerical simulation, and this is what was found in the calculations reported below.
3.0 NUMERICAL IMPLEMENTATION
The governing equations given above are solved with the Finite Volume Method (FVM) on a three-dimensional unstructured grid of polygonal cells via OpenFOAM, which is open-source software for CFD. For the present study, the software package OpenFOAM V4.1 is employed, which includes a variety of density-based and pressure-based solvers for compressible flows. As this research is mainly concerned with supersonic and hypersonic flows, starting from the density-based compressible solver, rhoCentralFoam, a new solver with the three different optional turbulent kinetic energy terms in Table 1 has been implemented. In the numerical simulation of high-speed flows, the representation of the convection terms is crucial for accuracy. In this study, the AUSMPW+ scheme is used to calculate the convection term for all the variables and the second-order central scheme is used for diffusion terms. To keep the solution bounded, and achieve high resolution, the second-order Total Variation Diminishing (TVD) interpolation with the van Albada limiter is used for the interface calculations. Local time marching is used to accelerate convergence to steady state, and the local Courant number is set to be less than 5.0 for the purpose of retaining numerical stability.
Like the new density solver, the CLS non-linear turbulence model has also been implemented into OpenFOAM V4.1, in a similar way to how the Launder–Sharma k–
$ \varepsilon $ model is coded. The major difference between the two is the use of the non-linear constitutive relation of Equation (5) in the former, using the model coefficients from Ref. (Reference Craft, Iacovides and Yoon23).
Since the total temperature of all the simulations in this paper is below 1800K, the working medium, dry air, is assumed to be a thermally and calorically perfect gas. For wall boundary conditions on the turbulence variables, 
$k_{wall}$ is set to zero. In the LSY model and CLS models 
$\tilde{ \varepsilon }_{wall}=0$, whilst in the SST model 
$\omega _{1}=\frac{6 \mu _{wall}}{ \rho _{wall} \beta _{1} \left( \Delta d_{1} \right) ^{2}_{\phantom{f}}}$, where 
$\omega _{1}$ is the 
$\omega$ value at the centre of the first cell above the wall, 
$\beta _{1}=0.075$, and 
$\Delta d_{1}$ is the distance from the centre of the first cell to the nearest wall. To generate inlet conditions for the SWBLI cases, separate simulations were performed of developing flat plate boundary layers, and profiles of mean and turbulence quantities were extracted at a suitable location to match the reported experimental conditions. The freestream turbulence levels in these boundary layers were set as:
LSY and CLS models:
\begin{equation}\frac{k_{\infty}^{\frac{1}{2}}}{U_{\infty}}=2\%, \ \ \ \ \frac{ \mu _{t}}{ \mu _{\infty}}=\frac{C_{ \mu } \cdot \left( k^{2}/\tilde{ \varepsilon } \right) }{ \mu _{\infty}/ \rho _{\infty}}=10\end{equation}
SST model (following the recommended values of Menter(Reference Menter25)):
\begin{equation}k_{\infty}=\frac{0.1U_{\infty}^{2}}{Re_{L}},\ \ \ \ \omega _{\infty}=5\frac{U_{\infty}}{L},\end{equation}
where L is the length of the computational domain.
4.0 RESULTS AND DISCUSSION
4.1 Zero-pressure-gradient hypersonic flat plate
The 2D zero-pressure-gradient flat plate boundary layer is a fundamental benchmark case which is widely used for the verification and validation of RANS models as well as numerical solvers on compressible turbulent boundary layer predictions. The data used for comparison here are obtained directly from the online database named ‘Turbulence Modelling Resource hosted by Langley Research Center. The boundary conditions are given in Table 2. A mesh of 128
$\times$128 cells is used to run all the simulations in this section. The size of the first cell off the wall is set to be 
$2.0 \times 10^{-6}$ m, which maintains the 
$y^+$ at the wall for all the simulations to be less than 1.0.
Boundary conditions of the 2D zero-pressure-gradient boundary layer flow
Comparison of wall skin friction of the LSY, SST and CLS models for the Ma = 5 flat plate case.
Comparison of velocity profiles at 
$\mathrm{Re}_{ \theta }$ =10,000 of the LSY, SST and CLS models for the Ma = 5 flat plate case.
In Fig. 1, profiles of the wall skin friction predicted by the LSY, SST and CLS models, versus 
$\mathrm{Re}_{ \theta }$, the Reynolds number based on momentum thickness, between 4,000 and 12,000, are compared. The velocity profiles which result from the use of the same models, at 
$\mathrm{Re}_{ \theta }$ of 10,000, are compared in Fig. 2. Firstly, the predicted 
$C_f$, wall friction coefficient, and velocity profiles from all the models match both the theoretical results and each other within a reasonable error tolerance for all three turbulence models, which confirms that the new density-based compressible solver and the non-linear CLS model have been implemented properly within the OpenFOAM V4.1 framework.
Secondly, the results above indicate that including or neglecting the turbulent kinetic energy terms (a), (b) or (c) has only a limited effect on the solution of the RANS equations in this case, implying that the assumptions 
$\bar{ \rho }k \ll \bar{p}$ and 
$k \ll \,\tilde{\!E}$ do hold in the case of zero-pressure-gradient hypersonic boundary layer flow. To confirm this, the ratios of 
$ \frac{2}{3}\bar{ \rho }k/{\bar{p}} $ and 
${k}/\,{\tilde{\!E}}$ predicted by the CLS model, at 
$\mathrm{Re}_{ \theta }$ of 10,000, are presented in Fig. 3. The ratio of 
$ \frac{2}{3}\bar{ \rho }k/{\bar{p}} $ is less than 6%, and the ratio of 
${k}/\,{\tilde{\!E}}$ is even lower, with a maximum value of less than 2%. So, for the simple zero-pressure-gradient boundary layer flows, the introduction of the k terms results in little difference and does not affect the validation of the performance of turbulence models.
Comparison of the ratios 
$ \frac{\text{2}}{\text{3}}\bar{ \rho }\textit{k}/{\bar{\textit{p}}} $ and 
${\textit{k}}/\,{\tilde{\!\textit{E}}}$ predicted by the CLS model, at 
$ \text{Re}_{ \theta } $ = 10,000 for the Ma = 5 flat plate case.
Schematic of 2D impinging shock flow.
4.2 Impinging shock
A 2D/axisymmetric impinging shock (Fig. 4) occurs when an externally generated oblique shock impinges on a flat/cylinder surface boundary layer. The interaction between an impinging oblique shock and turbulent boundary layer embodies the problems of compressibility and, depending on the strength of the shock, can lead to flow separation and turbulence amplification. The supersonic DNS case by Pirozzoli and Grasso(Reference Pirozzoli and Grasso15) at Mach number of 2.25 and the hypersonic experiment of Kussoy and Horstman(Reference Kussoy and Horstman26) at Mach number of 8.18 are tested here. For the hypersonic case, the data are obtained from the database of Marvin et al.(Reference Marvin, Brown and Gnoffo9)
4.2.1 Ma $ \boldsymbol{=} \textbf{2.25} $ impinging shock wave interaction
The freestream condition, as well as the boundary-layer information at the nominal impingement point, of the DNS are listed in Table 3. A two-dimensional structured grid, shown in Fig. 5, uniform in the streamwise direction and expanding exponentially in the direction normal to the wall, is used for the RANS simulations. Similar to the DNS calculation, the lower part of the inlet boundary is obtained from a separate flat plate boundary layer simulation which guarantees the momentum thickness 
$\theta_{0}$ of the boundary layer matches that of the DNS data at the impingement point. At the upper part, the shock is artificially applied by applying fixed values of the flow field satisfying the oblique shock wave relation with a shock angle of 
$33.2^{\circ}$ (the corresponding turn angle is about 
$8.1^{\circ}$). Zero-gradient boundary conditions for pressure and temperature and no-slip boundary condition for velocity are specified at the wall.
DNS flow conditions for the Ma = 2.25 impinging shock case by Pirozzoli and Grasso(Reference Pirozzoli and Grasso15)
Grid for the Mach = 2.25 impinging shock case.
Comparison of the predicted wall pressure (left) and skin friction coefficients (right) by three sets of grids for the Mach = 2.25 impinging shock case.
Three sets of grids, consisting of 
$80 \times 65\ \mathrm{cells}$, 
$160 \times 130\ \mathrm{cells}$ and 
$240 \times 195\ \mathrm{cells}$, maintaining the value of 
$y^+$ at all the near-wall nodes less than 0.1, were used to carry out a grid independence test. The profiles of the wall pressure and skin coefficient, resulting from the use of the non-linear CLS model with the k-related terms included in all transport equations (combination C3), are shown in Fig. 6. As is seen from the figures, little difference in the distribution of these two wall quantities exists between the results of the medium and the fine grids, but the separation length predicted by the coarse grid is shorter than that of the other two, which may be caused by the stronger numerical dissipation induced by larger grid spacing. The minor discrepancy between the results by the medium and fine grids shows that the grid density of 
$160 \times 130\ \mathrm{cells}$ is sufficient to achieve a high level of grid independence, thus the medium grid was selected to run all the simulations. Similar tests were also carried out for the other cases reported here, to ensure that the results presented are grid independent. A summary of the grid densities used in them is given in the Appendix at the end of this paper.
In Figs 7 and 8, the wall pressure and friction coefficient predicted by the three turbulence models, with different combinations of turbulent kinetic energy terms included in the transport equations, are compared with those of the DNS data (the nominal impingement point is at x = 22.04cm). In the interaction region, the LSY prediction of the wall friction coefficient is in the closest agreement with the DNS data, while the CLS model returns the best pressure distribution near the impingement point. Use of the SST model results in the prediction of a much larger separation and lower wall friction downstream of the impingement point. As also observed in the zero-pressure-gradient flat plate case, the presence of different k-term combinations in the transport equations results in only minor variations in the predicted wall quantities for all three models. Figure 9 shows the iso-lines of mean pressure superimposed on contours of the ratios of 
$\frac{2}{3}{\bar{ \rho }k}/{\bar{p}}$ and 
${k}/\,{\tilde{\!E}}$ predicted by the CLS model with combination C3 around the interaction region. The peak values of these ratios are less than 
$5.5\%$ and 3%, respectively, which explains why only minor difference can be found in the predicted wall quantities. This phenomenon again illustrates that the k terms of (a)–(c) in Equations (7) and (8) have little effect at this Mach number of 2.25 for the SWBLI case if the main focus is on the wall quantities.
Comparison of wall pressure, by the LSY, SST and CLS models, with the DNS data for the Ma = 2.25 impinging shock case.
Comparison of wall friction coefficients, by the LSY, SST and CLS models, with the DNS data for the Ma = 2.25 impinging shock case.
There is, however, one noticeable difference caused by the inclusion of the k terms in the transport equations, in the mean pressure contour patterns in the near-wall region behind the impingement point. Figure 10 displays the iso-lines of mean pressure predicted by the CLS model with different k-term combinations in the transport equations (the other two models return predictions similar to those of the CLS model), as well as the corresponding DNS iso-lines. The DNS iso-lines show an S-shaped feature downstream of the separation bubble, highlighted by the red box, while RANS predictions resulting from combination C1 (no k terms in the transport equations) fail to capture this pattern, showing iso-lines that approach the wall almost vertically. By contrast, the RANS predictions from the corresponding combinations C2 and C3, both of which include the term labelled (a) in the momentum and energy transport equations, capture the S-shaped feature in the pressure iso-lines. This finding implies that, although the inclusion of the turbulent kinetic energy terms in the transport equations for the supersonic impinging shock case has little effect on the prediction of the wall parameters, even the relatively small contribution of up to around 5% of the 
$ \frac{2}{3}\bar{ \rho }k/{\bar{p}} $ term to the effective pressure seen by the solver is sufficient to alter the shape of the pressure distribution in this region.
4.2.2 Ma $ \boldsymbol{=} \textbf{8.18} $ impinging shock wave interaction
The experimental configuration of Kussoy and Horstman(Reference Kussoy and Horstman26) is shown in Fig. 11, and the wedge angle (turn angle) of 
$10^{\circ}$ is considered here. Initially, without the shock generator present, the authors reported an undisturbed boundary layer profile, described in Table 4, at a distance of 187cm from the leading edge of the flat plate (47cm downstream of where the wedge leading edge was subsequently placed). The inlet to the computational domain was, therefore, taken as 5cm ahead of the leading edge of the wedge, and a separate flat plate boundary layer calculation was carried out to provide appropriate boundary layer profiles to impose at this location. The domain extended to 19cm behind the trailing edge of the wedge. A structured grid of 
$560 \times 180$ cells (see Fig. 12) is employed, clustered at the leading edge of the wedge and uniform in the interaction region, with the 
$y^+$ at the near-wall nodes kept at a value of less than 0.1 on the test bed and less than 0.5 on the wedge. Zero-gradient condition is applied for all variables at the outlet boundary, and at the wall the no-slip condition is applied with temperature fixed to the experimental value and zero gradient for pressure.
Undisturbed boundary layer parameters of the Ma = 8.18 impinging shock case
Iso-lines of mean pressure superimposed on the contours of the ratios of 
$ \frac{\text{2}}{\text{3}}\bar{ \rho }\textit{k}/{\bar{\textit{p}}} $ (left) and 
${\textit{k}}/{\tilde{\!\textit{E}}}$ (right) predicted by the CLS model with combination C3 for the Ma = 2.25 impinging shock case.
Comparison of mean pressure, predicted by the CLS model with combinations C1, C2 and C3, with the DNS data of Pirozzoli and Grasso(Reference Wu and Martin16), for the Ma = 2.25 impinging shock case.
Experimental configuration of Kussoy and Horstman(Reference Kussoy and Horstman26).
Figures 13 and 14 show comparisons of measured and predicted wall pressure and heat flux (the projection of the wedge leading edge onto the flat plate is located at x = 0cm, and the nominal impingement point is around x = 36.89cm). The prediction of the wall static pressure is not very sensitive to changes in the turbulence model; all predictions are in close agreement with the data, with the only exception being the prediction of the start of the flow separation indicated by rise in wall pressure, which erroneously occurs later with the linear k–
$ \varepsilon $ model. The local wall heat flux comparisons, on the other hand, show that all the EVMs tested over-estimate the peak levels in the interaction region and there are now greater differences between the model predictions, with the linear k–
$ \varepsilon $ producing the greatest over-estimate and the non-linear version the lowest (but still significant) over-estimate. The predicted profiles of the local wall heat flux also show a noticeable dependence on the presence, or absence, of the k terms in the transport equations, with the predictions that result when these terms are neglected (option C1) being furthest from the experimental data.
Grid for the Ma = 8.18 impinging shock case.
Comparison of wall pressure, by the LSY, SST and CLS models, with the experimental data for the Ma = 8.18 impinging shock case.
Comparison of wall heat flux, by the LSY, SST and CLS models, with the experimental data for the Ma = 8.18 impinging shock case.
Further information on the predictive characteristics of the different modelling strategies is presented in Fig. 15, in the form of profiles of the wall friction coefficient predicted by the three models, with various treatments of turbulent kinetic energy terms in the transport equations. As seen from these figures, unlike the supersonic impinging shock case, more significant differences now appear in the prediction of the locations of flow separation. For all three turbulence models, combination C1 with no k terms in the transport equations predicts that flow separation starts the furthest downstream, while combination C2, where only the k term identified as (a) is present in the transport Equations (7) and (8), by contrast, predicts the earliest onset of flow separation. Combination C3, with all k terms included, results in predictions of the locations of the start of flow separation somewhere in the middle between those of combinations C1 and C2. The predicted location of the reattachment point, however, appears to be almost unaffected by the inclusion of the k terms in the transport equations.
The separation lengths predicted by the three different models with three treatments of the k terms in the transport equations are listed in Table 5. In the last column of Table 5, for each model the separation lengths are shown normalised with the value predicted when all k terms are included, C3. Combination C1, ignoring all three k terms, tends to underpredict the separation by at least 23%. Combination C2, however, which retains the k term (a) in the transport equations, overpredicts the separation by at least 29%. This suggests that the k terms in the transport equations exert a strong influence on the conservation of the mean momentum and energy equations, and the assumptions of 
$ \bar{ \rho }k \ll \bar{p} $ and 
$ k \ll \,\tilde{\!E} $ no longer hold in the interaction region of this hypersonic SWBLI case.
Comparison of separation lengths of the hypersonic impinging shock case
Comparison of wall friction coefficient, by the LSY, SST and CLS models, for the Mach = 8.18 impinging shock case.
Since the three turbulence models return similar variation patterns, the ratios of 
$ \frac{2}{3}\bar{ \rho }k/{\bar{p}} $ and 
$ {k}/\,{\tilde{\!E}} $ in the interaction region produced only by the CLS model are presented in Fig. 16. Near the impingement point, the ratio of 
$ \frac{2}{3}\bar{ \rho }k/{\bar{p}} $ can reach up to 32%, at its peak, which would cause a large imbalance in the momentum equation if ignored. The k, turbulent kinetic energy, in the same region can make a contribution of up to 13% of the total mean energy, and neglecting it can also lead to significant loss of energy conservation. This explains why the hypersonic case results in a more significant variation in the prediction of the flow separation, as well as the skin friction, due to the treatment of the k-related terms, whereas it did not show a strong influence in the earlier, lower Mach number, case.
4.3 Compression corner
The compression corner (see Fig. 17) is a benchmark case which is widely used and accepted to test the ability of turbulence models to predict a high-speed flow field with shock wave-induced separation. The DNS case by Wu and Martin(Reference Wu and Martin16) at Mach number of 2.9 and the experiments of Kussoy and Horstman(Reference Kussoy27) at Mach number of 7.05 serve as the test cases reported here. The statistical data of the DNS case have been extracted from the published paper of Priebe et al.(Reference Priebe, Wu and Martin28), while the experimental data for the hypersonic case have again been obtained from the database of Marvin et al.(Reference Marvin, Brown and Gnoffo9)
4.3.1 Ma $ \boldsymbol{=} \textbf{2.9} $ compression corner
The incoming flow conditions of this DNS case, where a ramp angle of 
$24^{\circ}$ was considered, are listed in Table 6. A structured grid of 
$ 200 \times 180 $ cells, uniform in the streamwise direction and expanding exponentially in the wall-normal direction, with 
$y^+$ along the near-wall nodes less than 0.1, shown in Fig. 18, is used for the RANS simulations. The inlet boundary conditions are obtained through a preliminary computation of the flow over an undisturbed flat plate which produces a turbulent velocity profile of the same displacement thickness, 
$ \delta _{0}^{\ast} $, as that provided by the DNS study at a specified location upstream of the compression corner. The wall boundary is set to be no-slip for velocity, with zero gradient for pressure and fixed temperature. The ‘extra’ flat part of wall, after the ramp surface, is added to guarantee a supersonic outlet boundary. Zero-gradient conditions are applied for all variables at the outlet boundaries.
DNS flow conditions for the Ma = 2.9 compression corner case by Wu and Martin(Reference Wu and Martin16)
Iso-lines of mean pressure superimposed on contours of the ratios of 
$ \frac{\text{2}}{\text{3}}\bar{ \rho }\textit{k}/{\bar{\textit{p}}} $ (left) and 
$ {\textit{k}}/\,{\tilde{\!\textit{E}}} $ (right) predicted by the CLS model with combination C3 for the Ma = 8.18 impinging shock case.
Schematic of 2D hypersonic compression corner flow.
Figures 19 and 20 display the wall pressure and skin friction coefficients (the corner is located at 
$x=0$cm). In the interaction region, the LSY model returns the best overall agreement with the DNS data in terms of wall pressure and skin friction coefficients, while the SST model again returns a much longer separation bubble than that of the DNS. Similar comparisons have been reported in the paper of Georgiadis and Yoder(Reference Georgiadis and Yoder29), in which they discussed the possible reason and recalibrated the SST model for flows with SWBLIs. The CLS model also over-estimates the separation length, but to a far lesser extent than the SST model, and in the downstream region the profiles of the wall pressure and skin coefficient returned by the CLS model are in close agreement with those of the DNS data.
The separation lengths predicted by the three models are listed in Table 7. Compared with the Mach = 2.25 impinging shock case, the presence of the turbulent kinetic energy terms in the mean flow transport equations exerts a stronger influence on the prediction of flow separation. The separation lengths predicted by each model with combination C1 are about 10% shorter than those predicted with combination C3, while the results with combination C2 are approximately 10% longer. As also seen in the comparison of Figs 18 and 19, the linear LSY model predicts separation length closest to that of the DNS.
Comparison of separation lengths of the Ma = 2.9 compression corner case
Grid for the Ma = 2.9 compression corner case.
Comparison of wall pressure, by the LSY, SST and CLS models, with the experimental data for the Ma = 2.9 compression corner case.
Comparison of skin friction coefficients, by the LSY, SST and CLS models, with the experimental data for the Ma = 2.9 compression corner case.
The stronger influence of the inclusion of the k terms in the mean flow transport equations on the prediction of flow separation in this case implies that the ratios 
$ \frac{2}{3}\bar{ \rho }k/{\bar{p}} $ and 
$ {k}/\,{\tilde{\!E}} $ must be higher than those of the supersonic impinging shock case. Contours of these two k-related parameter ratios, resulting from the use of the CLS model, are plotted in Fig. 21. The peaks of 
$ \frac{2}{3}\bar{ \rho }k/{\bar{p}} $ and 
$ {k}/\,{\tilde{\!E}} $, in the interaction region, are around 13% and 6%, respectively, being more than twice as high as those of the corresponding ratios in the supersonic impinging shock case. Since these two cases are both in the supersonic regime, the much higher values suggest that it may not simply be the somewhat higher Mach number in the corner case that leads to the stronger influence of the k terms, but that it could also be due to differences in the strength of the shock and interaction. To illustrate this, Fig. 22 shows ratios of the local to freestream static pressure in the two cases, where the significantly denser iso-lines, and the greater pressure increase across the shock, in the compression corner case can be clearly seen. Some of this increase in shock strength is due to the higher Mach number (2.9 in the corner case, compared with 2.25 in the impinging one), but some is also due to the greater angle the flow is turned through in the corner case. To quantify this, if the Mach number of the impinging case were increased to 2.9, the shock angle in that case would then be 
$26.4^{\circ}$, whereas the one seen in the compression corner case (Fig. 21) is actually around 
$28.9^{\circ}$. The above results show that even at moderate Mach numbers the omission of the k-related terms from the mean transport equation can have an effect when the shock interactions become more complex.
Iso-lines of mean pressure superimposed on the contours of the ratios of 
$ {\frac{\text{2}}{\text{3}}\bar{ \rho }\textit{k}}/{\bar{\textit{p}}} $ (left) and 
$ {\textit{k}}/\,{\tilde{\!\textit{E}}} $ (right) predicted by the CLS model with combination C3 for the Ma = 2.9 compression corner case.
Comparison of ratios of static pressure to freestream static pressure by the CLS model with combination C3 for the Ma = 2.25 (left) and Ma = 2.9 compression corner case (right).
4.3.2 Ma $ \boldsymbol{=} \textbf{7.05} $ compression corner
The test model of the experiment is shown in Fig. 23, for which a flare angle 
$\theta$ of 35o is considered. The nominal freestream conditions are given in Table 8. The inlet boundary condition was obtained by first calculating an undisturbed cylinder boundary layer flow and matching the measured and computed displacement thickness at 6cm ahead of the corner. A no-slip condition is applied at the wall, with temperature fixed to the experimental value and zero gradient for pressure. Zero-gradient condition is applied for all variables at the outlet boundaries. Since the test model is axisymmetric, the ‘wedge’ boundary condition of OpenFOAM is used to run axisymmetric simulations using one cell in the circumferential direction, with a wedge angle of 
$1^{\circ}$. The computations are performed using a 
$200 \times 200$ cell structured grid (see Fig. 24), uniform in the streamwise direction and expanding exponentially in the wall normal direction. The 
$y^+$ at the first grid point off the wall is kept at a value of less than 0.1. The lengths of cylinder and flare are both set to be 20cm, which are long enough to cover the whole interaction region.
Incoming boundary conditions for the Ma = 7.05 compression corner case
The test model of the experiment of Kussoy and Horstman(Reference Kussoy27).
Grid for the Ma = 7.05 compression corner case.
Comparison of wall pressure, by the LSY, SST and CLS models, with the experimental data for the Ma = 7.05 compression corner case.
Comparison of wall heat flux, by the LSY, SST and CLS models, with the experimental data for the Ma = 7.05 compression corner case.
Comparison of wall friction coefficient, by the LSY, SST and CLS models, for the Ma = 7.05 compression corner case.
The predicted surface pressure and wall heat flux with different combinations of the k terms in the transport equations, from all three models, are compared in Figs 25 and 26 (x is in the streamwise direction along the test model, and the corner is located at 
$x = 0$cm). The non-linear CLS model returns the best prediction, compared with the measured wall pressure profiles, with respect to both the peak value and the pressure plateau caused by the flow separation. The LSY model under-estimates the separation, while the SST model, by contrast, over-estimates it. The predicted peak values of the wall pressure show less sensitivity to the inclusion of the k terms in the transport equations than the corresponding peaks in the wall heat flux, although the locations of the peaks in both are significantly influenced by the inclusion of the terms. All three models overpredict the peak of wall heat flux severely around the reattachment region. Again, the predictions of the linear LSY model are in closer overall agreement with the experimental data, though as far as the pressure variation is concerned the non-linear CLS model is in closer agreement. Predicted profiles obtained without any k terms in the transport equations appear to be furthest from the experimental ones for both the CLS and LSY models.
The wall friction coefficients predicted by the three models tested, with different combinations of the k terms in the transport equations, are shown in Fig. 27, and the corresponding separation lengths are listed in Table 9. The separation location measured was at 
$x = -6.3$cm, while no data for the attachment location were reported. For each model, the predicted lengths of the flow separation are strongly influenced by the inclusion of the k terms in the mean flow transport equations. Combination C1, with no k terms in the transport equations, returns the smallest separation bubble, and compared with the results for combination C3, with all k terms included, the bubble is shorter by over 35%. In contrast, combination C2 results in a separation bubble almost 40% longer than that for combination C3. Like in the hypersonic impinging shock case, the separation point locations are strongly influenced by the variations in the treatment of the k terms, but now the influence also extends to the prediction of the location of the reattachment point. One possible reason for this phenomenon is that, for the impinging shock flow, the nominal impingement point is determined by the turn angle of the shock generator, so the high-pressure region downstream of the reattachment point is then also largely determined by this angle itself. As a result, the reattachment point of the flow is largely determined by the angle of the incident shock (see Fig. 4), and the treatment of the k terms, which directly or indirectly change the pressure gradient through the mean transport equations, would therefore mainly affect the separation point. However, for the compression corner flow, at hypersonic speeds, the impinging shock is indeed the ‘separation shock’ (see Fig. 17), so the impingement point would change as the separation point moves upstream or downstream. Thus, the different combinations of the k-related terms could cause a great impact on both the separation and reattachment points. Finally, accepting that the C3 versions of the mean flow equations result in a most representative solution of each model, the CLS model produces the most reliable prediction of the location of the separation point.
The significant differences in the separation length, again, indicate that 
$ \frac{2}{3}\bar{ \rho }k $ and k must be a considerable proportion of 
$ \bar{p} $ and 
$ \,\tilde{\!E} $, respectively. The comparison of the ratios of 
$ \frac{2}{3}\bar{ \rho }k/{\bar{p}} $ and 
$ {k}/\,{\tilde{\!E}} $ by the CLS model, with combination C3, is shown in Fig. 28. As seen from the figures, the peak values of the ratios 
$ \frac{2}{3}\bar{ \rho }k/{\bar{p}} $ and 
$ {k}/\,{\tilde{\!E}} $ reach up to 45% and 18%, respectively. These comparisons suggest that the assumptions 
$ \bar{ \rho }k \ll \bar{p} $ and 
$ k \ll \,\tilde{\!E} $ no longer hold in the interaction region of the hypersonic flows involving SWBLIs, especially with flow separation. The k terms of (a), (b) and (c), therefore, should not be neglected when computing such flows.
Comparison of separation lengths of the Ma = 7.05 compression corner case
Iso-lines of mean pressure superimposed on the contours of the ratios of 
$ \frac{\text{2}}{\text{3}}\bar{ \rho }\textit{k}/{\bar{\textit{p}}} $ (left) and 
$ {\textit{k}}/\,{\tilde{\!\textit{E}}} $ (right) predicted by the CLS model with combination C3 for the Ma = 7.05 compression corner case.
5.0 CONCLUDING REMARKS
The role of the terms which involve the turbulent kinetic energy, which are often neglected from the compressible RANS equations, has been investigated by examining typical supersonic and hypersonic benchmark cases. Predictions resulting from three different levels of inclusion of these terms in the mean momentum and energy equations, ranging from all of them being left out to all included, have been compared.
In simple zero-pressure-gradient hypersonic boundary layer flow, the numerical results showed that neglecting the three k terms, as is standard practice in incompressible flow computations, has only limited impact on the predicted flow field and surface variables. For supersonic flows with SWBLI, the impact of the inclusion of the k terms in the mean momentum and energy transport equations is greater than that for the zero-pressure-gradient boundary layer case, but the effects on wall quantities, such as wall pressure and friction coefficients, may not be very significant. However, with greater shock wave strength, the discrepancy caused by ignoring all the k terms or part of them, taking the separation length as an example in the compression corner case, could be around 10
$\%$. It should, also, be noted that, although the influence on the wall variables may not be substantial for some engineering applications, if the spatial distribution of flow field variables, such as mean pressure, is the focal point, neglecting the turbulent kinetic energy terms from the mean momentum equations could lead to the incorrect prediction of the shape of iso-line patterns.
For hypersonic flows with SWBLI, the inclusion of turbulent kinetic energy terms plays a very important role in the conservation of mean momentum and energy equations. Compared with including all the k terms in the RANS equations, ignoring all of them could delay the predicted flow separation, leading to an influence of up to 40% on the prediction of the separation length. In contrast, neglecting only the terms in the mean total energy definition and the molecular diffusion and turbulent transport could hasten the separation significantly, with an overprediction ratio for the separation length of up to 1.4, compared with the results when including all the k terms. The different treatments of k-related terms mainly affect the prediction of separation point for the impinging shock flow but would cause a strong influence on the prediction of both the separation and reattachment points for the compression corner flow. While the full inclusion of these terms does not always result in predictions that agree better with DNS/experimental data, this is likely caused by the fact that their exclusion cancels out effects of other flaws in the RANS models employed. Since the k terms have such a significant influence on the numerical solution of the RANS equations, to get a fair, uniform and reliable validation of different turbulence models for hypersonic flows, it is necessary to fully include these terms in the mean flow transport equations. Moreover, although only axisymmetric/2D cases were tested here, it can be expected that, with the increase of complexity, the impact of including, or excluding, the k-related terms might have an even greater impact on flows with three-dimensional SWBLIs.
From the above findings, it is recommended that all three of these terms be included when running hypersonic, or even supersonic, turbulent flow simulations, especially for flows with shock wave-induced separations. As a consequence, the use of simple turbulence models that do not provide turbulent kinetic energy explicitly, and are therefore unable to represent these terms when simulating supersonic and hypersonic turbulent flows with SWBLIs, is questionable.
ACKNOWLEDGEMENTS
H. Zhang is grateful for the support of China Aerodynamic Research and Development Center throughout this work.
SUPPLEMENTARY MATERIAL
To view supplementary material for this article, please visit https://doi.org/10.1017/aer.2020.93
APPENDIX
The mesh used for the grid independence study of each SWLBI case is listed below.
Mesh used for grid independence study
