An Effective Simulation Scheme for Predicting the Aerodynamic Heat of a Scramjet-Propelled Vehicle

: Currently, aerothermal research into scramjet-propelled vehicles characterized by a wedge-shaped section is relatively sparse. Based on the Mach number, grid strategy, and numerical method, an effective simulation scheme for predicting the aerodynamic heat of a scramjet-propelled vehicle during ﬂight is proposed in this paper. At different Mach numbers, the appropriate grid strategy and numerical method were determined by validation tests. Two-dimensional external ﬂow ﬁeld models based on wedge sections were established and, unlike in blunt bodies, the tests showed that at the high supersonic stage, the ideal cell Reynolds number should be no larger than 16. At the hypersonic stage, the ideal cell Reynolds number and aspect ratio of wall cells near the shock should be no larger than 40, and the AUSM+ ﬂux type performs better than Roe’s FDS ﬂux type at the above stages. The aerothermal prediction indicates that during a ﬂight time of about 34 s, the temperature change reaches about 1913.35 ◦ C, and the maximum average temperature change rate reaches 115 ◦ C/s.


Introduction
Aerodynamic heat prediction plays an important role in the design and optimization of hypersonic vehicles [1,2]. In general, there are three technical approaches to aerodynamic heat prediction: the engineering estimation [3], the wind tunnel experiment [4], and computational fluid dynamics (CFD) [5,6]. With the development of computer technology, CFD has made great progress; its results offer higher accuracy than engineering estimation, and, unlike the wind tunnel experiment, its simulation ability can cover the entire flight trajectory. The CFD technology has become an indispensable means in the study of aerothermal predictions. Its results are influenced by many factors, such as Mach number, gird strategy, and numerical methods, which include turbulence models, upwind scheme, limiter, and so on. Domestic and foreign scholars have carried out research into these factors. The aerodynamic characteristics of different Mach numbers are different, and the corresponding grid strategy and numerical method are also different [7]. In terms of grid strategy, the height of the first layer of wall cells has a great influence on the numerical calculation results, and the mesh needs to be adjusted several times according to the CFD results, which is a time-consuming part of the CFD process. Ref. [8] discusses the effect of the grid on aerothermal prediction by taking an external flow field simulation model based on a blunt body as an example. The results show that the deviation of the results could reach more than 20% when the Reynolds number of the wall grid is doubled while other conditions are unchanged. Ref. [9] considers the effect of cell Reynolds numbers on aerothermal prediction by applying a cylinder as an example. The results show that it is necessary to keep the cell Reynolds number less than 10 in the aerothermal calculation, and Refs. [10,11] presents similar research cases. Ref. [12] shows that the property of shock stability influences the accuracy of heating prediction. Ref. [13] reveals that the aspect ratio Pa, temperature is 222.54 K, air density is 0.0339 kg/m 3 and dynamic viscosity is 1.454e-05pa.s. Cruising phase: the flight altitude is about 30 km, Mach number is 6.5, pressure is 1197 Pa, temperature is 226.51 K, air density is 0.0181 kg/m 3 and dynamic viscosity is 1.475e-05pa.s.
According to the different Mach numbers, the vehicle is in different external flow field environments [7,19]. At the transonic speed stage (when the Mach number is between 0.6 and 1), compressibility effects such as flow choking becoming important. At the supersonic speed stage (Mach number between 1 and 3), depending on the specific shape of the aircraft, a shock wave may be produced in the airflow. At the high supersonic (when the Mach number is between 3 and 5) and hypersonic speed stage (when the Mach number is bigger than 5), aerodynamic heating becomes very important. The flight path of the RBCC scramjet-propelled vehicle spans all of the above stages; therefore, a specific piecewise analysis is needed. According to the Mach number of the flight path, the simulation scheme is divided into four stages. To ensure the reliability of the prediction results, the simulation parameters should be confirmed by validation tests at each stage.
For flows involving compressibility and heat transfer, three conservation equations for mass, momentum, and energy need to be solved. When the flow is turbulent, additional transport equations are also solved, adopting the one-equation RANS Spalart-Allmaras model, which is the most widely applied [20], and the two-equation RANS SST k-omega turbulence model [21]. The control-volume-based technique is used to discretize the above equations. Firstly, the external flow field model is divided into discrete control volumes using a computational grid. Secondly, the integration of equations on the individual control volumes to algebraic equations for the discrete dependent variables and conserved scalars is accomplished using an upwind scheme. Next, the above-discretized equations are linearized, and the resultant linear equation system is solved to yield variables. Two numerical solvers, a pressure-based solver and a density-based solver, were applied in this study. More details about these equations and the finite volume method could be found in Ref. [22]. According to the different Mach numbers, the vehicle is in different external flow field environments [7,19]. At the transonic speed stage (when the Mach number is between 0.6 and 1), compressibility effects such as flow choking becoming important. At the supersonic speed stage (Mach number between 1 and 3), depending on the specific shape of the aircraft, a shock wave may be produced in the airflow. At the high supersonic (when the Mach number is between 3 and 5) and hypersonic speed stage (when the Mach number is bigger than 5), aerodynamic heating becomes very important. The flight path of the RBCC scramjet-propelled vehicle spans all of the above stages; therefore, a specific piecewise analysis is needed. According to the Mach number of the flight path, the simulation scheme is divided into four stages. To ensure the reliability of the prediction results, the simulation parameters should be confirmed by validation tests at each stage.
For flows involving compressibility and heat transfer, three conservation equations for mass, momentum, and energy need to be solved. When the flow is turbulent, additional transport equations are also solved, adopting the one-equation RANS Spalart-Allmaras model, which is the most widely applied [20], and the two-equation RANS SST k-omega turbulence model [21]. The control-volume-based technique is used to discretize the above equations. Firstly, the external flow field model is divided into discrete control volumes using a computational grid. Secondly, the integration of equations on the individual control volumes to algebraic equations for the discrete dependent variables and conserved scalars is accomplished using an upwind scheme. Next, the above-discretized equations are linearized, and the resultant linear equation system is solved to yield variables. Two numerical solvers, a pressure-based solver and a density-based solver, were applied in this study. More details about these equations and the finite volume method could be found in Ref. [22].

Grid Strategy
The experimental data in Ref. [23] are used to validate the simulation parameters at this stage, and the simulation external flow field model is as shown in Figure 2. The bold Appl. Sci. 2021, 11, 9344 4 of 26 black dot is the (0,0) point of the coordinate system. Line f is the wedge section, line g is the flat plate.

Grid Strategy
The experimental data in Ref. [23] are used to validate the simulation parameters at this stage, and the simulation external flow field model is as shown in Figure 2. The bold black dot is the (0,0) point of the coordinate system. Line f is the wedge section, line g is the flat plate.
where L is the characteristic length of the wedge, μ is the dynamic viscosity of the air, ρ is the air density, and U is the velocity of freestream. Considering that the simulation model is wedge-shaped, an empirical correlation formulation to estimate the skin friction coefficient ( f C ) is adopted, which is for the fully developed turbulent flow over a flat plate: Having computed the skin friction coefficient, the wall shear stress ( W τ ) is calcu- After that, we can calculate the friction velocity ( τ u ) from the wall shear stress ( W τ ): The equation for + y is: Turbulent flows are affected by the presence of walls; therefore, the accurate representation of the flow in the near-wall region determines the successful prediction of wall-bounded turbulent flows, and it is necessary to ensure the first layer cells of the boundary layer are in the sub-viscosity layer. The calculation of the height of the first layer cells of the boundary layer begins with the flow Reynolds number (R e ): where L is the characteristic length of the wedge, µ is the dynamic viscosity of the air, ρ is the air density, and U is the velocity of freestream. Considering that the simulation model is wedge-shaped, an empirical correlation formulation to estimate the skin friction coefficient (C f ) is adopted, which is for the fully developed turbulent flow over a flat plate: Having computed the skin friction coefficient, the wall shear stress (τ W ) is calculated: After that, we can calculate the friction velocity (u τ ) from the wall shear stress (τ W ): The equation for y + is: Therefore, the height of the wall-adjacent cell centroid from the wall (y P ) is: Finally, the first height of the boundary layer cell (y H ) is double y P : Appl. Sci. 2021, 11, 9344

of 26
The result of the first height of the boundary layer cell is an estimate and needs to be updated using information from the CFD analysis. In addition to calculating the y H , the growth ratio (r) normal to the wall and inflation layers (N) may also be required by mesh generation, which allows the layers to inflate normally to the wall. For aerodynamic flows using RANS turbulence modeling, CFD often aims for N = 15-30 through the thickness of the boundary layer. By adopting this criterion, we can calculate r normally to the wall, which is started by applying the formulation to estimate the boundary layer thickness, where R e is bigger than 5 × 10 5 : For a given number of r and N, we can calculate the total thickness of the inflation layers. Next, we match this to the boundary layer thickness δ to allow the inflation layers to cover the boundary layer thickness. Assume that the inflation layer has two layers, then the total thickness (y T ) is: Similarly, the total thickness of four inflation layers is: Therefore, the total thickness of N inflation layers is: Equation (11) is the formula of a geometric series, which could be rewritten using the identity for convenience: Hence, the total thickness of the inflation layer is: Now, the error between the total thickness of the inflation layer y T and the thickness of the boundary layer δ is: Equation (14) is a function of the growth ratio r and we need to find a value of r to make the error equal to 0, which is the total height of the inflation layer equal to the boundary layer thickness. The Newton-Raphson root-finding algorithm is adopted, and once the algorithm has converged, the solution r is the maximum growth ratio. Because, in mesh generators, the r is always greater than 1 and less than 2, the initial guess of 1 and 2 could be used to bound the root-finding algorithm to accelerate the convergence.
Three different Mach numbers, 0.768, 0.817. and 0.86 were applied, respectively, in the validation tests at the subsonic stage. The static temperature of freestream T was 300 K and the dynamic viscosity µ was 1.7894 × 10 −5 Pa·s; the density of air was ideal-gas, the specific heat at a constant air pressure was 1006.43 j/kg-k, and the characteristic length c was 0.59696 m. Initially, y + was 5 to ensure the first layer cells of the boundary layer were in the sub-viscosity layer, and y H was 1.53 × 10 −5 , according to Equations (1)- (7). Next, the CFD analysis showed that the maximum value of y + at the WALL surface reached 6.5. The simulation tests showed that the y + value of 1 could guarantee that the first layer cells of the boundary layer are in the sub-viscosity layer during the simulation process. The turbulent viscosity ratio (TVR) was adopted to evaluate whether the boundary layer profile was captured well. Figure 3 shows the corresponding plots of the TVR for the mesh Appl. Sci. 2021, 11, 9344 6 of 26 with y + = 1 and N = 20, since the y + value is 1, which means presence within the laminar sub-layer, and absence of turbulent viscosity. As the mesh moved through the buffer region and into the logarithmic region, the TVR value increased, and the maximum value occurred near the middle of the boundary layer. Subsequently, the turbulence gradually dissipated as the mesh approached the free stream condition, which was expected. Therefore, the values of y + and N were sufficient, so we could calculate the maximum growth ratio r by Equation (14). The first layer height y H and the maximum growth ratio were (1.49 × 10 −6 , 1.42), (1.38 × 10 −6 , 1.42) and (8.42 × 10 −6 , 1.34) respectively.
Next, the CFD analysis showed that the maximum value of + y at the WALL surface reached 6.5. The simulation tests showed that the + y value of 1 could guarantee that the first layer cells of the boundary layer are in the sub-viscosity layer during the simulation process. The turbulent viscosity ratio (TVR) was adopted to evaluate whether the boundary layer profile was captured well. Figure 3 shows the corresponding plots of the TVR for the mesh with + y = 1 and N = 20, since the + y value is 1, which means presence within the laminar sub-layer, and absence of turbulent viscosity. As the mesh moved through the buffer region and into the logarithmic region, the TVR value increased, and the maximum value occurred near the middle of the boundary layer. Subsequently, the turbulence gradually dissipated as the mesh approached the free stream condition, which was expected. Therefore, the values of + y and N were sufficient, so we could calculate the maximum growth ratio r by Equation (14).  To determine the best domain size, computational region analysis was performed. For this purpose, three sizes were chosen, as shown in Figure 4a [24]. The INLET boundary contained line a, the OUTPUT boundary contained line d, the FARDIELD boundary contained lines h, i, and j, the SYMMETRY boundary contained line e, and the WALL boundary contained lines f and g. Figure 4b offers a close view of the grids near the WALL region. Taking the numerical simulation adopting freestream Mach number of 0.768 and the Spalart-Allmaras model as an example, the values of the local Mach number m, the pressure coefficient (cp), and the temperature are as shown in Table 1, where x is the Xaxis location, and the computational domain with 8c was chosen for low errors. Further, to be sure of grid independency, three levels of different mesh were applied. As can be seen in Table 2, it is clear that the values of m at 480,000 were better; therefore, the mesh with 480,000 cells was selected. To determine the best domain size, computational region analysis was performed. For this purpose, three sizes were chosen, as shown in Figure 4a [24]. The INLET boundary contained line a, the OUTPUT boundary contained line d, the FARDIELD boundary contained lines h, i, and j, the SYMMETRY boundary contained line e, and the WALL boundary contained lines f and g. Figure 4b offers a close view of the grids near the WALL region. Taking the numerical simulation adopting freestream Mach number of 0.768 and the Spalart-Allmaras model as an example, the values of the local Mach number m, the pressure coefficient (c p ), and the temperature are as shown in Table 1, where x is the X-axis location, and the computational domain with 8c was chosen for low errors. Further, to be sure of grid independency, three levels of different mesh were applied. As can be seen in Table 2, it is clear that the values of m at 480,000 were better; therefore, the mesh with 480,000 cells was selected.   At the transonic stage, the pressure-based solver was used with a coupled algorithm [25]. The second-order upwind interpolation scheme was selected for the density [26]. The gradients were computed according to the least-squares cell-based method. For preventing spurious oscillations, a gradient limiter was adopted, which used the Minmod function [17] to limit and clip the reconstructed solution overshoots and undershoots on the cell faces. The second-order scheme was used for interpolating the pressure values at the faces, which reconstructed the face pressure using a central differencing scheme. The validation parameters of the convergence numerical results were as shown in Table 3.  Figure 5 demonstrates the distributions of m under three different Mach numbers of freestream M, where the X-axis was the ratio between x and c and the Y-axis was the value of m. EXP represents the experimental data, SST and SA represent the simulation results based on SST k-omega and Spalart-Allmaras, respectively. The overall tracking effect of the simulation data on the experiment data was good. At M = 0.768, the average error rate of SST k-omega was 9.69%, and that of Spalart-Allmaras was 8.34%. The variance of SST k-omega error was 0.000368, and that of Spalart-Allmaras is 0.000122. At M = 0.817, the average error rate of SST k-omega was 9.52%, and that of Spalart-Allmaras was 8.18%. The variance of SST k-omega error was 0.00035, and that of Spalart-Allmaras was 0.000152. At M = 0.86, the average error rate of SST k-omega was 10.14%, and that of Spalart-Allmaras was 8.52%. The variance of SST k-omega error was 0.000359, and that of Spalart-Allmaras was 0.000284. The total mean error rates were 9.78% for SST and 8.35% for Spalart-Allmaras. In general, the trend of the simulation results was in good agreement with the experimental data, and of the two models, the Spalart-Allmaras simulation results were better than those of SST k-omega at the transonic stage. In conclusion, the mesh with the value of + y as 1, = the turbulence model applying Spalart-Allmaras and the pressure-based solver, adopting a coupled algorithm, is preferred at the transonic stage.  In conclusion, the mesh with the value of y + as 1, = the turbulence model applying Spalart-Allmaras and the pressure-based solver, adopting a coupled algorithm, is preferred at the transonic stage.

Grid Strategy
The experimental data in Ref. [27] were used to validate the simulation parameters at this stage. Given the existence of the shock wave in the airflow, the value of the shock wave angle was estimated firstly. The Mach number M was 2.85 and the wedge angle α is 30 • , which meets the Equation (15), where γ = 1.4: The shock angle s was 54 • , which was calculated by Equation (16): Next, we established and divided the external flow field model for the simulation, as shown in Figure 6. The bold black dot is the (0,0) point of the coordinate system. Lines o and q are the flat plate and line p is the wedge section. The bold black square point is the origin of the wedge section. In conclusion, the mesh with the value of + y as 1, = the turbulence model applying Spalart-Allmaras and the pressure-based solver, adopting a coupled algorithm, is preferred at the transonic stage.

Grid Strategy
The experimental data in Ref. [27] were used to validate the simulation parameters at this stage. Given the existence of the shock wave in the airflow, the value of the shock wave angle was estimated firstly. The Mach number M was 2.85 and the wedge angle α is 30°, which meets the Equation (15), where γ = 1.4: The shock angle s was 54°, which was calculated by Equation (16): Next, we established and divided the external flow field model for the simulation, as shown in Figure 6. The bold black dot is the (0,0) point of the coordinate system. Lines o and q are the flat plate and line p is the wedge section. The bold black square point is the origin of the wedge section.  At the supersonic stage, the total temperature of the freestream T t was 270 K and the total pressure of the freestream p t was 1.7 atm. The static temperature of the freestream T was calculated by Equation (17) and T was 102.88 K: The static pressure of the freestream p was calculated by Equation (18) and p was 5881 Pa: The speed of sound C air was calculated by Equation (19) and C air was 203.36348 m/s: The Mach number M was 2.85 and the velocity of the freestream U was 579.5859 m/s. The corresponding unit Reynolds number of the freestream was 18 × 10 6 , the characteristic length L was 1.08 m, the density of air was ideal-gas, the specific heat at constant air pressure was 1006.43 j/kg-k, and the viscosity adopted three-coefficient Sutherland Law. Initially, we referred to the value of y + at the transonic stage to estimate y H 3.25 × 10 −6 , while the test results showed that the maximum y + of the wall surface exceeded 1 during the simulation. As described in the previous section, considering the maximum y + of wall surface during the CFD simulation and the boundary layer profile through TVR, the y + value was 0.8 and the y H value was 2.5 × 10 −6 .
To determine the best domain size, a computational region analysis is performed. For this purpose, three sizes were chosen, as shown in Figure 7a [24]. The INLET boundary contains line a and h, the OUTPUT boundary contains lines d and k, the FARDIELD boundary contains lines l, m, and n and the WALL boundary contains lines o, p, and q. Because of the shock wave, in order to obtain more accurate results, more attention was needed for the grids near the shock wave regions. Figure 7b demonstrates the close view of the grids near the shock wave region. Taking the SST k-omega model as an example, the values of the ratio between the wall surface pressure p and p t , and temperature are as shown in Table 4, where x is the distance to the origin point of the wedge section, and the computational domain with 8c was chosen for low ratio errors. Next, three levels of different mesh were applied to ensure grid independency. From Table 5, the mesh with 520,000 cells was applied at this stage.

p
The speed of sound air C was calculated by Equation (19) and air C was 203.36348 m/s: The Mach number M was 2.85 and the velocity of the freestream U was 579.5859 m/s. The corresponding unit Reynolds number of the freestream was 18e + 06, the characteristic length L was 1.08 m, the density of air was ideal-gas, the specific heat at constant air pressure was 1006.43 j/kg-k, and the viscosity adopted three-coefficient Sutherland Law. Initially, we referred to the value of + y at the transonic stage to estimate H y 3.25e-6, while the test results showed that the maximum + y of the wall surface exceeded 1 during the simulation. As described in the previous section, considering the maximum + y of wall surface during the CFD simulation and the boundary layer profile through TVR, the + y value was 0.8 and the H y value was 2.5e-6.
To determine the best domain size, a computational region analysis is performed. For this purpose, three sizes were chosen, as shown in Figure 7a [24]. The INLET boundary contains line a and h, the OUTPUT boundary contains lines d and k, the FARDIELD boundary contains lines l, m, and n and the WALL boundary contains lines o, p, and q. Because of the shock wave, in order to obtain more accurate results, more attention was needed for the grids near the shock wave regions. Figure 7b demonstrates the close view of the grids near the shock wave region. Taking the SST k-omega model as an example, the values of the ratio between the wall surface pressure p and t p , and temperature are as shown in Table 4, where x is the distance to the origin point of the wedge section, and the computational domain with 8c was chosen for low ratio errors. Next, three levels of different mesh were applied to ensure grid independency. From Table 5, the mesh with 520,000 cells was applied at this stage.

Numerical Method
At the supersonic stage, the numerical method is similar to the transonic stage. The validation parameters of the convergence numerical results are shown in Table 6.  Figure 8 shows the wall surface pressure distributions from the experiment and the simulation, where the X-axis x is the distance to the original point of the wedge section and the Y-axis is the ratio between the wall surface pressure and the total pressure of the freestream. The overall tracking effect of the simulation data on the experiment data was acceptable. The simulation errors of the two models were very close at the upper stream of the wedge section. At the wedge section, where x was bigger than 0, the mean error rate of SST k-omega data was 4.76%, and that of Spalart-Allmaras was 5.19%. The variance of SST k-omega simulation error was 0.001043, and that of Spalart-Allmaras was 0.003091. Figure 9 indicates the density profiles at x = 0 from the experiment and the simulation, the X-axis is the ratio between local density and the freestream density, and the Y-axis is the y location. The average error rate of the SST k-omega data was 7.92%, and that of Spalart-Allmaras was 10.04%. The variance of the SST k-omega simulation error was 0.00307, and that of Spalart-Allmaras was 0.003417. In general, the trend of the simulation results was in good agreement with the experimental data, and the SST k-omega simulation results were better than those of Spalart-Allmaras at the supersonic stage. Figure 9 indicates the density profiles at x = 0 from the experiment and the simulation, the X-axis is the ratio between local density and the freestream density, and the Y-axis is the y location. The average error rate of the SST k-omega data was 7.92%, and that of Spalart-Allmaras was 10.04%. The variance of the SST k-omega simulation error was 0.00307, and that of Spalart-Allmaras was 0.003417. In general, the trend of the simulation results was in good agreement with the experimental data, and the SST k-omega simulation results were better than those of Spalart-Allmaras at the supersonic stage.  Figure 9 indicates the density profiles at x = 0 from the experiment and the simulation, the X-axis is the ratio between local density and the freestream density, and the Y-axis is the y location. The average error rate of the SST k-omega data was 7.92%, and that of Spalart-Allmaras was 10.04%. The variance of the SST k-omega simulation error was 0.00307, and that of Spalart-Allmaras was 0.003417. In general, the trend of the simulation results was in good agreement with the experimental data, and the SST k-omega simulation results were better than those of Spalart-Allmaras at the supersonic stage.

Grid Strategy
The experimental data in Ref. [28] were used to validate the simulation parameters at this stage. The Mach number M was 5 and the wedge angle α was 15 • . The shock angle s was 24.32 • , which was calculated by Equation (16). Next, we established and divided the external flow field model for the simulation, as shown in Figure 10. The bold black dot is the (0,0) point of the coordinate system. Line m is the flat plate and line n is the wedge section. The bold black square point is the origin of the wedge section. L is the distance between the flat plate's leading edge and the wedge section origin, which was 0.25 m.
at this stage. The Mach number M was 5 and the wedge angle α was 15°. The shock angle s was 24.32°, which was calculated by Equation (16). Next, we established and divided the external flow field model for the simulation, as shown in Figure 10. The bold black dot is the (0,0) point of the coordinate system. Line m is the flat plate and line n is the wedge section. The bold black square point is the origin of the wedge section. L is the distance between the flat plate's leading edge and the wedge section origin, which was 0.25 m.  (17), and T was 79.17 K; the static pressure of the freestream p was calculated by Equation (18), and p is 850.5173 Pa; the speed of sound air C was calculated by Equation (19), and air C was 178.3999 m/s. The Mach number M was 5 and the velocity of the freestream U was 891.99969 m/s. The density of air was idealgas, the specific heat at constant air pressure was 1006.43 j/kg-k, and the viscosity adopted three-coefficient Sutherland Law. The corresponding unit Reynolds number of the freestream was 6e+06 and the characteristic length c was about 0.34659 m. At this stage, shock waves occurred and play important roles that affected the simulation. As in the previous stage, the computational region analysis was performed and three sizes were chosen, as shown in Figure 11a [24]. The INLET boundary contains lines a, e and o, the OUTPUT boundary contains lines d and h, the FARDIELD boundary contains lines i, j, k, and p, and the WALL boundary contains lines m, n, and q. Figure 11b shows the close view of the refined mesh near the shock wave region. Taking the SST k-omega model at W T = 290 K as an example, the values of the wall pressure coefficient and wall Stanton number (st) were as shown in Table 7, where x is the X-axis location, and the computational domain with 8c was chosen for low ratio errors. Further, three levels of different mesh were applied to be sure of grid independency and, from Table 8, the mesh with 540,000 cells was preferred. At the high supersonic stage, the total temperature of the freestream T t was 475 K and the total pressure of the freestream p t was 4.5 × 10 5 Pa. The thermal conditions at the wall were T w = 100 K leading to a ratio T W /T r = 0.24 and T w = 290 K leading to a ratio T W /T r = 0.72, where T r is the recovery temperature. The static temperature of the freestream T was calculated by Equation (17), and T was 79.17 K; the static pressure of the freestream p was calculated by Equation (18), and p is 850.5173 Pa; the speed of sound C air was calculated by Equation (19), and C air was 178.3999 m/s. The Mach number M was 5 and the velocity of the freestream U was 891.99969 m/s. The density of air was ideal-gas, the specific heat at constant air pressure was 1006.43 j/kg-k, and the viscosity adopted three-coefficient Sutherland Law. The corresponding unit Reynolds number of the freestream was 6 × 10 6 and the characteristic length c was about 0.34659 m. At this stage, shock waves occurred and play important roles that affected the simulation. As in the previous stage, the computational region analysis was performed and three sizes were chosen, as shown in Figure 11a [24]. The INLET boundary contains lines a, e and o, the OUTPUT boundary contains lines d and h, the FARDIELD boundary contains lines i, j, k, and p, and the WALL boundary contains lines m, n, and q. Figure 11b shows the close view of the refined mesh near the shock wave region. Taking the SST k-omega model at T W = 290 K as an example, the values of the wall pressure coefficient and wall Stanton number (st) were as shown in Table 7, where x is the X-axis location, and the computational domain with 8c was chosen for low ratio errors. Further, three levels of different mesh were applied to be sure of grid independency and, from Table 8, the mesh with 540,000 cells was preferred.     Ref. [9] suggests that for predicting the aerodynamic heat of flow around a 2D cylinder, the grid with the cell Reynolds number that represents the minimum normal grid distance to the wall should be less than 8. Therefore, keeping the total number of mesh cells unchanged, meshes with y + = 0.3, 0.15 and 0.08 were applied when looking for the suitable cell Reynolds number for the wedge section at this stage. The corresponding y H and the cell Reynolds number were about (2.46 × 10 −6 , 16), (1.23 × 10 −6 , 8) and (6.56 × 10 −7 , 4), respectively. Ref. [13] demonstrates that the aspect ratio of the wall cells near the shock is a major factor that influences the simulation performance, and we changed the aspect ratio of the wall cells near the shock to confirm an appropriate value under the following conditions: (1) The total number of mesh cells is unchanged. (2) The cell Reynolds number is unchanged. (3) The changes of aspect ratio are operated in small regions of the wall near shock.

Numerical Method
At the high supersonic stage, the density-based solver employing time-derivative preconditioning was used [29]. Ref. [16] concludes that the numerical results of AUSM+ [30,31] and Roe's FDS flux type are closer to the experiment data; these two types were adopted at the high supersonic stage. The validation parameters of the convergence numerical results were as shown in Table 9.

Simulation Parameters
Based on the SST k-omega and Spalart-Allmaras turbulence models, combined with two flux types, we adopted three different cell Reynolds numbers to carry out multiple groups of simulation at T W = 100 K. Figure 12a presents the wall pressure coefficient distributions from the experiment and the simulation adopting the SST k-omega model, where X-axis is the ratio between the location of x and L, and the Y-axis is the wall pressure coefficient. The overall distribution trend of the simulation data matches well with the experiment data and the simulation results are very close to the experiment data at the upper stream of the wedge section, where x/L is between 0 and 1. At the wedge section of interest, where x/L is greater than 1, the error ratio distribution of numerical data with the experimental data is illustrated in Figure 12b. The results based on cell Reynolds number of 16 and the AUSM+ flux type were the closest to the experimental data and the average error ratio was about 10.73%. Further, another four simulations with different aspect ratios based on a cell Reynolds number of 16 and the AUSM+ flux type were performed. The aspect ratios 400, 200, 100, and 50 were adopted; Figure 13a reveals the distribution diagram of the error ratio of the pressure coefficient. As the aspect ratio decreased, the error ratio decreased. When the aspect ratio was equal to 100, the optimal simulation results could be obtained, and the mean error, error variance, and mean error ratio were 0.022, 5.47 × 10 −5 , and 9.1%, respectively. Subsequently, as the aspect ratio decreased, the error ratio increased.  Figure 13b is a comparison of the wall pressure coefficient distribution from the experiment and the simulation adopting the Spalart-Allmaras model at T W = 100 K. The overall distribution trend of the simulation data matched well with the experiment data and the simulation results were very close to the experiment data at the upper stream of the wedge section. The error ratio distribution of the numerical data with the experiment data at the wedge section is displayed in Figure 14a. The results based on a cell Reynolds number of 8 and the AUSM+ flux type were the closest to the experimental data and the average error ratio was about 11.18%. The aspect ratios 400, 200, 100, and 50 were adopted; Figure 14b shows the distribution diagram of the error ratio of the pressure coefficient. As the aspect ratio decreased, the error ratio decreased. When the aspect ratio was equal to 200, the optimal simulation results could be obtained, and the mean error, error variance, and mean error ratio were 0.027, 8.95 × 10 −5 , and 10.9%, respectively. Subsequently, as the aspect ratio decreased, the error ratio increased.
To sum up, at T W = 100 K, a grid strategy with a cell Reynolds number of 16 and aspect ratio of 100, combined with a numerical method, adopting the SST k-omega turbulence model, the density-based solver, and the AUSM+ flux type, obtained an optimal simulation.
Similarly, other multiple groups of simulation were performed at T W = 290 K. The overall distribution trend of the simulation data adopting the SST k-omega model matched well with the experiment data, as can be seen in Figure 15a. Figure 15b describes the error ratio distribution of numerical data with the experiment data at the wedge section. The results based on a cell Reynolds number of 16 and an AUSM+ flux type were the closest to the experimental data and the average error ratio was about 8.27%. Further, using the same number of cells, another four simulations with different aspect ratios based on a cell Reynolds number of 16 and an AUSM+ flux type were performed. The aspect ratios 400, 200, 100, and 50 were adopted; Figure 16a demonstrates the distribution diagram of the error ratio of the pressure coefficient. As the aspect ratio decreased, the error ratio decreased. When the aspect ratio was equal to 200, the optimal simulation results could be obtained, and the mean error, error variance and mean error ratio were 0.018, 8.16 × 10 −5 , and 7.45%, respectively. Subsequently, as the aspect ratio decreased, the error ratio increased. Figure 16b sums up the wall pressure coefficient distributions from the experiment and the simulation at T W = 290 K adopting the Spalart-Allmaras model. The overall distribution trend of the simulation data matched well with the experiment data and the simulation results were very close to the experiment data at the upper stream of the wedge section. Figure 17a shows the error ratio distribution of the wall pressure coefficient at the wedge section. The results based on a cell Reynolds number of 16 and an AUSM+ flux type were the closest to the experimental data and the average error ratio was about 11.18%. The aspect ratios 400, 200, 100, and 50 were adopted; Figure 17b reveals the distribution diagram of the error ratio of the pressure coefficient. As the aspect ratio decreased, the error ratio decreased. When the aspect ratio was equal to 200, the optimal simulation results could be obtained, and the mean error, error variance. and mean error ratio were 0.027, 8.95 × 10 −5 , and 10.9%, respectively. Subsequently, as the aspect ratio decreased, the error ratio increased.
Primarily, at T W = 290 K, the optimal results were achieved based on a grid strategy applying a cell Reynolds number of 16 and an aspect ratio of 200, and a numerical method accepting the SST k-omega model, the density-based solver, and the AUSM+ flux type.
closest to the experimental data and the average error ratio was about 8.27%. Further, using the same number of cells, another four simulations with different aspect ratios based on a cell Reynolds number of 16 and an AUSM+ flux type were performed. The aspect ratios 400, 200, 100, and 50 were adopted; Figure 16a demonstrates the distribution diagram of the error ratio of the pressure coefficient. As the aspect ratio decreased, the error ratio decreased. When the aspect ratio was equal to 200, the optimal simulation results could be obtained, and the mean error, error variance and mean error ratio were 0.018, 8.16e-05, and 7.45%, respectively. Subsequently, as the aspect ratio decreased, the error ratio increased. Figure 16 (b) sums up the wall pressure coefficient distributions from the experiment and the simulation at W T = 290 K adopting the Spalart-Allmaras model. The overall distribution trend of the simulation data matched well with the experiment data and the simulation results were very close to the experiment data at the upper stream of the wedge section. Figure 17a shows the error ratio distribution of the wall pressure coefficient at the wedge section. The results based on a cell Reynolds number of 16 and an AUSM+ flux type were the closest to the experimental data and the average error ratio was about 11.18%. The aspect ratios 400, 200, 100, and 50 were adopted; Figure 17b reveals the distribution diagram of the error ratio of the pressure coefficient. As the aspect ratio decreased, the error ratio decreased. When the aspect ratio was equal to 200, the optimal simulation results could be obtained, and the mean error, error variance. and mean error ratio were 0.027, 8.95e-05, and 10.9%, respectively. Subsequently, as the aspect ratio decreased, the error ratio increased.
Primarily, at W T = 290 K, the optimal results were achieved based on a grid strategy applying a cell Reynolds number of 16 and an aspect ratio of 200, and a numerical method accepting the SST k-omega model, the density-based solver, and the AUSM+ flux type.

Grid Strategy
The experimental data in Ref. [32] were used to validate the simulation parameters at this stage. The Mach number M was 7.05 and the wedge angle α was 20°. The shock angle s was 27.24°, which was calculated by Equation (16). We then established and divided the external flow field model for the simulation, as shown in Figure 18. Line d is the flat plate and line g is the wedge section. The bold black dot is the (0,0) point of the coordinate system. The bold black square point is the origin of the wedge section (or the interaction of the plate and the wedge). At the hypersonic stage, the static temperature of the  Figure 17. (a) error ratio of wall pressure coefficient distribution diagram at the wedge section adopting Spalart-Allmaras at T W = 290 K; (b) error ratio of wall pressure coefficient distribution diagram at the wedge section adopting four different aspect ratios.
In conclusion, the performance of the SST k-omega model was better than that of the Spalart-Allmaras model for the wedge section at the high supersonic stage. When T W /T r is less than 0.3, the mesh with y + is 0.3, the cell Reynolds number is 40 and the aspect ratio is 100, the turbulence model applying Spalart-Allmaras and the density-based solver applying the AUSM+ flux type is preferred. When T W /T r is greater than 0.3, the mesh with y + is 0.3, the cell Reynolds number is 16, and the aspect ratio is 200, the turbulence model applying SST k-omega and the density-based solver applying the AUSM+ flux type is preferred.

Grid Strategy
The experimental data in Ref. [32] were used to validate the simulation parameters at this stage. The Mach number M was 7.05 and the wedge angle α was 20 • . The shock angle s was 27.24 • , which was calculated by Equation (16). We then established and divided the external flow field model for the simulation, as shown in Figure 18. Line d is the flat plate and line g is the wedge section. The bold black dot is the (0,0) point of the coordinate system. The bold black square point is the origin of the wedge section (or the interaction of the plate and the wedge). At the hypersonic stage, the static temperature of the freestream T was 81.2K, the wall temperature T W was 311 K, and the static pressure of the freestream p was 576 Pa. M was 7.05 and the velocity of the freestream U was 1274 m/s. The corresponding unit Reynolds number of the freestream was 5.8 × 10 6 , the characteristic length c was 2 m, the air density s was ideal-gas, the specific heat at constant air pressure was 1006.43 j/kg-k, and the viscosity adopted three-coefficient Sutherland Law. m/s. The corresponding unit Reynolds number of the freestream was 5.8e + 06, the characteristic length c was 2 m, the air density s was ideal-gas, the specific heat at constant air pressure was 1006.43 j/kg-k, and the viscosity adopted three-coefficient Sutherland Law. As in the previous stage, the computational region analysis was performed and three sizes were chosen, as shown in Figure 19a [24]. The INLET boundary contains lines a and h, the OUTPUT boundary contains lines c and j, the FARDIELD boundary contains lines k and l and the WALL boundary contains lines d and g. Figure 19b shows the close view As in the previous stage, the computational region analysis was performed and three sizes were chosen, as shown in Figure 19a [24]. The INLET boundary contains lines a and h, the OUTPUT boundary contains lines c and j, the FARDIELD boundary contains lines k and l and the WALL boundary contains lines d and g. Figure 19b shows the close view of the refined mesh near the shock wave region. Taking the Spalart-Allmaras model as an example, the values of T/T I NF and the values of pressure at S = 0.155 m were as shown in Table 10. S is the distance to the origin point of the wedge along the wedge surface, y is the distance normal to the wedge surface, and T I NF is the freestream's static temperature ahead of interaction, which was 81.2 K. From Table 10, the computational domain with 8c was chosen for low ratio errors. Three levels of different mesh were tested to be sure of grid independency and from Table 11, the mesh with 650,000 cells was preferred. As in the previous stage, the computational region analysis was performed and three sizes were chosen, as shown in Figure 19a [24]. The INLET boundary contains lines a and h, the OUTPUT boundary contains lines c and j, the FARDIELD boundary contains lines k and l and the WALL boundary contains lines d and g. Figure 19b shows the close view of the refined mesh near the shock wave region. Taking the Spalart-Allmaras model as an example, the values of INF T T and the values of pressure at S = 0.155 m were as shown in Table 10. S is the distance to the origin point of the wedge along the wedge surface, y is the distance normal to the wedge surface, and INF T is the freestream's static temperature ahead of interaction, which was 81.2 K. From Table 10, the computational domain with 8c was chosen for low ratio errors. Three levels of different mesh were tested to be sure of grid independency and from Table 11, the mesh with 650,000 cells was preferred.   As in the high supersonic stage, keeping the total number of mesh cells unchanged, meshes with y + = 0.7, 0.35, and 0.18 were applied to looking for the suitable cell Reynolds number for the wedge section at this stage, according to Ref. [14]. The corresponding y H and the cell Reynolds number were about (6.8 × 10 −6 , 40), (3.4 × 10 −6 , 20), and (1.75 × 10 −6 , 10) respectively. Next, based on the appropriate cell Reynolds number, four different aspect ratios were adopted to find the proper aspect ratio for the wedge section under conditions in which the total mesh number and the cell Reynolds number were unchanged.

Numerical Method
As in the high supersonic stage, the density-based solver employing time-derivative preconditioning was used. Roe's FDS and AUSM+ flux types were applied. The validation parameters of the convergence numerical results are shown in Table 12.

Simulation Parameters
Based on the SST k-omega and Spalart-Allmaras turbulence models, combined with two flux types, we adopted three different cell Reynolds numbers to carry out multiple groups of simulations, and the simulation tests showed that a cell Reynolds number of 40 was the appropriate value for the wedge section. S is the distance to the original point of the wedge along the wedge surface and Y is the normal distance to the wedge surface. T INF is the freestream's static temperature ahead of interaction, which was 81.2 K. U INF is the freestream's velocity ahead of interaction, which was 1274 m/s. T and U are the local freestream temperature and velocity, respectively. Figure 20 shows distribution diagrams of T/T INF and U/U INF versus Y from the experiment and the simulation, respectively, at S = 0.055 m. The identical distribution diagrams at S = 0.105 m and S = 0.155 m are depicted in Figures 21 and 22, respectively. The overall distribution trend of the simulation data matched well with the experiment data. The results based on the Spalart-Allmaras turbulence model, the cell Reynolds number of 40 and AUSM+ flux type were the closest to the experimental data, and the average error ratio of T/T INF and U/U INF were about (7.41% 3.1%), (5.46% 2.14%), and (6.32% 2.22%), respectively. Further, another four simulations with different aspect ratios based on Spalart-Allmaras, a cell Reynolds number of 40 and an AUSM+ flux type were performed. The aspect ratios 80, 40, 10, and 2 were adopted, and Figures 23-25 represent the distribution diagrams of the error ratios of T and U. As depicted in Figure 23, at S = 0.055 m, the mean errors of T and U corresponding to the four aspect ratio conditions were (7.41%, 3.1%), (6.04%, 2.07%), (5.97%, 2.05%), and (5.43%, 1.93%), and the optimal simulation results were obtained when the aspect ratio was 2. In Figures 24  and 25, we can observe the same consequences. At S = 0.105 m, the mean errors of T and U corresponding to the four aspect ratio conditions were (5.46%, 2.18%), (4.96%, 2.17%), (4.96%, 2.17%), and (4.89%, 2.16%). At S = 0.155 m, the mean errors of T and U corresponding to the four aspect ratio conditions were (6.32%, 2.59%), (6.04%, 2.56%), (6.03%, 2.56%), and (6.01%, 2.54%). When the aspect ratio decreased from 80 to 40, the error ratios of T and U were reduced to a certain extent. However, with the further reduction of the aspect ratio, the optimization effect of the numerical calculation decreased. At S equals 0.055 m, the maximum changes of error ratio of T and U were 0.54% and 0.12%, respectively; at S equals 0.105 m and 0.155 m, the error ratio of T and U was almost unchanged. To sum up, the aspect ratio for the wedge section at the hypersonic stage was no larger than 40. almost unchanged. To sum up, the aspect ratio for the wedge section at the hyperso stage was no larger than 40.
In conclusion, the performance of the Spalart-Allmaras model was better than that the SST k-omega model for the wedge section at the hypersonic stage. The mesh with is 0.7, the cell Reynolds number is 40, and the aspect ratio is 40, the turbulence mod applying Spalart-Allmaras and the density-based solver applying the AUSM+ flux type preferred.

Simulation Scheme and Aerothermal Prediction
Regardless of the influence of fuel and engine performance, it was assumed that the climbing phase would be uniformly accelerated and that the acceleration process would be completed instantly. According to the flight path of the RBCC vehicle in Figure 1b, the climbing phase was divided into 34 sub-phases in seconds. Each sub-phase remained at a constant speed and the acceleration of the next sub-phase occurred instantly. Every substage was classified into the corresponding stage described in Section 2, based on Mach

Simulation Scheme and Aerothermal Prediction
Regardless of the influence of fuel and engine performance, it was assumed that the climbing phase would be uniformly accelerated and that the acceleration process would be completed instantly. According to the flight path of the RBCC vehicle in Figure 1b, the climbing phase was divided into 34 sub-phases in seconds. Each sub-phase remained at a constant speed and the acceleration of the next sub-phase occurred instantly. Every substage was classified into the corresponding stage described in Section 2, based on Mach  In conclusion, the performance of the Spalart-Allmaras model was better than that of the SST k-omega model for the wedge section at the hypersonic stage. The mesh with y + is 0.7, the cell Reynolds number is 40, and the aspect ratio is 40, the turbulence model applying Spalart-Allmaras and the density-based solver applying the AUSM+ flux type is preferred.

Simulation Scheme and Aerothermal Prediction
Regardless of the influence of fuel and engine performance, it was assumed that the climbing phase would be uniformly accelerated and that the acceleration process would be completed instantly. According to the flight path of the RBCC vehicle in Figure 1b, the climbing phase was divided into 34 sub-phases in seconds. Each sub-phase remained at a constant speed and the acceleration of the next sub-phase occurred instantly. Every sub-stage was classified into the corresponding stage described in Section 2, based on Mach number. An effective simulation scheme for predicting the aerodynamic heat of a scramjet-propelled vehicle is shown in Table 13. The external flow field model for the simulation was established based on a wedge-shaped head and the wedge angle α was 20 • , as shown in Figure 26.     Figure 27 presents the aerodynamic heat of the RBCC scramjet-propelled vehicle in the flight path, as shown in Figure 1b. During the ejection phase, which lasted for about 10 s, the vehicle rose from 11 km to 16 km in altitude, the velocity accelerated from 0.8 Ma to 3 Ma, and the aerodynamic heat increased from an initial temperature of 216.65 K (−56.35 • C) to 543 K (270 • C), with a temperature change of 326.35 • C and an average temperature change rate of 32.635 • C/s. During the sub-combustion stamping stage, which lasted for about 20 s, the vehicle rose from 16 km to 26 km in altitude, the velocity accelerated from 3 Ma to 6 Ma, and the aerodynamic heat increased from a temperature of 543 K (270 • C) to 1550 K (1277 • C), with a temperature change of 1007 • C and an average temperature change rate of 50.35 • C/s. During the super-combustion stamping stage, which lasted for about 4 s, the vehicle rose from 26 km to 30 km in altitude, the velocity accelerated from 6 Ma to 6.5 Ma, and the aerodynamics heat increased from a temperature of 1550 K (1277 • C) to 2010 K (1737 • C), with a temperature change of 460 • C and an average temperature change rate of 115 • C/s. During the cruising phase, the vehicle maintained a constant speed and altitude of 6.5 Ma and 30 km, and the aerodynamic heat finally increased to 2130 K (1857 • C). With the increase of velocity, the average temperature change rate also increased. During a flight time of about 34 s, the aerodynamic heat increased from 216.65 K to 2130 K, and the temperature change reached about 1913.35 • C, which was consistent with the descriptions in Refs. [33][34][35], and the maximum average temperature change rate reached 115 • C/s. The RBCC hypersonic vehicle experienced severe aerodynamic thermal problems during the flight.  Figure 27 presents the aerodynamic heat of the RBCC scramjet-propelled vehicle in the flight path, as shown in Figure 1b. During the ejection phase, which lasted for about 10 s, the vehicle rose from 11 km to 16 km in altitude, the velocity accelerated from 0.8 Ma to 3 Ma, and the aerodynamic heat increased from an initial temperature of 216.65 K (−56.35 °C) to 543 K (270 °C), with a temperature change of 326.35 °C and an average temperature change rate of 32.635 °C/s. During the sub-combustion stamping stage, which lasted for about 20 s, the vehicle rose from 16 km to 26 km in altitude, the velocity accelerated from 3 Ma to 6 Ma, and the aerodynamic heat increased from a temperature of 543 K (270 °C) to 1550 K (1277 °C), with a temperature change of 1007 °C and an average temperature change rate of 50.35 °C/s. During the super-combustion stamping stage, which lasted for about 4 s, the vehicle rose from 26 km to 30 km in altitude, the velocity accelerated from 6 Ma to 6.5 Ma, and the aerodynamics heat increased from a temperature of 1550 K (1277 °C) to 2010 K (1737 °C), with a temperature change of 460 °C and an average temperature change rate of 115 °C/s. During the cruising phase, the vehicle maintained a constant speed and altitude of 6.5 Ma and 30 km, and the aerodynamic heat finally increased to 2130 K (1857 °C). With the increase of velocity, the average temperature change rate also increased. During a flight time of about 34 s, the aerodynamic heat increased from 216.65 K to 2130 K, and the temperature change reached about 1913.35 °C, which was consistent with the descriptions in Ref. [33][34][35], and the maximum average temperature change rate reached 115 °C/s. The RBCC hypersonic vehicle experienced severe aerodynamic thermal problems during the flight.

Conclusions
In this paper, an effective simulation scheme for predicting the aerodynamic heat of a scramjet-propelled vehicle during flight is proposed. Each phase adopts a different grid strategy and numerical method, which were confirmed by the validation tests described in Section 3. In particular, the flux type dependencies, the cell Reynolds number dependencies, and the aspect ratio dependencies are discussed in Sections 3.3 and 3.4. The aerothermal simulation based on a wedge-shaped head of the RBCC scramjet-propelled vehicle during the climbing phase and cruising phase were performed. Through systematic analysis, several conclusions can be drawn, as follows: 1.
For the wedge section in the scramjet-propelled vehicle heating prediction, the suggested value of cell Reynolds number and aspect ratio of wall cells near the shock are different from those for blunt bodies. At the high supersonic stage, the ideal cell Reynolds number should be no larger than 16. At the hypersonic stage, neither the ideal cell Reynolds number, nor the aspect ratio of the wall cells near the shock should be larger than 40.

2.
For the wedge section in the scramjet-propelled vehicle heating prediction, with the same grid strategy, upwind scheme, solver, and so on, the AUSM+ flux type performs better than Roe's FDS flux type. The Spalart-Allmaras turbulence model performs better at the transonic stage and the hypersonic stage, while the SST k-omega performs better at the supersonic stage and the high supersonic stage. 3.
Smaller cell Reynolds numbers are not necessarily better. On the one hand, this may cause the quality of the mesh to decrease. On the other hand, as described in Sections 3.3 and 3.4, their improvement of calculation accuracy is limited; they even reduce the accuracy. The value of the aspect ratio is similar to the cell Reynolds number. It is necessary to confirm the appropriate values through CFD analysis based on numerical results and TVR plots.