key: cord-0068736-uknip3jr authors: Qu, Feng; Sun, Di; Liu, Qingsong; Bai, Junqiang title: A Review of Riemann Solvers for Hypersonic Flows date: 2021-10-21 journal: Arch Comput Methods Eng DOI: 10.1007/s11831-021-09655-x sha: 6fe778f107bd6abc8cdac915090db103f266b154 doc_id: 68736 cord_uid: uknip3jr In the design of the hypersonic airliner which can greatly shorten the flight time and conduct space travel, it is of great significance for a Riemann solver to accurately simulate hypersonic flows. In this survey, the research process on the Riemann solver for hypersonic flows is reviewed, including the constructions of the traditional Riemann solvers, improvements of the traditional Riemann solvers for the shock anomaly, and the importance of the all-speed Riemann solver for hypersonic flows. Moreover, constructions and applications of the all-speed Riemann solvers are presented. It is obvious that the Riemann solver should be capable of obtaining physical solutions at low speeds and be robust against shock anomaly near strong shocks at the same time in simulating hypersonic flows. With the explosive development of the aerospace technology in the last decades, it is common for an airliner to conduct cruising flight at transonic speeds [1, 2] , which would still cost as long as 12 h for an international flight from New York to London. Moreover, people's desire for space travelling can not be met by such traditional airliners [3, 4] . Therefore, it is in an urgent demand to design a new generation of airliners, which is capable of hypersonic flight and performing space travel. In terms of the long time cruising flight at hypersonic speeds, numerous studies have been conducted on this topic to indicate that the new generation airliner's thermal protection system will encounter severe heating loads in these cases [5] . In order to overcome this issue, large-scale thinwall structures and reusable lightweight materials, which have not been widely used in aerospace engineering, should be adopted to build the thermal protection system [6] [7] [8] [9] . Therefore, it is of great significance to accurately predict the aerodynamic force and heating flux of hypersonic airliners. In general, engineers are inclined to employ four different ways to predict the aerodynamic force and heating loads of the hypersonic airliners, including flight test [10] , empirical estimation methods [11] , wind tunnel experiments [12, 13] , and Computational Fluid Dynamics (CFD) method [14] [15] [16] [17] . Among these technologies, the empirical estimation method can not be applied to complex flows. The wind tunnel experiment method is incapable of reproducing real flight environment due to the scale limitation of the wind tunnel. The flight test can reveal the real flight condition by avoiding the limitation of the wind tunnel, and many flight tests for verification vehicles, such as X-43A [18] , X-51A [19, 20] , Hifire series [21] [22] [23] , HYFLEX [24] , and X-37B [25] , have been conducted to help researchers study the aerodynamic loads of the aircrafts at hypersonic speeds. However, these tests are so expensive that few researchers can adopt such methods to study the aerodynamic property of the hypersonic airliner in details. In comparison, the CFD technology is capable of conducting simulations covering the whole flight range of the hypersonic airliner and avoiding the scale limitation of the wind tunnel. Moreover, it can be conveniently applied to simulate three dimensional complex flows accurately. Consequently, CFD has been attracting more and more attentions in the past few years. For example, in order to promote the usage of CFD in the design process of the hypersonic vehicle, numerous studies have been carried out and show that the performance of Riemann solvers is of significance in the prediction of the hypersonic aerodynamic heating loads [26] [27] [28] . Nowadays, the widely used Riemann solvers in both academic and engineering areas are built upon the Godunov's idea to solve the compressible Euler equations [29] . However, the mathematical properties of the partial difference equations depicting the flows vary at different speeds [30] . In case the Mach number is subsonic, the equations satisfy the elliptic property. In contrast, the equations are with the parabolic property at transonic speeds and the hyperbolic property at supersonic/hypersonic speeds. Therefore, the traditional Riemann solvers are incapable of generating physical solutions at subsonic speeds in theory [31, 32] . Through complex mathematical manipulation, Guillard verified such conclusions based upon the asymptotic analyses of the compressible Euler equations [33, 34] . In order to make these traditional Riemann solvers applicable to simulate the low-speed flows, most of the CFD engineers and scholars employ preconditioning matrix technologies [35] [36] [37] [38] . Numerical experiments show that these methods can help traditional Riemann solvers avoid unphysical solutions at low speeds. However, they will encounter the global cut-off problem in supersonic and hypersonic problems. In other words, if hypersonic flows around an airliner are simulated, these preconditioning matrix methods have little influence on the traditional Riemann solvers and the low-speed flows in the boundary layers can not be accurately captured. As is well known, the airliner will encounter severe force and heating loads at hypersonic speeds, which are relative to the pressure, velocity, and temperature distributions in the boundary layers. Therefore, these methods may be incapable of obtaining physical solutions in studying the hypersonic airliners' aerodynamic properties. On the other hand, Roe indicated that a Riemann solver should avoid the appearance of the unphysical shock anomaly phenomenon to simulate the hypersonic heating load accurately [28] . Unfortunately, most widely used Riemann solvers, such as the Roe type schemes, the AUSM type schemes, and the HLL type schemes, tend to generate unphysical shock anomalies at hypersonic speeds [39] [40] [41] [42] [43] . The research development of the widely used Riemann solvers, especially those which are capable of obtaining physical solutions at low speeds and being robust against the unphysical shock anomaly at hypersonic speeds, will be reviewed in this paper. This manuscript is organized as follows. In the second section, the Gudunov method for the compressible Euler equations will be briefly presented. The constructions of the traditional Riemann solvers widely used today and their improved versions for the shock anomaly will be introduced in the third section and the fourth section, respectively. The fifth section will discuss the importance of the all-speed Riemann solver for hypersonic flows. After this, detailed survey on the constructions and applications of the all-speed Riemann solvers will be shown in Sect. 6. Remarks and conclusions will be presented in the last section. In general, the Euler equations are used to depict motions of inviscid flows [44] . In the two-dimensional case, they are written as follows: where F(q) and G(q) are the inviscid flux vectors, q is the conservative variable vector, and ρ is density, u is velocity, e is internal energy, and p is static pressure [45] . In conservative form, the Godunov method for the Euler equations can be written as follows in the finite volume form: where n indicates the time level and F i+1∕2 , are the numerical fluxes at the cell interfaces. Assume the data q n i at time level n is known, the solution at time level n + 1 will be advanced via Eq. (5) after computing the numerical fluxes. To solve this problem, Godunov adopted the direction splitting technology [29] and assumed two states interact with each other at the cell interface. Take the interface labeled i + 1/2 as an example, the initial conditions at time level n are Then the flux at the cell interface i + 1/2 can be obtained by solving the general Riemann problem with the following possible wave patterns as in Fig. 1 . More details on these wave patterns and the construction of the Godunov method can be found in Refs [46, 47] . As can be seen in Sect. 2, the Godunov scheme considers all of ten possible wave patterns as in Fig. 1 According to the homogeneity property of the Euler equations, Eq. (1) can be rewritten as follows: The indices A and B indicate the Jacobian matrices. In the two dimensional case, where indicates the ratio of specific heats. The Jacobian matrix B can be given in a similar manner. Based on the mathematical analysis, the compressible Euler equations satisfy the hyperbolic property and where K indicates the right eigenvector matrix of the Jacobian matrix A, is the characteristic matrix of the Jacobian matrix A, and K −1 is the left eigenvector matrix of the Jacobian matrix A. For the two-dimensional case, these matrices can be written as follows: The possible wave patterns for the evaluation of the Godunov flux in Eq. (5) [29] where a indicates the speed of sound and H is the total specific enthalpy. In order to avoid the expensive nonlinear operations of the Godunov scheme, Roe linearized the Euler equations by replacing the Jacobian matrix of the flux vector with the linearized Jacobian matrix as follows [48] : where the symbol ' ~ ' indicates the Roe average and Numerical cases in Ref [48] indicate that the Roe scheme is with a high resolution in capturing shock discontinuities and linear contact discontinuities. In addition to the Roe scheme, the HLL type schemes are also well-known FDS type Riemann solvers [49] [50] [51] [52] [53] [54] [55] . Among them, the HLLE scheme proposed by Einfeldt has drawn great attention because it can capture shocks sharply and avoid the appearance of the unphysical expansion shock. In the two-dimensional case, the HLLE scheme can be written as follows: where As can be seen in Eqs. (15) and (16) , HLLE is built on the two-wave model and ignores the middle linear degenerate field. Therefore, it can not accurately resolve contact discontinuities and shear waves. In order to address this issue, Einfeldt proposed the following HLLEM scheme by adding anti-diffusion terms to the HLLE scheme in the linear degenerate fields [56] . Compared with the FDS type Riemann solvers, the FVS type solvers were constructed by considering the salient features of the physical phenomena modeled by the Euler equations. Also, they are with the identifications of being more efficient than the FDS type solvers [57] . Nowadays, widely used FVS type Riemann solvers are the Steger-Warming Splitting scheme [58] and the van Leer's FVS scheme [59] . The constructions of these two schemes will be presented in the following sections. In order to follow the upwind property of hyperbolic Euler equations, Steger and Warming split the flux vector F into with where the eigenvalues satisfy the condition In order to satisfy Eq. (25), Steger and Warming proposed a splitting of the eigenvalues with definitions Borrowing Steger and Warming's idea, van Leer also split the Euler equations as Eq. (21) . However, he improved the Steger and Warming's scheme by proposing the following extra desirable properties: All eigenvalues of the matrix dF + dq are larger or equal to zero, and all eigenvalues of the matrix dF − ∕dq are smaller or equal to zero. 3. The functions of F ± should be continuous, The function dF ± dq should be continuous. 6. When |Ma| < 1 , one of the eigenvalues of the dF ± dq should disappear. 7. The flux vector F ± should be with the same order of the Mach number as the flux vector F , and the order should be as low as possible. Through mathematical analyses, the van Leer's FVS scheme which satisfies these properties aforementioned was constructed in the following two-dimensional form: The FDS type Riemann solvers and the FVS type Riemann solvers have been widely used in both scholar and engineering areas. However, the FDS solvers are much more time consuming than FVS, while FVS can not accurately capture contact discontinuities in comparison with the FDS type solvers [60] . In order to combine the merits of the FDS type Riemann solvers with the FVS type Riemann solvers, Liou and Steffen proposed the AUSM type schemes to split the flux vector F into a convective part and a pressure part in two-dimensional case as follows [61, 62] : where The AUSM Scheme In Ref [63] , Liou and Steffen proposed the first version of the AUSM-type scheme, which is called AUSM, based on the Liou-Steffen flux splitting procedure. In general, it can be written in the following form: In order to improve the AUSM scheme's robustness, Liou proposed the AUSM + scheme by unifying the definition of the speed of sound [64, 65] . This scheme is almost the same as the AUSM scheme except In Ref [28] , Roe and his co-workers found that it is challenging to compute hypersonic flows accurately. A good Riemann solver for hypersonic flows should be with a high boundary-layer resolution and robust against the unphysical shock anomaly phenomenon. Also, they conduct extensive numerical tests to show that the FVS type Riemann solvers are more robust against the shock anomaly than the others. However, they can not capture boundary layers accurately due to their intrinsic quality of Boltzmann. In contrast, both the FDS type Riemann solvers and the Flux Splitting type Riemann solvers are with high resolutions in capturing boundary layers. For these, most studies are conducted on the improvements of the FDS type Riemann solvers and the Flux Splitting type Riemann solvers for the shock anomaly. In Ref [39] , Peery initially indicated that the Roe scheme is inclined to show unphysical shock anomaly phenomenon. This phenomenon can be clearly seen in Fig. 2a . In order to overcome this defect, some entropy fix methods were proposed [66, 67] . However, these methods contain problemdependent coefficients. Improper coefficient's definition may deteriorate the resolution [68] . On the other hand, Quirk conjectured that any scheme that does not survive the odd-even decoupling is incapable to be free from the appearance of the shock anomaly [69] [70] [71] . Following Quirk's idea, Kim conducted analyses on the odd-even decoupling and argued that the feeding rate of the pressure perturbation and the damping rate of the density perturbation make the Roe scheme shock unstable simultaneously [72] . By controling these perturbations near shocks, he proposed the RoeM scheme and the RoeM2 scheme in two-dimensional cases as follows: RoeM where (49) As can be seen in Fig. 2 which shows the schemes' pressure contours and shock profiles along the stagnation streamline, both RoeM and RoeM2 improve the Roe scheme's robustness against the shock anomaly remarkably. For the HLLEM scheme presented in Sect. 3, Xie and his coworkers conducted numerical analyses based on the odd-even decoupling case and indicated that the HLLEM scheme can not ensure the consistency of mass flux across the normal shock, which makes it easy to produce unphysical shock anomalies [73, 74] . Similar studies can also be found in Refs [75, 76] . Based on such conclusion, Xie proposed the HLLEMS scheme by controling the dissipation of the shear waves around the shock waves with the following two-dimensional form Fig. 3 Density Contours of the Odd-even grid perturbation problem in the primary grid (Ms = 6.0) [77] Figures 3 and 4 compare the HLLEM scheme and the HLLEMS scheme's performances in capturing shocks. As can be clearly seen, HLLEMS scheme does improve the HLLEM scheme's robustness against the shock anomaly. Solvers for the Shock Anomaly By analyzing the hypersonic flow over the blunt body case, Kim found that the AUSM + scheme tends to show unphysical pressure overshoots near strong shocks and in the zones where the convective velocity is small (Fig. 5) . Also, he attributed these phenomena to the fact that the AUSM + scheme only considers the density at the left state when simulates the mass flux. In order to relieve this defect, he proposed the AUSMPW scheme by adding a weight to the AUSM + scheme to consider the states both at the left and the right [78] . In two-dimensional case, it is written as follows: where The unphysical pressure overshoots AUSM + encounters [78] (66) In Ref [79] , Kim analyzed the ASUMPW scheme and indicated that the functions f L,R in Eq. (65) are time consuming due to their complex forms. In order to improve this drawback, he proposed the AUSMPW + scheme by modifying these functions to The definitions of the symbols p 1,L , p 1,R , p 2,L , p 2,R are shown in Fig. 6 . Also, to dispel the appearance of the unphysical expansion shock wave, the numerical speed of sound is constructed as follows: where and the definitions of the V L and V R are shown in Fig. 7 . Figure 8 displays the schemes' pressure contours in hypersonic flows around a blunt wedge case (Mach number 15). As can be clearly seen, AUSMPW + avoids the unphysical pressure oscillations which are encountered by AUSM + . In Ref [28] , Roe conducted the hypersonic flow over a blunt body case to compare the Riemann solvers' performances in hypersonic heating predictions. The initial freestream conditions and the computational mesh are illustrated in the Sect. 4 of this study. Figures 9 and 10 compare the traditional Riemann solvers' profiles of heat transfer rates over the blunt body. As can be clearly seen, RoeM and AUSMPW + , which are highly robust against the shock anomaly, are shown to be much more stable and more accurate than the other methods on both grids. Inspired by the Liou-Steffen splitting procedure, scholars proposed two more procedures to split the Euler flux vector into a convective part and a pressure part. These two procedures are the Toro splitting procedure [80] [81] [82] and the Zha-Bilgen splitting procedure [83] [84] [85] . Borrowing Liou and Steffen's idea, Toro proposed a different idea to split the flux vector F into a convective part and a pressure part as follows [80] : where The Zha-Bilgen Splitting Procedure In Ref [83] , Zha and Bilgen proposed a distinct idea to split the flux vector F into a convective part and a pressure part in the following two-dimensional case: Schemes' profiles of heat transfer rates over the blunt body (coarse grid) [28] Schemes' profiles of heat transfer rates over the blunt body (dense grid) [28] where By conducting characteristic analyses of these splitting procedures, Qu found that the Zha-Bilgen splitting procedure is the best one in terms of consistency with theoretical physical properties. Based on such conclusion, Qu combined the Zha-Bilgen splitting procedure with the AUSMPW + scheme's idea and proposed the E-AUSMPW scheme [86] . In two-dimensional case, it is constructed as follows: It should be reminded that the definitions of the speed of sound is the same as the AUSMPW + scheme's. Figure 11 illustrates the schemes' pressure, temperature, and density distributions in the inverse receding flow case where the initial conditions are ( , p, u) L = (1.0, 1.0, −0.8c L ) , ( , p, u) R = (1.0, 1.0, 0.8c R ) . It is reminded that the results were obtained with the CFL number 0.5, the iteration count 50, the number of grid points 300, and the first order accuracy in both space and time. As can be seen, the E-AUSMPW is much more robust against the unphysical 'overheating' phenomenon than the other schemes by avoiding the appearance of the unphysical temperature overshoot. In general, it is challenging for existing schemes to produce satisfactory heating profiles on sphere [87] . Therefore, in order to show the E-AUSMPW scheme's performance in hypersonic flow's computations, Qu conducted hypersonic flow's computation over a blunt cone [88] . Figure 12 shows the model and the computational mesh. The model's length is 447 mm while the radius r = 27.94 mm. The free stream conditions are M ∞ = 10.6, T ∞ = 47.3 K, Re ∞ /m = 3.937 × 10 6 . Also, two different meshes are adopted: 77 (circumferential) × 50 (wall normal) × 90 (flow direction) called primary1 grid, 77 (circumferential) × 80 (wall normal) × 90 (flow direction) called primary2 grid. More details can be found in Ref [86] . The schemes' residual histories in Fig. 13 indicate that EAUSMPW performs much better in terms of the convergence property than the others. Figures 14 and 15 compare the schemes' heating profiles. As can be clearly seen, EAUSMPW is the only one that can give a smooth and symmetry heating profile. Also, its peak heating value agrees well with the stagnation point which is consistent with the theoretical analysis as well. Solver for Hypersonic Flows In order to study the traditional Riemann solvers' performance at low speeds, Qu conducted numerical cases in Ref [89] . Figure 16 shows the schemes' pressure contours of inviscid flows around the NACA0012 airfoil in case the Mach number is 0.05 and the 1st order reconstruction technology is adopted. All these traditional Riemann solvers are incapable of obtaining physical solutions in this lowspeed case. Figures 17 and 18 display the solvers' pressure contours in combination with the 2nd order reconstruction technology [90] and the 3rd order reconstruction technology [91] . As can be seen, higher order reconstruction method may generate non-physical solution. Therefore, it seems to be necessary to adopt an all-speed Riemann solver to simulate low-speed flows accurately. Nowadays, the existing RANS-LES hybrid methods, including DES [92] [93] [94] , SAS [95, 96] , and PANS [97, 98] , etc., are constructed by recovering the RANS methods in the boundary layer and using the LES elsewhere. Therefore, it is in a high demand to analyze the low-speed turbulent flow's accuracy using RANS methods, and Refs [99, 100] conducted studies on this issue with several test cases. For example, Fig. 19 compares different Riemann solvers' surface pressure distributions in the turbulent flow past a NACA0012 airfoil case (M ∞ = 0.15, Re ∞ = 6 × 10 6 ) [101] [102] [103] . It should be reminded that the AUSMPWM scheme and the F-Roe scheme are all-speed Riemann solvers built upon the AUSMPW + scheme and the Roe scheme, respectively. More details on these schemes will be introduced in Sect. 6. Similarly, Fig. 20 compares different Riemann solvers' surface pressure distributions in the turbulent flow past a NACA4412 airfoil case (M ∞ = 0.09, α = 13.87°, Fig. 11 Schemes' pressure, temperature, and density distributions at the 50th step in the mesh of 300 grid points (first order both in space and time, a AUSM+, b AUSMPW+, c E-AUSMPW, d Roe's FDS) [86] and Re ∞ = 1.52 × 10 6 ) [104] . Results in these figures indicate that a high level of accuracy at low speeds is of significance to the steady RANS simulations. As well known, both the force and the heating loads are relative to the pressure, velocity, and temperature distributions in the boundary layers. Therefore, theoretically these methods are incapable of obtaining physical solutions in studying the hypersonic airliners' aerodynamic properties, and the Riemann solver should be free from the unphysical solutions at low speeds if employed in the design of the hypersonic airliner. In order to verify this assumption, Qu adopted different Riemann solvers to conduct studies on this issue [26, 27] . For example, the three dimensional, viscous, and hypersonic flow over a blunt cone case was conducted [105, 106] . In this case, the model's length is 447 mm, the radius r is 27.94 mm, and the grid type is C type. The free stream conditions are M ∞ = 10.6, T ∞ = 47.3 K, and Re ∞ /m = 3.937 × 10 6 . Moreover, the following different cell Reynolds numbers based on the minimum cell size [107] [108] [109] were used with the number of the grids unchangeable: 2, 4, 6, 8, 10, 20, 40, 60, 80, and 100. Figure 21 compares the schemes' performances with different cell Reynolds numbers. It should be noticed that the Roe scheme and the F-Roe scheme are combined with the Muller's entropy fix to improve the robustness in high-speed flows. As can be seen, by improving the accuracy at low speeds, both F-Roe and AUSMPWM are more accurate than their original counterparts in hypersonic heating computations. Moreover, Both F-Roe and AUSMPWM show much less dependencies on the cell Reynolds number than their original versions due to their high levels of accuracy at low speeds. As introduced before, the widely used Riemann solvers in both scholar and engineering areas are built upon the Godunov's idea to solve the compressible Euler equations. They are incapable of obtaining physical solutions at subsonic speeds in theory [110, 111] . In order to overcome this defect, preconditioning matrix methods are widely used. However, these methods will encounter the global cut-off problem in hypersonic problems. In other words, if hypersonic flows past an airliner are simulated, these preconditioning matrix methods have little influence on the traditional Riemann solvers and the low-speed flows in the boundary layers can not be accurately captured as shown in Sect. 5. Therefore, these methods are incapable of obtaining physical solutions in studying the hypersonic airliners' aerodynamic properties, and the Riemann solver should be improved to be free from the unphysical solutions at low speeds if employed in the design of the hypersonic airliner. In this section, asymptotic analyses of the compressible Euler equations on low speeds' issues will be reviewed. Also, the constructions and the applications of the all-speed Riemann solvers will be presented. As introduced before, the mathematical properties of the Euler equations at different speeds are different. If the Mach number is low and the flow is nearly incompressible, the equations are with the elliptic property. If the Mach number is supersonic, the equations are with the hyperbolic property. Therefore, the traditional Riemann solvers built upon the compressible Euler equations should be incapable of obtaining physical solutions at subsonic speeds in theory. Through complex mathematical manipulation, Guillard conducted asymptotic analyses of the compressible Euler equations at low speeds and verified such conclusion [33, 34] . In this section, the singular limit of the compressible Euler equations with the Mach number approaching zero will be reviewed briefly. As the first step, the compressible Euler equations are non-dimensionalized with the variables * = max ( ) , [86] where x * indicates an arbitrary length scale. By combining the two-dimensional compressible Euler equations in Sect. 2 with these non-dimensional variables, the Euler equations can be non-dimensionalized as follows: Continuity Equation Equation where ̃ = (u, v) T is the velocity vector and M * = √ u * ⋅u * +v * ⋅v * c * is the reference Mach number. In the second step, the non-dimensional variables in Eq. (83) are expanded into powers of the reference Mach number aysmptotically such as In the third step, by combining the Eqs. (87) (88) (89) (90) (91) with the non-dimensionalized Eqs. (84) (85) (86) and collecting the terms with the same power of M * , the following equations can be obtained (the superscript ∧ is dropped for convennience): Therefore, with the Mach number approaching zero, the pressure should vary with the square of the Mach number in space as follows: Following Guillard's idea, Li conducted asymptotic analyses of the Roe scheme on low speeds' issues in Ref [112] and found that the original Roe scheme is incapable of obtaining physical solutions at low speeds. Similarly, Qu studied the HLLE type schemes' performance theoretically and indicated that both HLLEM and HLLEMS are incapable of obtaining physical solutions [77] . In this section, the constructions of the FDS type Riemann solvers for all speeds will be introduced briefly. The A-Roe Scheme In Ref [113] , Li unified the Low-speed Roe scheme and the shock-capturing Roe scheme with a function of the local Mach number to construct a new version Roe scheme. In the two-dimensional case, it can be written as follows: [99] where The T-Roe Scheme Different form Li's idea, Thornber indicated that the modification of the right eigenvector matrix can also improve the Roe scheme's performance at low speeds [114, 115] . In the two-dimensional case, it is with the following form: where where a' is defined in Eq. (99) . The F-Roe Scheme Asymptotic analyses of the Roe scheme indicated that the accuracy problem of the Roe scheme is attributed to the dissipation corresponding to the term Δu [113] . Based on such conclusion, Fillion proposed the F-Roe scheme to improve the dissipation of the term Δu with the following two-dimensional form [116, 117] : where The RoeMAS Scheme In Ref [118] , all improvements of the Roe scheme aforementioned are numerically studied. Results show that these improvements perform identically in low-speed problems. However, they are built upon the original Roe scheme and encounter the unphysical shock anomaly in hypersonic cases [119] . By modifying the maximum/minimum wave speeds and improving the density perturbation's damping rates based on KIM's analyses, Qu proposed a new Roe-type scheme called RoeMAS for all speeds [120] . In two-dimensional case, the RoeMAS scheme can be written as follows: Stagnation heating values obtained by different schemes within meshes of different cell Reynolds numbers [26] where = (0, 1, 0 ,ũ) T all speeds, several numerical cases were conducted. The first case is the two-dimensional inviscid flow over the NACA0012 airfoil. The Mach numbers are 0.001, 0.01, and 0.1. Figure 22 shows different schemes' pressure fluxtuations Ind(p) = (P max -P min )/P max VS the Mach number on different grids. As can be seen, only RoeMAS and F-Roe are capable of agreeing with the theoretical predictions in these low-speed cases. Also, the entropy fix, which is widely adopted to improve the Roe type schemes' robustness at high [120] speeds, is shown to deteriorate the F-Roe's high resolution at low speeds. The second case is the odd even decoupling case. The initial conditions are given as ( , u, v Figure 23 shows the schemes' density contours on grids of 20 (Y) × 200 (X). As can be clearly seen, RoeMAS is more robust in terms of shock anomaly at high speeds. As analyzed in Sec 4, Xie proposed the HLLEMS scheme to improve robustness against the shock anomaly. In terms of the low-speed issues, Qu carried out asymptotic analysis of the HLLEMS scheme's behavior at low speeds in Ref [77] . Results indicate that both HLLEM and HLLEMS are incapable of obtaining physical solutions at low speeds [73, 121] . Based on such analyses, the HLLEMS-AS scheme for all speeds is proposed to control the numerical dissipation in the momentum equations at low speeds. In two-dimensional case, it can be written as follows: where The definitions of the other symbols can be found in Sects. 3.1.2 and Sect. 4.1.2. Figure 24 compares the schemes' pressure fluctuations Ind(p) = (P max -P min )/P max VS the inlet Mach number on different grids in the low Mach number inviscid flow around the NACA0012 airfoil case. As can be seen, HLLEMS-AS improves the HLLEMS scheme's accuracy at low speeds remarkably. On the other hand, the density contours of the Odd-even grid perturbation problem in Fig. 25 indicate that HLLEMS-AS is as robust against the shock anomaly as HLLEMS at high speeds. Flux Splitting Riemann Solvers In order to achieve high levels of accuracy at all speeds, Liou reviewed the AUSM + scheme and proposed the following improvements in his studies [122, 123] : Scale the speed of sound to control the numerical dissipation at low speeds. Introduce the concept of pressure dissipation to the convective flux. Add the dissipation of velocity to the pressure flux. Numerical tests indicate that the AUSM + UP scheme with these improvements of the AUSM + scheme performs well at both low speeds and high speeds. In general, it is written as follows: where K u = 0.75, K p = 0.25, = 1.0 , and M ∞ indicates the freestream Mach number. Different from Guillard's analysis method, Eiji Shima and his co-workers conducted dimensional analyses on the dissipation terms of the AUSM + scheme. They found that the numerical dissipations of the pressure flux in Eq. (37) is too large in low speeds' simulations. Based on such conclusion, (120) they proposed the SLAU scheme by controlling the dissipation in the pressure flux [124] . In two-dimensional cases, they can be written as follows: Figure 26 shows the SLAU schemes' pressure contours of inviscid flows around NACA0012 in cases the Mach number is 0.05 and three different reconstruction technologies are adopted. Compared these results with Figs. 16, 17, and 18, it can be clearly seen that SLAU is much more accurate than traditional Riemann solvers. Also, the all-speed Riemann (137) solver is of significance in simulating low-speed flows in comparison with the reconstruction technologies. Inspired by the construction of the SLAU scheme, Qu adopted the controlling function in Eq. (137) to control the AUSMPW + scheme's numerical dissipation in the pressure flux at low speeds as in Ref [31] . Extensive numerical test cases demonstrate that AUSMPWM performs better than SLAU at both low speeds and high speeds. However, by theoretical asymptotic analyses, Qu found that AUSMPWM is still incapable of obtaining physical solutions in the low speed limit. To overcome this defect, Qu conducted asymptotic analyses of the AUSMPWM scheme and proposed the AUSMPWM2 scheme [125] . In mathematical form, the mass flux of the AUSMPWM2 scheme is the same as that of the AUSMPW + scheme, and the pressure flux is modified into the following form: where Figure 27 displays the schemes' pressure fluctuations Ind(p) = (P max − P min )/P max VS the Mach number in case they are combined with different reconstruction methods to simulate the 2d inviscid NACA0012 airfoil case. As analyzed in Sect. 6.1, the pressure fluctuations should (140) f 2 = min (|M|, 1) , a 1st order, b 2nd order, c 3rd order) [89] scale with the square of the Mach number in this continuous case. Therefore, the results in this figure indicates that AUSMPWM2 is the only one that can obtain physical solutions at low speeds. In addition, it should be reminded that the higher order reconstruction methods provide limited help to the Riemann solvers' performances at low speeds. In order to show the AUSMPWM2 scheme's accuracy in predicting hypersonic heating loads' distributions, Ref [125] presents its application to the hypersonic viscous flow over the X33 aerocraft case. The X33 aircraft is a reusable manned space verification vehicle proposed by USA, and the NASA Langley Research Centre has conducted numerous wind tunnel tests to study the aerodynamic property of this vehicle [126, 127] . The three-view drawing and the multiblock structured mesh adopted for the X33 aerocraft are shown in Fig. 28 . The free stream conditions are M ∞ = 5.99, T ∞ = 47.3 K, = 0.0628 kg/m 3 , and the angle of attack is 30. More details on this case can be found in Ref [125] . Figure 29 shows the Mach number contours in the symmetry plane and pressure distributions at the wall obtained by the AUSMPWM2 scheme. The bow shock, the expansion wave, the separation zone at the tail and the shear layer , a 1st order reconstruction, b 2nd order reconstruction, c 3rd order reconstruction) [125] are captured accurately. Table 1 lists the schemes' computational errors at different locations along the meridian. Obviously, the AUSMPWM2 scheme is capable of obtaining much better heat distributions than SLAU and Roe. As can be seen in Sect. 6.3.2, the SLAU scheme adopts the Roe scheme's mass flux, which is not robust against the shock anomaly as analyzed in Sect. 4. Also, its pressure flux still can not obtain physical low-speed results as shown in Sec 6.3.3. Inspired by the constructions of the RoeMAS scheme and the AUSMPWM2 scheme, Qu proposed a shock stable Riemann solver for all speeds to accurately predict the aerodynamic force and heating loads of hypersonic airliners [128] . This solver, which is called SLAU-M1 and built upon the SLAU scheme, constructs the mass flux as RoeMAS and calculates the pressure flux as AUSMPWM2. In twodimensional mathematical form, it can be written as follows: where the mass flux ṁ SLAU -M1 is the same as the RoeMAS scheme's as in Sect. 6.2, and the pressure flux F SLAU -M1 p,i+1∕2 is More details can be found in Ref [128] . In order to show the SLAU-M1 scheme's accuracy in simulating hypersonic flows, Qu conducted two cases. The first case is the hypersonic viscous flow over a 2D double ellipsoid [129] , and the second is the hypersonic viscous flow over the hypersonic flight experiment (Hyflex) aerocraft [130] [131] [132] . In terms of the first case, the freestream conditions are M ∞ = 8.15, T ∞ = 56 K, Re ∞ /m = 7.68 × 10 6 , and P ∞ = 370.7 Pa. Figure 30 compares the schemes' Stanton number distributions along the meridians of the windward and leeward sides with the experimental data. The all-speed SLAU-M1 scheme and SLAU scheme agree with the experimental data much better than the other traditional Riemann solvers. Especially, SLAU-M1 is more accurate than SLAU. As the second case, the Hyflex aircraft, which is a reusable manned space verification vehicle proposed by the Japanese Aerospace Laboratory and National Space Development Agency of Japan, is investigated. The free stream conditions are M∞ = 12.23, T∞ = 129.41 K, = 0.0033678 kg/m 3 , and the angle of attack is 49 0 . (144) P = (0, p, 0, 0) T Table 2 illustrates the heating loads' distributions obtained along the meridian at the windward by different schemes. It also shows that SLAU-M1 is more accurate than the others (Fig. 31 ). Accurate prediction of aerodynamic force and heating flux is critical for the design of the new generation hypersonic airliner, and the selection of the Riemann solver is with a decisive influence on the accuracy. In this survey, the research progresses on the Riemann solvers for hypersonic flows, including the Godunov method for the compressible Euler equations, the traditional Riemann solvers for compressible flows, the improvements of traditional Riemann solvers for the shock anomaly, and the all-speed Riemann solvers, have been reviewed in details. The results described in these sections lead to the following concluding remarks: The constructions and the applications of all-speed Riemann solvers for hypersonic flows have been drawing increasing attentions. The existing all-speed solvers are all built to improve the traditional solvers' robustness in terms of the shock anomaly and accuracy at low speeds. Extensive low-speed numerical cases and theoretical analyses indicate that traditional Riemann solvers initially built for compressible flows are incapable of obtaining physical solutions in the low Mach number limit. In addition, the high order reconstruction technologies are shown to be of Systematic numerical cases, including the low-speed NACA0012 airfoil case, the odd-even decoupling case, the hypersonic flow over the X33 aerocraft case, the hypersonic flow over the 2D double ellipsoid case, and the hypersonic flow over the Hyflex aerocraft case, are illustrated. Results indicate that the all-speed Riemann solvers do improve the traditional solvers' accuracy in predicting the aerodynamic loads of the vehicles at hypersonic speeds. Till now, the study of the all-speed Riemann solvers are still based on the direction splitting procedure, which is widely used in most commercial CFD softwares. More work on the construction and application of the genuinely multidimensional all-speed Riemann solver should be conducted to further improve the accuracy in predicting hypersonic complex flows. Establishment methodology of comfort parameters series for civil aircraft cabin Aerodynamic optimization of civil aircraft with wing-mounted engine jet based on adjoint method Concept definition of new-generation multi-purpose manned spacecraft Hypersonic aerodynamic interference investigation for a two-stage-to-orbit model Multi-objective design optimization of blunt body with spike and aerodisk in hypersonic flow Fluid-thermal analysis of aerodynamic heating over spiked blunt body configurations Mahulikarb, aero-thermal analysis of lifting body configurations in hypersonic flow Modeling of aerodynamic heat flux and thermo elastic behavior of nose caps of hypersonic vehicles Drag and heat flux reduction mechanism induced by the spike and its combinations in supersonic flows: a review HIFiRE-5 flight vehicle design Advances in computational heat transfer Supersonic mixing in airbreathing propulsion systems for hypersonic flights Experimental investigation on aero-heating of rudder shaft within laminar/turbulent hypersonic boundary layers 31 Experiment for the HYFLEX vehicle in the ONERAF4 wind tunnel and the Mach number contours in the symmetry plane obtained by SLAU-M1 An effective flux scheme for hypersonic heating prediction of re-entry vehicles A new flux splitting scheme for the Euler equations II: E-AUSMPWAS for all speeds Numerical investigation of the supersonic stabilizing parachute's heating loads Aerothermoelastic-acoustics simulation of flight vehicles The hyper -X launch vehicle: challenges and design considerations for hypersonic flight testing The X-51A Scramjet engine flight demonstration program Performance of a generic X-51 waverider-thrust and drag computed using the MASIV reduced order model. AIAA propulsion and energy forum and exposition HIFiRE-5 boundarylayer transition measured in a Mach-6 quiet tunnel with infrared thermography Quiet tunnel measurements of HIFiRE-5 boundary-layer transition HIFiRE-1 preliminary aerothermodynamic experiments Aerodynamic characteristics evaluation of hypersonic flight experiment vehicle based on flight data The development of the X-37 Re-entry vehicle Investigation into the influences of the low speed's accuracy on the hypersonic heating computations A study of upwind schemes on the laminar hypersonic heating predictions for the reusable space vehicle Evaluation of Euler fluxes for hypersonic heating computations Riemann solvers and numerical methods for fluid dynamics Partial differential equation A Parameter-free scheme for all speeds' simulations Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number On the behaviour of upwind schemes in the low mach number limit On the behaviour of upwind schemes in the low mach number limit: II. Godunov type schemes Preconditioning applied to variable and constant density flows Preconditioning technique in computational fluid dynamics Viscous airfoil computations using local preconditioning Preconditioning method and engineering application of large Eddy simulation Blunt-body flow simulations Numerical instabilities inupwind methods: analysis and cures for the "Carbuncle" phenomenon Mechanismderived shock instability elimination for Riemann-solver-based shock-capturing scheme Numerical experiments on anomalies from stationary, slowly moving, and fast-moving shocks A cure for numerical shock instability in HLLC Riemann solver using antidiffusion control Fundamentals of aerodynamics Computational fluid dynamics: principles and applications, 1st edn A finite difference method for the numerical calculation of discontinuous solutions of hydrodynamic equations Study on the dissipative effect of approximate Riemann solver on hypersonic heat flux simulation Approximate Riemann solvers, parameter vectors and difference schemes A two-dimensional HLLC Riemann solver for conservation laws: application to Euler and magnetohydrodynamic flows A simple two-dimensional extension of the HLL Riema2nn solver for hyperbolic systems of conservation laws Robust HLL-type Riemann solver capable of resolving contact discontinuity A new formulation for twowave Riemann solver accurate at contact interfaces A HLLCtype Riemann solver and high-resolution Godunov method for a two-phase model of reactive flow with general equations of state On Godunov-type methods near low densities A low diffusion flux splitting method for inviscid compressible flows A new efficient formulation of the HLLEM Riemann solver for general conservative and nonconservative hyperbolic systems An implicit flux-vector splitting scheme for the computation of viscous hypersonic flow Flux vector splitting of the inviscid gas-dynamics equations with application to finite difference methods Flux vector splitting for the Euler equations, eighth international conference of numerical methods in fluid dynamics An accurate and robust flux splitting scheme for shock and contact discontinuities Ten years in the making-AUSM family Reflections on the early development of the "AUSM family" of Riemann solvers A new flux splitting scheme A sequal to AUSM: AUSM+ Progress towards an improved CFD method: AUSM+. AIAA Paper Modified entropy correction formula for the roe scheme Simple improvements of an upwind TVD scheme for hypersonic flow Artificial viscosity to cure the shock instability in high-order Godunov-type schemes Carbuncle phenomena and other shock anomalies in three dimensions A contribution to the great Riemann solver debate Mass flux schemes and connection to shock instability Cures for the shock instability, development of a shock-stable roe scheme An accurate and robust HLLC-type Riemann solver for the compressible Euler system at various Mach numbers On numerical instabilities of Godunov-type schemes for strong shocks Strategies to cure numerical shock instability in the HLLEM Riemann solver A simple cure for numerical shock instability in the HLLC Riemann solver A new all-speed flux scheme for the Euler equations An improvement of AUSM schemes by introducing the pressure-based weight functions Methods for the accurate computations of hypersonic flows I: AUSMPW+ scheme Flux splitting schemes for the Euler equations A flux splitting method for the Euler equations A robust flux splitting method with low dissipation for all-speed flows Numerical solutions of Euler equations by using a new flux vector splitting scheme A low-diffusion E-CUSP upwind scheme for transonic flows Low-diffusion efficient upwind scheme A new flux splitting scheme for the Euler equations Multidimensional, inviscid flux reconstruction for simulation of hypersonic heating on tetrahedral grids AIAA Paper Aero-thermal analysis of lifting body configurations in hypersonic flow A study of parameterfree shock capturing upwind schemes on low speeds' issues Towards the ultimate conservation difference scheme V: a second-order sequal to Godunov's method Accurate, efficient and monotonic numerical methods for multidimensional compressible flows. Part II: multi-dimensional limiting process Hybrid LES/RANS methods for the simulation of turbulent flows Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach A new version of detached-eddy simulation, resistant to ambiguous grid densties Development and application of SST-SAS turbulence model in the DESIDER project The scale-adaptive simulation method for unsteady turbulent flow predictions, Part 1: theory and model description Partially-averaged Navier-Stokes model for turbulence: a Reynolds-averaged Navier-Stokes to direct numerical simulation bridging method Partially-averaged Navier-Skokes method for turbulent flows: k-ω model implementation Investigation into the influences of the lowspeed flows' accuracy on RANS simulations Investigation of all-speed schemes for turbulent simulations with low-speed features A critical assessment of wind tunnel results for the NACA 0012 Airfoil Langley Research Center: turbulence modelling resource Effects of independent variation of mach and Reynolds numbers on the low-speed aerodynamic characteristics of the NACA 0012 airfoil section Flying-hot-wire study of flow past an NACA 4412 airfoil at maximum lift Improved k-w-r transition model by introducing the local effects of nose bluntness for hypersonic heat transfer Numerical investigation of blunt body's heating load reduction with combination of spike and opposing jet A grid strategy for predicting the space plane's hypersonic aerodynamic heating loads Grid criteria for numerical simulation of hypersonic aerothermodynamics in transitional regime A wall grid scale criterion for hypersonic aerodynamic heating calculation All-speed Roe scheme for the large eddy simulation of homogeneous decaying turbulence A low-Mach number method for the numerical simulation of complex flows The momentum interpolation method based on the time-marching algotithm for all-speed flows An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behavior An improved reconstruction method for compressible flows with low Mach number features Numerical dissipation of upwind schemes in low Mach flow A low-Mach number fix for Roe's approximate Riemann solver A new platform for core thermal-hydraulic studies Mechanism of Roe-type schemes for allspeed flows and its application Cures for expansion shock and shock instability of Roe scheme based on momentum interpolation mechanism A new Roe-type scheme for all speeds AUSM-like expression of HLLC and its all-speed extension An enhanced AUSM + -up scheme for high-speed compressible two-phase flows on hybrid grids A sequal to AUSM, Part II: AUSM+-up for all speeds Parameter-Free simple low-dissipation AUSM-family scheme for all speeds An improvement on the AUSMPWM scheme for hypersonic heating predictions X-33 Reusable launch vehicle structural technologies Computational/experimental aeroheating predictions for X-33 phase II vehicle Shock-stable flux scheme for predicting the hypersonic airliner's aerodynamic heating loads Attempt to evaluate the computations for test case 6.1: cold hypersonic flow past ellipsoidal shapes, hypersonic flows for reentry problems Aerodynamic characteristics of lifting body HYFLEX under high-temperature real gas condition. 55th AIAA aerospace sciences meeting Overview of the HYFLEX Project CFD analysis of aerodynamic heating for HYFLEX high enthalpy flow tests Acknowledgements This study was co-supported by the National Natural Science Foundation of China (No. 11972308 and No. 11902265) and the National Major Basic Research Project 1912. The first author would like to give his greatest appreciate to his father, who passed away on January 4th, 2021. In fact, the first draft was accomplished by the first author during the Spring Festival in 2020, in which most Chinese people kept themselves home due to the appearance of the COVID-19, and the first author and his father had each other. His father's attendance helped the first author a lot in finishing this paper. How the first author wishes he could see his father again! Wish him happiness in heaven! The authors declare these is no conflict of interest regarding the publication of this paper.