Efficient Method for Heat Flux Calculations within Multidisciplinary Analyses of Hypersonic Vehicles

Abstract

A large amount of heat flux from aerodynamic heating acts on reusable spacecraft; thus, an accurate heat flux prediction around spacecraft reentry is essential for developing a high-performance reusable spacecraft. Although the approximate convective-heating equations can calculate the heat flux with high efficiency and sufficient fidelity, the heat flux should be evaluated over a thousand times for the entire trajectory in multidisciplinary analyses. For these reasons, it is necessary to develop an efficient method for calculating the heat flux for multidisciplinary analysis. In this paper, an efficient method for heat flux calculation that is adoptable by multidisciplinary analyses for hypersonic vehicles, such as spacecraft, is developed. Approximate convective-heating equations were adopted to relieve the computational cost of estimating the heat flux, and an adaptive time step method for heat flux calculations was developed to reduce the number of heat flux calculations required across the entire flight trajectory. A dynamic factor was introduced that adjusts the time step between each instance of the heat flux calculation. Since the time step using this factor could increase under low heat flux conditions, the number of heat flux calculations decreases by approximately one-tenth with over 90% accuracy. Therefore, the efficiency was improved with high accuracy using the adaptively-determined time step according to this dynamic factor.

1. Introduction

Currently, private enterprises, such as Space Exploration Technologies Corp. (SpaceX), Virgin Galactic, and Blue Origin, are developing launch vehicles and spacecraft. SpaceX is developing launch vehicles and spacecraft to reduce space transportation costs and colonize Mars. Virgin Galactic is developing commercial spacecraft for space tourists, and the company has succeeded in spaceflight with a full crew that included the founder of the company. Moreover, Blue Origin focuses on developing rocket-powered vertical takeoff and landing for access to suborbital and orbital space. These operations imply that space development has shifted from governmental to commercial areas, based on increasing profits through reduced expenditures by employing reusable spacecraft and diversifying the purposes of the flights.

Various technologies are necessary to develop a high-performance reusable spacecraft. One of the most essential technologies involves accurate heat flux prediction during spacecraft reentry because a large amount of heat flux is generated by aerodynamic heating. Severe aerodynamic heating throughout reentry may destroy or damage a spacecraft. To ensure that a spacecraft can withstand this severe aerodynamic heating, the heat flux should be accurately predicted in advance.

Various methods have been proposed to predict heat fluxes including high-fidelity methods, such as the Reynolds-averaged Navier–Stokes (RANS) equation [1,2,3] and large eddy simulation (LES) [4]; and low-fidelity methods, such as the one-equation method [5], which estimates the entire heat flux on a vehicle with one equation, the flat plate equation [6], which is considered at the same velocity and incidence angle of the surface, and a correlation the friction coefficient [7] and Stanton number [8]. The RANS equation or LES can predict the heat flux on a spacecraft with high accuracy. However, their computational cost is too high since these simulations can require tens of hours. On the other hand, the computational cost of the one-equation method is lower. However, this method cannot estimate the heat flux at a specific point on the surface of a spacecraft because one equation is used to predict the heat flux on the entire vehicle. The approximate convective-heating equations proposed by E. V. Zoby et al. [9] can calculate the heat flux at a specific point on the surface with the sufficient accuracy and high efficiency using an axisymmetric analogy. These equations utilize a modified Reynolds analogy form based on ${Re}_{\theta }$ to relate the heat transfer to skin friction. Originally, these equations were developed to predict the heat flux on an axisymmetric body based on an axisymmetric analogy. The approximate convective-heating equations have been extended to three-dimensional problems by adopting streamline coordinates and metrics. [10] Although these equations are more efficient than high-fidelity methods and have sufficient accuracies for various shapes [10,11,12], they require over one hour of calculation. The reason for this phenomenon is that they are based on the inviscid flow field properties determined through the finite volume method, which is insufficient for application to an entire trajectory. Alternatively, J. S. Zhao et al. [13] combined the equations with Newtonian theory and linear fitting, reducing the computational time to several seconds.

Despite the significant improvement in the computational efficiency due to the aforementioned approaches, challenges remain in the estimation of heat fluxes. First, heat fluxes should be evaluated along the entire trajectory of a vehicle, resulting in over a thousand evaluation stages. Second, heat flux evaluations along the trajectory are usually simulated over a thousand times during vehicle shape optimization [14,15,16,17]. Thus, heat flux evaluations are performed over several million times in a multidisciplinary optimization process. For these reasons, although the computation time for the heat flux evaluation at each instant is very short, an enormous number of heat flux evaluations are required for multidisciplinary optimization. This process requires excessive computation time. Therefore, it is necessary to develop an efficient method for calculating the heat flux along the entire trajectory for multidisciplinary optimization.

A practical solution for reducing the computational cost is increasing the trajectory-analysis time step. Since the heat flux is calculated along the entire trajectory, decreasing the number of calculations for the trajectory by increasing the time step relieves the computational cost. Several adaptive time-step methods for time-dependent problems have been proposed, such as a method using a simple expression [18,19,20] and a method using a proportional-integral (PI) [21,22,23,24,25,26] or proportional-integral-derivative (PID) controller [27,28,29]. These methods control the time step based on the error (the rate of change in the solutions). Thus, these methods are suitable for steady-state problems because the error decreases along the time integral due to the convergence of solutions. These methods are not appropriate for trajectory analysis because of the substantial changes in the flow conditions during the analysis.

Another solution is splitting the time step for multidisciplinary analysis, including the heat flux evaluation along a trajectory, into individual time steps for individual analyses. The required time steps of the individual analyses for accurate computational processes are different. If the same time step is applied for all analyses in multidisciplinary analysis, the minimum required time step is selected for accurate results. However, if the time step differs for each analysis, the number of calculations decreases. As a result, the computational cost is reduced. Although splitting a time step into individual time steps has been attempted [1,30,31,32], the ratio of the time steps between individual analyses was constant. If the ratio is constant, the computational cost may not decrease drastically because each individual time step is not controllable using the change in the conditions for each analysis. Therefore, a method that adjusts each individual time step is needed.

In this paper, an efficient method for heat flux calculations that is adoptable in the multidisciplinary analysis of hypersonic vehicles is developed. To achieve this goal, approximate convective-heating equations are adopted to diminish the computational cost of computing the heat flux. In addition, to reduce the number of heat flux calculations along a trajectory, an adaptive time-step method is developed. A dynamic factor that adjusts the time step between each heat flux calculation is introduced. The time step adjusted by the factor decreases when a large amount of heat flux is generated. Moreover, the time step increases when a small amount of heat flux is generated. The time step of heat flux calculations with the developed dynamic factor is adaptively determined to increase the efficiency of multidisciplinary analysis, while detailed information on the heat fluxes in the high heat flux conditions is obtained.

2. Heat Flux Calculation

In this paper, inviscid flow on the vehicle surface was estimated using modified Newtonian theory [33] to calculate the heat flux on the surface. This theory can be used to compute the surface pressure. The surface velocity magnitude was obtained using the assumption of transforming the normal velocity on the surface into pressure. In addition, the direction of surface velocity was calculated as follows:

$$\frac{\mathit{V}}{\left|\mathit{V}\right|}=\frac{\mathit{n}\times {\mathit{V}}_{\mathrm{\infty }}\times \mathit{n}}{\left|\mathit{n}\times {\mathit{V}}_{\mathrm{\infty }}\times \mathit{n}\right|},$$

(1)

where $\mathit{V}={\left(u,v,w\right)}^{\mathrm{T}}$, ${\mathit{V}}_{\mathit{\infty }}$, and $\mathit{n}$ are the surface velocity, freestream velocity, and normal vector on the surface, respectively. Based on these flow properties, the heat flux on the surface was estimated using the approximate convective-heating equation [10,11,13].

2.1. Streamline Calculation

Since the approximate convective-heating equations compute the heat flux along a streamline, inviscid surface streamlines are needed. Regarding the streamline calculations, two integration methods exist: forward integration, which integrates the streamline segments from the stagnation point along the surface velocity direction, and backward integration, which integrates the segments from seeds (starting points of streamline calculation) to the opposite direction to that of the surface velocity. The results of both integrations are similar for simple geometries, such as a sphere, ellipsoid, and spherically blunted cone. However, the streamline using forward integration is not well distributed on complex geometries, such as wing–body configurations, because the differences in the well-distributed streamlines are overly small near the stagnation point. Thus, backward integration was adopted in this paper.

To calculate the streamline, a candidate point is chosen at a short distance from the endpoint of the streamline (i.e., the seed point) as the inverse direction of the surface velocity using Equation (1). If the candidate point falls within the same grid as the endpoint [34], it becomes the new streamline endpoint. Otherwise, the candidate point is projected onto the plane that includes other grids, and the projection is checked to determine if it falls within the same grid. If the projected point is within the grid, it becomes the new streamline endpoint. This process is repeated from the seed point to the stagnation point.

The seeds for calculating the streamline are selected from the midpoints of all grids that include the trailing edge of the wing or the aft of the body to distribute the streamline over the wing and body. As illustrated in Figure 1, this approach results in a well-distributed streamline on the wing and body.

2.2. Streamline Metrics

On the surface, the inviscid streamline coordinates $\xi$ and $\beta$ were employed as described in Figure 2. $\xi$ is a coordinate along the streamline, and $\beta$ is a coordinate perpendicular to the streamline [13].

The differentials of the arc lengths on the surface were $ds=h_{s}d\xi$ and $\mathrm{d}{s}_{\perp }=h\mathrm{d}\beta$, and the differential of the position vector ($\mathrm{d}\mathit{R})$ on the surface was written as follows:

$$\mathrm{d}\mathit{R}=h_{\mathrm{s}}\mathrm{d}\xi \widehat{{\mathit{e}}_{\mathrm{s}}}+h\mathrm{d}\beta \widehat{{\mathit{e}}_{\perp }},$$

(2)

where $h_{\mathrm{s}}$ and $h$ are the metric coefficients associated with $\xi$ and $\beta$, respectively; and $\widehat{{\mathit{e}}_{\mathrm{s}}}$ and $\widehat{{\mathit{e}}_{\perp }}$ are the unit vectors in the same direction and the perpendicular direction of the streamline on the surface, respectively. $h$ indicates the convergence or divergence of streamlines, and it is identical to the local radius of the body in an axisymmetric flow.

In the application of the axisymmetric analog [35], the streamline metric $h$ should be calculated to substitute it for the radius of an equivalent axisymmetric body. An efficient method was adopted that calculates the streamline metric by using only two independent variables in Cartesian coordinates ($x, y, z$) for the integration. If ($y, z$) were selected, $v=v(y, z)$, $w=w(y, z)$, and the streamline metric was calculated as follows [11]:

$$h=\frac{\left|\nabla \mathit{F}\right|}{{\mathit{F}}_x\left|\mathit{V}\right|}\left[w{\left(\frac{\partial y}{\partial \beta }\right)}_z\right],$$

(3)

where $\mathit{F}$ is the body surface described by $\mathit{F}\left(x, y, z\right)=0$ and $\nabla \mathit{F}={\left(F_x,F_y,F_z\right)}^{\mathrm{T}}$ is the normal vector to the surface. The differential equation of the partial derivative is as follows:

$$\frac{\mathrm{d}}{\mathrm{d}s}\left[w{\left(\frac{\partial y}{\partial \beta }\right)}_z\right]=\frac{1}{\left|\mathit{V}\right|}\left[w{\left(\frac{\partial y}{\partial \beta }\right)}_z\right]\left[\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right],$$

(4)

where $s$ is the surface distance along the streamline. The equations for using the independent variables $\left(x, y\right)$ and $(x, z)$ were derived in Reference [4]. Since the denominator in Equation (3) had a component of $\nabla \mathit{F}$, the variables of integration were selected to maximize the component in the denominator.

2.3. Heating Equations

The approximate convective-heating equations for calculating the heat flux on the surface were proposed by Zoby et al. [9]. Above an altitude of 50 km, the flow condition is predominantly laminar [36]. Since reentry vehicles spend the majority of their flight time above 50 km, the flow is mostly laminar [37,38]. Therefore, under these flight conditions, the flow is assumed to be laminar. For laminar flow,

$${\dot{q}}_{\mathrm{w},\mathrm{L}}=0.22{\left(R{e}_{\theta ,\mathrm{e}}\right)}^{-1}\left(\frac{{\rho }^{*}}{{\rho }_{\mathrm{e}}}\right)\left(\frac{{\mu }^{*}}{{\mu }_{\mathrm{e}}}\right){\rho }_{\mathrm{e}}{u}_{\mathrm{e}}\left(H_{\mathrm{a}\mathrm{w}}-H_{\mathrm{w}}\right){\left(P{r}_{\mathrm{w}}\right)}^{-0.6},$$

(5)

where ${\dot{q}}_{\mathrm{w},\mathrm{L}}$ is the heat flux on the surface; $R{e}_{\theta }$ and ${Pr}_{\mathrm{w}}$ are the momentum thickness Reynolds number and Prandtl number, respectively; $\rho$, $\mu$, $u$, and $H$ are the density, viscosity, velocity, and enthalpy, respectively; subscripts (_e), (_w), and (_aw) indicate estimations at the boundary layer edge, wall, and adiabatic wall, respectively; and superscript (*) indicates the evaluations by Eckert’s reference enthalpy relation [39] to consider compressible effects.

The laminar momentum thickness ${\theta }_{\mathrm{L}}$ was computed as follows:

$${\dot{q}}_{\mathrm{w},\mathrm{L}}=0.22{\left(R{e}_{\theta ,\mathrm{e}}\right)}^{-1}\left(\frac{{\rho }^{*}}{{\rho }_{\mathrm{e}}}\right)\left(\frac{{\mu }^{*}}{{\mu }_{\mathrm{e}}}\right){\rho }_{\mathrm{e}}{u}_{\mathrm{e}}\left(H_{\mathrm{a}\mathrm{w}}-H_{\mathrm{w}}\right){\left(P{r}_{\mathrm{w}}\right)}^{-0.6},$$

(6)

The following equation from Kemp et al. [40] was adopted to correct the momentum thickness to account for the approximate effect of the velocity gradient on laminar heating:

$${\overline{\theta }}_{\mathrm{L}}=\frac{{\theta }_{\mathrm{L}}}{1+0.09\sqrt{\overline{\beta }}}$$

(7)

Thus, ${\overline{\theta }}_{\mathrm{L}}$ in Equation (7) replaces ${\theta }_{\mathrm{L}}$ in Equation (6) before the calculation of laminar heating. $\overline{\beta }$ is a velocity gradient parameter defined as follows:

$$\overline{\beta }=\frac{2\zeta }{u_{\mathrm{e}}}\left[\left(\frac{\mathrm{d}{u}_{\mathrm{e}}}{\mathrm{d}s}\right)/\left(\frac{\mathrm{d}\zeta }{\mathrm{d}s}\right)\right],$$

(8)

where $\zeta$ is the parameter determined by the Lees-Dorodnisyn transformation [40]. This parameter can be expressed as follows:

$$\mathsf{\zeta }={\int }_0^s{\rho }_w{\mu }_w{u}_{\mathrm{e}}{h}^2\mathrm{d}s$$

(9)

2.4. Heat Flux near the Stagnation Point

The aforementioned approximate convective-heating equations could not compute the heat flux at the stagnation point because integration began at the stagnation point. The stagnation-point heat flux was calculated as follows [41]:

$${\dot{q}}_{\mathrm{w},\mathrm{s}}=0.767\sqrt{\frac{\mathrm{d}{u}_{\mathrm{e}}}{\mathrm{d}x}}{\left(\rho \mu \right)}^{0.5}\left(H_{\mathrm{s}}-H_{\mathrm{w}}\right){\left(P{r}_{\mathrm{w}}\right)}^{-0.6},$$

(10)

where subscript s indicates the properties at the stagnation point and $\mathrm{d}{u}_e/\mathrm{d}x$ is the gradient of velocity at the stagnation point. DeJarnette et al. proposed the following equation to calculate the gradient [42]:

$$\frac{\mathrm{d}{u}_{\mathrm{e}}}{\mathrm{d}x}=\frac{\left|{\mathit{V}}_{\infty }\right|}{R}\sqrt{1.85\frac{{\rho }_{\infty }}{{\rho }_{\mathrm{s}}}}$$

(11)

If the approximate convective-heating equations in the previous section were adopted for calculating the heat flux, a physically unrealistic pattern would appear near the stagnation point because of the velocity singularity at this point. To address this problem, a surface curve, named the ε-curve, perpendicular to the inviscid surface streamlines and surrounding the stagnation point was used [11].

The approximation of the similar $\left({\rho }^{*}{\mu }^{*}\right)$ and the linear $u_{\mathrm{e}}$ and $h$ along a streamline could be applied near the stagnation point. Based on this approximation, the integration for computing the momentum thickness ${\theta }_{\mathrm{L}}$ on the ε-curve was replaced as follows [11]:

$${\int }_0^{s_{\epsilon }}{\rho }^{*}{\mu }^{*}{u}_{\mathrm{e}}{h}^2\mathrm{d}s={\left({\rho }^{*}{\mu }^{*}\right)}_s{\left(u_e\right)}_{\epsilon }{h}_{\epsilon }^2\frac{S_{\epsilon }}{4}$$

(12)

The heat flux on the ε-curve was calculated without integration because the momentum thickness on the ε-curve was computed using Equation (12). The heat flux inside the ε-curve was computed by interpolating the heat flux at the stagnation point and the ε-curve.

2.5. Heat Flux Calculation Procedure

The first step in estimating the heat flux using approximate convective-heating equations is to calculate the streamline. Next, the heat flux can be predicted at each point along the streamline from the stagnation point to the endpoint (seed) of the streamline. If the calculation point is inside the ε-curve, the heat flux is obtained by interpolating the heat flux at the stagnation point using Equations (10) and (11) and on the ε-curve using Equation (12). If the calculation point is outside the ε-curve, the boundary layer edge properties are calculated based on the modified Newtonian theory. Then, the streamline metric and momentum thickness are computed from these properties using Equations (3) and (4) and (6)–(9). Finally, the heat flux at this point is estimated from the computed momentum thickness using Equation (5). This procedure is illustrated in Figure 3.

2.6. Heat Flux Calculation Validation

To validate the approximate convective-heating equations used in this paper, the heat flux on a sphere with a 0.0508 m radius was compared with experimental results [43] and a high-fidelity computational fluid dynamics (CFDs) method using the RANS solver. In this paper, ANSYS Fluent, a commercial CFD software package, was employed to evaluate the RANS equations with steady-state and laminar flows. The RANS solver used in this paper adopted the implicit AUSM+ flux [44], Green–Gauss node-based gradient [45], and blending of the first- and second-order upwind schemes. The freestream and wall conditions for the sphere are described in

Table 1. Freestream and wall conditions of the sphere.

	Value
$M_{\infty }$	9.74
$\gamma$	1.4 (perfect gas)
${\rho }_{\infty }$	0.004272 kg/m3
$T_{\infty }$	53.17 K
$T_{\mathrm{w}}$	300.0 K

. The heat flux that was obtained by using the approximate convective-heating equations was in good agreement with those found from the experiment and RANS, as plotted in Figure 4.

Additionally, the heat flux on an 8° sphere cone with a 0.06349 m radius at α = 0° was compared with the experimental result [43] and RANS result. The setting of RANS solver is the same as the calculation of the sphere. The freestream and wall conditions for the sphere are described in

Table 2. Freestream and wall conditions of the 8° sphere cone α = 0°.

	Value
$M_{\infty }$	5.00
$\gamma$	1.4 (perfect gas)
${\rho }_{\infty }$	0.0500 kg/m3
$T_{\infty }$	115.28 K
$T_{\mathrm{w}}$	300.0 K

. The heat flux obtained by using the approximate convective-heating equations was also in good agreement with those found from the experiment and RANS, as plotted in Figure 5.

M. Brchnelova et al. compared the heat flux on the space shuttle geometry between the calculation using the approximate convective-heating equations with modified Newtonian theory and an experiment [37]. The result shows the sufficient accuracy of the approximate convective-heating equations.

Furthermore, the heat fluxes on the wing–body configuration using the RANS and the approximate convective-heating equations were calculated for comparison with more complex geometry. As shown in Figure 6, Korea Aerospace Research Institute’s KSP-1, which is a three-ton class vehicle with a 7 m fuselage and 4 m span wing, was adopted. The freestream and wall conditions of KSP-1 are described in

Table 3. Freestream and wall conditions of KSP-1.

	Case 1	Case 2	Case 3	Case 4
$M_{\infty }$	21.39	16.68	12.20	8.01
$\gamma$	1.4	1.4	1.4	1.4
$\alpha$	40°	40°	40°	40°
${\rho }_{\infty }$	0.000439 kg/m3	0.000508 kg/m3	0.000645 kg/m3	0.001302 kg/m3
$T_{\infty }$	254.83 K	258.16 K	263.76 K	270.65 K
$T_{\mathrm{w}}$	300.0 K	300.0 K	300.0 K	300.0 K

. The computational mesh used for the RANS calculation is depicted in Figure 7. An unstructured hybrid mesh with 50 prism layers was utilized, and the number of node points was approximately 10,000,000. To calculate the heat flux using the approximate convective-heating equations, the body surface was discretized into 80 grids in the longitudinal direction and 58 grids to represent a cross-section. The wing surface was discretized into 20 grids in the span and 128 grids to represent the airfoil.

The results of the RANS and approximate convective-heating equation, which are shown in Figure 8, are similar to each other. For a detailed comparison, the heat fluxes along the wing section of 1.85 m are plotted and displayed in Figure 9. The peaks of the heat fluxes are in good agreement, and the trends of the heat fluxes with variations in the x-coordinates are similar. From the sphere and KSP-1 results, the approximate convective-heating equations used in this paper showed sufficient accuracy for calculating the heat flux.

3. Adaptive Time Step for the Heat Flux Calculation

The approximate convective-heating equations described in the previous section were coupled with trajectory analysis to evaluate the heat flux of the spacecraft. Moreover, the analysis was combined with weight, propulsion, and aerodynamic analyses [14]. In this paper, the KSP-1 trajectory was analyzed from mission orbit to landing using three-degree-of-freedom equations. The trajectory started at an altitude of 300 km with a speed of 7000 m/s, a flight path angle of 0°, and an incline angle of 80°. Along this trajectory, over one thousand heat flux evaluations should be performed to estimate the heat flux during flight.

To reduce the number of heat flux calculations, the time step of the calculations was increased by multiplying it by a factor $C$ as follows:

$$∆t_{\mathrm{h}}=C∆t_{\mathrm{t}},$$

(13)

where $∆t_{\mathrm{t}}$ and $∆t_{\mathrm{h}}$ are the time steps of the trajectory and heat flux calculation, respectively. $∆t_{\mathrm{t}}$ was identical to $∆t_{\mathrm{h}}$ if the heat flux was computed at all positions in the trajectory. If $C$ increased, the number of heat flux calculations decreased because the time step $∆t_{\mathrm{h}}$ increased. However, if $C$ decreased, the number of heat flux calculations increased because the time step $∆t_{\mathrm{h}}$ decreased.

3.1. Constant $C$

First, the constant $C$ was determined numerically. According to

Table 4. Total computational cost with the constant $C$.

$$\mathit{C}$$	Number of Computations	Maximum Stagnation Heat Flux (W/m2)	Total Computational Cost (s)
1 (original)	1197	2,133,778	1553.9525
2	600	2,133,778	832.1355
7	174	1,751,234	272.4919
11	111	1,700,364	211.6848
17	72	1,875,972	137.2899
23	54	634,614	114.4419
31	39	1,513,749	103.5552

, the total computational cost decreased with an increasing $C$ because a large $\Delta t_h$ reduced the number of computations. However, when using the constant C, the time at the original maximum stagnation heat flux was skipped due to the increased $\Delta t_h$, resulting in inaccurate maximum stagnation heat fluxes, as observed in Table 4 and Figure 10. The maximum heat flux was particularly important because the spacecraft should withstand the maximum value.

3.2. Dynamic $C$

In the previous section, if a constant value of $C$ was used, the heat flux calculation could not be performed for the maximum heat flux position in the trajectory. To address this problem, it is necessary to use a small $\Delta t_h$ at a high heat flux for accurate calculations, and a large $\Delta t_h$ at a low heat flux for efficient calculations. Thus, a dynamic $C$ was adopted using the heat flux at the stagnation point of ${\dot{q}}_{\mathrm{s}}$ as the control value of $C$ because the estimation of the heat flux at this point did not require a time-consuming streamline calculation. Furthermore, a stagnation heat flux could represent the characteristics of heat fluxes in certain conditions because the changes in the heat fluxes at the stagnation point could express the changes in the general heat fluxes. Notably, the stagnation heat flux may not be at its maximum because of the shape of the shock wave and the laminar–turbulent flow transition. Using the heat flux at the stagnation point ${\dot{q}}_{\mathrm{s}}$, the dynamic factor $C$ was defined as follows:

$$C=\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\left[{\left(\frac{{\dot{q}}_{\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}-{\dot{q}}_{\mathrm{s}}}{{\dot{q}}_{\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}}\right)}^b\times a\right]+1,$$

(14)

where $a$ and $b$ are the parameters that control the time step of the heat flux evaluation along the trajectory, and ${\dot{q}}_{\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum heat flux at the stagnation point throughout the trajectory. The dynamic factor $C$ in Equation (14) was determined based on the difference between the heat flux at an instant (${\dot{q}}_{\mathrm{s}}$) and the maximum heat flux at the stagnation point throughout the trajectory (${\dot{q}}_{\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}$). The maximum heat flux is an important property that needs to be analyzed throughout the trajectory of hypersonic vehicles to evaluate the possibility of destruction and damage. According to Equation (14), the time step decreased when the heat flux at the stagnation point was high, and the time step increased when the heat flux was low. This phenomenon indicated that the efficiency of the heat flux calculation within multidisciplinary analysis was improved with an elongated time step under low heat flux conditions. In addition, detailed information on the heat flux was obtained by shortening the time step under high heat flux conditions because periods with high heat flux are important due to the possibility of destruction and damage to the vehicles. The maximum time step increased with a larger value of $a$ because the maximum value of $C$ increased. If $b$ increased, the time step decreased at a low stagnation heat flux because the value of $C$ was low, as described in Figure 11.

To find the characteristics of $a$ and $b$ in Equation (14), the calculations were conducted a hundred times with $a$ from 10 to 100 at intervals of 10 and $b$ from 1 to 10 at intervals of 1. For brevity, only several results are presented in

Table 5. Total computational cost with the dynamic $C$.

$$\mathit{a}$$	$$\mathit{b}$$	Number of Calculations	Maximum Stagnation Heat Flux (W/m2)	Total Computational Cost (s)
Original		1197	2,133,778	1553.9525
20	1	75	2,133,778	125.2944
10	2	156	2,133,778	215.9298
20	2	93	2,133,778	139.2946
30	2	72	2,133,778	103.8202
40	2	57	2,133,778	81.8680
20	3	108	2,133,778	132.3144
30	3	81	2,133,778	98.2031
20	4	123	2,133,778	128.7518
30	4	96	2,133,778	85.8918

and Figure 12. The total computational cost was dramatically reduced by approximately one-tenth with Equation (14), as seen in Table 5. Furthermore, the maximum heat fluxes were identical to the original flux, as shown in Table 5 and Figure 12, because the time step at the high stagnation heat flux was reduced. The number of calculations decreased by increasing $a$ with the same $b$ because the maximum value of $C$ became large. In contrast, the number of calculations increased by increasing $b$ with the same $a$ because a low value of $C$ was maintained at a relatively low stagnation heat flux, as shown Figure 11.

To analyze the performance of Equation (4) according to $a$ and $b$, the efficiency and accuracy were defined as follows:

$$(\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{y})=\frac{\left(\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{calculation}\mathrm{with}\mathrm{original}\right)-\left(\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{calculation}\mathrm{using}\mathrm{Equation}\left(14\right)\right)}{\left(\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{calculation}\mathrm{with}\mathrm{original}\right)}$$

(15)

$$(\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y})=1-\frac{{\int }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{y}}^{}\left|\left(\mathrm{heat}\mathrm{flux}\mathrm{of}\mathrm{original}\right)-\left(\mathrm{heat}\mathrm{flux}\mathrm{of}\mathrm{using}\mathrm{Equation}\left(14\right)\right)\right|d\left(Mach\right)}{{\int }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{y}}^{}\left(\mathrm{heat}\mathrm{flux}\mathrm{of}\mathrm{original}\right)d\left(Mach\right)}$$

(16)

The efficiency indicates the reduction rate of the number of calculations by using dynamic $C$. The accuracy represents how close the heat flux with dynamic $C$ to the heat flux without dynamic $C$ is.

If $a$ decreased, the efficiency decreased because $C$ was small; if $a$ increased, the efficiency increased because $C$ was large. If $b$ decreased, the accuracy decreased because $C$ became constant; if $b$ increased, the accuracy increased because $C$ trended toward one. Figure 13 shows a plot of the efficiency and accuracy according to $a$ and $b$. The shadow region in Figure 13 is the zone with efficiencies and accuracies higher than 0.9. In

Table 6. Total computational cost to achieve high efficiency and accuracy.

$\mathit{a}$	$\mathit{b}$	Number of Calculations	Efficiency	Maximum Stagnation Heat Flux (W/m2)	Accuracy
40	2	57	0.95238	2,133,778	0.91737
20	1	75	0.93042	2,133,778	0.93484
20	3	108	0.90727	2,133,778	0.95610

, the three samples in the region (indicated by bold squares in Figure 13) are described. When $a$ and $b$ are 40 and 2, respectively, the efficiency is high, and the accuracy is relatively low. However, the accuracy is high with $a=20$ and $b=3$ despite the relatively low efficiency. The efficiency and accuracy of the other are compromised.

4. Conclusions

In this paper, an efficient method for heat flux calculations that could be adopted in multidisciplinary analyses of hypersonic vehicles was developed. The approximate convective-heating equations were adopted to diminish the computational cost of estimating the heat flux. Furthermore, an adaptive time step method for calculating the heat flux of spacecraft along a trajectory was developed to reduce the number of heat flux evaluations required along the entire flight trajectory. A dynamic factor that adjusts the time step between each heat flux evaluation was introduced. The dynamic factor was determined based on the difference between the heat flux at an instant and the maximum heat flux over the entire trajectory. With the developed adaptive time step method, detailed information on the heat flux was obtained by shortening the time step when the heat flux was high. The efficiency of the heat flux calculations for the entire trajectory was improved by increasing the time step under low heat flux conditions. As a result, the number of heat flux calculations was reduced from 1197 to approximately 100, while the same maximum heat flux is estimated. In other words, this method decreased the number of heat flux calculations by approximately 1/10 with over 90% accuracy. If the time step was adaptively determined according to this dynamic factor, not only could the efficiency be improved, but the accuracy could also be high. Therefore, this method was used to strengthen a multidisciplinary analysis by improving the efficiency with high accuracy. In particular, this technique can facilitate the rapid and iterative design of a spacecraft that can safely reenter the atmosphere.