A High-Order Hybrid Numerical Scheme for Hypersonic Flow Over A Blunt Body

A hybrid scheme is developed for direct numerical simulations of hypersonic flows over a blunt body. The scheme switches to the first-order AUSMPW+ scheme near the bow shock to provide sufficient dissipation to handle the carbuncle phenomenon. In the smooth part of the computational domain, a sixth-order central scheme with an eighth-order low-pass filter is adopted to provide high spatial accuracy to resolve turbulence. The hybrid scheme is shown to be able to obtain smooth and accurate predictions of laminar hypersonic flows over a blunt body. Using the hybrid scheme, a direct numerical simulation of a Mach 6 hypersonic flow over a circular cylinder is conducted. The result shows the turbulent structures in the near-wall region are well resolved by the hybrid scheme, and the bow shock is also captured without introducing any numerical oscillations. With the boundary layer transition on the cylinder’s surface, the simulation indicates that the heat flux peak shifts from the stagnation point to the transitional zone and its peak value is increased by 50%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\%$$\end{document}.


Introduction
An accurate predicton of hypersonic flow over a blunt body is one of the key problems in high-speed aerodynamics due to its close relation to the aerothermal characteristics of hypersonic vehicles. The temperature at the surface of a re-entry vehicle can exceed 10,000 K at hypersonic speed, and peaks of temperature and wall heat flux are usually observed on the windward side of a hypersonic blunt body (Anderson 2006). As a result, structural failure could happen, due to effects such as ablation unless a thermal protection system (TPS) is employed. Designing the TPS relies on an accurate prediction of the wall heat flux (Schneider 2004(Schneider , 2006Wang et al. , 2012. An under-prediction of wall heat flux may lead to structural failure due to ablation, and an over-prediction of the wall heat flux would increase the volume and weight of the TPS and lower the overall performance of a 1 3 hypersonic vehicle. It is commonly believed that the wall heat flux reaches its maximum at the stagnation point of a blunt body, which is generally true when the flow is in the laminar regime, and many experimental and numerical studies have been conducted for laminar flows (Lees 1956;Kemp et al. 1955;Peery and Imlay 1988;Kopriva 1993;Kitamura et al. 2010;Rodionov 2017;Chen et al. 2018). However, experiments on re-entry vehicle models presented by Hollis and Collier (2008); Hollis (2012) show that, when the Reynolds number is large enough, a laminar-turbulent transition may take place in the boundary layer downstream of the stagnation point, and the heat transfer in the transitional region could be equivalent or even stronger than that at the stagnation point. When calculating the wall heat transfer rate or skin friction in these cases, a pure laminar flow solution would be far from accurate. Therefore, a study of transitional and turbulent flows over a hypersonic blunt body is of a particular importance for both fundamental high-speed aerodynamics and engineering applications.
Direct numerical simulation (DNS), which resolves turbulence down to the Kolmogorov scale and provides comprehensive flow details and accurate statistics, is highly popular in turbulence research (Moin and Mahesh 1998). DNS has been widely applied to study turbulent boundary layer and laminar-turbulent transition. Zhong ( , 1998; Zhong and Tatineni (2003); Zhong and Wang (2012) conducted a series of DNS of hypersonic flows over parabolic shaped bodies with hypersonic incoming Mach numbers. A shockfitting method was used to capture shock-waves and a high-order compact upwind scheme was adopted to resolve the boundary layer. Their research focused on the receptivity and instability of hypersonic boundary layers, but the Reynolds numbers were relatively low (below 10, 000 Zhong 1998;Zhong and Tatineni 2003), and therefore only the early stage of the boundary layer transition was captured. Li et al. (2008) studied a M ∞ = 6 flow over a blunt cone using a DNS approach, and a complete transitional process on both windward and leeward sections were obtained by using a seventh-order WENO scheme. Very fine flow structures were resolved in their DNS. However, their computational domain was limited to a thin layer near the wall, and neither the nose of the cone or the shock-wave were included in their simulations. A range of DNS studies of turbulent boundary layer flow up to M ∞ = 12 were conducted by Martin (2007) and Duan et al. (2010Duan et al. ( , 2011, and the characteristics of hypersonic turbulent boundary layers were systematically discussed. However, their research was limited to a flat-plate configuration. Lagha et al. (2011) has also conducted a series of DNS studies of flat-plate boundary layers up to M ∞ = 20 , and they reported that the turbulence characteristics of hypersonic boundary layers are essentially incompressible. To the best of the authors' knowledge, a DNS study of hypersonic transition to a fully developed turbulent boundary layer over a blunt body with the shock-wave being resolved has not been reported in the literature.
For DNS of high-speed turbulent flows, high-order low-dissipative shock-capturing schemes are usually used to capture shocks and resolve small-scale fluctuations. Typical approaches include the seventh-order WENO schemes (Shu and Osher 1988;Martin et al. 2006;Capdeville 2008), seventh-order low-dissipative Monotonicity-Preserving (MP) scheme (Fang et al. 2013), and, more recently, high-order targeted ENO schemes (Fu et al. 2016(Fu et al. , 2017. High-order central schemes are also sometimes used for DNS of shock-free high-speed flows (Pirozzoli 2010;Pirozzoli and Bernardini 2011). Generally, high-order shock-capturing schemes have higher resolution and lower dissipation than some second-order and third-order upwind-biased schemes at the high-wavenumber end, so that they can efficiently resolve small-scale turbulent fluctuations. However, the low-dissipative property may cause problems in blunt body flows, especially when the Mach and Reynolds numbers are high. In fact, these problems have been realized and discussed for decades. Peery and Imlay (1988) first reported the notorious carbuncle phenomenon when studying a hypersonic blunt body problem using Roe's scheme (Roe 1981). Further studies (Robinet et al. 2000;Elling 2009;Kitamura et al. 2012) discovered that anomalous solutions, including bad shape of the shock-wave as well as oscillations of flow field and wall variables, may emerge when solving a strong bow shockwave in front of a blunt body. Quirk (1994) gave an explanation of this unphysical phenomenon by pointing out that the dissipation provided by a Riemann solver is insufficient to achieve a numerically stable solution. Sanders et al. (1998) proposed an entropy fix for Roe's scheme to cure the carbuncle phenomenon. Liou and Steffen (1993) proposed an advection upstream splitting method (AUSM) which is more robust than Roe's scheme near the strong shockwave by splitting the mass fluxes and pressure components separately. Recently, some new artificial viscosity methods have been proposed by Rodionov (2017) and Chen et al. (2018) to overcome the carbuncle phenomenon. Pandolfi and D'Ambrosio (2001) further confirmed that the carbuncle phenomenon is related to a variety of factors such as flow conditions, reconstruction scheme, density and aspect ratio of the mesh. Kitamura et al. (2010) gave a detailed evaluation of popular numerical methods used to simulate hypersonic blunt body flows, including the family of AUSM schemes (Liou and Steffen 1993;Kim et al. 2001), the modified Roe scheme (Kim et al. 2003), the Harten-Lax-van Leer-Einfeldt (HLLE) scheme (Einfeldt 1988) and the Godunov methods. Among all the tested schemes, the AUSMPW+ scheme proposed by Kim et al. (2001) presents a superior performance against others in terms of the prediction of wall heat flux, although the study was limited to first-order and second-order schemes. Nishikawa and Kitamura (2008) combined two different Riemann solvers with a rotated Riemann solver approach, and they demonstrated that the combined Riemann solver was able to provide carbuncle-free results with a good boundary layer resolution due to the reduced numerical dissipation. Their analysis was also limited to the first-order and secondorder schemes.
As it has been identified that numerical dissipation plays a key role in solving highspeed blunt body flow, it is necessary to check if a high-order scheme, which normally has a low level of dissipation, can give a stable solution for this type of flow. In the present study, high-order shock-capturing schemes are first evaluated for hypersonic flow over a blunt body. Two classical and widely accepted flux splitting methods, the Steger-Warming scheme (Steger and Warming 1981), which was once reported to be free of the carbuncle phenomenon (Pandolfi and D'Ambrosio 2001), and the AUSMPW+ scheme (Kim et al. 2001), which presented the best performance in the study of Kitamura et al. (2010), are evaluated in the present study by combining with high-order flux reconstruction methods. The tested high-order schemes are unable to give a satisfactory solution at high Reynolds numbers. Therefore, we propose a hybrid method that use a low-order shock-capturing scheme around the strong bow shock-waves and a high-order central scheme inside the boundary layer to accurately resolve near-wall turbulence. We will demonstrate that the new hybrid scheme can give stable and well converged solutions for laminar flows at all the studied Reynolds and Mach numbers. A DNS of hypersonic flow around a circular cylinder at M ∞ = 6 is then conducted by using the proposed hybrid scheme, and both the bow shock-wave and stagnation point were included in the simulation. The boundary layer transition is triggered by using wall blowing and suction at the windward side of the cylinder, and a fully developed turbulent boundary layer is achieved. Compared to a laminar flow, the peak wall heat flux shifts from the stagnation point to the transitional region, and its value is increased by 50%.

Governing Equations
The compressible Navier-Stokes equations in a general, time-invariant, curvilinear coordinate system are numerically solved using the finite difference method. The set of equations are written in a strong conservative form as, where J = | (x,y,z) ( , , ) | , and tr represents the transpose of the matrix. The Navier-Stokes equations are non-dimensionalised with the freestream velocity, u ∞ , the temperature, T ∞ , the density, ∞ , the viscosity, ∞ , and a suitable length scale, here the radius of the cylinder, r 0 . The reference Mach number M and the Reynolds number Re are , The pressure, p, is non-dimensionalised with ∞ u 2 ∞ , and is related to the temperature, T and the density, , via the ideal gas law, p = T M 2 . The total energy is calculated as E = The convection and diffusion flux vectors in Eq. 1 can be respectively written as, and For the convenience of representing matrices and vectors, (x 1 , x 2 , x 3 ) , (u 1 , u 2 , u 3 ) and ( 1 , 2 , 3 ) are set equivalent to (x, y, z), (u, v, w) and ( , , ) , respectively. The standard Einstein summation notation is used. In Eqs. (2) and (3), the grid transformation is used for all metric coefficients. The contravariant velocity components are written as The stress tensor, ij , and the heat flux vector, ̇q i , are respectively expressed as,

Numerical Method
For the calculation of the convection terms, Eq. 2, with an upwind-biased scheme, a flux vector splitting (FVS) procedure is required to split the convection flux into its upwind and downwind parts, as ̂ =̂ + +̂ − . Two different flux splitting methods, the Steger-Warming scheme (Steger and Warming 1981) and the AUSMPW+ scheme proposed by Kim et al. (2001), are studied. The Steger-Warming scheme is efficient and has proved to show good positivity preserving properties (Gressier et al. 1999;Witherden and Jameson 2018), and was once believed to be free of the carbuncle phenomenon (Pandolfi and D'Ambrosio 2001). The AUSMPW+ scheme is an important member of the AUSM family (Liou 1996(Liou , 2006(Liou , 2010, which was designed specifically for hypersonic flows. The AUSMPW+ scheme has drawn plenty of attention in the high-speed community and has shown good performance in overcoming the carbuncle phenomenon (Kitamura et al. 2010;Chen et al. 2018). The details of the Steger-Warming and AUSMPW+ schemes in the curvilinear coordinate system are given in the Appendices A and B, respectively.
To approximate the derivative, ̂ , on the i th node of a one dimensional grid, = 1, 2, ., i, .., N , a numerical reconstruction of ̂ at the interface locations on the two sides of the i th node is required, and the numerical difference on the node can be expressed as, are the reconstructed values at the interfaces between nodes i − 1 and i and between nodes i and i + 1 , respectively.
For a first-order reconstruction, the left and right variables at the interface i + 1 2 simply take the values from the i th and (i + 1) th nodes, respectively. For high-order reconstructions, however, the flux vectors at the interface are reconstructed with a high-order approach using multiple nodes. In the present paper, two high-order shock-capturing schemes, the seventh-order WENO (WENO7) scheme (Balsara and Shu 2000) and seventh-order MP (MP7) scheme (Suresh and Huynh 1997;Li and Jaberi 2011) are incorporated. Details of the implementation of the WENO7 and MP7 schemes can be found Appendices C and D, respectively. A local characteristic decomposition (Qiu and Shu 2002) is adopted for highorder schemes to improve the quality of the solutions, but this is not needed for the firstorder scheme.

3
The diffusion terms, Eq. 3, are solved using a sixth-order central scheme. The primitive variables u i and T are firstly differentiated, and the stress tensor as well as the heat flux vector are then formed at each node. The diffusion terms are then solved by differentiating the stress and heat flux with another application of the first-derivative solver. This method is more efficient than the direct calculation of second-order derivatives, although the later method can be more numerically stable. After all the spatial terms are solved, a three-step third-order total variation diminishing Runge-Kutta method is used for temporal integration (Gottlieb and Shu 1998). Close to the boundaries, to avoid the stencil of schemes go beyond the boundary, we shorten the width of the stencils by lowering the order of accuracy and using the biased stencil. At boundary nodes, the six-order scheme is reduced to a fourth-order biased scheme, and the WENO7 and MP7 schemes are replaced with a thirdorder biased scheme.
The flow solver used in the present study is ASTR, an open-source computational fluid dynamics code previously applied in DNS of various compressible turbulent flows (Fang et al. 2014(Fang et al. , 2015Ni et al. 2016;Fang et al. 2017;Ni et al. 2018Ni et al. , 2019Fang et al. 2020;Liu et al. 2021;Yang et al. 2022).

Laminar Hypersonic Flow Past a Circular Cylinder
Different FVS procedures and reconstruction schemes are first assessed in the case of a 2-D hypersonic flow past a circular cylinder in the laminar flow regime. A sketch of the flow is shown in Fig. 1a, and we can see that a strong bow shock appears in front of the body with a hypersonic incoming flow of M ∞ ≫ 1 . The grid, computational domain and . The grid is uniformly distributed in the circumferential direction, but concentrated towards the wall using the stretching function proposed by Vinokur (1983) , so that the grid spacing reaches the minimum, Δ min , at the wall. The stretching function is defined as, where i = 0, 1, ..., N , N is the total number of nodes in the wall-normal direction, d(i) is the normal distance of the i th node to the wall, D is the furthest wall-normal distance of the grid to the wall. The stretching factor is defined by, where Δd 1 is the distance of the first node away from the wall. Three cases are studied to assess the numerical schemes, whose flow conditions are listed in Table 1. The sets of grid for all the three cases have the same size of 251 × 251 , and the grid Reynolds numbers, Re grid = ∞ u ∞ Δ min ∞ , satisfy the criterion of Re grid ≤ 3 proposed by Klopfer and Yee (1988), as listed in Table 1. The flow condition of Case 1 are taken from the experiment of Tewfik and Giedt (1960). Case 2 has the same flow condition as the experiment of Holden et al. (1988), with higher Mach and Reynolds numbers than Case 1. To assess the effect of the Reynolds number, a third case (Case 3) is investigated with the same Mach number as Case 1 but the Reynolds number of Case 2. Case 1 and Case 2 have also been studied by Zhong (1998) and Peery and Imlay (1988) to assess the shock fitting method and to report the carbuncle phenomenon.
The wall heat transfer coefficient, C h , and pressure coefficient, C p , obtained with different combinations of FVS and reconstruction schemes for Case 1 and Case 2 are shown in Figs. 2 and 3, respectively, where C h and C p are defined as,  Figure 2 shows that, for Case 1, C h and C p obtained by the four evaluated schemes lie on top of each other. They all agree well with the former calculation of Kopriva (1993), although all the numerical simulations predict a lower level of heat transfer rate than the measurement. The results of Case is shown in Fig. 3, which we can observe the disparities of the different schemes near the stagnation point ( = 0 • ). The wall heat transfer rates from high-order reconstruction methods present asymmetric wiggles around the stagnation point, which is known as the carbuncle phenomenon. On the contrary, the AUSMPW+ scheme with the first-order reconstruction can give a smooth result, which is consistent with the result of Kitamura et al. (2010) The temperature fields of Case 1 and Case 2, obtained with the different schemes, are presented in Figs. 4 and 5. For Case 1, all the schemes exhibit smooth temperature fields, and the results from different schemes are similar with each other. However, for Case 2, obvious differences can be observed from the results of different schemes, as shown in Fig. 5. The AUSMPW+ scheme combined with both the MP7 and WENO7 schemes generates severe oscillations on the temperature field. The SW scheme, on the other hand, gives a smoother solution, but an abnormal temperature distribution can be seen near the stagnation line, resulting in an error in predicting the heat transfer rate near the stagnation region, as shown in Fig. 3   Clearly, all the studied schemes produce fairly good results for Case 1, but only the AUSMPW+ scheme with a first-order reconstruction is able to give a smooth result for Case 2. To clarify the reason for the difference, Case 3, whose Mach number is the same as Case 1 but Reynolds number is the same as Case2, is further adopted to identify the main factor for the carbuncle phenomenon. The temperature fields of Case 3 from different schemes are presented in Fig. 6, in which similar abnormalities and oscillations as for Case 2 are observed. Therefore, the carbuncle phenomenon tends to occur at a higher Reynolds  (Pandolfi and D'Ambrosio 2001), which indicates the critical role of dissipation. At a lower Reynolds number, the viscous dissipation provided by the flow is stronger, and the carbuncle phenomenon is less likely to happen. The numerical dissipation also play a similar role in avoiding numerical oscillations. For all the studied cases, the first-order upwind reconstruction gives smooth results, because of its higher numerical dissipation. The high-order reconstruction, on the other hand, produces some numerical wiggles. Even for the low Reynolds number case, oscillations can be observed behind the shock-wave, as shown in Fig. 4 (b), although they are limited to a small region. Based on the research of Kim and Kim (2005), a multidimensional limiting process may be able to further overcome the post-shock oscillations.

A Hybrid Scheme for Hypersonic Blunt Body Flow
Sect. 2 shows the critical importance of the numerical dissipation to avoid the carbuncle phenomenon for hypersonic flows over a blunt body. The first-order AUSWMPW+ scheme is able to give smooth results. However, to resolve fine turbulent structures in a DNS, excessive dissipation would destroy the small-scale turbulence structures and the simulation would become under-resolved (Mittal and Moin 1997;Larsson et al. 2007). A similar problem is encountered for DNS of shock-wave/turbulence interaction, where a high level of numerical dissipation is needed to preserve the monotonicity near the shock-wave , but a low-or none-dissipative scheme is required to resolve turbulence efficiently. A large number of hybrid schemes were therefore proposed to solve the problem, and a review of these schemes can be found in ref (Pirozzoli 2011). The basic idea of these schemes is to use a high-dissipative shock-capturing scheme near a shock-wave and a low-dissipative upwind scheme or none-dissipative central scheme in the smooth region. Examples of such a scheme are the hybrid upwind/WENO schemes (Adams and Shariff 1996 Ren et al. 2003;Chao et al. 2009;, and the hybrid central/shockcapturing schemes (Visbal and Gaitonde 2005;Touber and Sandham 2009;Sjögreen et al. 2019). A shock sensor is essential for a hybrid scheme to distinguish the flow field between smooth and non-smooth parts. Several shock sensors have been proposed, and have been compared in Pirozzoli (2011) and Zhao et al. (2020). The Ducros et al. (1999), which identifies discontinuities using the local pressure gradient and velocity divergence, is widely used by many hybrid schemes.
For the simulation of a hypersonic flow over a blunt body at a high Reynolds number, an even higher level of numerical dissipation is required around the strong bow shock to overcome the carbuncle phenomenon, compared to the simulations of oblique shock-waves. From the findings of Sect. 2, the first-order AUSMPW+ scheme can provide a non-oscillatory solution in a high Reynolds number laminar flow. However, the scheme is clearly too dissipative to use for DNS. An important characteristic of hypersonic blunt flow is that the bow shock is distant from the wall boundary layer where turbulence exists, and there is no direct interaction between the bow shock and turbulence, which would naturally fit a hybrid scheme.

A Hybrid Central-AUSMPW + Scheme
Considering the high level of numerical dissipation required around the bow shock and the shock being far from the boundary layer, a new hybrid central-AUSMPW+ scheme is proposed for DNS of hypersonic flows over a blunt body. The detail of its implementation is presented for a 1-D case, and the extension of the scheme to 2-D (resp. 3-D) is carried out by directional discretisation. The hybrid scheme is implemented according to the following steps: • The Ducros sensor is adopted to identify the shock-wave, shown as node i in Fig. 7. • The shock zone is created by including four guard nodes on each side of node i. • The convection flux vectors in the shock zone (from i − 5∕2 to i + 7∕2 shown in Fig. 7) are calculated with the first-order AUSMPW+ scheme to provide sufficient dissipation to control the carbuncle phenomenon. • In the shock free zone, the scheme switches to a sixth-order central scheme to reconstruct the convection flux at the interface vector as, • To remove the small-scale aliasing errors, an eighth-order low bypass filter is applied to the vector, , in the shock free zone at the end of each Runge-Kutta sub-steps as, The high-order filter limits the filtering operation only at high wavenumbers, which minimises the damping of the flow structures, and it is, therefore, popular in DNS and LES of compressible flows (Gaitonde et al. 2000;Visbal and Rizzetta 2002;Kawai and Fujii 2008;Wang et al. 2017). In the present study, the filter is applied to the flow variables in the computational coordinates ( , , ) . The dissipation property of the eighth-order low-pass filter can be seen in refs (Visbal and Gaitonde 1999;Yee and Sjögreen 2011;Hadjadj et al. 2012).
For the proposed hybrid scheme, there is no need to incorporate the local characteristic decomposition, as only the first-order scheme is used near discontinuities. Four guard nodes on each side of an identified discontinuity are adopted here, which ensures that the stencils of the sixth-order central scheme and eighth-order filter won't go across the identified discontinuities. The risk is that the first-order scheme may be activated in the smooth part of the flow and potentially damp the flow structures. For the cases studied in the present paper, however, the risk can be avoided because there is no direct interaction between the shock and turbulence. For a multi-dimensional case, once a node is identified in the shock zone, the first-order scheme is applied in all the directions. Note that the proposed hybrid scheme is not recommended for the flow with a direct interaction between shockwave and turbulence, as the first-order scheme might be triggered in the turbulent region and its strong dissipation could damp turbulent structures.

Assessment of the Hybrid Method
The hybrid scheme is first assessed for Case 2 (laminar case). Its flow conditions are listed in Table 1. The distribution of the Ducros sensor's function and the shock zone are shown Sixth-order central scheme for shock free zone in Fig. 8. We can see that the shock is well identified by the shock sensor, and the shock zone is far from the wall of the cylinder. The temperature field obtained by the hybrid scheme is compared with the MP7-SW and MP7-AUSMPW+ schemes in Fig. 9. The hybrid scheme generates a smooth solution, suggesting that the numerical instability has been avoided. The wall pressure coefficients and the wall heat transfer rate are further compared in Fig. 10. The pressure coefficients from different schemes agree with each other well, and also for the pressure fields shown in Fig. 11. For the wall heat transfer rate, both the MP7-SW and the MP7-AUSMPW+ schemes show some1 abnormal wiggles near the stagnation point, while the hybrid method gives a smooth result. This is due to the fact that the hybrid scheme returns to the first-order AUSMPW+ scheme near the shock-wave. Even if the dissipation becomes higher near the shock-wave, the high-order shock-capturing scheme is still not stable enough to overcome the numerical instability, as explained by Quirk (1994).
The Reynolds number of Case 2 is not high enough to trigger a laminar-turbulence transition. Therefore, a new case, named Case 4, is added, and its flow conditions are listed in Table 2. The Reynolds number is of one order of magnitude higher than that of Case 2. The temperature, T ∞ , is set to be the temperature of air at an altitude of 30km based on the standard atmosphere model.
The case is first assessed in a 2-D laminar flow to test the hybrid method. Two sets of grids are used to conduct a grid independence study. The grid sizes are 256 × 256 and 512 × 512 , and the Re grid for these two grids are 15 and 7.5, respectively. The temperature fields obtained with the two grids are shown in Fig. 12, and the result obtained by the MP7-AUSMPW+ scheme on the coarse grid is given for comparison. The MP7-AUSMPW+ scheme generates strong numerical oscillations, similar to the ones observed in Case 2. The hybrid scheme, however, presents smooth results for both grids, indicating its good numerical stability.
The wall heat transfer rate C h profiles are presented in Fig. 13. For the case of hypersonic flow, it is generally required Re grid ≤ 3 , in order to obtain an accurate wall heat   Klopfer and Yee (1988) However, by using the high-order central scheme in the near-wall region, the hybrid scheme ensures the prediction of wall heat flux at the stagnation point converges to the value given by the Fay-Riddell equation , C h,FR , with a much larger Re grid . The difference between the stagnation heat transfer rates from the two sets of grids is less than 1.5%. The time history of the density residual is presented in Fig. 14, where the residual is defined as the maximum difference of the density fields between two computational steps at n and n + 1 as, It shows that the hybrid scheme can reduce the residuals to the magnitude of 10 −14 , while the MP7-AUSMPW+ scheme gives a much higher level of residuals, due to the numerical oscillations.
The computing times cost per step by different schemes for Case 4 with the 256 × 256 mesh are compared in Table 3. The computing times are obtained using 8 cores of an AMD EPYC 7742 Processor. From Table 3 we can see the hybrid scheme costs the least time, compared with MP7 and WENO7 schemes. The WENO7 scheme is the most computationally expensive scheme due to the reconstruction of fluxes in multiple sub-stencils.

DNS of Hypersonic Flow Around a Circular Cylinder
The test conducted for a laminar flow over a cylinder at a high Reynolds number show that the proposed hybrid scheme is able to prevent numerical instability caused by the carbuncle phenomenon, and mesh independence can be achieved on a relatively coarse mesh, due to the high-order central scheme being used in the near-wall region. In this section, the hybrid scheme developed in Sect. 3 is applied to conduct a DNS of a hypersonic flow over a blunt body.
The flow configuration is the same as for Case 4 and a cross-section of the cylinder is schematically shown in Fig. 15. The computational domain and boundaries on the cross-section plane are the same as those laminar cases shown in Fig. 1. A periodic blowing and suction is applied on the cylinder's surface in the region 29 • ≤ ≤ 31 • to trigger boundary layer transition, as shown in Fig. 15. The wall-normal velocity of the blowing and suction, v bs , is given as, where A bs is the amplitude of the blowing and suction. The terms f bs , g bs , and h bs express the variations of v bs in the circumferential ( ), spanwise (z), and temperal (t) directions, respectively. They are given as: where Lz is the spanwise width of the computational domain, a and b are the circumferential locations of the beginning and the end of the blowing and suction region, and 0 , 1 and 2 are three random numbers ranging from 0 to 1. The parameters of the wall blowing and suction can be found in Table 4. Equation 15 was first proposed by Rai et al. (1995) to trigger a boundary layer transition for the DNS of a supersonic boundary layer and has been widely adopted to study compressible boundary layers (Gao et al. 2005;Sayadi et al. 2013;Fang et al. 2020;Di Renzo and Urzay 2021). It can be seen that the wall blowing and suction used in the present case have two modes in the streamwise direction, one mode in the spanwise direction, and two temporal modes. No extra mass flux is introduced into the flow field from the wall since the net flow rate of the disturbance is zero. The wall-parallel velocity components in the blowing and suction region are set to zero. The isothermal no-slip boundary condition is applied to the rest of the cylinder surface, and the wall temperature is three times of incoming freestream temperature. A periodic boundary condition is applied in the spanwise direction, and the spanwise length of the computational domain, Lz, is set to 0.1. A body-fitted mesh with 3256 × 320 × 512 nodes in the circumferential, wall-normal and spanwise directions, respectively, is used to mesh the domain (see Fig. 16). In the circumferential direction, 3, 000 nodes are distributed in the upper part of the computational domain from = 30 • to 85 • to resolve turbulence. The flow in the bottom part of the domain is supposed to remain laminar, and the mesh resolution is similar to the one used for the 2-D laminar cases. The mesh is hyperbolically stretched in the wall-normal direction, and the grid Reynolds number is Re grid = 10 . In the spanwise direction, the mesh is uniformly distributed. The mesh resolutions of the first node away from wall in the local wall units at = 60 • , where the turbulence is fully developed, is 19.00, 0.62 and 8.70 in the streamwise, wall-normal and spanwise directions, respectively. The ratio between the effective mesh resolution to the local Kolmogorov scale, Δ∕ , is 2.74 at the wall and 0.75 near the edge of the boundary layer, in which Δ is defined as the cube root of the mesh sizes in the three directions, and is estimation via the resolved dissipation in the DNS.
The flow field is initialised by a converged 2-D laminar flow field at the same flow condition, which is copied across the spanwise direction. A constant time step, Δt = 3.5 × 10 −5 , is used for the simulation, and the corresponding Courant-Friedrichs-Lewy (CFL) number is around 0.8. The temporal evolution of the pressure, temperature and x-velocity component on a sample point located at = 50 • , d = 0.0015 are shown in Fig. 17, in which d is the normal distance to the wall. The fluctuations grow after the initial stage at about t = 2 , indicating that the sample point is in a turbulent state. A statistically steady state is reached at about t = 6 , after which the data are collected to compute the statistics. The simulation ends at t = 19 . Figure 17 shows that the mean and high-order statistics are stable and independent from the number of samples near the end of the simulation. (15)

Instantaneous Flow Field
The DNS data is analysed to study the flow characteristics. Fig. 18 shows the instantaneous density schlieren extracted from the middle cross section, at z = 0.05 , where the schlieren is calculated according to the local gradient of density field as, Figure 18a shows the bow shock-wave in the inviscid region of the flow field and the turbulence in the thin layer near the wall. By zooming into the near-wall region (see Figure 18b, it is observed that the inner turbulence is sharply separated from the outside essentially irrotational flow. The turbulence structures are observed as bulges and rings, which is associated with the so-called hairpin or horseshoe structure of wall turbulence (Adrian et al. 2000). Similar structures were also observed for compressible simulations over flat plate boundary layers (Pirozzoli et al. 2008;Wu and Martín 2007;Duan et al. 2011). Figure 18 also confirms that the bow shock-wave is sufficient far away from the cylinder boundary layer, and there is no direct interaction between the shock-wave and turbulence.
The turbulence structures can be further observed through near-wall velocity fluctuations, known as streaks (Kline et al. 1967). The velocity vector is decomposed into the circumferential (also streamwise), wall-normal and spanwise components, respectively, and the streamwise velocity fluctuations on a near-wall slice are shown in Fig. 19. The slice is located at d = 0.5422 × 10 −4 , which corresponds to d + = 23.7 (there the wall unit is based on the wall variables at = 60 • ). Fig. 19a shows that the streamwise elongated velocity streaks are generated downstream of = 35 • . The streaks are mainly observed around = 45 • . At a further downstream location at about = 75 • (shown in Fig. 19b, the streaks become wider and their fluctuating intensity becomes weaker, which means the turbulence is weakening due to the effects of a favourable pressure gradient and curvature. The turbulent coherent structures in the transitional zone and fully developed zone are shown in Fig. 20. The coherent structures are visualised with the iso-surfaces of swirling strength, ci , defined as the imaginary part of the complex eigenvalue pair of the velocity gradient tensor (Zhou et al. 1999). In the transitional zone, diamond shaped vortex packages can be seen immediately downstream of the wall blowing and suction, indicating a typical bypass transition triggered by the wall blowing and suction (Durbin and Wu 2006), and the distribution of coherent structures present a clear intermittency. Hairpin-like structures with counter rotating legs and ring-shaped head can be observed as highlighted in Fig. 20 (a), which is consistent with the characteristics of coherent structures present in incompressible boundary layers (Zhou et al. 1999;Green 2007;Jeong et al. 1997). In the fully developed zone at about = 60 • , a forest of coherent structures with counter-rotating streamwise vortices are present in the boundary layer, which is also observed in the flow of flat plate boundary layers (Wu and Moin 2009;Pirozzoli et al. 2008;Ringuette et al. 2008). The transporting effect caused by the vortical motion of these structures is responsible for the high skin-friction and wall heat flux.

Mean Flow Field
Mean temperature, ⟨T⟩ , and pressure, p , fields are presented in Fig. 21, where the Reynolds-averaged mean pressure, p , is obtained by averaging p in the temporal and spanwise directions, and the Favre-averaged mean temperature, ⟨T⟩ , is calculated as, ⟨T⟩ = T∕ .

3
The fluctuations from the Reynolds-averageing and Favre-averaging are denoted with ′ and ′′ , respectively. It can be seen in Fig. 21 that the mean flow field is smooth and oscillationfree, which indicates that the proposed hybrid scheme works well for DNS of a hypersonic flow past a blunt body. Compared with the laminar result shown in Fig. 12, no obvious difference is observed for the mean flow field. This is expected because the turbulence only exists in a thin layer near the wall.
The mean and instantaneous coefficients of wall pressure, C p , heat transfer, C h , and skin friction, C f , are compared with the laminar results in Fig. 22, where C f is defined as,  in which n i and t j are the normal and tangential unit vectors of the cylinder's surface. The mean wall mean pressure does not show any obvious change when the boundary layer transitions from laminar to turbulent, as observed in Fig. 22a. For the wall heat transfer and skin friction, however, the transition from laminar to turbulent pose a (17) Fig. 14    huge impact to these two coefficients. The peak of C h moves from the stagnation point ( = 0 • ) to the transitional region at = 40 • , and its peak value is predicted to increase by about 50% , from C h = 0.0019 to C h = 0.0028 . The instantaneous wall heat transfer coefficient fluctuates intensively, demonstrating the impact of turbulence on the instantaneous heat transfer. It can be seen that the skin friction coefficient is increased by the laminar-to-turbulent boundary layer transition, but the location of the peak of C f is not strongly affected. However, as a consequence of the turbulent transition, the skin-friction is seen to fluctuate significantly. The comparisons of the mean streamwise velocity and mean temperature profiles at = −60 • and = 60 • , where the boundary layer states are respectively laminar and turbulent, are shown in Fig. 23. Due to the enhanced momentum and heat exchanges caused by turbulence, both the streamwise velocity and temperature at = 60 • have larger values in the near-wall region and smaller values in the outer part, compared to the corresponding laminar profiles. This leads to higher local skin friction and wall heat transfer rate. Outside the boundary layer ( d > 0.02 ), the laminar and turbulent profiles are close to each other, which explains why the mean flow field shown in Fig. 21  1 3 profile at = 60 • presents a law of the wall with a linear layer and a log layer, although the log layer level is higher than for the classical log law. The upper shift of the log layer is a common characteristic of a boundary layer with a favourable pressure gradient (Mcdonald 1969), or over a convex wall . The further study of turbulent boundary layer is conducted by analysing the turbulence kinetic energy (TKE) budget equations given as Adumitroaie et al. (1999), is the convection term, is the pro- ij is the dissipation term. For the fully resolved turbulence, the sum of the right-hand-side of Eq. 18 (i.e., the balance term) should be close to zero. Therefore, by checking the level of the balance term, we could identify if a DNS is fully resolved.
The profiles of the TKE budget terms at = 60 • are plot in Fig. 25, from which we can see that the characteristic of the TKE budget is very similar with that in a equilibrium turbulent boundary layer. The TKE is mainly produced in the buffer zone, transported to the near-wall region by the turbulent and viscous transport terms, and consumed by the dissipation term (Kim et al. 1987). In the outer part of the boundary layer, the budget shows a local balance between production and dissipation. The balance term is also shown in Fig. 25, and we can confirm that the term is close to zero, meaning TKE budget terms is well balanced and the turbulence is fully resolved by the adopted numerical scheme.

Summary and Conclusion
Several flux vector splitting and reconstruction schemes are first assessed in the paper for simulating 2-D hypersonic laminar flows past a circular cylinder. It is found that the tested high-order schemes are not able to give a satisfactory result because of the carbuncle phenomenon, especially when the Reynolds number is high.
To provide a stable high-order scheme for DNS of hypersonic flows past a blunt body, a hybrid central/AUSMPW+ scheme is introduced. The scheme switches to the first-order AUSMPW+ scheme in the shock zone to provide sufficient dissipation to overcome the carbuncle phenomenon. In the smooth part of the computational domain, a sixth-order central scheme with an eighth-order filter is used to provide a high resolution numerical approach to resolve fine turbulence structures. The assessment of the hybrid scheme for laminar flows shows that the scheme can give stable and non-oscillatory solutions. Mesh independence is observed for coarser grids than with the other schemes tested because of the high-order scheme being used in the near-wall region. Note that the hybrid scheme is proposed for the flow past a blunt body, and it is not recommended for the problems like shock-wave/boundary layer interaction, as the first-order scheme might be triggered in the turbulent zone and over-damp the structures of turbulence. The proposed hybrid scheme is used to perform a DNS of a M ∞ = 6 flow over a circular cylinder. By applying wall blowing and suction on the upper surface of the cylinder, a bypass transition to fully developed turbulence is obtained. The detached bow shockwave is well captured, and a non-oscillatory mean flow is obtained. The fine turbulent structures are resolved in the near-wall region, and the turbulence structures are similar to those observed in incompressible boundary layer simulations. In the transitional zone, diamond-shaped vortical packages and hairpin-like coherent structures are observed, and a large number of streamwise elongated vortices are shown in the fully-developed zone. We show that the wall heat flux and skin-friction are greatly increased because of the boundary layer transition. The peak of the wall heat flux is seen to shift from the stagnation point to the transition zone, and its peak value is also increased by 50% , showing it is necessary to account for the boundary layer state in the design of a thermal protection system. This DNS describes a complete transition process for the case of a hypersonic flow over a circular cylinder. More in-depth analysis of the flow characteristics, including Reynolds analogy, will be carried out as future work. x + 2 y + 2 z . The positive and negative flux vector, ̂ ± is given as,

Appendix A: The Steger-Warming Flux Splitting Scheme
(3− )(̃± 4 +̃± 5 )c 2 2( −1) , P = 2p( − 1)̃± 1̃x (̃yw −̃zv) , ̃x i = x i ∇ , the local speed of sound is c = √ T M , and ̃± i , is calculated by, with = 0.04 . The eigenvalues of the Jacobian matrix of the Euler flux vector, i , is given as, After the split of the Euler flux vector at each grid's node, the positive and negative fluxes at the interface location between the i th and i + 1 th nodes, ̂ + i+1∕2 and ̂ − i+1∕2 , can be reconstructed by using a upwind-biased scheme using values from nodes within its stencil. The total flux vector, ̂ i+1∕2 , is then obtained through Eq.(A1). Finally, the convectional term in Eq. (1) at the i th can be calculated as,

Appendix B: The AUSMPW+ Flux Splitting Scheme
For the AUSMPW+, the flow variables at left and right sides of a node's interface are first reconstructed using upwind-biased and downwind-biased schemes, respectively, and the final flux vector is then calculated using the AUSMPW+ scheme by using the left and right variables as, where c 1∕2 is the speed of sound at the node's interface, and L,R and the pressure vectors L,R are expressed as, The enthalpy, H, is defined as, The subscripts, L, R, stand for variables at the left and right sides of a node's interface, calculated by the upwind-biased and downwind-biased schemes, respectively. It is noted that the coordinate transformation matrix, i , and its Jacobian, J, are also reconstructed onto the node's interface to preserve the geometric conservation. (B14) c 1∕2 = c 2 s ∕max(|U L |, c s ), for 1 2 (U L + U R ) > 0, c 2 s ∕max(|U R |, c s ), for 1 2 (U L + U R ) < 0, 1 3

Appendix C: Reconstruction Scheme of WENO7
For the WENO7 reconstruction,the positive flux at the interface, f + i+1∕2 , is reconstructed from sub-stencils as, in which the reconstruction from each sub-stencils, q k , is give as, and the weight, k , is expressed as, In the expression of k , C k and IS k , are respectively, the optimal weights and smoothness estimators to ensure the WENO scheme avoids reconstruction using a sub-stencil containing a discontinuity, and returns to the optimal linear scheme in the smooth region. In Eq. C16, is set to 10 −10 to avoid a zero denominator. The smoothness estimators IS k are calculated as, and the optimal weights, C k , are give as, For a reconstruction in a fully smooth region, the WENO7 scheme will have k = C k , and the scheme returns to the standard seventh-order linear upwind scheme.