Abstract
The design of the thermal protection system requires high-precision and high-reliability CFD simulation for validation. To accurately predict the hypersonic aerodynamic heating, an overall simulation strategy based on mutual selection is proposed. Foremost, the grid criterion based on the wall cell Reynolds number is developed. Subsequently, the dependence of the turbulence model and the discretization scheme is considered. It is suggested that the appropriate value of wall cell Reynolds number is 1 through careful comparison between one another and with the available experimental data. The excessive number of cells is not recommended due to time-consuming computation. It can be seen from the results that the combination of the AUSM+ discretization scheme and the Spalart-Allmaras turbulence model has the highest accuracy. In this work, the heat flux error of the stagnation point is within 1%, and the overall average relative error is within 10%.
1. Introduction
The axisymmetric shape is widely used in hypersonic missiles, such as the Kinzhal missile and the AHW [1] program. Hypersonic vehicles achieve high speeds (Mach 5-20) via rocket propulsion or scramjet engines [2], high lift-drag design [3, 4], and the advantage of near-space. On the one hand, high-speed flight makes it difficult to detect and track the hypersonic vehicle. On the other hand, high-speed flight will cause violent friction between the surface and the air and generate a lot of heat. Ref. [5–7] illustrated that at Mach 6, the skin temperature may exceed 900 K, and the leading edge temperature will easily exceed 1000 K, which presents significant challenges for thermal protection. Thermal protection system (TPS) can provide hypersonic vehicles the thermal insulation from the aerodynamic heating [5, 8, 9]. This subsystem of a vehicle is used to ensure the safety of the aircraft, especially when there exists a large temperature difference between the cabin interior and the TPS structure. Therefore, the reliable prediction of surface heat fluxes is essential to the TPS, and its size has to be the minimum value able to guarantee suitable protection [10]. The following three methods [11] are commonly used to compute heat flux to guide the design of TPS: experiment/wind tunnel experiment, computational fluid dynamics (CFD) simulation, and engineering experience estimation method. Compared with the engineering empirical estimation method, CFD simulation is more accurate. Compared with the experimental/wind tunnel experiment, CFD simulation cycle is shorter, the model size does not need to be reduced, and the simulation parameters can be flexibly modified. Therefore, to better guide the design of TPS for hypersonic aircraft [12], how to accurately predict aerodynamic heating with CFD method has become a desirable concern.
In the whole process of CFD simulation, each step will affect the final numerical simulation accuracy, including mesh generation, turbulence model, and discretization scheme setting. The poor simulation strategy may lead to errors of 4 times [13] or more. As the heat transfer rate is determined by the temperature gradient [14], the mesh quality becomes the main factor. For this reason, a lot of research work has been carried out on the grid strategy, especially on the first layer grid [15, 16]. Wall cell Reynolds number is an important scale criterion. Its reasonable value range becomes a research hotspot. Ref. [17, 18] showed that a cell Reynolds number of about 20 in the first grid layer was reasonable. It is found that the heat flux converges when the wall cell Reynolds number is less than or equal to 3 [19]. However, the value range mentioned above needs to be more specific to meet the needs of high-precision simulation. In addition to grid research, there are many issues worthy of attention, such as turbulence models and discretization schemes. The results of hypersonic verification and evaluation show that different turbulence models and discretization schemes have great differences in predicting the accuracy of aerodynamic heating [13, 20]. The existing studies separate influencing factors independently, but CFD simulation is a whole process. It is limited to analyze single factors on the accuracy, which need to be considered together [21]. Besides, the detailed and quantitative strategy is desirable for the main factor grid.
In this study, a new grid strategy based on the wall cell Reynolds number is proposed by comparing the accuracy of aeroheating prediction. After getting the appropriate mesh in the preprocessing stage, the dependence between the turbulence model and the discretization scheme is found and characterized. Considering all factors comprehensively rather than individually, the overall CFD simulation strategy proposed in this paper can realize high-precision aerodynamic heating simulation, which provides a data basis for thermal protection design.
2. Methods
As shown in Figure 1, the entire process can be divided into four steps. Our approach requires consideration of the whole simulation to achieve high-accuracy prediction. Foremost, the geometry part is necessary. The focus of the first step is that the geometry model should be as simple as possible and as complicated as necessary. The second step is to divide the imported geometry model into discrete cells. The mesh quality has an important influence on the convergence and results. After importing the meshing document, the setup of a simulation can be specified. At the end of the procedure, the analysis of the results can give feedback to the simulation.
2.1. Geometry and Meshing
In this work, we take the blunt double cone as an example, regarding the wind tunnel test report of NASA TP-2334 [22]. The geometry model is created in UG NX 10.0, and its planar diagram is shown in Figure 2(a). The curvature radius of the head is 0.3835 cm, and the half cone angles are 12.84° and 7°, respectively, and the length of the model is 12.224 cm.
(a)
(b)
Since this research focuses on the flow of the boundary layer, the structured mesh is used to avoid the undesirable large discrete error. Equation (1) illustrates that the heat flux is determined by the temperature gradient, so the set of first layer height is essential. As shown in Figure 2(b), the ICEM-CFD software is used to mesh the whole computational domain into structured grids. Because the model is axisymmetric, it is divided into six blocks with half of the computational domain. Besides, the grid is encrypted in the boundary layer, especially near the head with severe heating.
The free flow conditions refer to the wind tunnel test conditions, as shown in Table 1. The boundary conditions are the far-field inlet, pressure outlet, and isothermal wall, respectively.
The wall cell Reynolds number [23] refers to the local Reynolds number based on the velocity in the calculation cell and the length scale of the cell, as shown in Equation (2). In which, denotes the Reynolds number in Table 1. The means the normal height of the boundary layer.
The quantitative grid strategy proposed in this work is based on the following two criteria: (1) On account of the prediction accuracy, the first cell height will be set to 2.8-7 m, 5.6-7 m, …,1.8-5 m, corresponding to the wall cell Reynolds number ranging from 0.5 to 32. (2) The normal expansion ratio of the grid is fixed at 1.1, so the total number of cells is only related to the density of the first layer. At any of the wall cell Reynolds numbers, the value of total elements ranges from 1 million to 5 million.
2.2. Turbulence Model
CFD simulation is essentially a process of giving equations and solving discretization solutions, while the turbulence model is used to describe the transport equation(s). Because direct numerical simulation (DNS) requires vast computing resources, the Reynolds average (RANS) method is widely used to process the time average, as shown in Equations (3) and (4), to solve four physical quantities: three velocity components and pressure .
Although the RANS method reduces the computational cost, it also introduces a nonlinear term, which leads to the unclosed N-S equation. Therefore, it is necessary to introduce additional turbulent flow variable(s) and differential equation(s). In this paper, the Spalart-Allmaras model and the Shear Stress Transport (SST) model, which are widely used and perform well in aerothermal simulation, are selected [24, 25].
Spalart-Allmaras model [26] is a one-equation model. The introduction of the variable is equivalent to turbulent kinematic viscosity. The transport equation of this model is as follows: where is the production of turbulent viscosity, is the destruction of turbulent viscosity, and is constant.
As a one-equation model, the SA model is naturally missing the term (turbulent kinetic energy). But this is not a major effect in thin shear flows, and the addition of to the diagonal elements of the stress tensor is approximate in any case. The SA model is the eddy viscosity model which has a robust numerical formulation and has shown promising results for a wide variety of flows. Unlike most two-equation models, the SA model does not require additional wall function when dealing with near-wall flow problems, so it is quite insensitive to the wall spacing. The surface properties varied by less than 2.5% for the flow when the wall spacing was varied between 0.01 and unity [27]. While stable for large values, the maximum for accurate solutions of the SA model should be roughly .
SST two-equation model [26] was proposed and developed by Menter. It combines the standard model and the standard model, which will use the model for the near-wall region and the model for the free surface. The SST model introduces two variables: the turbulence kinetic energy and the rate of dissipation of the eddies . The transport equation of this model is as follows:
where and are the production items of item and item , respectively, is the effective diffusivity, and is the dissipation of turbulence.
The SST model utilizes the same blending function () between the and models as the baseline (BSL) model; however, the former model also employs a modified form of the eddy viscosity definition which accounts for the transport of the Reynolds stress. The eddy viscosity is redefined in the following way: where is a blending function that guarantees that the model is satisfied whereas the original formulation is used for the rest of the flow. Based on the philosophy underlying the Johnson-King model, the SST model leads to a significant improvement for all flows involving adverse pressure gradients and should be the model of choice for aerodynamic applications. The SST model has greatly benefited from the strength of the underlying turbulence models. Compared with the other two-equation models, the SST model is more suitable for many industrial applications and most general-purpose commercial CFD codes.
2.3. Discretization Scheme
Discretization is the process of transferring continuous functions, models, variables, and equations into discretized counterparts. As shown in Ref. [20, 21, 28], the Roe-FDS scheme and ASUM+ scheme have been widely used and achieved good results. Therefore, the following work focuses on the impact of these two discretization schemes:
2.4. Roe-FDS [29]
Because the Godunov format requires a lot of computation, Roe linearizes the Jacobian matrix and proposes a new flux difference scheme on this basis. This scheme is a flux differential splitting scheme (FDS) that approximates the Riemann solution: where represents the Roe-averaged values and represents the diagonal matrix of characteristic velocity. and are left and right eigenvectors, respectively.
2.5. AUSM+ [30]
The AUSM scheme is an improvement of the Vanleer scheme in terms of scheme structure, and it is a composite scheme of FVS (flux vector splitting scheme) and FDS in terms of the dissipation term. The AUSM + scheme proposed by Liou, which divides the inviscid flux into the convection term and the pressure term : where the interface sound velocity and .
3. Results and Discussion
In this section, the simulation results are imported into CFD-post for intuitive expression and analysis. By comparing the results of different simulation parameter settings, the influence of each stage and mutual dependencies are analyzed. The inlet is set as the far-field boundary, and the parameters of free flow are shown in Table 1. The outlet is set as the pressure out boundary. The wall is the surface of the biconic model, and a no-slip isotherm wall condition is adopted. In order to save computing resources, half of the computational domain is used and symmetric boundary conditions are adopted.
3.1. Influence of Meshing Grid
We will discuss the influence of meshing, i.e., grid strategy, on CFD accuracy during hypersonic aerodynamic thermal numerical simulation in this section. According to Ref. [16, 18, 31], the number of grids and the height of the first layer play a significant role in the accuracy of simulation results and even determine whether they can converge. For this reason, the structured grid is divided in this paper. Figure 3(a) shows the computational domain of a blunt biconic, and Figure 3(b) shows the boundary layer and stagnation point details. Considering the angle of attack as 0°, 35 sets of grids are divided. The grid parameters are specified in Table 2, and the total number is in terms of .
(a)
(b)
In this case, the closer to the nose cone head, the more severe the heat. Therefore, the plane is selected to display the flow field pressure contour and temperature contour, as shown in Figure 4. The two different kinds of grid No. 2-4 and No. 3-1 (see in Table 2) are chosen to explain the impact of grid strategy clearly. As shown in Figure 5, the data on the solid black line in Figure 4 are selected for comparative analysis to clearly explain the influence of grid strategy. It can be seen that even with different wall cell Reynolds numbers and different densities of the grid, the pressure contours are similar. The difference between the temperature results in the two grids is so small that it is negligible. The conclusion we can draw from this is that the temperature and pressure are not particularly dependent on the grid, so there is no need for severe simulation conditions.
(a)
(b)
(c)
(d)
(a) Pressure profile at mm
(b) Temperature profile at mm
Figure 6 is the wall heat flux simulation results with two different grids. It can be seen that the overall distribution is quite different, especially at the stagnation point. In which, the maximum relative error can be 10%. We can conclude that the wall heat flux result is highly dependent on the grid, so it is desirable to show an appropriate grid strategy.
(a) Wall heart flux contour with grind No. 1-4
(b) Wall heart flux contour with grind No. 3-1
(c) Wall heart flux distribution with different grinds
As shown in Figure 7, after comparing the results obtained by CFD simulation with the wind tunnel test value , the following conclusions can be drawn: (1) The accuracy of the heat transfer rate at the stagnation point increases with the increase in the number of grids at first, but after the number of elements reaches 3 million, the improvement of accuracy is almost negligible. (2) Meanwhile, the accuracy relies on the height of the first layer grid. The accuracy is the highest when the wall cell Reynolds number is about 1. When the wall cell Reynolds number is greater than 16, the prediction is quite poor. (3) Compared with the previous studies, the simulation accuracy can be significantly improved (as shown in Table 3) by using an appropriate grid strategy, and the relative error is less than 1%. (4) As shown in Table 3, the method presented in this paper performs better than the theory of Fay and Riddell [32] that commonly used in engineering. (5) Two sets of grids are used in the subsequent analysis, of which 3 million elements are used with Reynolds numbers of wall elements set at 1 and 2 (grid No. 3-1 and grid No. 3-2).
3.2. Influence of Turbulence Model and Discretization Scheme
In this section, the analysis of the influence of the turbulence model and discretization scheme on the simulation accuracy will be carried out through careful comparison of experimental results and theory. There are several criteria for verifying the quality of different turbulence models: (1) Turbulence model sensitivities (grid wall spacing and turbulence level); (2) convergence of time and step iteration; (3) quantitative precision of the numerical simulation when compared with experimental data. According to the analysis in the previous section, although the SA model is insensitive to mesh refinement and wall spacing, in general, the heat flux simulation error is larger when is larger. Therefore, the grid strategy proposed in this paper can significantly improve the simulation accuracy. Besides, in our cases, uniform iteration step size and iteration time are adopted to ensure that the residual of each run is within 1-3.
In this paper, the accuracy of numerical simulation and robustness of turbulence model are analyzed through careful comparison between one another and with the available experimental data. When the wall cell Reynolds number is large (), the transition is easy to occur in the boundary layer, and turbulence is more common in the flow. Compared with laminar flow, turbulent flow is more chaotic and irregular. In this work, two sets of grids (grid No. 3-1 and grid No. 3-2) are used for comparative analysis on the simulation accuracy for different turbulence models. The heat flux on the wall surface that we obtained is normalized (divided by ) for quantitative comparison. Because the heat flux at the stagnation point is too large, this work only shows the results of the second half, as shown in Figure 8.
In addition to simulation accuracy, the calculation time is also an important factor that must be measured in practical engineering. The overall mean errors of the selected points on the wall surface with different models are also compared, as shown in Table 4. From the above discussion, the following conclusions can be drawn: (1) In our case, the solution of the boundary mainly focuses on the viscous sublayer. Although it is based on the accurate and robust near-wall formulation of the Wilcox model, the additional near-wall functions are still somewhat unwieldy, so the performance of the SA model is still better than that of the SST model. (2) Near the outer edge of the boundary layer and in shear layers, the SST model blends into a transformed version of the formulation, thus providing the same good predictions for temperature and pressure in the flow fields as the SA model (see in Figures 4–6). (3) The main additional complexity in the SST model formulation compared to standard models lies in the necessity to compute the distance from the wall, which results in a longer calculation time than the SA model.
For the transient N-S equation based on density solver, the discretization scheme is an important factor that affects the accuracy and stability of the solution. Compared with the central scheme, the upwind scheme has its dissipation characteristic, which is closer to the real flow field. Therefore, in this work, the influence on the accuracy of CFD simulation with the AUSM+ scheme and Roe-FDS scheme is discussed.
Similarly, the simulation data on the surface of the blunt double cone are normalized. Because the heat flux at the stagnation point is too large, this work only shows the results of the second half. The results with the two grids (grid No. 3-1 and grid No. 3-2) are shown in Figure 9. It can be seen that the results of the AUSM+ scheme and Roe-FDS scheme are almost consistent, which can predict the aerodynamic heating. The prediction accuracy of heat flux in our work is within 1%, as shown in Table 5.
3.3. The Dependence between Turbulence Model and Discretization Scheme
The above is the independent influence analysis of the turbulence model and the discretization scheme, which is also the main research theme of most current work. In this work, the dependence between the two factors is discussed. Similarly, the first pretreatment is normalization. The results of the wall heat flux are shown in Figure 10. The following conclusions can be drawn: (1) The maximum relative error is within 15%. (2) As for the SST turbulence model, the AUSM+ scheme and Roe-FDS scheme are desirable. (3) As for the SA turbulence model, the AUSM+ scheme has performed far better than the Roe-FDS scheme.
4. Conclusions
It has been widely concerned to achieve accurate prediction of the heating load under the condition that the hypersonic flow is strongly compressed and frictional. In this work, an overall simulation strategy based on a quantitative grid criterion and the interdependence of different stages is proposed. Through careful comparisons between one another and with the available experimental data, the conclusions are as follows: (1)Although the SA model and SST model are quite insensitive to the wall spacing when compared with other turbulence models. The results show that the heat flux is strongly dependent on the grid quality. Based on the analysis of the heat flux numerical results at stagnation point with different meshing grids, the wall cell Reynolds number is suggested to be 1(2)The dependence between the turbulence model and the discretization scheme cannot be ignored. It seems hard to analyze these two factors independently. Indeed, there is neither an absolutely good turbulence model nor an absolutely good discretization scheme. For the SST model, Roe-FDS is more suitable than AUSM+; conversely, AUSM+ is a better choice for the SA model. Different collocation will result in a great accuracy difference. In this paper, the average relative error of surface heat flux is reduced from 14.3% to 9.3% by comprehensive consideration(3)The leading edge region has the highest temperature and heat flux and is most likely to be damaged. Therefore, accurately predicting the aerodynamic heating loads of the hypersonic vehicle is of great significance. In this work, the relative error of stagnation point heat flux is less than 1% with the overall CFD simulation strategy. Accurately predicting heat flux can be used to provide some reference for the design of thermal protection systems to meet the requirements of the severe flight conditions and lightweight
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This research was funded by the National Natural Science Foundation of China (No. 61572307).