Fluid / thermal / chemical non-equilibrium simulation of hypersonic reentry vehicles

A multi-physics frame work has been setup for the simulation of surface heat flux for nonablating hypersonic reentry vehicles and presented in this paper. The main goal of this work was to set up a simple approach for the heat flux prediction during the reentry of the vehicle. The vehicle considered in the calculation is an axisymmetric vehicle flying at zero degree angle of attack. Chemical nonequilibrium in the flowfield is simulated by implementing a set of finite rate equations in the laminar finite rate model in FLUENT. The frame work set up was validated with the results available in the literature. Good correlation was observed between the results from the commercial code with the implemented equations and the results from the literature.


INTRODUCTION
The hypersonic reentry vehicle is immersed in a hot gas as the kinetic energy of the flow is dissipated by viscous effects.The detached shock wave in front of the vehicle also increases the gas temperature.Thermal Protection System (TPS) is used to insulate the hypersonic reentry vehicle from the high temperature during the reentry [1].Numerical simulation of the heat transfer during the reentry flight and the material thermal response of the TPS can be challenging due to the multiphysics interactions like chemically reacting flow, radiation and heat conduction.A literature survey was conducted to examine the fluid thermal coupling procedures for hypersonic reentry vehicles.It was identified that aero thermal heat flux prediction and material thermal response is evaluated using in house codes like SACCARA (Sandia Advanced Code for Compressible Aerothermodynamics Research and Analysis) [2], GIANTS (Gauss-Siedel Implicit Aerothermodynamic Navier-Stokes code with thermo chemical surface conditions) [3], COYOTE II [2], FIAT (Fully Implicit Ablation and Thermal response code) [3] and so on.With the advances in commercial Computational Fluid Dynamics (CFD) and Finite Element Analysis (FEA) codes, the prediction of aero thermal heat flux and the material thermal response of the TPS using these codes are highly feasible.The primary issue of this work is to identify and understand the capability of these codes for predicting the chemically reacting hypersonic flow.FLUENT is used as the CFD code as it is capable of modeling the species transport apart from the basic fluid flow [4].The chemically reacting hypersonic flow for an axisymmetric vehicle at zero degree angle of attack is solved using FLUENT and in the current work.

GOVERNING CFD EQUATIONS
The computational fluid dynamics calculations are performed using the numerical code FLUENT.The fluid has been modeled as a reacting gas in thermal equilibrium and chemical non-equilibrium.The flow is assumed to be laminar.The governing equations are as follows: The ideal gas equation of state is used

CHEMICAL NON-EQUILIBRIUM MODEL
As mentioned earlier the fluid which is air surrounding the reentry vehicle is at a high temperature due to kinetic energy dissipation and shock wave.The high temperature causes air to dissociate and even ionize.The temperature in the nose area of Apollo reentry was about 11,000 K at a Mach number of 35 [5].The constant specific heat assumption becomes invalid at those temperatures.Chemical nonequilibrium assumption says that the characteristic chemical reaction time is the same as the characteristic time of fluid motion.The laminar finite-rate model in FLUENT computes the chemical species production rate using modified Arrhenius equation.The net source term of chemical species i is computed as [4] (10) The net rate of creation of species i in reaction r is given by (11) The forward rate constant for reaction r is modeled using the Arrhenius expression (12) The backward rate constant for reaction r is computed from the forward rate using (13) The two-temperature model of Park [2] consisting of five species (N 2 , O 2 , NO, N, O) is used.The Park model uses the average of translational and vibrational temperatures to calculate the rate constants [5].The reaction rate coefficients and the characteristic temperature of dissociation are given in Table 1.In [2] the Park model is used with the assumption of thermo-chemical nonequilibrium whereas in this work it is used with the assumption of thermal equilibrium and therefore it is assumed that all the internal energy modes are in equilibrium at temperature T.

TRANSPORT PROPERTIES
The transport properties of the species are obtained from the kinetic theory of gases [4].The fluid viscosity is defined using kinetic theory as (14) The Lennard-Jones parameters used were obtained from ref [5] and [6] and are given in Table 2.
The thermal conductivity of a pure gas is defined using kinetic theory as The mass diffusion coefficient is defined using kinetic theory as (16) The mass diffusion coefficient of species i in the mixture is defined as The mixture values of µ and k for the chemically reacting gas is defined using Wilke's rule Where (20)

THERMAL EQUILIBRIUM MODEL
When the fluid temperature increases the internal energy modes get excited and the specific heat ratio is not constant.Thermodynamic properties of high temperature air in the form of polynomial curve fits exist which can used to implement this model.This model assumes that the thermal relaxation time is much smaller than characteristic time of fluid motion.
The values of specific heat of all the species are required for the case of thermal equilibrium.The curve fits provided in [7] are used.Originally they are curve fit with the following polynomial with five ranges 300 Since FLUENT allows only three ranges for the temperature dependent specific heat and also because of the limitation that specific heat values cannot be modified in FLUENT [4], the first and second ranges and the fourth and fifth ranges are combined and curve fit with seventh degree polynomials.The specific heat of the mixture is calculated as a mass fraction average of the pure species heat capacities [4] (22)

BOUNDARY CONDITIONS
Pressure far-field conditions are used in FLUENT to model a free-stream condition at infinity, with free-stream Mach number and static conditions being specified [4].The freestream conditions for the calculations are given in Table 3 [2].The free stream static pressure, Mach number and static temperature are given as input to the pressure far-field.
Wall boundary condition is used to bind the solid and fluid region in FLUENT [4].For nonablating case, a no-slip condition is used.It indicates that the fluid sticks to the wall.An isothermal wall is used for the thermal boundary condition.A constant wall temperature of 294.4 K was assumed for the initial simulation.Fully catalytic wall with species concentrations equal to the freestream composition (77% N 2 and 23% O 2 by weight) is assumed.The Axis boundary condition is specified at the centerline of the axisymmetric geometry.

NUMERICAL TECHNIQUE
The density-based solver in FLUENT solves the governing equations of continuity, momentum, energy and species transport simultaneously The inviscid flux vector F in Eqn 23 is computed using a flux-vector splitting scheme Advection Upstream Splitting Method (AUSM) [4].The spatial discretization is specified using a second order upwind scheme.Time discretization is accomplished using an implicit method.
Computational domain used in the CFD.

RESULTS AND DISCUSSION
3.1.COMPUTATIONAL DOMAIN An axisymmetric IRV-2 vehicle [2] at zero degree angle of attack is considered for the simulations (see Figure 1).The vehicle, which is a sphere-biconic-cylinder, has a nose radius of 0.01905 m and a total length of 1.3866 m.The biconic angles are 8.42°and 6.10°with a break at 0.1488m.The flow field grid for the CFD simulation is a structured grid created using Gridgen [8] and is shown in Figure 2. The first cell spacing near the wall is chosen as 1 × 10 -6 m. Figure 3 shows the surface heat flux evaluated by FLUENT for the trajectory point 1 in Table 3.
The surface heat flux is plotted versus the y-coordinate which corresponds to the radial direction.The stagnation point heat flux is 601.99415W/cm 2 for the grid size 64 × 64 cells in axial and radial directions.The grid was refined to 128 × 128 to study the effect on the convective heat transfer.The stagnation point heat flux value did not change appreciably compared to the coarser 64 × 64 grid which is illustrated in figure 3. Two cases are considered for the species at the wall:

Non catalytic wall
In this case the mass fraction gradients at the wall are assumed to be zero [9] which translates to no reactions at the wall that is,

Fully catalytic wall
In this case all the dissociated atoms are recombined at the wall.The wall material, i.e., the TPS material, in this case, catalyzes chemical reactions at the surface [5].The species concentration at the surface is assumed to be of freestream composition Figure 4 compares the surface heat flux in the two cases.The stagnation point heat flux in the first case is 452.32 W/cm 2 whereas in the second case it is 601.99415W/cm 2 due to the recombination of dissociated species which is an exothermic reaction.It is observed that the stagnation point heat flux in the case of non catalytic wall is about 25% less than that of the flux determined for fully catalytic wall.The contours of static temperature and pressure for the first trajectory point are shown in figures 5 and 6.The surface heat flux distributions for the first three trajectory points are shown in figure 7 and the respective heat flux distributions are shown in figure 8.After the third trajectory point it became difficult to determine the CFD solution using second order numerical scheme with the grid 128 × 128.Nevertheless using first order numerical scheme enabled to achieve a solution with a compromise in the accuracy.At this time, the CFD mesh is refined further and being investigated with second order scheme and hence the solution is not presented for further trajectory points.

VALIDATION
Figure 9 compares the total surface heat flux determined using the present method with the heat flux obtained from [2] where the authors performed the computational fluid dynamics calculations using SACCARA (Sandia Advanced Code for Compressible Aerothermodynamics Research and Analysis).The authors have assumed the flow to be in thermo-chemical nonequilibrium.The heat flux data in the form of graph from [2] is extracted using xy Extract, software which extracts data from a graphic contained in a bitmap file.
Three different diffusion models were implemented in the CFD code, they were: 1. Constant Lewis Number = 1.4 [10] Lewis number is the ratio of thermal diffusivity to mass diffusivity [11] (26) 2. Constant Schmidt Number = 0.5 [2] Schmidt number is the ratio between momentum diffusivity and mass diffusivity [11] (27)

Maxwell-Stefan Equations
The multicomponent diffusion coefficient D im is arrived by simplifying the Maxwell-Stefan equations [20] (28) The multicomponent diffusion model is computationally expensive as it involves the calculation of binary mass diffusion coefficient.The influence of the different diffusion models on the surface heat flux is shown in Figure 10.The one with constant Schmidt number is found to be a good approximation to the exact diffusion model.

VALIDATION OF CONTINUUM ASSUMPTION
One of the aspects of hypersonic flow is the low density at the reentry conditions.Typically the densities at the reentry conditions are very small and hence the mean free path which is the average distance between molecular collisions is comparable to the characteristic length.Therefore the flowfield cannot be considered as continuum and hence the Navier-Stokes equations become non applicable.A non dimensional number called Knudsen number Kn is used to determine the nature of the fluid and it is defined as [12] (29) where λ the mean free path and L is the characteristic flow field dimension.The mean free path is calculated from kinetic theory using [5] (30)    where k is the Boltzmann constant which is equal to 1.38 × 10 −23 J/K and πd 2 is called collision cross section.The mean free path for the reentry conditions of IRV-2 vehicle is calculated with the characteristic molecular diameter of air as 3.711 × 10 −10 m [6] λ = 6.2846e − 4 m.The Knudsen number for the reentry conditions of IRV-2 vehicle is calculated with the characteristic [4].The system of governing equations solved by the FLUENT solver in vector form is given by[4]

Figure 3 2 )Figure 4
Figure 3 Grid refinement study for surface heat flux.

Figure 6
Figure 6 Contours of static pressure for trajectory Point 1.

Figure 5
Figure 5 Contours of static temperature for trajectory point 1.

Table 1
Reaction rate constants M vib = N 2 , O 2 , NO M atom = N, O

Table 3
Freestream conditions A r , β r , E r constants in the Arrhenius expression of the rate constant C pi specific heat of species i, J/kg.K C p specific heat of the mixture, J/kg⋅K C H convective heat transfer coefficient, W/m 2 K C i molar concentration of species of i, kgmol/m 3 D ij binary mass diffusion coefficient of species i in species j, m2/s D im mass diffusion coefficient of species i in the mixture, m2/s