Numerical Study of High-Temperature Nonequilibrium Flow around Reentry
Vehicle Coupled with Thermal Radiation

Accurate aerodynamic heating prediction is of great significance to current manned space flight and deep space exploration missions. The temperature in the shock layer surrounding the reentry vehicle can reach up to 10,000 K and result in remarkable thermochemical nonequilibrium, as well as considerable radiative heat transfer. In general, high-temperature flow simulations coupled with thermal radiation require appropriate numerical schemes and physical models. In this paper, the equations governing hypersonic nonequilibrium flow, based on a three-temperature model combined with a thermal radiation solving approach, are used to investigate the radiation effects in the reentry shock layer. An axisymmetric spherical case shows that coupling the flow-field simulation with radiation has a scarce influence on the convective heating prediction, but has some impact on the radiative heating calculation. In particular, for the Apollo capsule reentry, both the absorption coefficient and incident radiation are remarkable inside the shock layer. The radiative heating maximum reaches nearly 38% of that of the convective heating making a considerable contribution to the total aerodynamic heating. These results indicate that in the hypersonic regime, in order to account for the total heating, it is necessary to simulate the high-temperature thermochemical nonequilibrium flows coupled with thermal radiation.


Introduction
When a reentry vehicle enters the Earth atmosphere at hypervelocity (about 10 km/s), the surrounding high-temperature shock layer resulting in strong aerodynamic heating seriously threatens the safety of the vehicle. The aerodynamic thermal environment prediction and thermal protection design have become the key technical barriers that have to face in current manned space flight and deep space exploration missions [1]. the stagnation point of FIRE II capsule using HARA code combined with the viscous shock layer technique, finding that the radiation heat flux obtained from the flow simulation coupled with thermal radiation was reduced by about 30% compared with that from the uncoupled simulation. Palmer et al. [7] directly integrated the NEQAIR code into the flowfield solver DPLR to simulate the flow over FIRE II coupled with or without thermal radiation in a certain reentry state, finding that the coupled stagnation heat transfer was 15.1% lower than the uncoupled value. Martin [8] used a CFD code NH7AIR including multi-species and multi-temperature models with SPRADIAN to solve the reentry flows coupled with radiation loosely, in which the radiative transfer was calculated using the tangential slab or finite volume methods. Feldick [9] applied the one-dimensional tangential slab approximation, Monte-Carlo simulation, and the P 1 -Monte Carlo hybrid methods in combination with DPLR, respectively, to simulate the hypersonic flows surrounding the MPCV capsule coupled thermal radiation. It was found that the radiative heat flux of the tangential slab approximation was overestimated by nearly 30% than that of the Monte Carlo simulation.
However, in order to solve the radiative transfer in the high-temperature reentry flowfield, some studies have adopted a simple one-dimensional tangential slab integration model, which obviously has insufficient accuracy, while some studies have directly adopted direct numerical simulation techniques such as the Monte Carlo method, which causes very heavy computational loads. Therefore, Wang et al. [10] proposed a thermal radiation solving method library consisting of the discrete coordinate method, P 1 approximation, and optically limiting approximation, which can achieve the effective coupling simulation with the equations governing the thermochemical nonequilibrium flows around reentry vehicles.
In this paper, the hypersonic reentry flows coupled with thermal radiation were simulated using the thermochemical nonequilibrium flow equations combined with the thermal radiation solving method library. The high-temperature air radiation properties were calculated by the gray gas model [11] or the two-step model [12]. An axisymmetric spherical case under reentry conditions was first solved to investigate the effects of coupled thermal radiation simulations on the stagnation convective heating and radiative heating. Then, the three-dimensional reentry flowfield of Apollo at Ma = 31.5 was simulated to discuss the aerothermodynamic characteristics in the high-temperature thermochemical flowfield coupled with thermal radiation.

Flowfield Governing Equations Based on Three-Temperature Model
Three-temperature model is employed to describe the thermodynamic nonequilibrium in the hightemperature flowfield surrounding the reentry vehicle [13], assuming that the molecular rotational mode is fully excited and is in equilibrium with the heavy particle translational mode, corresponding to a translational-rotational temperature, T; the molecular vibrational mode corresponds to a vibrational temperature, T v ; the electronic energy and the electron translation energy correspond to an electronicelectron temperature, T e . Under the three-temperature assumption, the flowfield governing equations include the species density, total momentum, total energy, total vibrational energy, and total electronicelectron energy equations. The details are as follows: @qe e @t þ @qe e u j @x j ¼ À @q e;j @x j À @ @x j where N s is the total number of air species; the index 'dia.' denotes a diatomic molecule or an diatomic ion species; ρ s and ρ are the densities of the species s and the mixture, respectively; u and p are the flow velocity and pressure; τ is the shear stress tensor; e and h are the specific total energy and specific total enthalpy of the mixture; e v and e e are the specific vibrational energy and the specific electronic-electron energy of the mixture; h s , e v,s and h e,s are the specific enthalpy, specific vibrational energy, and specific electronicelectron enthalpy of the species s; J s,j is the mass diffusion flux of the species s in the j-direction; q j , q v,j and q e,j are the total heat flux, vibrational heat flux, and electronic-electron heat flux in the j-direction, respectively; ω s is the chemical reaction source term; ω v and ω e represents the vibrational energy and electronic-electron energy source terms.

Thermochemical Nonequilibrium Models
The translational and rotational energies of each air species in Eqs. (1) to (5) are directly calculated based on the equipartition theorem in statistical physics. The vibrational energy of the diatomic species follows a harmonic oscillator formula as follows: where R s represents the gas constant of the species s; θ v,s is the characteristic vibrational temperature of the species s. The high bound electron energy levels deviate far from the Boltzmann distribution, and make very small contributions to the total electronic-electron energy. Therefore, only the ground state and the first excited state of the electronic energy are considered in this paper. The specific electronic-electron energy is calculated using the expressions as follows: where θ e,s is the characteristic electron temperature; g 0 and g 1 represent the degeneracy of the ground state and the first excited state of the electron energy, respectively.
The chemical reactions are computed using the Gupta eleven-species twenty-reaction model [14]. Most data of the Gupta model are derived from the Dunn-Kang model [15], and the dissociation reaction data are taken from Blottner [16]. The back reaction rate in the Gupta model is calculated using the Arrhenius formula. There is a coupling effect between the chemical reactions and thermodynamic nonequilibrium. For example, the molecular vibrational energy excitation significantly affects the dissociation reaction rate, while the molecular dissociation directly causes the loss of vibrational energy. The controlling temperature is employed to take into account the effects of the thermodynamic nonequilibrium on the chemical reactions [17].
The source terms of the total vibrational and electronic-electron energy equations in the governing Eqs. (4) and (5) can be further decomposed as follows: where the translational-vibrational energy transportation term ω t−v is expressed by the Landau-Teller model [18], and the vibrational relaxation time is calculated by the Millikan-White relations [19], and the Park's high-temperature correction is included [20]; The vibration-electronic-electron energy transportation term ω v−e still follows the Landau−Teller form, and the vibration-electronic-electron relaxation time is given by Kim et al. [17]; The translation-electronic-electron energy transportation term ω t−e needs to calculate the effective collision cross-section between electrons and heavy particles, and in this paper, the electronneutral particle and electron-ion collision cross-sections are given by Gupta et al. [14] and Kim et al. [17], respectively; ω chem,v and ω chem,e represent the effects of chemical reactions on vibrational energy and electronic-electron energy, respectively, and are calculated using non-preferential models [20]; The electron pressure term ω epg models the work done by the charge-induced electric field on the electrons, and ω eii represents the electron impact ionization term, the expressions of which are taken from Gnoffo et al. [13] and Lee [21].

Numerical Schemes for Flow Simulation
The finite volume method is used to solve the above-mentioned hypersonic thermochemical nonequilibrium governing equations. The modified Steger-Warming scheme [22] is employed to discretize the convective flux, and the MUSCL method [23] is used to improve it to second-order accuracy. The viscous flux is calculated using the second-order central difference. The line relaxation algorithm [24] is used for time marching.

Thermal Radiation Solver
The radiative transfer is calculated using the thermal radiation solving method library based on the discrete coordinate method, P 1 approximation and optically limiting approximation, which is proposed by Wang et al. [10] and can achieve the effective coupling simulation with the governing equations of the thermochemical nonequilibrium flows around reentry vehicles. The radiation absorption coefficients of the high-temperature air are calculated using Wang's gray gas model [11] and Anderson's two-step model [12]. The Wang's gray gas model is expressed as follows: where the empirical parameters are a = 1.68 × 10 −19 m 2 /(kg·K 5 ), b = 1, c = 5. For Anderson's two-step model, when the radiation wave length λ ranging from 0 to 1100 Å, the radiation absorption coefficient is calculated by the formula as follows: where ρ and T are the air density and temperature, respectively; the referential density ρ 0 = 1.225 kg/m 3 ; when λ > 1100Å, the radiation absorption coefficient is calculated as follows: where I b,λ is the blackbody radiative intensity given by the Planck function. The parameters in Eq. (12) are list in Tab. 1. For the thermochemical nonequilibrium flow in this paper, the translational-rotational temperature is used as the temperature T in the two-step model.

Grid Independence and Validation
In this section, a benchmark case of FIRE II with the flight test data of radiative transfer and total heat flux is employed to show the reliability of the present numerical solver. Only the axisymmetric flowfield is simulated surrounding the blunt FIRE II forebody as shown in Fig. 1a. The freestream parameters are taken as the velocity u ∞ = 11.31 km/s, the air density ρ ∞ = 8.57 × 10 −5 kg/m 3 and the temperature T ∞ = 210 K. The air is composed of the nitrogen and oxygen with the mass fraction Y N2 = 0.77 and Y O2 = 0.23. The FIRE II surface is set to be a fully non-catalytic wall with the temperature T w = 810 K and emissivity ε = 1.0.

Grid Independence
The stagnation convective heat flux is chosen as the quantity for the grid refinement study, which has been proved very sensitive to the computational grid [25]. Additionally, this study focuses on the grid effect in the normal direction, because the numerical results are much more sensitive to the grid refinement in that direction. Five different meshes are used with 155 (tangential) × 70 (normal), 155 × 90, 155 × 110, 155 × 130, and 155 × 150, and the corresponding initial grid spacing Δn are 5 × 10 −5 m, 1 × 10 −5 m, 5 × 10 −6 m, 1 × 10 −6 m and 5 × 10 −7 m, respectively. The 155 × 130 grids are shown in Fig. 1b. Fig. 2 presents the effect of the grid refinement on the predicted convective heating. The 155 × 130 grids are refined enough to obtain a convergent value of the stagnation convective heating. In fact, when Δn = 1 × 10 −6 m, the Reynolds number based on the initial grid spacing, Re w , has been less than 1, which also follows the requirement rule of Re w~1 [26] for getting a convergent heat flux. Therefore, for all the following simulations in this paper, the initial grid spacing near the wall is refined small enough to ensure the grid independence with Re w~1 .

Validation
The 155 × 130 grids are used to simulate the thermochemical nonequilibrium flow around FIRE II coupled with thermal radiation, in which both the convective and radiative heatings are calculated for validation. Fig. 3 compares the present results of the convective heat flux with those given by two other popular CFD solvers. The present convective heating profile is consistent with those of LAURA and DPLR. Tab. 2 further lists both the stagnation convective and radiative heat flux predicted by the present codes, LAURA/NEQAIR and DPLR/NEQAIR, with the total stagnation heat flux measured in the flight test, which also shows that the present solver can give the reliable simulation results.

Sphere
The freestream conditions of the axisymmetric spherical case with R = 1.0 m is taken as the density, ρ ∞ = 3.1 × 10 −4 kg/m 3 , the temperature, T ∞ = 247 K, the air composition, Y N2 = 0.77, Y O2 = 0.23, which actually corresponds to 60 km High-altitude atmospheric parameters. The wall temperature is set to be T w = 1000 K and the fully catalytic wall is used. The wall emissivity is taken as ε = 1.0. In this case, both the uncoupled and coupled simulations with thermal radiation are performed for Ma = 15, 20, 25, 30, and 35, respectively. Due to the low computational loads of this case, the flow simulation coupled with radiation can be closely solved using the high-precision discrete ordinates method in the thermal radiation solving method library.  Fig. 4 also shows whether or not coupling the flowfield simulation with thermal radiation has little effect on the convective heating prediction. Fig. 5 shows the variations of the radiative heat flux at the sphere stagnation point with Mach number predicted by the simulations coupled with or without radiation. The comparison between Wang's gray model and Anderson's two-step model is also included. It needs to note that the vertical coordinate in Fig. 5 is logarithmic with a base of 10. As Ma increases, the radiative heat flux increases exponentially. The   radiative heat transfer given by Wang's gray gas model is much greater than that of Anderson's two-step model, which indicates that the absorption coefficient model has a significant influence on the radiative heating prediction. It is also seen in Fig. 5 that the uncoupled radiative heat transfer is greater than the coupled. Although the difference shown in Fig. 5 is not significant, it is actually remarkable when converted from the logarithmic value to the original value. According to the data at Ma = 35 shown in Tab. 3, the uncoupled and coupled convective heat fluxes are very close with only 2.59% difference, while the uncoupled radiative heat flux is 0.8 MW/m 2 by nearly 30% higher than the coupled.

Apollo
During reentry into the Earth atmosphere at hypervelocity, the temperature in the shock layer surrounding Apollo spaceship can exceed 10,000 K with the considerable radiation effects, which requires the flowfield simulation closely coupled with thermal radiation. The Apollo geometry used in this section is shown in Fig. 6a. However, only a part of the Apollo afterbody is included in the simulation. The flight Mach number of Apollo is set to be Ma = 31.5 at the angle of attack, α = 24.5°. The freestream density is taken as ρ ∞ = 3.96 × 10 −4 kg/m 3 , and the temperature is T ∞ = 252.5 K. The air compositions are Y N2 = 0.77 and Y O2 = 0.23. The wall temperature is set to be T w = 800 K, with the fully non-catalytic wall condition, and the wall emissivity is ε = 1.0. The total number of domain grids is about 655,000, as shown in Fig. 6b. The discrete ordinate method in the thermal radiation solving method library [10] was used to solve the radiative transfer equation, and the Anderson's two-step model was used to calculate the absorption coefficients of high-temperature air. Fig. 7 shows the Apollo's flowfield results with a typical high-temperature shock layer. The bow shock detached distance along the stagnation line is only about 1/30 of the forebody radius. Due to flight at an angle of attack, the pressure distribution on the Apollo forebody surface is eccentrically annular. The flow expands greatly over the Apollo shoulder to its 33°conical afterbody. Immediately behind the strong shock, the translational-rotational temperature in the flowfield can reach up to 17,000 K resulting in severe thermochemical nonequilibrium, and then the translational-rotational temperature decreases rapidly. Near  the stagnation region in the shock layer, both the translational-rotational and vibrational temperature maintains at a level of slightly above 10,000 K. It is noticeable that a high-concentration electron sheath surrounds Apollo with the number density above 10 19 m −3 in most areas and even up to 10 22 m −3 around the stagnation region, which certainly leads to communication 'black barrier' during the Apollo reentry. Fig. 8 shows the radiation characteristics in Apollo flowfield. Corresponding to the high temperature, both the average absorption coefficient and the incident radiation are remarkable in the shock layer, while the distribution of the latter is more diffused. Both peaks of the average absorption coefficient and the incident radiation occur in the stagnation region. Fig. 9 shows the aerodynamic heating distributions on the Apollo surface. The peak of the convective heat flux (q cw ) is located at the Apollo lower shoulder, rather than at the stagnation point. This is the result of flying with an angle of attack and a much smaller shoulder radius than the blunt forebody radius. The distribution of the radiative heat flux (q rw ) reflects the integral effect of thermal radiation. The radiative heating maximum still appears near the stagnation point with the highest temperature and density forming the strongest radiative transfer. The peak of the convective heat transfer is 2.29 MW/m 2 , while the maximum of the radiative heat flux is 0.872 MW/m 2 , which shows that, in such flight condition, the radiative heating has been in a considerable level of about 38% of the convective heating.

Conclusion
The axisymmetric spherical case shows that whether or not coupling the flowfield simulation with thermal radiation has a little influence on the convective heating prediction, but has some impacts on the radiative heating calculation. At Ma = 35, the radiative heat flux predicted by the flow simulation uncoupled with radiation is almost 30% higher than that coupled. The simulation results of Apollo reentry closely coupled with radiation show that the thermal radiation distributions mainly correspond to the temperature distribution. Both the average absorption coefficient and incident radiation are remarkable in the shock layer, while the distribution of the latter is more diffused. The peak of the radiative heating reaches nearly 38% of the maximum of the convective heating on the Apollo surface and makes a considerable contribution to the total aerodynamic heating. Therefore, under hypervelocity conditions, it is necessary to simulate the high-temperature thermochemical nonequilibrium flows around reentry vehicles coupled with thermal radiation, and take into account the radiation effects in the aerodynamic heating.