An engineering method of aerodynamic heating prediction for hypersonic blunt body vehicles

Abstract

This paper develops a rapid engineering method for predicting aerodynamic heating of hypersonic blunt body vehicles and implements a complete computational program in C++. This engineering method is used to reduce the time required to predict the aerodynamic heating for vehicles. It utilizes the Kemp-Riddell formula to calculate stagnation point heat flux, while the downstream region heat flux is determined using the reference enthalpy method. Additionally, a simplified streamline tracing approach is proposed to calculate inviscid surface streamlines, achieving a 90% improvement in computational efficiency. A mean filtering method is also introduced to triangular surface meshes to effectively improve the smoothness of the local Reynolds number on low-density surface meshes. The engineering method is validated against experimental data from a spherically blunted cone and an Orbiter Vehicle model, showing good agreement in wall heat flux predictions for small angles of attack, with a relative error of less than 15% in the stagnation and non-expanded downstream regions. A comparison between the perfect gas model and the chemical equilibrium gas model indicates that, while their results are generally similar, the perfect gas model reduces computational time by 70% in related calculation processes, making it suitable for most conditions.

Introduction

The prediction of hypersonic aerodynamic heating forms the foundation for designing and developing thermal protection systems intended for use in hypersonic vehicles. It carries substantial significance during the preliminary design stage of reentry vehicles like spacecraft return capsules, space shuttles, and missiles. The techniques employed to predict aerodynamic heating mainly encompass numerical simulation, experiment, along with engineering approximation. Although ground wind tunnel tests can yield dependable data regarding aerodynamic heating, their applicability across diverse real-flight scenarios or ability to generate extensive amounts of aero-thermodynamics environmental information swiftly is limited due to equipment constraints coupled with high economic costs associated with conducting these experiments over a short timeframe. Although there have been notable advancements in computational fluid dynamics (CFD) techniques and computer performance, numerical methods continue to demand substantial computational time for solving aerodynamic heating problems. Consequently, they are not well-suited for rapid computations or optimization studies. In contrast, engineering approximation methods offer expedited calculation efficiency. Hence, employing engineering approximation methodologies to predict aerodynamic heating assumes a crucial role during the preliminary design phase of hypersonic aircraft.^1–3

The majority of engineering methods used to calculate aerodynamic heating are based on Prandtl’s boundary layer theory. These methods consider flow outside the boundary layer as inviscid while introducing viscous correction inside the boundary layer. The inviscid flow can be typically determined by solving Euler equations with numerical method. Alternatively, independent panel methods grounded in Newtonian theory can be employed to compute surface pressures. ⁴ Other gas parameters can be estimated assuming an isentropic flow condition. Viscous corrections inside the boundary layer can be categorized into simplified approximation methods and the axisymmetric analogy. Simplified approximation methods approximate different components of an aircraft as simple two-dimensional models and estimate aerodynamic heating on the aircraft’s surface based on analytical solutions derived from these two-dimensional models. ² On the other hand, Cooke first proposed the axisymmetric analogy, ⁵ which simplified three-dimensional boundary layer equations into the same form as for the axisymmetric flow. Several researchers have studied and improved the axisymmetric analogy method,^6–9 which enhanced its computational accuracy and speed. Nevertheless, determining continuous inviscid surface streamlines and calculating streamline metrics remain intricate tasks within this framework. In preliminary aircraft design stage, there is a preference for employing simpler calculation methods. One such method that combines Reynolds analogy ¹⁰ with reference enthalpy theory, ¹¹ provides a simplified approach for calculating convective heat transfer within boundary layers as compared to utilizing an axisymmetric analogy. Li and Gao ¹² developed an engineering combination method specifically designed for predicting aerodynamic heating in complex reentry vehicle shapes. They utilized the Fay-Riddell method ¹³ at stagnation points, employed approximate plate methods for downstream regions, and adopted empirical hemispherical methods ¹⁴ for the junction zone.

This paper presents an engineering method for predicting aerodynamic heating in blunt-body reentry vehicles based on the approach developed by Li and Gao. To enhance computational efficiency and simplify programming, this study introduces a simplified approach for calculating inviscid surface streamlines. Points along these streamlines are positioned at the centroids of corresponding surface elements. In this study, independent panel methods grounded in Newtonian theory are employed to predict surface pressure while estimating other gas parameters under assumed conditions of isentropic flow. The heating rate at the stagnation point is computed using the Kemp-Riddell method, ¹⁵ while approximate flat plates methods based on Reynolds analogy and Eckert’s reference enthalpy theory are utilized for downstream regions. Additionally, surface heating rate distributions are computed for spherically blunted cone and the OV model (orbital vehicles) to verify the accuracy of this engineering method.

Aerodynamics properties evaluationon the edge of boundary layer

The pressure on surface

For high Reynolds number flows with an attached boundary layer on the aircraft surface, the pressure gradient along the surface normal within the boundary layer is zero. Therefore, the pressure on the edge of this boundary layer equals the surface pressure predicted by Newtonian theory.

The surface pressure $p_e$ can be calculated using the definition of the compressible flow pressure coefficient $C_p$ , as expressed by

$$C_p=(\frac{p_e}{p_{\infty }}-1)\frac{2}{\gamma {\mathrm{Ma}}_{\infty }^2}$$

(1)

where $\gamma$ represents the specific heat ratio. $M{a}_{\infty }$ and $p_{\infty }$ represent the Mach number and pressure of the freestream flow. The modified Newtonian theory proposed by Lees ¹⁶ is adopted to calculate the pressure coefficient on the windward side, given by

$$C_p=C_{p_{\max }}\mathrm{si}{n}^2\theta .$$

(2)

Here, $C_{p_{\max }}$ is the maximum value of pressure coefficient, determined based on the stagnation pressure behind normal shock as

$$\begin{array}{c}{C}_{p_{\max }}=\frac{2}{\gamma {\mathrm{Ma}}_{\infty }^2}\{{\left[\frac{{(\gamma +1)}^2{\mathrm{Ma}}_{\infty }^2}{4\gamma {\mathrm{Ma}}_{\infty }^2-2(\gamma -1)}\right]}^{\gamma /(\gamma -1)}\hfill \\ \multicolumn{1}{c}{\left[\frac{1-\gamma +2\gamma {\mathrm{Ma}}_{\infty }^2}{\gamma +1}\right]-1\}.}\end{array}$$

(3)

The pressure coefficient on the leeward side is calculated using the Prandtl-Meyer formula, ¹⁷ given by

$$C_p=-\frac{\gamma +1}{2}{\theta }^2\{\sqrt{1+{\left[\frac{4}{(\gamma +1)M{a}_{\infty }\theta }\right]}^2}-1\}.$$

(4)

In the above equations, $\theta$ represents impact angle of the surface element, defined as

$$\theta =\arccos \left(\frac{\overrightarrow{\mathit{n}}•\overrightarrow{\mathit{V}}}{\left|\overrightarrow{\mathit{n}}\right|•\left|\overrightarrow{\mathit{V}}\right|}\right)-\frac{\pi }{2}$$

(5)

where $\overrightarrow{\mathit{n}}$ is the outward normal vector of the surface element. $\overrightarrow{\mathit{V}}$ is the velocity vector of the freestream. A positive $\theta$ corresponds to the windward side, while a negative $\theta$ corresponds to the leeward.

Other gas properties

Once the pressure distribution on the surface has been obtained, additional thermodynamic parameters on the edge of the boundary layer can be determined by combining assumptions of isentropic flow and a gas model (used for non-perfect gases).

In case of a perfect gas, it becomes necessary to compute various parameters including post-bow shock values such as gas pressure $p_2$ , density ${\rho }_2$ , and Mach number $M{a}_2$ . These parameters can be calculated using the following isentropic flow oblique shock relations ¹⁷ :

$$p_2=p_{\infty }[\frac{2\gamma {\mathrm{Ma}}_{\infty }^2\mathrm{si}{n}^2\beta }{1+\gamma }-\frac{\gamma -1}{\gamma +1}]$$

(6)

$${\rho }_2={\rho }_{\infty }\left[\frac{(\gamma +1){\mathrm{Ma}}_{\infty }^2\mathrm{si}{n}^2\beta }{2+(\gamma -1){\mathrm{Ma}}_{\infty }^2\mathrm{si}{n}^2\beta }\right]$$

(7)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{Ma}}_2^2=\frac{{\mathrm{Ma}}_{\infty }^2+\frac{2}{\gamma -1}}{\frac{2\gamma }{\gamma -1}{\mathrm{Ma}}_{\infty }^2\mathrm{si}{n}^2\beta -1}}\\ \multicolumn{1}{c}{+\frac{\frac{2}{\gamma -1}{\mathrm{Ma}}_{\infty }^2\mathrm{co}{s}^2\beta }{{\mathrm{Ma}}_{\infty }^2\mathrm{si}{n}^2\beta +\frac{2}{\gamma -1}}}\end{array}$$

(8)

The parameter $\beta$ denotes the shock wave angle, while ${\rho }_{\infty }$ represent the density of the freestream flow ahead of the shock wave.

However, determining the precise local shock wave angle on the surface is challenging. In engineering estimations for blunt bodies, it is common to assume a value of $\beta ={90}^{°}$ . Under this assumption, thermodynamic parameters ${\rho }_e$ , $T_e$ , and $h_e$ on the edge of the boundary layer can be calculated, respectively, by $p_e=p_2{p}_e^{\gamma }/{\rho }_2^{\gamma }$ , ${\rho }_e=T_{e}R/p_e$ , and $h_e=\gamma R{T}_e/(\gamma -1)\times {10}^{-3}$ . Here, ideal gas constant $R=287\mathrm{J}/\mathrm{kg}·\mathrm{K}$ . The unit of $h_e$ takes $\mathrm{kJ}/\mathrm{kg}$ . Then, Mach number $M{a}_e$ and velocity $u_e$ on the edge of the boundary layer can be calculated using the following relationships

$$M{a}_e=\sqrt{[{\left(\frac{p_{02}}{p_e}\right)}^{\frac{\gamma -1}{\gamma }}-1]\frac{2}{\gamma -1}}$$

(9)

$$u_e=M{a}_e\sqrt{\gamma R{T}_e},$$

(10)

where $p_{02}$ is the total pressure behind the shock wave.

Equation (9) is derived from the relationship between total pressure and static pressure for an isentropic gas, while equation (10) represents the relationship between Mach number and velocity for a perfect gas. Additionally, the viscosity coefficient ${\mu }_e$ on the edge of the boundary layer can be calculated using the Sutherland formula. When dealing with non-perfect gases, it’s typically necessary to utilize appropriate gas models or fitting formulas in order to calculate thermodynamic and transport parameters at the edge of the boundary layer. In engineering calculations, a common assumption made places the starting position of inviscid streamlines on a blunt body surface at its stagnation point while assuming an isentropic flow along these streamlines. Therefore, it can be assumed that entropy on the edge of the boundary layer equals that observed at the stagnation point. Consequently, calculating entropy at stagnation point serves as the first step.

For high temperature chemically equilibrium gases, the following group of fitting formulas of thermodynamic and transport parameters are employed. ¹⁸ The entropy can be calculated by

$$S=\{\begin{array}{c}3.07311-0.28711\ln \left(\frac{p}{p_0}\right)+1.0048\ln \left(\frac{h}{h_0}\right),\hfill \\ \multicolumn{1}{c}{167.5<h\le 1161.8(\mathrm{kJ}/\mathrm{kg})}\\ \multicolumn{1}{c}{[5.98796-0.1957\ln \left(\frac{p}{p_0}\right)]\exp \left(\frac{{\left(\frac{h}{418.68}\right)}^{0.275}}{3.548}\right),}\\ \multicolumn{1}{c}{1161.8<h<62383(\mathrm{kJ}/\mathrm{kg})}\end{array}$$

(11)

where $p_0=101325\mathrm{Pa}$ and $h_0=4.1868\mathrm{kJ}/\mathrm{kg}$ . By substituting $h_s$ and $P_s$ into equation (11), one can determine the entropy at the stagnation point as well as on the edge of the boundary layer. If the entropy $S_e$ and pressure $P_e$ are known, the enthalpy $h_e$ on the edge of the boundary layer can be determined using the bisection method and equation (11). After obtaining the enthalpy $h_e$ on the edge of the boundary layer, the velocity $u_e$ on the edge of the boundary layer can be determined by applying the energy conservation equation, given by

$$u_e=\sqrt{2(h_s-h_e)\times {10}^3}$$

(12)

where the enthalpy of stagnation $h_s$ can be calculated by ¹⁴ $h_s=1.0084T_{\infty }+\frac{V_{\infty }^2}{2}\times {10}^{-3}$ . Density ${\rho }_e$ on the edge of the boundary layer can be calculated by the following fitting equations ¹⁹ :

$$\rho =0.213833\left(\frac{p}{p_0}\right){\left(\frac{h}{1755.52}\right)}^{-b}$$

(13)

$$b=\{\begin{array}{c}0.972,167.5<h\le 1755.52(\mathrm{kJ}/\mathrm{kg})\hfill \\ \multicolumn{1}{c}{0.715+1.38974\times {10}^{-2}\ln \left(\frac{p}{p_0}\right),}\\ \multicolumn{1}{c}{1755.52<h<35026(\mathrm{kJ}/\mathrm{kg})}\end{array}$$

(14)

Temperature $T_e$ can be evaluated by ¹²

$$T=\{\begin{array}{c}1.2431h^{0.9614},167.5<h\le 1881.1(\mathrm{kJ}/\mathrm{kg})\hfill \\ \multicolumn{1}{c}{1.748\times {10}^3{(0.31469\ln \left(\frac{h}{h_0}\right)-0.9221)}^c,}\\ \multicolumn{1}{c}{1.881.1<h<46055(\mathrm{kJ}/\mathrm{kg})}\end{array}$$

(15)

$$c=2.41+7.09637\times {10}^{-2}\ln \left(\frac{p}{p_0}\right).$$

(16)

The dynamic viscosity ${\mu }_e$ is given as ¹²

$$\mu =\{\begin{array}{c}0.296109\times {10}^{-6}{h}^{0.7259},167.5<h\le 711.7(\mathrm{kJ}/\mathrm{kg})\hfill \\ \multicolumn{1}{c}{1.05311\times {10}^{-6}{h}^{0.5325},711.7<h\le 1881.1(\mathrm{kJ}/\mathrm{kg})}\\ \multicolumn{1}{c}{58.4182\times {10}^{-6}{[0.31469\ln \left(\frac{h}{h_0}\right)-0.9221]}^d,}\\ \multicolumn{1}{c}{1881.1<h<51497(\mathrm{kJ}/\mathrm{kg})}\end{array}$$

(17)

$$d=1.75+3.99551\times {10}^{-2}\ln \left(\frac{p}{p_0}\right).$$

(18)

The local speed of sound $a_e$ is calculated using the following fitting equation ¹²

$$a=\frac{687.8\sqrt{\frac{h}{1755.52}}}{\sqrt{{\left(\frac{h}{1755.52}\right)}^{1-b}-0.2698b}}$$

(19)

where $b$ is given by equation (14).

Approximate method for heating rate

Heating rate in stagnation region

Assuming that the maximum and minimum radius of principal curvature at a stagnation point on the three-dimensional surface are $R_1$ and $R_2$ respectively, the formula for calculating the three-dimensional stagnation point heating rate ${\left(q_{\mathrm{sw}}\right)}_{3D}$ is given by ²⁰

$${\left(q_{\mathrm{sw}}\right)}_{3D}=\sqrt{\frac{1+k}{2}}{\left(q_{\mathrm{sw}}\right)}_{2D,\mathrm{Axis}}$$

(20)

where ${\left(q_{\mathrm{sw}}\right)}_{2D,\mathrm{Axis}}$ represents the axisymmetric stagnation point heating rate with $R_1$ as the radius, and $k$ is the three-dimensional correction factor defined as $k=R_1/R_2$ . The Kemp-Riddell method, ¹⁵ widely recognized in engineering, is used to calculate the heating rate at the stagnation point in the form of

$$q_{\mathrm{sw}}=\frac{131884.2}{\sqrt{R_n}}{\left(\frac{{\rho }_{\infty }}{{\rho }_0}\right)}^{0.5}{\left(\frac{v_{\infty }}{v_0}\right)}^{3.25}(1-\frac{h_{\mathrm{sw}}}{h_s})$$

(21)

where ${\rho }_0=1.225\mathrm{kg}/{\mathrm{m}}^3$ , $v_0$  = 7900 m/s. $R_n$ is radius of the stagnation. The wall enthalpy $h_{\mathrm{sw}}$ is calculated according to the following relations ¹² :

$$h_{\mathrm{sw}}=\{\begin{array}{c}{R}^{3.5}{T}_w,T_w\le 170K\hfill \\ \multicolumn{1}{c}{0.796329T_w^{1.041},170.0K<T_w<1748.0K}\\ \multicolumn{1}{c}{78.4187\exp (3.178{\left(\frac{T_w}{1748}\right)}^f),T_w\ge 1748K}\end{array}$$

(22)

$$f=\frac{1.0}{2.41+7.09637\times {10}^{-2}\ln \left(\frac{p_{\mathrm{se}}}{101325}\right)}$$

(23)

In equation (20), the calculation of heating rate requires the radius of curvature at the stagnation point. However, determining the curvature radius at different angles of attack and for various configurations can be challenging for surfaces defined by triangular grids. To address this issue, a quadratic surface is fitted as an approximation to a subset of nodes surrounding the stagnation point. The primary curvature at the stagnation point is estimated by utilizing the primary curvature at the point closest to the stagnation point on the fitted quadratic surface. The quadratic surface is fitted using an equation

$$\begin{array}{c}F(x,y,z)=a_1{x}^2+a_2{y}^2+a_3{z}^2+\hfill \\ \multicolumn{1}{c}{a_4\mathrm{xy}+a_5\mathrm{yz}+a_6\mathrm{xz}+}\\ \multicolumn{1}{c}{a_{7}x+a_{8}y+a_{9}z+a_{10}}\end{array}$$

(24)

where $a_1$ , $a_2$ , $\cdots$ , $a_{10}$ represent fitting coefficients. When the coordinates $(0,0,0)$ lie on the surface being fitted, $a_{10}$ is equal to 0; otherwise, $a_{10}$ is equal to 1. Least squares method is employed to solve for remaining nine fitting coefficients using no fewer than nine coordinate points surrounding the stagnation point.

Kindlmann’s method ²¹ can be utilized to calculate the principal curvatures at a specific point on the quadratic surface determined. To apply this method, begin by determining the Hessian matrix $\mathit{Hess}\mathit{F}$ of $F(x,y,z)$ , which is given by ²²

$$\mathit{Hess}\mathit{F}=\left(\begin{array}{ccc}{F}_{\mathrm{xx}}\hfill & F_{\mathrm{xy}}\hfill & F_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{F_{\mathrm{yx}}}& F_{\mathrm{yy}}\hfill & F_{\mathrm{yz}}\hfill \\ \multicolumn{1}{c}{F_{\mathrm{zx}}}& F_{\mathrm{zy}}\hfill & F_{\mathrm{zz}}\hfill \end{array}\right).$$

(25)

Here, $F_{\mathrm{xx}}$ , $F_{\mathrm{xy}}$ , $\cdots$ , and $F_{\mathrm{zz}}$ represent second-order derivatives of $F(x,y,z)$ . With this information, the mean curvature $H$ and Gaussian curvature $K_G$ on the quadratic surface can be calculated by ²³

$$H=\frac{\mathrm{tr}\left(\mathrm{adj}\left({\mathit{H}}_1\right)\right)-\det \left(\mathit{Hess}\mathit{F}\right)}{2{|\mathbf{\nabla }\mathbf{F}|}^3}$$

(26)

$$K_G=-\frac{\det {\mathit{H}}_1}{{|\mathbf{\nabla }\mathit{F}|}^4}$$

(27)

where ${\mathit{H}}_1$ is defined as

$${\mathit{H}}_1=\left(\begin{array}{cc}\mathit{Hess}\mathit{F}\hfill & \mathbf{\nabla }{\mathbf{F}}^{\mathbf{T}}\hfill \\ \multicolumn{1}{c}{\mathbf{\nabla }\mathbf{F}}& 0\hfill \end{array}\right)$$

(28)

$\mathbf{\nabla }\mathbf{F}$ is the gradient vector, given by

$$\mathbf{\nabla }\mathbf{F}=\left(\begin{array}{ccc}{F}_x\hfill & F_y\hfill & F_z\hfill \end{array}\right)$$

(29)

Then, the primary curvatures on the quadratic surface $k_1$ and $k_2$ can be calculated by ²²

$$\begin{array}{c}{k}_1=H+\sqrt{H^2-K_G}\hfill \\ \multicolumn{1}{c}{k_2=H-\sqrt{H^2-K_G}}\end{array}$$

(30)

Once the heating rate at the stagnation point and in the downstream region is obtained, the heating rate on the elements of the stagnation region, excluding the stagnation element, can be calculated through interpolation.

Heating rate in downstream region

Convective heat transfer calculations in high velocity flows typically depend on Eckert’s reference enthalpy theory. ¹¹ The reference enthalpy $h^{*}$ is usually expressed as

$$h^{*}=\alpha h_e+(1-\alpha )h_w+\beta r\left(\frac{\gamma -1}{2}\right){\mathrm{Ma}}_e^2{h}_e$$

(31)

where $r$ denotes the recovery factor. For laminar flow, $r$ is equal to $\sqrt{\mathrm{Pr}}$ , while for turbulent flow, $r$ is equal to $\sqrt[3]{\mathrm{Pr}}$ . The Prandtl number $\mathrm{Pr}$ , which is 0.72 for air, is also used in the calculation. Ecket’s proposed form of reference enthalpy uses $\alpha =0.5$ and $\beta =0.22$ . Additionally, Meador and Smart ²⁴ presented an alternative reference enthalpy with $\alpha =0.45$ , $\beta =0.16$ for laminar flow and $\alpha =0.5$ , $\beta =0.16$ for turbulent flow. Once the reference enthalpy is determined, using equations (13) and (17), the reference density ${\rho }^{*}$ and reference viscosity ${\mu }^{*}$ can be calculated accordingly. Subsequently, the heating rate $q_w$ in the downstream region can be directly computed for laminar flows as ¹¹

$$q_w=0.332P{r}^{-\frac{2}{3}}{\rho }_e{u}_e^{\sqrt{\frac{{\rho }^{*}{\mu }^{*}}{{\rho }_e{\mu }_e}}}{\mathrm{Re}}_s^{-0.5}(h_r-h_w)$$

(32)

and for turbulent flows

$$q_w=\frac{0.0296{\rho }_{e{u}_e}{\left(\frac{{\rho }^{*}}{{\rho }_e}\right)}^{0.8}{\left(\frac{{\mu }^{*}}{{\mu }_e}\right)}^{0.2}\left(h_{r-h_w}\right)}{P{r}^{\frac{2}{3}}{\mathrm{Re}}_s^{0.2}}$$

(33)

Here, $h_r$ represents the recovery enthalpy given by $h_r=h_e+{\mathrm{ru}}_e^2/2$ , and $h_w$ represents wall enthalpy that can be calculated using equation (24). To distinguish between laminar and turbulent flows in the downstream region, the transition local Reynolds number is given by ²⁵

$$R{e}_{\mathrm{trL}}=1.8\times {10}^5{e}^{0.394M{a}_e}$$

(34)

Inviscid surface streamlines

Determining the local Reynolds number on the aircraft’s surface is typically necessary to calculate convective heat transfer in high enthalpy flows. The local Reynolds number, defined as

$$R{e}_s={\int }_0^s\frac{{\rho }_e{u}_{e}ds}{{\mu }_e},$$

(35)

is calculated by integrating along a given streamline starting from the stagnation point, with $s$ representing the length. Calculating this integral requires determining the coordinates of points along the streamline, as well as the density ${\rho }_e$ , velocity $u_e$ , and viscosity coefficient ${\mu }_e$ at these points. Efficient determination of the streamlines relies on the Newtonian steepest descent concept employed in surface streamline tracing. In accordance with the Newtonian steepest descent concept, the direction of the velocity vector $\overrightarrow{\mathit{u}}/\left|\overrightarrow{\mathit{u}}\right|$ on a surface element can be expressed as ²⁶

$$\frac{\overrightarrow{\mathit{u}}}{\left|\overrightarrow{\mathit{u}}\right|}=\frac{\overrightarrow{\mathit{n}}\times \overrightarrow{{\mathit{V}}_{\infty }}\times \overrightarrow{\mathit{n}}}{\Vert \overrightarrow{\mathit{n}}\times \overrightarrow{{\mathit{V}}_{\infty }}\times \overrightarrow{\mathit{n}}\Vert }$$

(36)

where $\overrightarrow{{\mathit{V}}_{\infty }}$ represents freestream velocity vector and $\overrightarrow{\mathit{n}}$ denotes outward normal direction vector of the surface element. The relationship among these vectors is illustrated Figure 1. Streamline traces are obtained by identifying the points of intersection between the edges of the element and a line constructed based on the velocity direction. The process of tracing streamlines is extensively described in Deters, ²⁶ and Figure 2(a) shows the streamline tracing process on the surface of a double ellipsoid model^27,28 using this method.

Figure 1.

Direction of velocity on the surface element.

Figure 2.

The process of streamlines tracing using two different methods on a double ellipsoid model: (a) points along streamline located on the edges and (b) points along streamline located at the centroid.

However, the process of streamline tracing may encounter special cases where conflicting directions occur at element boundaries and streamline directions at element nodes, ²⁶ as illustrated in Figure 3. This complexity adds challenges to the programming aspect. Moreover, the parameters related to the density ${\rho }_e$ and the viscosity ${\mu }_e$ at the points along these streamlines must be obtained through interpolation, thereby augmenting the computational effort required. In this paper, to avoid the special case depicted in Figure 3 and reduce computational complexity, the points along the streamline are positioned at the centroid of elements. As depicted in Figure 4, for a downstream element A, its upstream elements are chosen from the cells that share at least one node with it. The selection process depends on the angle , which is formed by the line connecting the centroid of A and the centroid of the surrounding elements, and the inverse direction of A’s velocity. The element that minimizes the angle $\alpha$ is selected as A’s closest upstream element. Typically, the starting point of the streamline is the stagnation point, and the process of tracing the streamline can be completed by reaching the stagnation point using the aforementioned method. However, due to the limited number of elements, the streamline tracing may not precisely reach the stagnation point but instead encounters a specific stagnation region, which signifies the completion of the streamline tracing process. Finally, the stagnation point is added to determine the entire streamline. The relationship between the stagnation point and the stagnation region, as well as two different approaches of streamline tracing, is depicted in Figure 2(b) (the element density on the model has been decreased to enhance image clarity). Typically, a blunt body vehicle’s surface flow exhibits only one stagnation point. In this study, the stagnation element is identified as the element with the greatest angle between the free stream vector and the outer normal vector. Its centroid is considered the stagnation point. Additionally, the calculation program established in this study incorporates a variable angle parameter that adjusts the boundary of the stagnation region to accommodate the heat flux calculations for different shaped aircraft. When the angle between the normal vector of a given surface element and that of the stagnation elements is less than this angle, the surface element is included in the stagnation region.

Figure 3.

Special cases in streamline tracing: (a) conflicting directions at element boundaries and (b) streamline direction at element node.

Figure 4.

Determination method of upstream surface elements.

Figure 5 illustrates the CPU time consumed by both the simplified stream tracing method and the traditional stream tracing method for computing streamlines over the surface of a double ellipsoid at a 0-degree angle of attack, using seven different grid configurations. The computations were conducted on an Intel Core i7-9750H CPU utilizing a single thread. The results indicate that as the number of grids increases, both methods experience an increase in computation time; however, the rate of increase for the simplified stream tracing method is significantly lower than that of the traditional method. This results in approximately a 90% enhancement in computational speed compared to the traditional approach, with the improvement becoming more pronounced as the grid resolution increases. Therefore, when the number of grids is substantial and the computational task is heavy, the time savings offered by the simplified stream tracing method can be considerable.

Figure 5.

The performance of the two streamline tracing methods across different grid resolutions.

However, the local Reynolds number computed based on streamlines obtained through this simplified method is susceptible to influence from the surface mesh. Figure 6(a) demonstrates that an inhomogeneous, low-density mesh results in a less smooth computed local Reynolds number. However, increasing the mesh density to reduce non-smoothness comes at the cost of significantly increased computation time, contradicting the purpose of this engineering research.

Figure 6.

Effect of filtering on local Reynolds number of OV model: (a) without filter and (b) with filter.

To smooth the local Reynolds numbers on the surface, a mean filter was implemented in this study. The mean filter is a widely used and straightforward technique in image processing. ²⁹ As shown in Figure 7(a), a 3 × 3 filter is illustrated, which replaces the value at the filter’s center with the average of the values within the 3 × 3 surrounding units, where each unit has a weight of 1/(3 × 3). To apply the mean filter to the surface triangular mesh, this study uses neighboring triangular elements as the region over which the mean filter averages values. This neighborhood can consist of elements that share common vertices with the central element, as illustrated in Figure 7(b). When the mean filter traverses and acts on every element on the surface, a complete mean filtering of the surface is achieved. The degree of smoothing by the mean filter can be adjusted by iterating the filtering process multiple times, modifying the weight of the central element, or altering the size of the filter. As shown in Figure 7(c), the neighborhood of the central element includes two layers of neighboring elements, with the second layer composed of elements sharing common vertices with the first layer. The size of the filter, and thus the degree of smoothing, can be controlled by varying the number of these layers. Additionally, to preserve important features during filtering, the filtering region typically avoids areas with abrupt changes or stationary points on the surface. Figure 6(b) illustrates the Reynolds number contour of the OV model after two of mean filtering using the filter shown in Figure 7(b). It is evident that the contour is significantly smoother than before filtering.

Figure 7.

Different mean filters: (a) the mean filter for image processing, (b) the mean filter applied to triangular surface meshes, and (c) the larger-sized mean filter applied to triangular surface meshes.

Flow of the aerodynamic heating calculation program

The aerodynamic heating calculation program used in this study comprises three main steps: input, calculation, and output, as illustrated in the flowchart in Figure 8. The input is divided into two parts: one consists of the necessary physical calculation conditions and calculation settings for the computations. The physical calculation conditions include the free-stream Mach number $M{a}_{\infty }$ , temperature $T_{\infty }$ , pressure $P_{\infty }$ , density ${\rho }_{\infty }$ , and wall temperature $T_w$ . The calculation settings primarily refer to the selection of the gas model, the Reynolds number filter settings, the division of the stagnation region, and output specifications. This part of the input is written in a specified format and order in a TXT file to facilitate modification and program reading. The second part of the input includes an ASCII-encoded STL mesh file that represents the geometrical shape of the aircraft. Most widely-used CFD mesh generation software can export this simple and universal STL file for the aircraft surface mesh. However, the STL file contains only the vertices of each triangular surface and their normal vectors. Therefore, before initiating calculations, it is necessary to determine the adjacency relationships between surface elements to assist in subsequent streamline tracing while ensuring that the normal vectors of all surface elements point outward from the aircraft.

Figure 8.

Program calculation process.

After completing the surface pressure calculations, the program will open two threads to calculate the gas parameters on the edge of boundary layer and wall heat flux for both the stagnation region and the downstream region, thereby improving computational efficiency. However, in the threading for calculating heat flux in the downstream region, to ensure the completeness of the output data, streamline tracing and local Reynolds number calculations will also be conducted for the stagnation region units, even though local Reynolds number is not required for the stagnation region heat flux calculations. Once the wall heat flux calculations are completed, the program will export the calculated data in the specified format according to the output settings and provide it to post-processing software.

Examples and results

Example 1: Spherically blunted cone

The spherically blunted model and its wind tunnel experiment data are sourced from the NASA TN D-5450 report. ³⁰ The geometric dimensions of model selected for this study, along with the grid used for the calculations, are shown in Figure 9. The Mesh refinement is applied near the nose of model. The calculation conditions for the calculations are listed in Table 1. When calculating the Reynolds number, the mean filter shown in Figure 7(a) was applied iteratively twice to the local Reynolds number for smoothing. The angle parameter for adjusting the stagnation region was set to 1.5°.

Figure 9.

The spherically blunted cone model: (a) geometry of the spherically blunted cone and (b) the surface mesh of the spherically blunted cone.

Table 1.

The calculation conditions for heating prediction of the spherically blunted cone.

Condition	$M{a}_{\infty }$	$P_{\infty }$ ( $\mathrm{Pa}$ )	${\rho }_{\infty }$ ( $\mathrm{kg}/{\mathrm{m}}^3$ )	$T_{\infty }$ ( $\mathrm{K}$ )	$T_w$ ( $\mathrm{K}$ )	AOA
1	10.6	300.8826	0.02426	43.20	294.4	0
2	10.6	300.8826	0.02426	43.20	294.4	5
3	12.1	300.8826	0.02426	43.20	294.4	0
4	13.6	300.8826	0.02426	43.20	294.4	0

Figure 10 illustrates the heat flux along the centerline of different cross-sections at angles $\phi$ , as defined in Figure 9(a), for the blunt cone under conditions 1 and 2. The results indicate that there is a 12% discrepancy between the calculated and experimental values at the stagnation point for both conditions. In the downstream region, the heat flux results for both gas models show a relative error of less than 10%, closely matching the experimental data at a 0-degree angle of attack. However, as shown in Figure 10(b), a significant deviation appears along the centerline of the windward side as the angle of attack increases, with the error in the downstream region growing to over 26%. Figure 10(c) and (d) demonstrate that the heat flux at other positions on the surface aligns well with the experimental data, with an error of less than 18%. The approximate method of heating rate prediction for plates can be used to predict the wall heat flux in the downstream region at small angles of attack.

Figure 10.

The heat flux along the centerline of different cross-sections under conditions 1 and 2: (a) condition 1, AOA 0, $\phi$  =  $0^{°}$ , (b) condition 2, AOA 5, $\phi$  =  ${180}^{°}$ , (c) condition 2, AOA 5, $\phi$  =  $0^{°}$ , and (d) condition 2, AOA 5, $\phi$  =  ${90}^{°}$ .

Figure 11 presents the heat flux contours calculated using a chemical equilibrium gas model. It can be observed that in operating conditions 1 and 2, the wall heat flux in the downstream region is relatively stable, indicating laminar flow. However, as the incoming Mach number increases, conditions 3 and 4 exhibit a abrupt rise in wall heat flux in the downstream area, with the flow transition being the primary cause of this abrupt change. As illustrated, the transition points in conditions 3 and 4 occur at nearly the same location, and the wall heat flux in the turbulent regime is approximately four times greater than that in the laminar regime. Figure 12 shows that the laminar heat flux in condition 4 is approximately 40% higher than in condition 3, while the wall heat flux in the turbulent region is approximately 60% greater. This indicates that as the Mach number increases, the rate of heat flux growth in turbulent flow typically surpasses that in laminar flow.

Figure 11.

Heat flux contours of the spherically blunted cone under different calculating conditions (chemical equilibrium gas model): (a) AOA 0, $\mathrm{Ma}$  = 10.6, (b) AOA 5, $\mathrm{Ma}$  = 10.6, (c) AOA 0, $\mathrm{Ma}$  = 12.1, and (d) AOA 0, $\mathrm{Ma}$  = 13.6.

Figure 12.

The heat flux along the centerline of the windward surface at different Mach numbers.

Figure 13 illustrates the CPU time consumption of two gas models during two calculation processes related to the gas models under varying grid densities. The figure reveals that the primary difference in computation time between the two models occurs in the calculation of gas parameters at the boundary layer edge. The chemical equilibrium gas model requires methods such as bisection or other iterative approaches to determine certain gas parameters, while the perfect gas model allows for direct computation. This difference is the main reason for the observed time discrepancies. Compared to the chemical equilibrium gas model, the perfect gas model can save 70%–75% of the time in computing gas parameters at the boundary layer edge and reduce the time for heat flux calculations in the downstream region by 5%–20%, resulting in an overall time reduction of 65%–70% for these two processes combined.

Figure 13.

Performance of the two gas models in gas model-related calculation processes under different grid resolutions.

Example 2: OV model

The Orbiter Vehicle (OV) model is derived from Greene et al. ³¹ The dimensions of the model and the grid used for calculations are shown in Figure 14, while Table 2 presents the freestream conditions used for calculations. In calculating the local Reynolds number, a mean filter from Figure 7(b) was applied iteratively twice to smooth the Reynolds number, with an angular parameter set to 0.5° to adjust the size of the stagnation region.

Figure 14.

The OV model: (a) the geometry of OV model and (b) the surface mesh of OV model.

Table 2.

The calculation conditions for heating prediction of OV model.

Condition	$M{a}_{\infty }$	$P_{\infty }$ ( $\mathrm{Pa}$ )	${\rho }_{\infty }$ ( $\mathrm{kg}/{\mathrm{m}}^3$ )	$T_{\infty }$ ( $\mathrm{K}$ )	$T_w$ ( $\mathrm{K}$ )	AOA
1	9.84	124.1056	0.0083753	52.0	300.0	8
2	9.84	124.1056	0.0083753	52.0	300.0	16

Figure 15 illustrates the heat flux contours of the OV model under conditions 1 and 2, calculated with two different gas models. The heat flux peak area around the stagnation point is larger when the chemical equilibrium gas model is used than when the perfect gas model is applied. The gas temperature around the stagnation point is generally higher, which may result in air dissociation. Consequently, the perfect gas model might fail to accurately predict the gas properties at the boundary layer edge, leading to a gradient of heat flux along the streamwise direction over relatively smooth surfaces. In contrast, the calculations using the chemical equilibrium gas model appear to better align with the actual physical situation.

Figure 15.

Heat flux contours of the OV model under different calculating conditions: (a) condition 1, perfect gas model, (b) condition 1, chemical equilibrium model, (c) condition 2, perfect gas model, and (d) condition 2, chemical equilibrium model.

Additionally, Figure 15 reveals that aside from the peak heat flux area at the stagnation point, a crescent-shaped peak heat flux region also forms at the edge of the OV model’s nose. The heat flux in this region is comparable to the stagnation heat flux under calculation condition 1. When the angle of attack increases to 18°, the extent of the peak heat flux region slightly decreases, although the wall heat flux experiences a significant increase. The rapid increase in gas velocity along the streamline direction at the model’s nose edge is likely the direct cause of this peak heat flux region. As shown in Figure 16(a), the pressure coefficient contours for calculation condition 1 indicate that gas pressure rapidly decreases along the streamline direction at the windward surface’s edge, suggesting that the gas is undergoing rapid expansion in this region. This leads to substantial increases in gas velocity, as illustrated in Figure 16(b). Although the gas enthalpy and density decrease, they remain relatively high compared to the leeward surface’s nose edge due to proximity to the stagnation point, as shown in Figure 16(c) and (d). Gas acceleration under higher enthalpy and density conditions contributes more significantly to the wall heat flux than under lower enthalpy and density conditions. Therefore, the peak heat flux occurs at the windward surface’s nose edge, but not at the corresponding edge on the leeward side.

Figure 16.

Gas parameters at the edge of the boundary layer on the OV model under condition 1: (a) pressure coefficient contours, (b) the velocity at the edge of boundary layer along the centerline, (c) the density at the edge of boundary layer along the centerline, and (d) the enthalpy at the edge of boundary layer along the centerline.

Figure 17 presents the dimensionless heat flux along the centerline of the windward surface under conditions 1 and 2, with experimental data sourced from Greene et al. ³¹ It is clear that the heat flux in the cylindrical section of the model closely matches the experimental measurements for both conditions, with a relative error of less than 15%. However, significant discrepancies are observed in the downstream expansion section of the OV model under both conditions. For condition 2, Greene et al. ³¹ indicates that flow transitions occurred in the downstream expansion section; however, the transition criterion used in this study failed to accurately predict the transition, which may be a primary reason for the observed errors.

Figure 17.

Dimensionless heat flux along the centerline of the windward side of the OV model: (a) condition 1 and (b) condition 2.

Additionally, in the surface extension section of the downstream region, a substantial amount of inviscid flow mass is entrained into the boundary layer. The surface streamlines originating from the stagnation region are engulfed by the boundary layer, leading to a disparity between the entropy at the boundary layer edge and the entropy at the stagnation point. This disparity could be one potential explanation for the observed discrepancies, as the computational method assumed in this paper may not be applicable in such scenarios.

Conclusion

This paper provides a detailed and comprehensive introduction to the aerodynamic heating calculation engineering method for hypersonic blunt body vehicles. By simplifying the streamline tracing approach, the paper enhances computational efficiency while validating the effectiveness of the engineering method through two examples. The conclusions drawn from this research are as follows:

1. The simplified streamline tracing method constructed in this study effectively enhances computational efficiency, achieving a 90% reduction in the time required for streamline calculations compared to traditional methods.

2. Validation results using a spherically blunted cone and OV models indicate that the method effectively predicts wall heat flux for blunt-body vehicles at small angles of attack. In the stagnation and non-expanded downstream regions, the method achieves a relative error of less than 15%.

3. In downstream expansion regions and at high angles of attack, the method shows larger errors in wall heat flux prediction. Incorporating a more accurate flow transition criterion and applying entropy corrections may help reduce these errors.

4. While the difference between results obtained using a perfect gas model and those from a chemical equilibrium gas model is generally small, using the perfect gas model can reduce computational time by 65%–70%. For most conditions, the perfect gas model is appropriate; however, for high-enthalpy flows and regions near stagnation points, the chemical equilibrium gas model provides more reasonable results.

In conclusion, the method presented in this paper strikes a practical balance between accuracy and efficiency, making it well-suited for the engineering design of hypersonic vehicles, particularly during the early design stages.