# Computationally Efficient Assessments of the Effects of Radiative Transfer, Turbulence Radiation Interactions, and Finite Rate Chemistry in the Mach 20 Reentry F Flight Vehicle

Authors: Gautham Krishnamoorthy; Lauren Elizabeth Clarke
Journal: Journal of Computational Engineering
 (2016)

---

## Abstract

Effects of finite rate chemistry, radiative heat transfer, and turbulence radiation interactions (TRI) are assessed in a fully coupled manner in simulations of the Mach 20 Reentry F flight vehicle. Add-on functions were employed to compute a Planck mean absorption coefficient and the temperature self-correlation term (for TRI effects) in the optically thin shock layer. Transition onset was induced by specifying a wall roughness height at the experimentally observed transition location. The chemistry was modeled employing eight elementary reactions and an equilibrium approach allowing species to relax towards their chemical equilibrium values over the process characteristic time scale. The wall heat fluxes in the turbulent region, density, and velocity profiles compared reasonably well against measurements as well as similar calculations reported previously. The density predictions were more sensitive to the choice of modeling options than the velocities. The radiative source term magnitude agreed closely with its measurements deduced from shock tube experiments. The TRI model predicted a 60% enhancement in emission due to temperature fluctuations in the turbulent boundary layer. While the variations in density and velocity predictions among the models diminished along the length of the body, the O and NO prediction variations extended well into the turbulent boundary layer.


---

## Body

## 1. Introduction

Planetary entry vehicles are subjected to significant aerothermal heating as they dissipate their kinetic energy into the atmosphere of their destination planet. The very high velocity of an earth bound reentry vehicle causes viscous dissipation and the formation of a shock wave resulting in very high air-temperatures surrounding the vehicle. Intense thermochemical processes occur in the region between the body surface and the shock. The molecules in the surrounding air dissociate and might be far from their equilibrium state [1, 2]. These reactions may be further catalyzed by the products resulting from the decomposition and ablation of the vehicle’s thermal protection system (TPS). For earth reentry, the concentrations of N and O and their recombination to produce NO molecules surrounding the vehicle in particular have been determined to be important towards determining the shock and boundary layer characteristics and for providing more reliable TPS sizing requirements [3, 4]. Further, at certain reentry conditions radiative transfer can contribute both to the surface heat fluxes and to the cooling of the shock and boundary layers surrounding the vehicle depending on the spectral emission and absorption properties of the surrounding molecules [5, 6]. Computational Fluid Dynamics (CFD) through a coupled modeling of nonequilibrium, reacting flows with radiative heat transfer can yield valuable insights into the flow characteristics within these hypersonic flight environments [7, 8]. However, due to the fact that radiative heat transfer calculations are inherently expensive as a result of its strong spatial, directional, and wave-length dependencies [9, 10], the radiation field is often computed in a 1D domain. Further, the coupling between the radiation and the thermal field is accomplished in a loose coupling/decoupled manner [11]. These approaches are justified by facts that, at high reentry velocities, the radiative heat flux to the stagnation point is significant. The shock stand-off distance is small compared with the nose radius of vehicles, and a spatial change of physical properties is much more severe in the direction normal to the shock than in the tangential direction. Since it is well known from simulations of combustion phenomena that a strong coupling between the fluid dynamics and radiative heat transfer is necessary to accurately predict the concentrations of minor species such as O, N, and NO [12, 13], this may hold true in hypersonic flows as well. Similarly, combustion simulations have also examined and highlighted the importance of turbulence-chemistry and turbulence radiation interaction (TRI) models on these minor specie concentrations [14, 15]. Therefore, the primary goals of this paper are as follows:(1)
To demonstrate a simple methodology for providing initial assessments of radiative transfer and TRI effects in a computationally efficient manner in hypersonic flow simulations for an optically thin shock layer.
(2)
To provide an initial assessment of fully taking into account radiative transfer and TRI effects on the density, velocity, and O and NO predictions surrounding a hypersonic vehicle.
(3)
To assess the differences between employing finite rate chemistry and equilibrium based thermochemical models on the density, velocity, and O and NO predictions surrounding the vehicle [16, 17].The hypersonic flow scenario specifically chosen for the purpose is the Reentry F flight vehicle which is a cone with a 5-degree half angle with a small spherical nose tip (Figure1) [18]. This has been a valuable dataset for evaluating the predictions of several CFD codes. The criteria that make this a good test case for this investigation are as follows:(1)
The gas temperature was found to be more than 6000 K at the nose tip with dissociation of oxygen and nitrogen occurring. Further dissociation of oxygen in the boundary layer was observed as a result of viscous dissipation resulting in temperatures around 3000 K [19]. The high temperatures make it a good test case to investigate the sensitivity of the predictions to equilibrium chemistry, finite rate chemistry, and radiation models.
(2)
At an altitude of 24.4 km (which is the conditions investigated in this study), the laminar to turbulence transition occurred at about 62.5% the total wall length translating to 2.5 m along the length of the simulated vehicle. This implies a significant interplay between turbulence, chemistry, and radiative transfer spanning the remaining 1.5 m of the simulated cone providing sufficient length of the domain to investigate the effects of TRI on the flow field.
(3)
Local thermodynamic equilibrium (LTE) can be assumed at this altitude. When shock passes through air, the energy goes initially into the translational energy of the molecules and is later redistributed into the rotational and vibrational energy levels. At 24.4 km this energy redistribution rapidly results from collisions but can be much slower at higher altitudes. However, LTE is a good approximation at altitudes between 15 km and 45 km [20]. Therefore, this alleviates the need to solve separate transport equations for the rotational and vibrational temperatures and subsequently use a geometric average of those temperatures to estimate the reaction rates, radiative transfer.
(4)
In general, the pathlines (the path followed by a fluid element) after the strongly dissociating shock layer will determine the chemical composition at any given location [2]. Therefore, at large distances downstream of the bow shock the chemistry is assumed to relax towards equilibrium. Considering the large length of the vehicle (4 m), the differences between the finite rate and equilibrium predictions can be quantified and the validity of this assumption can be assessed at a number of locations along the body including the turbulent boundary layer.
(5)
The availability of numerical results encompassing this scenario from previous investigators is needed to assess the validity of our simulation techniques.Figure 1
(a) Geometry of the reentry F vehicle. A 5° half-angle cone (represented as “A”) with a spherical nose tip of radius 0.00254 m (represented as “B”) and a total length of 4.0 m; (b) the mesh associated with the vehicle geometry employed in this study; (c) the entire computational domain along with the lines where the results were monitored.
(a)
(b)
(c)Previous simulations of this reentry condition have been carried out by several researchers [21–26] with a focus mainly on predicting the convective heat transfer rates on the surface. Since the contribution of surface radiative fluxes at these reentry conditions was expected to be negligible, these previous investigations assumed equilibrium chemistry and neglected radiative transfer. Consequently, finite rate chemistry, radiative transfer, and TRI effects on the temperature and species concentration have not been rigorously assessed. For instance, Wood et al. [21] provided three equilibrium air numerical solutions employing thin-layer Navier Stokes equation methods and viscous shock layer solution methods employing equilibrium air chemistry and neglecting the effects of radiative transfer. The viscous shock layer method was also determined to work well in the study by Thompson et al. [22]. Hamilton et al. [23] employed a variable entropy boundary layer code in conjunction with equilibrium chemistry. Recent studies have emphasized the fact that the prediction of transitional heating from this experiment is beyond the current state of the art [27]. Therefore, as a result of judicious a posteriori modeling of the transitional region (i.e., beginning and end of the transition region), the numerical predictions of Sutton et al. [24] were in excellent agreement with the transitional heating measurements. Further, by carrying out several detailed predictions of the Reentry F flight data, Wurster et al. [25] and more recently Barnhardt and Candler [26] have emphasized the sensitivities of nose bluntness and the angle of attack considered in the simulations on the heating rate predictions.


## 2. Methods

The simulations in this study were carried out employing the density based solver framework of the commercial code ANSYS FLUENT (version 12) [30]. Further, user-defined functions (UDFs) were developed to compute a temperature and pressure dependent absorption coefficient and model the TRI. In order to assess the influence of different models on the boundary and shock layer predictions, four scenarios were investigated:(1)
Finite rate chemistry without radiative transfer.
(2)
Fully coupled finite rate chemistry calculations with radiative transfer.
(3)
Fully coupled finite rate chemistry and radiative transfer calculations accounting for TRI.
(4)
Chemical equilibrium calculations without radiative transfer/TRI where the species were allowed to relax to their chemical equilibrium state depending on the characteristic (convective/diffusive + chemical) time scale associated with the process.


### 2.1. Geometry and Meshing

In order to get an initial assessment of the impact of different models, 2D axisymmetric flow conditions were assumed as in previous studies [19] by neglecting the experimentally reported angle of attack of 0.14 degrees. Further, the ablation of the nose tip was neglected. The nose tip radius (B, in Figure 1) was found to increase from 0.00254 m to 0.00343 m during the experiments. However, neglecting this has a minimal impact on the monitor locations along the body. The entire domain was meshed employing 280,000 cells employing a boundary layer meshing strategy. Comparing this against the 83,200 cells employed in a previously published simulation of these conditions [19] indicates the adequacy of this mesh resolution to obtain grid independent results.


### 2.2. Radiative Transfer and TRI Modeling

Since the shock layer surrounding the vehicle was determined to be optically thin (k
L
≪
1), where k is the absorption coefficient (in m−1) and L the path-length (in m), an appropriately determined Planck mean absorption coefficient was computed based on the local pressure and temperature [6] and implemented as a user-defined function (UDF). This is similar to the approach employed in the combustion modeling community where Planck mean absorption coefficients are routinely employed to couple radiation with the flow field in optically thin flames particularly when radiation is not the dominant energy transfer mechanism [31]. We also like to point out that the goal here was to provide an initial assessment on radiative transfer on the temperature and species concentration and not on a highly accurate radiative transfer calculation. A highly accurate spectral 1D calculation in itself would be 1 to 4 orders of magnitude longer than the time to compute a fluid solve [7].To fully account for the effects of TRI within the context of time-averaged flow simulations (such as Reynolds Averaged Navier Stokes Equation Models) requires time-averaging of the radiative transport equation (RTE) written in terms of the radiative intensity (I). In the absence of scattering, this may be written as(1)
∇
·
I
-
r
→
,
s
^
s
^
=
-
k
I
-
r
→
,
s
^
+
k
I
b
-
r
→
,
s
^
.In (1) “r
→” and “s
^” are the position and directional vectors, respectively, and I
b is the blackbody emissive power. Equation (1) results in several “closure terms” that require some form of modeling to fully account for TRI, mainly the temperature self-correlation (TSC) and absorption-temperature correlation (second term on the right hand side of (1)) and the absorption-intensity correlation (first term on the right hand side of (1)) [32]. However, for an optically thin shock layer, the thin eddy approximation may be invoked, which neglects the absorption-intensity correlation. With the additional assumption of computing the absorption coefficients based on mean species and temperature values we get (2)
∇
·
I
-
r
→
,
s
^
s
^
=
-
k
I
-
r
→
,
s
^
+
k
I
b
-
r
→
,
s
^
.The validity of these approximations has also been shown in recent studies of hypersonic turbulent boundary layers where the effects of absorption TRI due to radiation emitted in the shock layer were minimal [33]. Further, calculations in this study indeed showed that magnitude of emission was more than three orders of magnitude larger than the absorption term. Therefore, only the TSC term (the time-averaged I
b
- term in (2)) was employed to account for TRI effects. The time-averaged blackbody emissive power may be expanded in Taylor’s series about the mean temperature as [32](3)
I
b
T
=
∑
n
=
0
∞
T
′
n
n
!
d
n
d
T
n
I
b
T
T
¯
,where T
′ is the fluctuation in the temperature about the mean. Time-averaging equation (3) results in an expression for the time-averaged emissive power as(4)
I
b
¯
T
=
∑
n
=
0
∞
T
′
n
¯
n
!
d
n
d
T
n
I
b
T
T
¯
.Expanding (4) results in the following expression:(5)
T
4
¯
T
¯
4
=
1
+
6
T
′
2
¯
T
¯
2
+
4
T
′
3
¯
T
¯
3
+
T
′
4
¯
T
¯
4
.If the probability density functions (pdfs) of the temperature fluctuations can be assumed to be Gaussian, then odd-order terms can be assumed to be small or zero as suggested by data from the literature [34]. Further, experimental evidence points to the fact that [35](6)
T
4
¯
T
¯
4
≫
T
′
4
¯
T
¯
4
.Therefore, the TSC term in this study was modeled as(7)
T
4
¯
T
¯
4
=
1
+
6
T
′
2
¯
T
¯
2
.The turbulence was modeled in this study using the realizable k
-
ε turbulence model. A transport equation for the temperature variance (T
v
a
r) was solved [30] employing the computed values of turbulent viscosity (μ
t), turbulent kinetic energy (k), and the turbulence dissipation rates (ε) as (8)
∂
ρ
T
v
a
r
∂
t
+
∇
·
ρ
v
→
T
v
a
r
=
∇
μ
t
σ
t
∇
T
v
a
r
+
C
g
μ
t
∇
T
2
-
C
d
ρ
ε
k
T
v
a
r
.Constants σ
t, C
g, and C
d in (4) take the values 0.85, 2.86, and 2.0, respectively [30].


### 2.3. Chemistry Modeling

To account for finite rate chemistry effects, five species were considered: N2, O2, NO, N, and O. The temperatures encountered in this study were below the ionization temperatures (>9000 K) of the species. In this study, the kinetic rate constants corresponding to the two-temperature model of Park [1] were evaluated employing the translational temperature. In the Park model, an average of the translational and vibrational temperatures needs to be employed to evaluate the kinetic rates. However, since LTE can be assumed at the altitudes considered in this study, the need to solve a separate transport equation for the vibrational temperature is alleviated [20]. The rate constants employed in this study are listed in Table 1. The species transport properties were estimated employing the kinetic theory of gases employing Lennard-Jones parameters [30]. The specific heat of each gas species was computed from polynomial curve-fits encompassing the temperature range 300 K to 15000 K listed in Sockalingam and Tabiei [28, 29].Table 1
Reactions and rate constants employed in this study [28].
Reaction
Preexponential factor
Activation energy (J/kmol)
Temperature exponent
M = O2, N2, NO
M = N, O
M = O2, N2, NO
M = N, O
M = O2, N2, NO
M = N, O
N2 + M → 2N + M
7.0E + 18
3.0E + 19
9.4114E + 08
9.4114E + 08
−1.6
−1.6
O2 + M → 2O + M
2.0E + 18
1.0E + 19
4.9468E + 08
4.9468E + 08
−1.5
−1.5
NO + M → N + O + M
5.0E + 12
1.1E + 14
6.2771E + 08
6.2771E + 08
0.0
0.0
N2 + O → NO + N
6.4E + 14
3.1926E + 08
−1.0
NO + O → O2 + N
8.4E + 09
1.6171E + 08
0.0The general conservation equation for the mass fraction of specie “i” (Y
i) can be written as(9)
∂
∂
t
ρ
Y
i
+
∇
·
ρ
v
→
Y
i
=
-
∇
·
J
→
i
+
R
i
,where “ρ” denotes the density and “J
→” and “R” represent the diffusive fluxes and reaction source term of the species “i,” respectively. Based on an established solution obtained from solution to (9) employing the finite rate chemistry model, the deviation of the species compositions and temperature from their thermochemical equilibrium state was determined by allowing the compositions to relax to their chemical equilibrium values [30]. This was undertaken by allowing the reaction source term in (9) to be modeled as(10)
R
i
=
ρ
Y
i
e
q
-
Y
i
τ
c
h
a
r
,where Y
i is the mass fraction of species i and Y
i
e
q is the equilibrium mass fraction of species i (determined as a function of the local species, temperature, and pressure) and a characteristic time scale (τ
c
h
a
r). τ
c
h
a
r in turn is modeled as the sum of a convective/diffusion time scale in the computational cell and a chemical time scale determined from the chemical kinetic mechanism. Equation (10) basically implies that the species reaction towards their chemical equilibrium state over the τ
c
h
a
r period associated with the process.


### 2.4. Turbulence Modeling and Initiating the Onset of Transition

Predicting the onset of transition in this hypersonic flow regime has proven to be really challenging. However, fairly accurate predictions of the surface heat flux profiles have been obtained by judicious specification of the location for the onset of transition [19, 24, 27] with typical prediction accuracies ranging from 20 to 25% [27].The turbulence was modeled in this study employing the realizablek
-
ε turbulence model and the wall temperature of the vehicle was assumed to be constant at 500 K [19]. Since there is a factor of 3 variation among current empirical correlations for the extent and onset of transition (as a function of the edge Mach number) [36], the onset of transition at a distance of 2.5 m along the wall was artificially induced in this study by specifying a roughness height of 5 × 10−4 m along the wall from 2.5 m to 4 m. Wall roughness is one of the driving mechanisms by which disturbances can be induced in the boundary layer to initiate transition [27]. The predicted wall heat flux profiles in the turbulent boundary layer region are shown in Figure 2. It is worth noting that although the prediction accuracies at the peak heating locations are within the range reported in previous studies [19], the goal here is not a rigorous assessment of turbulence or transition prediction methods but on assessing the impacts of chemistry and radiative transfer models. We can also see in Figure 2 that as far as the surface heat fluxes are concerned, results from the different chemistry models (relaxation to chemical equilibrium versus finite rate chemistry) as well as the inclusion of radiative transfer/TRI all yield identical predictions. A complete summary of the modeling options employed in this study is listed in Table 2.Table 2
Summary of physical models employed in the simulations.
Physical models
Primary modeling option
Density
Ideal gas
Specific heat
Piecewise polynomial [29]
Thermal conductivity and mass diffusivity
Kinetic theory
Viscosity
Sutherland law
Wall boundary condition
No slip (fluid flow), constant temperature of 500 K (energy equation)
Far field boundary conditions
Gauge pressure: 2761 pascal
Mach number: 19.97
Turbulence
Realizablek
-
ε
Turbulence initiation
Roughness height of 5 × 10−4 m specified in the wall region 2.5 m to 4 m
Turbulence boundary conditions
Intensity: 0.01%
Viscosity ratio: 2.29E − 09
Radiative transport equation solver
Finite volume radiation model (angular resolution, theta × phi: 3 × 3)
Planck mean absorption coefficient (Kp) [6]∗
3.7516
×
10
-
6
·
(
P
)
1.31
·
exp
⁡
(
5.18
×
10
-
4
T
-
7.13
×
10
-
9
T
2
)
Gas phase kinetics
Table1
Relaxation to chemical equilibrium (kinetic mechanism for chemical time scale determination)
Table1
Turbulence-radiation interactions∗
Equation (7)
∗These models were implemented as user-defined functions (UDFs) in ANSYS FLUENT [30].Figure 2
Wall heat flux profiles in the turbulent boundary layer predicted by the realizablek
-
ε turbulent model (with a specified roughness height).


## 2.1. Geometry and Meshing

In order to get an initial assessment of the impact of different models, 2D axisymmetric flow conditions were assumed as in previous studies [19] by neglecting the experimentally reported angle of attack of 0.14 degrees. Further, the ablation of the nose tip was neglected. The nose tip radius (B, in Figure 1) was found to increase from 0.00254 m to 0.00343 m during the experiments. However, neglecting this has a minimal impact on the monitor locations along the body. The entire domain was meshed employing 280,000 cells employing a boundary layer meshing strategy. Comparing this against the 83,200 cells employed in a previously published simulation of these conditions [19] indicates the adequacy of this mesh resolution to obtain grid independent results.


## 2.2. Radiative Transfer and TRI Modeling

Since the shock layer surrounding the vehicle was determined to be optically thin (k
L
≪
1), where k is the absorption coefficient (in m−1) and L the path-length (in m), an appropriately determined Planck mean absorption coefficient was computed based on the local pressure and temperature [6] and implemented as a user-defined function (UDF). This is similar to the approach employed in the combustion modeling community where Planck mean absorption coefficients are routinely employed to couple radiation with the flow field in optically thin flames particularly when radiation is not the dominant energy transfer mechanism [31]. We also like to point out that the goal here was to provide an initial assessment on radiative transfer on the temperature and species concentration and not on a highly accurate radiative transfer calculation. A highly accurate spectral 1D calculation in itself would be 1 to 4 orders of magnitude longer than the time to compute a fluid solve [7].To fully account for the effects of TRI within the context of time-averaged flow simulations (such as Reynolds Averaged Navier Stokes Equation Models) requires time-averaging of the radiative transport equation (RTE) written in terms of the radiative intensity (I). In the absence of scattering, this may be written as(1)
∇
·
I
-
r
→
,
s
^
s
^
=
-
k
I
-
r
→
,
s
^
+
k
I
b
-
r
→
,
s
^
.In (1) “r
→” and “s
^” are the position and directional vectors, respectively, and I
b is the blackbody emissive power. Equation (1) results in several “closure terms” that require some form of modeling to fully account for TRI, mainly the temperature self-correlation (TSC) and absorption-temperature correlation (second term on the right hand side of (1)) and the absorption-intensity correlation (first term on the right hand side of (1)) [32]. However, for an optically thin shock layer, the thin eddy approximation may be invoked, which neglects the absorption-intensity correlation. With the additional assumption of computing the absorption coefficients based on mean species and temperature values we get (2)
∇
·
I
-
r
→
,
s
^
s
^
=
-
k
I
-
r
→
,
s
^
+
k
I
b
-
r
→
,
s
^
.The validity of these approximations has also been shown in recent studies of hypersonic turbulent boundary layers where the effects of absorption TRI due to radiation emitted in the shock layer were minimal [33]. Further, calculations in this study indeed showed that magnitude of emission was more than three orders of magnitude larger than the absorption term. Therefore, only the TSC term (the time-averaged I
b
- term in (2)) was employed to account for TRI effects. The time-averaged blackbody emissive power may be expanded in Taylor’s series about the mean temperature as [32](3)
I
b
T
=
∑
n
=
0
∞
T
′
n
n
!
d
n
d
T
n
I
b
T
T
¯
,where T
′ is the fluctuation in the temperature about the mean. Time-averaging equation (3) results in an expression for the time-averaged emissive power as(4)
I
b
¯
T
=
∑
n
=
0
∞
T
′
n
¯
n
!
d
n
d
T
n
I
b
T
T
¯
.Expanding (4) results in the following expression:(5)
T
4
¯
T
¯
4
=
1
+
6
T
′
2
¯
T
¯
2
+
4
T
′
3
¯
T
¯
3
+
T
′
4
¯
T
¯
4
.If the probability density functions (pdfs) of the temperature fluctuations can be assumed to be Gaussian, then odd-order terms can be assumed to be small or zero as suggested by data from the literature [34]. Further, experimental evidence points to the fact that [35](6)
T
4
¯
T
¯
4
≫
T
′
4
¯
T
¯
4
.Therefore, the TSC term in this study was modeled as(7)
T
4
¯
T
¯
4
=
1
+
6
T
′
2
¯
T
¯
2
.The turbulence was modeled in this study using the realizable k
-
ε turbulence model. A transport equation for the temperature variance (T
v
a
r) was solved [30] employing the computed values of turbulent viscosity (μ
t), turbulent kinetic energy (k), and the turbulence dissipation rates (ε) as (8)
∂
ρ
T
v
a
r
∂
t
+
∇
·
ρ
v
→
T
v
a
r
=
∇
μ
t
σ
t
∇
T
v
a
r
+
C
g
μ
t
∇
T
2
-
C
d
ρ
ε
k
T
v
a
r
.Constants σ
t, C
g, and C
d in (4) take the values 0.85, 2.86, and 2.0, respectively [30].


## 2.3. Chemistry Modeling

To account for finite rate chemistry effects, five species were considered: N2, O2, NO, N, and O. The temperatures encountered in this study were below the ionization temperatures (>9000 K) of the species. In this study, the kinetic rate constants corresponding to the two-temperature model of Park [1] were evaluated employing the translational temperature. In the Park model, an average of the translational and vibrational temperatures needs to be employed to evaluate the kinetic rates. However, since LTE can be assumed at the altitudes considered in this study, the need to solve a separate transport equation for the vibrational temperature is alleviated [20]. The rate constants employed in this study are listed in Table 1. The species transport properties were estimated employing the kinetic theory of gases employing Lennard-Jones parameters [30]. The specific heat of each gas species was computed from polynomial curve-fits encompassing the temperature range 300 K to 15000 K listed in Sockalingam and Tabiei [28, 29].Table 1
Reactions and rate constants employed in this study [28].
Reaction
Preexponential factor
Activation energy (J/kmol)
Temperature exponent
M = O2, N2, NO
M = N, O
M = O2, N2, NO
M = N, O
M = O2, N2, NO
M = N, O
N2 + M → 2N + M
7.0E + 18
3.0E + 19
9.4114E + 08
9.4114E + 08
−1.6
−1.6
O2 + M → 2O + M
2.0E + 18
1.0E + 19
4.9468E + 08
4.9468E + 08
−1.5
−1.5
NO + M → N + O + M
5.0E + 12
1.1E + 14
6.2771E + 08
6.2771E + 08
0.0
0.0
N2 + O → NO + N
6.4E + 14
3.1926E + 08
−1.0
NO + O → O2 + N
8.4E + 09
1.6171E + 08
0.0The general conservation equation for the mass fraction of specie “i” (Y
i) can be written as(9)
∂
∂
t
ρ
Y
i
+
∇
·
ρ
v
→
Y
i
=
-
∇
·
J
→
i
+
R
i
,where “ρ” denotes the density and “J
→” and “R” represent the diffusive fluxes and reaction source term of the species “i,” respectively. Based on an established solution obtained from solution to (9) employing the finite rate chemistry model, the deviation of the species compositions and temperature from their thermochemical equilibrium state was determined by allowing the compositions to relax to their chemical equilibrium values [30]. This was undertaken by allowing the reaction source term in (9) to be modeled as(10)
R
i
=
ρ
Y
i
e
q
-
Y
i
τ
c
h
a
r
,where Y
i is the mass fraction of species i and Y
i
e
q is the equilibrium mass fraction of species i (determined as a function of the local species, temperature, and pressure) and a characteristic time scale (τ
c
h
a
r). τ
c
h
a
r in turn is modeled as the sum of a convective/diffusion time scale in the computational cell and a chemical time scale determined from the chemical kinetic mechanism. Equation (10) basically implies that the species reaction towards their chemical equilibrium state over the τ
c
h
a
r period associated with the process.


## 2.4. Turbulence Modeling and Initiating the Onset of Transition

Predicting the onset of transition in this hypersonic flow regime has proven to be really challenging. However, fairly accurate predictions of the surface heat flux profiles have been obtained by judicious specification of the location for the onset of transition [19, 24, 27] with typical prediction accuracies ranging from 20 to 25% [27].The turbulence was modeled in this study employing the realizablek
-
ε turbulence model and the wall temperature of the vehicle was assumed to be constant at 500 K [19]. Since there is a factor of 3 variation among current empirical correlations for the extent and onset of transition (as a function of the edge Mach number) [36], the onset of transition at a distance of 2.5 m along the wall was artificially induced in this study by specifying a roughness height of 5 × 10−4 m along the wall from 2.5 m to 4 m. Wall roughness is one of the driving mechanisms by which disturbances can be induced in the boundary layer to initiate transition [27]. The predicted wall heat flux profiles in the turbulent boundary layer region are shown in Figure 2. It is worth noting that although the prediction accuracies at the peak heating locations are within the range reported in previous studies [19], the goal here is not a rigorous assessment of turbulence or transition prediction methods but on assessing the impacts of chemistry and radiative transfer models. We can also see in Figure 2 that as far as the surface heat fluxes are concerned, results from the different chemistry models (relaxation to chemical equilibrium versus finite rate chemistry) as well as the inclusion of radiative transfer/TRI all yield identical predictions. A complete summary of the modeling options employed in this study is listed in Table 2.Table 2
Summary of physical models employed in the simulations.
Physical models
Primary modeling option
Density
Ideal gas
Specific heat
Piecewise polynomial [29]
Thermal conductivity and mass diffusivity
Kinetic theory
Viscosity
Sutherland law
Wall boundary condition
No slip (fluid flow), constant temperature of 500 K (energy equation)
Far field boundary conditions
Gauge pressure: 2761 pascal
Mach number: 19.97
Turbulence
Realizablek
-
ε
Turbulence initiation
Roughness height of 5 × 10−4 m specified in the wall region 2.5 m to 4 m
Turbulence boundary conditions
Intensity: 0.01%
Viscosity ratio: 2.29E − 09
Radiative transport equation solver
Finite volume radiation model (angular resolution, theta × phi: 3 × 3)
Planck mean absorption coefficient (Kp) [6]∗
3.7516
×
10
-
6
·
(
P
)
1.31
·
exp
⁡
(
5.18
×
10
-
4
T
-
7.13
×
10
-
9
T
2
)
Gas phase kinetics
Table1
Relaxation to chemical equilibrium (kinetic mechanism for chemical time scale determination)
Table1
Turbulence-radiation interactions∗
Equation (7)
∗These models were implemented as user-defined functions (UDFs) in ANSYS FLUENT [30].Figure 2
Wall heat flux profiles in the turbulent boundary layer predicted by the realizablek
-
ε turbulent model (with a specified roughness height).


## 3. Results and Discussion

The results reported in this study correspond to the flow conditions at an altitude of 24.4 km that has been frequently investigated by other studies. The free-stream conditions corresponding to this altitude are adapted from [19] and are summarized in Table 3. Since the goal of this paper is to assess the impact of different chemistry and radiative transfer modeling options on the temperature and specie fields in the laminar and turbulent boundary and shock layers, our initial goal was to establish a shock and boundary layer field across the body that replicated the experimental scenario reasonably well. To ascertain this, we employed experimental measurements as well as simulations of this scenario carried out by previous investigators.Table 3
Free stream conditions (adapted from [19]).
Mach number (M
∞)
19.97
Angle of attack
0 degrees
Density (ρ
∞)
0.043523 kg/m3
Temperature (T
∞)
221.034 K
Pressure (P
∞)
2761.41 N/m2
Turbulent intensity
0.01%
Turbulent viscosity ratio
2.29E − 09


### 3.1. Shock Layer Profiles: Comparison against Previously Reported Solutions [21]

Since there are no previously reported solutions of temperatures and species concentrations across the boundary and shock layers for this geometry, an initial assessment of the validity of our modeling methodology was made by comparing our velocity and density predictions against those reported in the literature. The only published profiles of these variables at this altitude (24.4 km) in the open literature to the author’s best knowledge is the equilibrium air solutions provided by (Langley Aerothermodynamic Upwind Relaxation Algorithm) LAURA [21]. However, a few differences in the modeling methodologies employed by Wood et al. [21] and this study are to be noted:(1)
The calculations employing LAURA [21] were performed with a spherically blunted nose of 3.43 mm, whereas our calculations were performed with a blunted nose of 2.54 mm radius (as shown in Figure 1).
(2)
LAURA solutions in [21] take into consideration the angle of attack of 0.14 degrees, whereas this was neglected in this study. The impact of including this small angle of attack on the numerical predictions has been discussed by Barnhardt and Candler [26].
(3)
Wood et al. [21] assumed a wall temperature of 467.94 K as opposed to 500 K that was considered in this study.It is also worth noting that the calculations employing LAURA were reported only up to a distance of 2.0 m along the plate where the flow conditions were laminar [21].The relative density (ρ/ρ
∞) profiles predicted by the different modeling options across the shock layer in the laminar region are shown in Figure 3 at different axial locations from the nose cone. First, we note clear differences between the modeling scenarios on the relative density profiles. However, these differences reduce downstream of the bow shock along the vehicle body. Second, there are some differences between our equilibrium calculations and the equilibrium solutions reported in [21]. These may be attributed to the differences in the modeling methodology employed by Wood et al. [21] that have been listed above along with the differences in the thermochemical database employed in the calculations. Wood et al. [21] note that 7–10 percent discrepancies in the density profiles can arise due to differences in the thermochemical database alone. However, Figure 3 shows that the agreement between our equilibrium predictions and those employing LAURA improves at larger distances (2.0 m) as the influence of the differences in the nose bluntness goes away.Figure 3
Relative density (ρ/ρ
∞) profiles in the shock layer at different axial locations from the nose cone.The relative tangential velocity (V
t
a
n
g
e
n
t
i
a
l/C
∞) profiles (whereC
∞ is 297.3 m/s) across the shock layer at different axial locations from the nose cone are shown in Figure 4. First, the differences in the modeling scenarios appear to have had a greater impact on the density profiles (cf. Figure 3) than the velocity profiles. A similar observation was reported in Wood et al. [21], where differences in the thermochemical models employed to compute the chemical equilibrium influenced the density profiles more strongly than the velocities. Again, the agreement between our predictions and the predictions reported in [21] is only seen to improve at larger distances away from the nose as the influence of nose bluntness goes away.Figure 4
Relative velocity (V
t
a
n
g
e
n
t
i
a
l/C
∞) profiles in the shock layer at different axial locations from the nose cone.


### 3.2. Boundary Layer Profiles: Comparison against Previously Reported Solutions [21]

Next, the impact of the different modeling scenarios in the boundary layer is investigated to examine the prediction differences more closely. Figures5 and 6 compare the relative density and velocity profiles across the laminar boundary layer at different locations along the body. Similar to the observations in the shock layer, the differences between the modeling scenarios as well as the solutions from LAURA [21] are more apparently close to the nose and diminish at larger distances along the vehicle body. However, the inclusion of radiative transfer in the calculations appears to have a greater impact with an increase in axial distances downstream of the bow shock. This indicates that the magnitude of radiative emission as well as its region of influence grows with axial distance along the vehicle body. These are assessed next.Figure 5
Relative density (ρ/ρ
∞) profiles in the boundary layer at different axial locations from the nose cone.Figure 6
Relative velocity (V
t
a
n
g
e
n
t
i
a
l/C
∞) profiles in the boundary layer at different axial locations from the nose cone.


### 3.3. Radiative Transfer Calculations and the Effect of TRI

The variables of interest in most radiative transfer analysis are the distributions of radiative heat flux vectors (q
(
r
→
)) and the radiative source terms (-
∇
·
q
(
r
→
)). Distributions of the radiative heat flux vectors are necessary to perform accurate energy balances on any solid/liquid interfaces within the solution domain and at the boundaries. The radiative source term describes the conservation of radiative energy within a control volume and goes into the energy balance equation associated with the fluid flow. Therefore, these variables couple radiation with the other physical processes that occur in a multiphysics simulation. The radiative source term represents the radiative energy balance between absorption and emission within a control volume approximated using the same Planck mean absorption coefficient for the absorption and emission term as(11)
-
∇
·
q
r
→
=
k
G
r
→
-
4
π
I
b
r
→
,where I
b represents blackbody emission and G is the incident radiation (representing absorption) computed by integrating the radiative intensity (I) across all directions as(12)
G
r
→
=
∫
4
π
I
r
→
,
s
^
d
Ω
.Previous studies have highlighted the importance of accounting for the self-absorption on the shock layer even when the wall heat transfer is dominated by convection [37–39]. First, radiative emission and self-absorption can decrease the magnitude of the convective fluxes [38] and influence the shock stand-off distance when coupled strongly with ablating flows [39]. In this study the finite volume radiation model (which is a conservative variant of the discrete ordinates method) was employed to estimate (-
∇
·
q).Plots of the radiative source term are shown in Figure7 at two locations each along the laminar and turbulent boundary layers. The magnitude of (-
∇
·
q) at different temperatures and densities has been estimated from shock tube measurements and reported previously by Kivel [40]. At the densities and temperatures encountered in this study, -
∇
·
q was measured to range between 4,000 and 10,000 W/m3 which compares well with the values shown in Figure 7. We also note that the magnitude of -
∇
·
q as well as its “zone of influence” increases along the length of the vehicle. If Γ represents the ratio of radiative heat fluxes to the convective heat fluxes, Γ was found to be 0.001 in the laminar region and 0.003 in the turbulent region.Figure 7
Radiative source (-
∇
·
q) term profiles in the boundary layer (finite rate chemistry and radiation) at different axial locations from the nose cone.Figure8 investigates the effects of TRI at two axial locations along the turbulent boundary layer. The emission enhancement factor is the right hand side of  (7) and quantifies the increase in emission due to the turbulent fluctuations in temperature. Up to a 60% increase in emission is observed at some locations along the turbulent boundary layer. However, Figure 8 shows that the magnitude of the radiative source term and the impact of TRI are not large enough to have any bearing on the turbulent boundary layer.Figure 8
Temperature profiles (a) and the enhancement in emission (b) at two axial locations from the nose cone across the turbulent boundary layer.
(a)
(b)


### 3.4. Species Profiles across the Boundary Layer

The concentrations of the species along the boundary layer can significantly impact the absorption characteristics and shock stand-off distances. Further, accurate characterization of the ablation phenomena requires accurate determination of the species concentrations. It is generally perceived that at large distances from the bow shock, the chemistry within a fluid parcel tends towards their equilibrium concentrations. Temperatures over 4000 K result in significant dissociation of the O molecule. Figure9 shows the O mole fraction profiles in the boundary layer. At a V
∞ of 6 km/s which is close to these reentry conditions, the equilibrium concentration of O downstream of the shock has been reported to be about 0.3 [20] which compares well with the values shown in Figure 9. It is interesting to note the significant differences in the O mole fraction predictions between the two chemistry models with the finite rate chemistry predicting superequilibrium concentrations. The differences between the chemistry models however seem to improve with increase in distance downstream of the bow shock. The NO mole fraction predictions from the two chemistry models are shown in Figure 10. Contrary to the O concentrations, the NO mole fractions do not attain their equilibrium values even well into the turbulent boundary layer.Figure 9
O mole fraction profiles in the boundary layer at different axial locations from the nose cone.Figure 10
NO mole fraction profiles in the boundary layer at different axial locations from the nose cone.


## 3.1. Shock Layer Profiles: Comparison against Previously Reported Solutions [21]

Since there are no previously reported solutions of temperatures and species concentrations across the boundary and shock layers for this geometry, an initial assessment of the validity of our modeling methodology was made by comparing our velocity and density predictions against those reported in the literature. The only published profiles of these variables at this altitude (24.4 km) in the open literature to the author’s best knowledge is the equilibrium air solutions provided by (Langley Aerothermodynamic Upwind Relaxation Algorithm) LAURA [21]. However, a few differences in the modeling methodologies employed by Wood et al. [21] and this study are to be noted:(1)
The calculations employing LAURA [21] were performed with a spherically blunted nose of 3.43 mm, whereas our calculations were performed with a blunted nose of 2.54 mm radius (as shown in Figure 1).
(2)
LAURA solutions in [21] take into consideration the angle of attack of 0.14 degrees, whereas this was neglected in this study. The impact of including this small angle of attack on the numerical predictions has been discussed by Barnhardt and Candler [26].
(3)
Wood et al. [21] assumed a wall temperature of 467.94 K as opposed to 500 K that was considered in this study.It is also worth noting that the calculations employing LAURA were reported only up to a distance of 2.0 m along the plate where the flow conditions were laminar [21].The relative density (ρ/ρ
∞) profiles predicted by the different modeling options across the shock layer in the laminar region are shown in Figure 3 at different axial locations from the nose cone. First, we note clear differences between the modeling scenarios on the relative density profiles. However, these differences reduce downstream of the bow shock along the vehicle body. Second, there are some differences between our equilibrium calculations and the equilibrium solutions reported in [21]. These may be attributed to the differences in the modeling methodology employed by Wood et al. [21] that have been listed above along with the differences in the thermochemical database employed in the calculations. Wood et al. [21] note that 7–10 percent discrepancies in the density profiles can arise due to differences in the thermochemical database alone. However, Figure 3 shows that the agreement between our equilibrium predictions and those employing LAURA improves at larger distances (2.0 m) as the influence of the differences in the nose bluntness goes away.Figure 3
Relative density (ρ/ρ
∞) profiles in the shock layer at different axial locations from the nose cone.The relative tangential velocity (V
t
a
n
g
e
n
t
i
a
l/C
∞) profiles (whereC
∞ is 297.3 m/s) across the shock layer at different axial locations from the nose cone are shown in Figure 4. First, the differences in the modeling scenarios appear to have had a greater impact on the density profiles (cf. Figure 3) than the velocity profiles. A similar observation was reported in Wood et al. [21], where differences in the thermochemical models employed to compute the chemical equilibrium influenced the density profiles more strongly than the velocities. Again, the agreement between our predictions and the predictions reported in [21] is only seen to improve at larger distances away from the nose as the influence of nose bluntness goes away.Figure 4
Relative velocity (V
t
a
n
g
e
n
t
i
a
l/C
∞) profiles in the shock layer at different axial locations from the nose cone.


## 3.2. Boundary Layer Profiles: Comparison against Previously Reported Solutions [21]

Next, the impact of the different modeling scenarios in the boundary layer is investigated to examine the prediction differences more closely. Figures5 and 6 compare the relative density and velocity profiles across the laminar boundary layer at different locations along the body. Similar to the observations in the shock layer, the differences between the modeling scenarios as well as the solutions from LAURA [21] are more apparently close to the nose and diminish at larger distances along the vehicle body. However, the inclusion of radiative transfer in the calculations appears to have a greater impact with an increase in axial distances downstream of the bow shock. This indicates that the magnitude of radiative emission as well as its region of influence grows with axial distance along the vehicle body. These are assessed next.Figure 5
Relative density (ρ/ρ
∞) profiles in the boundary layer at different axial locations from the nose cone.Figure 6
Relative velocity (V
t
a
n
g
e
n
t
i
a
l/C
∞) profiles in the boundary layer at different axial locations from the nose cone.


## 3.3. Radiative Transfer Calculations and the Effect of TRI

The variables of interest in most radiative transfer analysis are the distributions of radiative heat flux vectors (q
(
r
→
)) and the radiative source terms (-
∇
·
q
(
r
→
)). Distributions of the radiative heat flux vectors are necessary to perform accurate energy balances on any solid/liquid interfaces within the solution domain and at the boundaries. The radiative source term describes the conservation of radiative energy within a control volume and goes into the energy balance equation associated with the fluid flow. Therefore, these variables couple radiation with the other physical processes that occur in a multiphysics simulation. The radiative source term represents the radiative energy balance between absorption and emission within a control volume approximated using the same Planck mean absorption coefficient for the absorption and emission term as(11)
-
∇
·
q
r
→
=
k
G
r
→
-
4
π
I
b
r
→
,where I
b represents blackbody emission and G is the incident radiation (representing absorption) computed by integrating the radiative intensity (I) across all directions as(12)
G
r
→
=
∫
4
π
I
r
→
,
s
^
d
Ω
.Previous studies have highlighted the importance of accounting for the self-absorption on the shock layer even when the wall heat transfer is dominated by convection [37–39]. First, radiative emission and self-absorption can decrease the magnitude of the convective fluxes [38] and influence the shock stand-off distance when coupled strongly with ablating flows [39]. In this study the finite volume radiation model (which is a conservative variant of the discrete ordinates method) was employed to estimate (-
∇
·
q).Plots of the radiative source term are shown in Figure7 at two locations each along the laminar and turbulent boundary layers. The magnitude of (-
∇
·
q) at different temperatures and densities has been estimated from shock tube measurements and reported previously by Kivel [40]. At the densities and temperatures encountered in this study, -
∇
·
q was measured to range between 4,000 and 10,000 W/m3 which compares well with the values shown in Figure 7. We also note that the magnitude of -
∇
·
q as well as its “zone of influence” increases along the length of the vehicle. If Γ represents the ratio of radiative heat fluxes to the convective heat fluxes, Γ was found to be 0.001 in the laminar region and 0.003 in the turbulent region.Figure 7
Radiative source (-
∇
·
q) term profiles in the boundary layer (finite rate chemistry and radiation) at different axial locations from the nose cone.Figure8 investigates the effects of TRI at two axial locations along the turbulent boundary layer. The emission enhancement factor is the right hand side of  (7) and quantifies the increase in emission due to the turbulent fluctuations in temperature. Up to a 60% increase in emission is observed at some locations along the turbulent boundary layer. However, Figure 8 shows that the magnitude of the radiative source term and the impact of TRI are not large enough to have any bearing on the turbulent boundary layer.Figure 8
Temperature profiles (a) and the enhancement in emission (b) at two axial locations from the nose cone across the turbulent boundary layer.
(a)
(b)


## 3.4. Species Profiles across the Boundary Layer

The concentrations of the species along the boundary layer can significantly impact the absorption characteristics and shock stand-off distances. Further, accurate characterization of the ablation phenomena requires accurate determination of the species concentrations. It is generally perceived that at large distances from the bow shock, the chemistry within a fluid parcel tends towards their equilibrium concentrations. Temperatures over 4000 K result in significant dissociation of the O molecule. Figure9 shows the O mole fraction profiles in the boundary layer. At a V
∞ of 6 km/s which is close to these reentry conditions, the equilibrium concentration of O downstream of the shock has been reported to be about 0.3 [20] which compares well with the values shown in Figure 9. It is interesting to note the significant differences in the O mole fraction predictions between the two chemistry models with the finite rate chemistry predicting superequilibrium concentrations. The differences between the chemistry models however seem to improve with increase in distance downstream of the bow shock. The NO mole fraction predictions from the two chemistry models are shown in Figure 10. Contrary to the O concentrations, the NO mole fractions do not attain their equilibrium values even well into the turbulent boundary layer.Figure 9
O mole fraction profiles in the boundary layer at different axial locations from the nose cone.Figure 10
NO mole fraction profiles in the boundary layer at different axial locations from the nose cone.


## 4. Conclusions

A methodology for performing fully coupled, multidimensional calculations of finite rate chemistry, radiative heat transfer, and turbulence radiation interactions (TRI) in hypersonic reentry scenarios in a computationally efficient manner was demonstrated. By considering the Mach 20 Reentry F vehicle as a case study, the impact of different modeling options on density, velocity, temperature, and O and NO predictions in the shock and boundary layers was assessed. Radiative transfer calculations were performed using the finite volume radiation model employing add-on functions to compute a Planck mean absorption coefficient and the temperature self-correlation term (for TRI effects) in the optically thin shock layer. The onset of transition was induced by specifying a wall roughness height beyond the experimentally observed transition location (2.5 m). The chemistry was modeled employing 8 elementary reactions and an equilibrium approach where the species were allowed to relax towards their chemical equilibrium values over the characteristic time scale associated with the process. From the results of this study, the following conclusions may be drawn:(1)
The wall heat fluxes in the turbulent region compared reasonably well against measurements and were not significantly impacted by the choice of modeling options investigated in this study.
(2)
The density and velocity profiles across the boundary and shock layers compared well against similar calculations reported in previous studies [21]. Density predictions were more sensitive to the choice of modeling options than the velocities.
(3)
The magnitude of the radiative source term employing the Planck mean absorption coefficient compared well against its measurements deduced from shock tube experiments [40]. The emission contribution to the radiative source term was more than three orders of magnitude greater than the absorption term, further justifying the usage of Planck mean absorption coefficient and modeling only the temperature self-correlation term to account for TRI.
(4)
The TRI model predicted up to a 60% enhancement in emission due to temperature fluctuations in the turbulent boundary layer. Further, the value ofΓ (the ratio of radiative heat fluxes to the convective heat fluxes) increased from 0.001 in the laminar region to 0.003 in the turbulent region. However, including the effects of radiation or TRI did not have any bearing on the temperature predictions in the boundary layer for this reentry scenario.
(5)
While the variations in density and velocities between the different chemistry modeling options diminished along the length of the body, the O and NO prediction variations extended well into the turbulent boundary layer.



---

*Source: 8401249-2016-06-09.xml*