Three-Temperature Model Applied to Thermochemical Non-Equilibrium Reentry Flows in 2D—Seven Species ()

Edisson Sávio de Góes Maciel^{}

Department of Mathematics, University of Pernambuco, Recife, Brazil.

**DOI: **10.4236/ojfd.2022.121001
PDF
HTML XML
188
Downloads
816
Views
Citations

Department of Mathematics, University of Pernambuco, Recife, Brazil.

In this work, a study involving the fully coupled Euler and Navier-Stokes reactive equations is performed. These equations, in conservative and finite volume contexts, employing structured spatial discretization, on a condition of thermochemical non-equilibrium, are analyzed. High-order studies are accomplished using the Spectral method of Streett, Zang, and Hussaini. The high enthalpy hypersonic flows around a circumference, around a reentry capsule, along a blunt body, and along a double ellipse in two-dimensions are simulated. The Van Leer, Liou and Steffen Jr., and Steger and Warming flux vector splitting algorithms are applied to execute the numerical experiments. Three temperatures, which are the translational-rotational temperature, the vibrational temperature, and the electron temperature, are used to accomplish the numerical comparisons. Excellent results were obtained with minimum errors inferior to 6.0%. The key contribution of this work is the correct implementation of a three temperature model coupled with the implementation of three algorithms to perform the numerical simulations, as well the description of energy exchange mechanisms to perform more realistic simulations.

Keywords

Thermochemical Non-Equilibrium, Three-Temperature Model, Van Leer Scheme, Liou and Steffen Jr. Scheme, Steger and Warming Scheme

Share and Cite:

de Góes Maciel, E. (2022) Three-Temperature Model Applied to Thermochemical Non-Equilibrium Reentry Flows in 2D—Seven Species. *Open Journal of Fluid Dynamics*, **12**, 1-35. doi: 10.4236/ojfd.2022.121001.

1. Introduction

Several conceptual designs for vehicles that would fly in the atmosphere at hypersonic speeds have been developed recently [1]. For typical flight conditions, the air that envelops these vehicles is chemically reacted, vibrationally excited, and ionized. These reactions and excitation processes occur at rates similar to the rate of fluid motion, which results in a state of thermochemical non-equilibrium. The thermochemical state of the gas influences its dynamics and affects the aerodynamic forces and heat transfer acting on the vehicle.

The motivation for multi-temperature model studies often originates for ionized flows owing to the wide disparities in the masses of the constituent species [2]. In this regard, the accuracy of macroscopic temperature models, such as the Landau-Teller translation-vibration (T-V) model, can lead to questioning of the underlying approximations. Such questioning in the past invariably has lead to the development of gas kinetic schemes, such as the Bhatnagar-Gross-Krook scheme [3], to describe the energy exchanges. In addition, the extended Navier-Stokes equations [4] are often used for high-altitude hypersonic flows to accurately capture the shock and provide more information than the classical Navier-Stokes equations can provide. However, the subjective question of the effect of the various energy exchanges using approximate macroscopic temperature models on the shock structure of hypersonic flows can be answered with greater ease if the analysis is made for flows within the continuum regime.

Spectral methods are considered a class of solution techniques using sets of known functions to solve differential equations [5]. Such methods are generally considered high order and capable of obtaining solutions with a high resolution. Unlike finite-difference and finite-element methods, spectral methods utilize an expansion in terms of global, rather than local, basis functions to represent the solution of a differential equation. When properly applied, these techniques accurately resolve phenomena on the scale of the mesh spacing. The order of truncation error decay with mesh refinement is also higher than which can be achieved with finite-difference and finite-element methods. For problems with smooth solutions, it is possible to produce spectral method whose truncation error goes to zero as faster than any finite power of the mesh spacing (exponential convergence).

In this work, a study involving the fully coupled Euler and Navier-Stokes reactive equations is performed. These equations, in conservative and finite volume contexts, employing structured spatial discretization, on a condition of thermochemical non-equilibrium, are analyzed. High-order studies are accomplished using the spectral method of [6]. The high enthalpy hypersonic flows around a circumference, around a reentry capsule, around a blunt body, and around a double ellipse in two-dimensions are simulated. The [7] [8] [9] flux vector splitting algorithms are applied to execute the numerical experiments. The Euler back-ward integration method is employed to march the schemes in time. The convergence process is accelerated to steady state condition through a spatially variable time step procedure, which has proved effective gains in terms of computational acceleration (see [10] [11] ). The reactive simulations involve Earth atmosphere chemical model of seven species and eighteen reactions, based on the [12] model. Three temperatures, which are the translational-rotational temperature, the vibrational temperature, and the electron temperature, are used to accomplish the numerical comparisons. Excellent results were obtained with minimum errors inferior to 6.0%. Moreover, the implementation of a three temperature model has a great impact to the CFD community in terms of numerical tools available to study thermochemical non-equilibrium flows in two-dimensions more realistically.

The main contribution and originality of this paper is concerned with the three-temperature model studied in it. As [21] has explained, the use of a three temperature model is more appropriated to describe the flow field in the case of numerical simulations of fluid dynamics. The translational and vibrational temperatures, case of two temperature models, are not sufficient to highlight the main features of the flowfield and the use of a third temperature, the electron temperature, is fundamental to describe more correctly the flow physics for the cases of weakly ionization, as is the case. On the other hand, this paper also incorporates the implementation of three upwind numerical algorithms on a context of finite volumes discretization and reactive equations and describes them with details. Although the numerical algorithms were not new for the CFD community, their implementations for such discretization and for reactive formulation is also a original contribution of this work. Moreover, the description of the most important energy exchanges between molecules and atoms, as well electrons, are also pointed out as key mechanisms to be considered in a more realistic scenario of fluid flow. The key contribution of this work is the correct implementation of a three-temperature model coupled with the implementation of three algorithms to perform the numerical simulations, as well the description of energy exchange mechanisms to perform more realistic simulations.

Small number of papers related to the three-temperature and multi-temperature models is available in the CFD literature. With it in mind, this paper aims to fill this space and to provide efficient tools to perform such numerical simulations. For example, the majority of papers lead to one- and two-temperature models and with a minor number of exchange mechanisms to describe the energy change between atoms, molecules, and electrons. The three-temperature model is more involved with the physics and chemistry of the reactive Euler and Navier-Stokes equations, obtaining better results than the one- and two-temperature counterparts. Moreover, the three algorithms studied in this paper have description for perfect gas formulation, being rare the implementation to models of one- and two-temperatures. The present manuscript is an original and detailed description of their implementations for reactive gas flow. The use of a spectral method for the reactive Euler and Navier-Stokes equations is also the extension of the perfect gas formulation to a reactive gas flow one. All of these requisites point out to the original contribution of this work to the CFD state-of-art in terms of non-equilibrium thermal and chemical modeling and reentry flows.

2. Spectral Method

On this current work, a spectral method is applied to assure that high order solutions are obtained by the numerical algorithms. Two classes of techniques for spectral discretization are referred to as tau and collocation methods [6] [13]. The latter technique is used here. In this scheme, the approximation series is determined by satisfying the differential equation exactly at a set of distinct collocation points. The locations of these points in the domain are linked to the choice of basis function. In this study, arbitrary collocation points are implemented. The collocation method is used here since enforcement of boundary conditions and evaluations of nonlinear terms are straightforward. Additionally, some accuracy advantage is seen in the collocation method over the tau method for a number of problems [6] [13]. The series expansion for a function *Q*(*x*) may be defined by

${Q}_{N}\left(x\right)={\displaystyle {\sum}_{n=0}^{\infty}{\stackrel{^}{Q}}_{n}{B}_{n}\left(x\right)}$, (1)

where *B _{n}*(

${B}_{n}\left(x\right)={T}_{n}\left(x\right)=2x{P}_{n-1}\left(x\right)-{P}_{n-2}\left(x\right),\text{\hspace{0.17em}}n\ge 2$, (2)

with: *P*_{0}(*x*) = 1 and *P*_{1}(*x*) = *x*. The Chebyshev-Gauss-Lobatto standard collocation points are defined by:

${x}_{l}=\mathrm{cos}\left(\frac{\pi l}{N}\right)$, $l=0,1,\cdots ,N$. (3)

The Chebyshev collocation points result from a simple change of variables, which relates the Chebyshev interpolation series to a Fourier cosine series [6, 13]. To evaluate the ${\stackrel{^}{Q}}_{n}$, the inverse relation is required. This is

${\stackrel{^}{Q}}_{n}={\stackrel{^}{c}}_{n}{\displaystyle {\sum}_{l=0}^{N}{w}_{l}{B}_{n}\left({x}_{l}\right){Q}_{i,j}}$, $n=0,1,\cdots $, (4)

with *w _{l}* being a normalized weighting function and
${\stackrel{^}{c}}_{n}$ a constant. These variables assume the following expressions to a Chebyshev-Gauss-Lobatto interpolation:

${\stackrel{^}{c}}_{n}=\frac{2}{N{\stackrel{\xaf}{c}}_{n}}$, where: ${\stackrel{\xaf}{c}}_{n}=\{\begin{array}{l}2,\text{\hspace{0.17em}}n=0\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{or}\text{\hspace{0.17em}}N\\ 1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}0<n<N\end{array}$ and ${w}_{l}=\frac{1}{{\stackrel{\xaf}{c}}_{l}}$. (5)

Legendre collocation is based on using Legendre polynomials as the basis function in Equation (1), e.g.,

${B}_{n}\left(x\right)=\frac{\left[\left(2n-1\right)x{P}_{n-1}\left(x\right)-\left(n-1\right){P}_{n-2}\left(x\right)\right]}{n},\text{\hspace{0.17em}}n\ge 2$, (6)

where: *P*_{0}(*x*) = 1 and *P*_{1}(*x*) = *x*. Interpolation via Legendre series cannot easily be related to trigonometric interpolation, so there is no simple expression to evaluate the
${\stackrel{^}{Q}}_{n}$ coefficients. Appeal must be made to the theory of numerical quadrature to form an approximation to the integrals which result from analytic Legendre interpolation [14]. Considering Equation (4), the normalized weights and constant of the Legendre-Gauss-Lobatto collocation points are

${w}_{l}=\frac{1}{N\left(N+1\right){B}_{N}^{2}\left({x}_{l}\right)}$ and ${\stackrel{^}{c}}_{n}=\{\begin{array}{l}2n+1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}n=0,1,\cdots ,\left(N-1\right)\\ N,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}n=N\end{array}$. (7)

In this current work, it was assumed that the Legendre-Gauss-Lobatto collocation points are the same as the Chebyshev-Gauss-Lobatto ones. It was also adopted the following collocation points and normalized weight for the Chebyshev-Gauss-Radau interpolation, based on the work of [15]:

${x}_{l}=\mathrm{cos}\left(\frac{2\pi l}{2N+1}\right)$, (8)

${w}_{l}=\{\begin{array}{l}\frac{N}{2N+1},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}l=0\\ \frac{N}{N+1},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{elsewhere}\end{array}$. (9)

For the Legendre-Gauss-Radau interpolation, also based in [15], the collocation points are defined by Equation (8) and the normalized weights are described by:

${w}_{l}=\{\begin{array}{l}\frac{1}{{\left(N+1\right)}^{2}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}l=0\\ \frac{1}{2{\left(N+1\right)}^{2}}\times \frac{1-{x}_{l}}{{B}_{N}^{2}\left({x}_{l}\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{elsewhere}\end{array}$. (10)

The same application to the vector of conserved variables *Q *is applied to the vector of flux *C*, to be defined in section 5. The proposed method has two collocation point options and two normalized weight functions to be considered by the Chebyshev and Legendre methods: Chebyshev-Gauss-Radau, Chebyshev-Gauss-Lobatto, Legendre-Gauss-Radau and Legendre-Gauss-Lobatto.

3. Reactive Navier-Stokes in 2D

As the Navier-Stokes equations tend to the Euler equations when high Reynolds number are employed, only the former equations are presented. The reactive Navier-Stokes equations in thermochemical non-equilibrium, where the rotational and vibrational contributions are considered, were implemented on conservative and finite volume contexts, in the two-dimensional space. In this case, these equations in integral and conservative forms can be expressed by:

$\frac{\partial}{\partial t}{\displaystyle {\int}_{V}Q\text{d}V}+{\displaystyle {\int}_{S}F\cdot n\text{d}S}={\displaystyle {\int}_{V}{S}_{CVE}\text{d}V}$, with: $F=\left({E}_{e}-{E}_{v}\right)i+\left({F}_{e}-{F}_{v}\right)j$, (11)

where: *Q *is the vector of conserved variables, *V* is the volume of a computational cell,
$F$ is the complete flux vector,
$n$ is the unity vector normal to the flux face, *S* is the flux area, *S _{CVE}* is the chemical, vibrational, and electron source term,

$Q=\left\{\begin{array}{c}\rho \\ \rho u\\ \rho v\\ e\\ {\rho}_{1}\\ {\rho}_{2}\\ {\rho}_{4}\\ {\rho}_{5}\\ {\rho}_{6}\\ {\rho}_{7}\\ \rho {e}_{V}\\ {\rho}_{e}{e}_{e}\end{array}\right\}$, ${E}_{e}=\left\{\begin{array}{c}\rho u\\ \rho {u}^{2}+p\\ \rho uv\\ \rho Hu\\ {\rho}_{1}u\\ {\rho}_{2}u\\ {\rho}_{4}u\\ {\rho}_{5}u\\ {\rho}_{6}u\\ {\rho}_{7}u\\ \rho {e}_{V}u\\ \rho {e}_{e}u\end{array}\right\}$, ${F}_{e}=\left\{\begin{array}{c}\rho v\\ \rho uv\\ \rho {v}^{2}+p\\ \rho Hv\\ {\rho}_{1}v\\ {\rho}_{2}v\\ {\rho}_{4}v\\ {\rho}_{5}v\\ {\rho}_{6}v\\ {\rho}_{7}v\\ \rho {e}_{V}v\\ \rho {e}_{e}v\end{array}\right\}$ ; (12)

in which: *ρ* is the mixture density; *ρ _{e}* is the electron’s density;

${E}_{v}=\frac{1}{Re}\left\{\begin{array}{c}0\\ {\tau}_{xx}\\ {\tau}_{xy}\\ {f}_{x}-{\varphi}_{x}\\ -{\rho}_{1}{v}_{1x}\\ -{\rho}_{2}{v}_{2x}\\ -{\rho}_{4}{v}_{4x}\\ -{\rho}_{5}{v}_{5x}\\ -{\rho}_{6}{v}_{6x}\\ -{\rho}_{7}{v}_{7x}\\ {q}_{v,x}-{\varphi}_{v,x}\\ {f}_{e,x}-{\varphi}_{e,x}\end{array}\right\}$, ${F}_{v}=\frac{1}{Re}\left\{\begin{array}{c}0\\ {\tau}_{xy}\\ {\tau}_{yy}\\ {f}_{y}-{\varphi}_{y}\\ -{\rho}_{1}{v}_{1y}\\ -{\rho}_{2}{v}_{2y}\\ -{\rho}_{4}{v}_{4y}\\ -{\rho}_{5}{v}_{5y}\\ -{\rho}_{6}{v}_{6y}\\ -{\rho}_{7}{v}_{7y}\\ {q}_{v,y}-{\varphi}_{v,y}\\ {f}_{e,y}-{\varphi}_{e,y}\end{array}\right\}$ ; (13)

*Q _{T}*

${S}_{CVE}=\left\{\begin{array}{c}0\\ 0\\ 0\\ 0\\ {\stackrel{\dot{}}{\omega}}_{1}\\ {\stackrel{\dot{}}{\omega}}_{2}\\ {\stackrel{\dot{}}{\omega}}_{4}\\ {\stackrel{\dot{}}{\omega}}_{5}\\ {\stackrel{\dot{}}{\omega}}_{6}\\ {\stackrel{\dot{}}{\omega}}_{7}\\ {Q}_{T\text{-}V}+{\displaystyle {\sum}_{s=mol}{\stackrel{\dot{}}{\omega}}_{s}{e}_{v,s}}\\ {S}_{e,pg}+{S}_{e,trans}+{S}_{e,inel}+{S}_{e,chem}\end{array}\right\}$. (14)

The viscous stresses, in N/m^{2}, are determined, according to a Newtonian fluid model, by:

$\begin{array}{l}{\tau}_{xx}=2{\mu}_{m}\frac{\partial u}{\partial x}-\frac{2}{3}{\mu}_{m}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right);\\ {\tau}_{xy}={\mu}_{m}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right);\\ {\tau}_{yy}=2{\mu}_{m}\frac{\partial v}{\partial y}-\frac{2}{3}{\mu}_{m}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right).\end{array}$ (15)

where *µ _{m}* is the molecular viscosity. Expressions to

${f}_{x}={\tau}_{xx}u+{\tau}_{xy}v+{q}_{x}+{q}_{e,x}+{q}_{v,x}$ ; (16)

${f}_{y}={\tau}_{xy}u+{\tau}_{yy}v+{q}_{y}+{q}_{e,y}+{q}_{v,y}$, (17)

where *q _{x}* and

${q}_{x}=k\frac{\partial T}{\partial x}$ and ${q}_{y}=k\frac{\partial T}{\partial y}$, (18)

where: *k* is the thermal conductivity due to translation and rotation. The *q _{v}*

${q}_{v,x}={k}_{v}\frac{\partial {T}_{v}}{\partial x}$ and ${q}_{v,y}={k}_{v}\frac{\partial {T}_{v}}{\partial y}$, (19)

with *k _{v}* being the mixture-vibrational-thermal conductivity, and

${q}_{e,x}={k}_{e}\frac{\partial {T}_{e}}{\partial x}$ and ${q}_{e,y}={k}_{e}\frac{\partial {T}_{e}}{\partial y}$, (20)

where *T _{e}* is the electron’s temperature and

${\rho}_{s}{v}_{sx}=-\rho {D}_{s}\frac{\partial {Y}_{MF,s}}{\partial x}$ and ${\rho}_{s}{v}_{sy}=-\rho {D}_{s}\frac{\partial {Y}_{MF,s}}{\partial y}$, (21)

with “*s*” referent to a given species, *Y _{MF}*

${Y}_{MF,s}=\frac{\frac{{\rho}_{s}}{{M}_{s}}}{{{\displaystyle \sum}}_{k=1}^{ns}\frac{{\rho}_{k}}{{M}_{k}}}$ (22)

and *D _{s}* is the species-effective-diffusion coefficient. “

${\varphi}_{x}={{\displaystyle \sum}}_{s=1}^{ns}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\rho}_{s}{v}_{sx}{h}_{s}$ and ${\varphi}_{y}={{\displaystyle \sum}}_{s=1}^{ns}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\rho}_{s}{v}_{sy}{h}_{s}$, (23)

being *h _{s}* the specific enthalpy (sensible) of the chemical species “

${\varphi}_{v,x}={\displaystyle {\sum}_{s=mol}{\rho}_{s}{v}_{sx}{h}_{v,s}}$ and ${\varphi}_{v,y}={\displaystyle {\sum}_{s=mol}{\rho}_{s}{v}_{sy}{h}_{v,s}}$, (24)

with *h _{v}*

$Re=\frac{{\rho}_{char}{V}_{initial}{L}_{REF}}{{\mu}_{m,char}}$, (25)

with “*char*” related to characteristic or freestream variables, *V _{initial}* is the flow initial velocity, and

For the electron energy equation, the following expressions are defined to *f _{e}*

${f}_{e,x}={\tau}_{xx,e}u+{\tau}_{xy,e}v+{q}_{e,x}$ and ${f}_{e,y}={\tau}_{xy,e}u+{\tau}_{yy,e}v+{q}_{e,y}$, (26)

where the electron viscous stresses are calculated with the Equation (15), being the electron’s molecular viscosity used in place of the mixture molecular viscosity. The electron’s diffusion flux is given by:

${\varphi}_{e,x}={\rho}_{e}{v}_{ex}{e}_{e}$ and ${\varphi}_{e,y}={\rho}_{e}{v}_{ey}{e}_{e}$. (27)

The caloric equation of state for the mass-averaged gas is derived subtracting the vibrational, kinetic, electron, and chemical energies from the total energy to yield the energy in the translational and rotational modes [20]. Assuming that the rotational energy modes are in equilibrium with the translational modes, it is possible to determine the translational-rotational temperature *T*,

$\sum}_{r\ne e}{\rho}_{r}{c}_{v,r}T}=e-\frac{1}{2}{\displaystyle {\sum}_{r\ne e}{\rho}_{r}\left({u}^{2}+{v}^{2}\right)}-{\displaystyle {\sum}_{s=mol}{\rho}_{s}{e}_{v,s}}-{\rho}_{e}{e}_{e}-{\displaystyle {\sum}_{r}{\rho}_{r}{h}_{r}^{0$, (28)

where *h*^{0} is the formation enthalpy. The pressure is the sum of the partial pressures,

$p={\displaystyle {\sum}_{r\ne e}{\rho}_{r}{R}_{r}T}+{p}_{e}$, (29)

being *p _{e}* the electron’s pressure. The equations of state for the electrons are

${\rho}_{e}{c}_{v,e}{T}_{e}={\rho}_{e}{e}_{e}-\frac{1}{2}{\rho}_{e}\left({u}^{2}+{v}^{2}\right)$ and ${p}_{e}={\rho}_{e}{R}_{e}{T}_{e}$. (30)

It was assumed that the energy contained in the excited electronic states of the molecules is negligible relative to the energy contained in the other modes.

4. Energy Exchange Mechanisms

The energy exchange mechanisms that appear on the *S _{CVE}* term must be modeled. The models are described below, based on [20] work.

4.1. Translation-Vibration Energy Exchanges

The rates of energy exchange between the translational and vibrational modes *Q _{T}*

${Q}_{T\text{-}Vs}={\rho}_{s}\frac{{e}_{vs}^{*}\left(T\right)-{e}_{vs}}{\langle {\tau}_{s}\rangle}$, (31)

with *T* the translational-rotational temperature. The Landau-Teller relaxation time is given by an expression from [21],

$\langle {\tau}_{s}\rangle =\frac{{{\displaystyle \sum}}_{r}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{X}_{r}}{{{\displaystyle \sum}}_{r}\frac{{X}_{r}}{{\tau}_{sr}}}$, for $r\ne e$, (32)

where *X _{r}* is the mole fraction and

${\tau}_{sr}=\frac{B}{{p}_{s}}{\text{e}}^{\left[A\left({T}^{-1/3}-0.015{\mu}_{sr}^{0.25}\right)-18.42\right]}$, (33)

with *B* = 1.013 × 10^{5} Pa.s, *p _{s}* is the species pressure in Pa,

4.2. Electron Pressure Gradient, *S _{e}*

This term is an approximation to the work done on electrons by the electric field induced by the electron pressure gradient [22]. Its expression is given by:

${S}_{e,pg}=-{p}_{e}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)$. (34)

4.3. Translational-Electron Relaxation, *S _{e}*

Energy exchange between translational and electron energies [20] [22]. This term is defined by:

${S}_{e,trans}=3{\rho}_{e}{R}_{univ}\left(T-{T}_{e}\right)\sqrt{\frac{8{R}_{e}{T}_{e}}{\pi}}{\displaystyle {\sum}_{r\ne e}\frac{{\rho}_{r}{N}_{Av}}{{M}_{r}^{2}}{\sigma}_{er}}$ (35)

where *σ _{er}*,

${\sigma}_{er}={a}_{r}+{b}_{r}{T}_{e}+{c}_{r}{T}_{e}^{2}$, (36)

with model constants defined in Table 1. For the case of electron-ion interactions, the effective Coulomb cross section is determined by:

${\sigma}_{e,ions}=\frac{8\pi}{27}\frac{{e}^{4}}{{\left(k{T}_{e}\right)}^{2}}\mathrm{ln}\left[1+\frac{9{\left(k{T}_{e}\right)}^{3}}{4\pi {N}_{e}{e}^{6}}\right]$, (37)

with “*e*” being the electron charge (1.6022 × 10^{−19} C), *N _{e}* is the electron’s density number, and

4.4. Rotational-Electron Relaxation, *S _{inel}*

This term is related to the inelastic energy exchange between electrons and molecules [22]. The energy relaxation between the rotational and electron energies should be considered in the electron energy equation because of the electron interactions with molecular multipoles. In order to simplify the relaxation term, an energy transfer rate factor, *g _{rot}*

${S}_{e,rot}=3{\rho}_{e}{R}_{univ}\left(T-{T}_{e}\right)\sqrt{\frac{8{R}_{e}{T}_{e}}{\pi}}{\displaystyle {\sum}_{r\ne e}{g}_{rot,r}\frac{{\rho}_{r}{N}_{Av}}{{M}_{r}^{2}}{\sigma}_{er}}$. (38)

For neutral species, the energy transfer rate factor is listed in Table 2. For molecular ions, we assume they have the same rate factor as their neutral molecules.

4.5. Electron Energy Due to Chemical Reactions

The chemistry related term gives added or removed electron energy by chemical reactions and can be expressed as:

${S}_{e,chem}={\stackrel{\dot{}}{\omega}}_{e}{e}_{e}$. (39)

Table 1. The coefficients of the electron-neutral collisional cross section.

Table 2. The transfer energy rate factors of rotational-electron energy relaxation.

5. Numerical Algorithms

Considering the two-dimensional and structured case, the flux vector splitting algorithms follow that described in [7] [8] [9] [18] [19].

5.1. Van Leer and Liou and Steffen Jr. Schemes

The system is solved in three parts separately, according to [23]. The first part takes into account the dynamic part, which considers the Navier-Stokes equations plus the electron energy equation, the second one takes into account the chemical part involving the chemical contributions, and finally, the third part considers only the vibrational contribution. Hence, the discrete-dynamic-convective flux, which solves the dynamic part, is given by:

$\begin{array}{c}{R}_{i+1/2,j}^{Dyn}={\left|S\right|}_{i+1/2,j}\{\frac{1}{2}{M}_{i+1/2,j}\left[{\left(\begin{array}{c}\rho a\\ \rho au\\ \rho av\\ \rho aH\\ \rho a{e}_{e}\end{array}\right)}_{L}+{\left(\begin{array}{c}\rho a\\ \rho au\\ \rho av\\ \rho aH\\ \rho a{e}_{e}\end{array}\right)}_{R}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.05em}}-\frac{1}{2}{\varphi}_{i+1/2,j}\left[{\left(\begin{array}{c}\rho a\\ \rho au\\ \rho av\\ \rho aH\\ \rho a{e}_{e}\end{array}\right)}_{R}-{\left(\begin{array}{c}\rho a\\ \rho au\\ \rho av\\ \rho aH\\ \rho a{e}_{e}\end{array}\right)}_{L}\right]\}+{\left(\begin{array}{c}0\\ {S}_{x}p\\ {S}_{y}p\\ 0\\ 0\end{array}\right)}_{i+1/2,j}\end{array}$ (40)

the discrete-chemical-convective flux is defined by:

$\begin{array}{c}{R}_{i+1/2,j}^{Chem}={\left|S\right|}_{i+1/2,j}\{\frac{1}{2}{M}_{i+1/2,j}\left[{\left(\begin{array}{c}{\rho}_{1}a\\ {\rho}_{2}a\\ {\rho}_{4}a\\ {\rho}_{5}a\\ {\rho}_{6}a\\ {\rho}_{7}a\end{array}\right)}_{L}+{\left(\begin{array}{c}{\rho}_{1}a\\ {\rho}_{2}a\\ {\rho}_{4}a\\ {\rho}_{5}a\\ {\rho}_{6}a\\ {\rho}_{7}a\end{array}\right)}_{R}\right]\\ \text{\hspace{0.17em}}-\frac{1}{2}{\varphi}_{i+1/2,j}\left[{\left(\begin{array}{c}{\rho}_{1}a\\ {\rho}_{2}a\\ {\rho}_{4}a\\ {\rho}_{5}a\\ {\rho}_{6}a\\ {\rho}_{7}a\end{array}\right)}_{R}-{\left(\begin{array}{c}{\rho}_{1}a\\ {\rho}_{2}a\\ {\rho}_{4}a\\ {\rho}_{5}a\\ {\rho}_{6}a\\ {\rho}_{7}a\end{array}\right)}_{L}\right]\}\end{array}$ (41)

and finally the discrete-vibrational-convective flux is given by:

$\begin{array}{c}{R}_{i+1/2,j}^{Vib}={\left|S\right|}_{i+1/2,j}\{\frac{1}{2}{M}_{i+1/2,j}\left[{\left(\rho {e}_{v}a\right)}_{L}+{\left(\rho {e}_{v}a\right)}_{R}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.05em}}-\frac{1}{2}{\varphi}_{i+1/2,j}\left[{\left(\rho {e}_{v}a\right)}_{R}-{\left(\rho {e}_{v}a\right)}_{L}\right]\}\end{array}$ (42)

where: ${S}_{i+1/2,j}={\left[{S}_{x}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{S}_{y}\right]}_{i+1/2,j}^{\text{T}}$ defines the normal area vector for the surface ( $i+1/2,j$ ), and the speed of sound is defined by the following expression:

$a=\sqrt{\left(\beta +1\right)\frac{p}{\rho}}$, (43)

with *β* a parameter to be defined, calculated at each interaction. The normal area components *S _{x}* and

The same definitions presented in [7] [8] [18] [19] are valid to these algorithms. The definition of the dissipation term $\varphi $ determines the particular formulation of the convective fluxes. The choice below corresponds to the [7] scheme, according to [24]:

${\varphi}_{i+1/2,j}={\varphi}_{i+1/2,j}^{VL}=\{\begin{array}{l}\left|{M}_{i+1/2,j}\right|,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\left|{M}_{i+1/2,j}\right|\ge 1;\\ \left|{M}_{i+1/2,j}\right|+0.5{\left({M}_{R}-1\right)}^{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{\hspace{0.05em}}0\le {M}_{i+1/2,j}<1;\\ \left|{M}_{i+1/2,j}\right|+0.5{\left({M}_{L}+1\right)}^{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{\hspace{0.05em}}-1<{M}_{i+1/2,j}\le 0.\end{array}$ (44)

and the [8] scheme is obtained by, according to [24]:

${\varphi}_{i+1/2,j}={\varphi}_{i+1/2,j}^{LS}=\left|{M}_{i+1/2,j}\right|$, (45)

with *M* being the Mach number and *R* and *L *right and left states. Both schemes

Table 3. Values of *S _{x}* and

Figure 1. Computational cell.

are first-order accurate in space and in time. The high-order spatial accuracy is obtained, in the current study, by the spectral method.

5.2. Steger and Warming Scheme

The [9] scheme was also implemented in the present work. The *U* contravariant velocity defined by the [9] scheme is

$U=u{n}_{x}+v{n}_{y}$. (46)

The eigenvalues of the Euler equations are defined as

${\lambda}_{1}=U$, ${\lambda}_{2}=U+a$ and ${\lambda}_{3}=U-a$. (47)

The separated eigenvalues, based on the work of [9], which defines positive and negative contributions, can be expressed as

${\lambda}_{1}^{+}=\frac{1}{2}\left({\lambda}_{1}+\sqrt{{\lambda}_{1}^{2}+{\epsilon}^{2}}\right)$, ${\lambda}_{2}^{+}=\frac{1}{2}\left({\lambda}_{2}+\sqrt{{\lambda}_{2}^{2}+{\epsilon}^{2}}\right)$ ; (48)

${\lambda}_{3}^{+}=\frac{1}{2}\left({\lambda}_{3}+\sqrt{{\lambda}_{3}^{2}+{\epsilon}^{2}}\right)$, ${\lambda}_{1}^{-}=\frac{1}{2}\left({\lambda}_{1}-\sqrt{{\lambda}_{1}^{2}+{\epsilon}^{2}}\right)$ ; (49)

${\lambda}_{2}^{-}=\frac{1}{2}\left({\lambda}_{2}-\sqrt{{\lambda}_{2}^{2}+{\epsilon}^{2}}\right)$, ${\lambda}_{3}^{-}=\frac{1}{2}\left({\lambda}_{3}-\sqrt{{\lambda}_{3}^{2}+{\epsilon}^{2}}\right)$, (50)

where *ε* is a small parameter, adopted as 0.04.

The coefficients of the [9] scheme, which defines the separated numerical flux vectors, are determined in Table 4. Six coefficients were defined. These coefficients were prescribed in terms of thermodynamic and sound speed variables.

The numerical flux vector was decomposed in three components, based on the dynamic, chemical, and vibrational parcels of the flow. The flux vector to each parcel is separated in positive and negative contributions, based on the homogeneity property of the Euler flux vectors. For the dynamic parcel, where it was included the electron energy equation, one has:

Table 4. Coefficients of the Steger and Warming scheme.

${F}_{Dyn}^{\pm}=\left\{\begin{array}{c}{C}_{1}{\lambda}_{1}^{\pm}+{C}_{2}{\lambda}_{2}^{\pm}+{C}_{2}{\lambda}_{3}^{\pm}\\ {C}_{1}{\lambda}_{1}^{\pm}u+{C}_{2}{\lambda}_{2}^{\pm}\left(u+a{n}_{x}\right)+{C}_{2}{\lambda}_{3}^{\pm}\left(u-a{n}_{x}\right)\\ {C}_{1}{\lambda}_{1}^{\pm}v+{C}_{2}{\lambda}_{2}^{\pm}\left(v+a{n}_{y}\right)+{C}_{2}{\lambda}_{3}^{\pm}\left(v-a{n}_{y}\right)\\ {C}_{3}{\lambda}_{1}^{\pm}+{C}_{4}{\lambda}_{2}^{\pm}+{C}_{5}{\lambda}_{3}^{\pm}\\ \left({C}_{1}{\lambda}_{1}^{\pm}+{C}_{6}{\lambda}_{2}^{\pm}+{C}_{6}{\lambda}_{3}^{\pm}\right){e}_{e}\end{array}\right\}$ ; (51)

For the chemical parcel, one has for the six equations, the following numerical flux vector, where *c _{s}* is the mass fraction and is numbered from 1 to 7, except the number 3, which is associated to the N

${F}_{Chem}^{\pm}=\left\{\begin{array}{c}\left({C}_{1}{\lambda}_{1}^{\pm}+{C}_{2}{\lambda}_{2}^{\pm}+{C}_{2}{\lambda}_{3}^{\pm}\right){c}_{s,1}\\ \left({C}_{1}{\lambda}_{1}^{\pm}+{C}_{2}{\lambda}_{2}^{\pm}+{C}_{2}{\lambda}_{3}^{\pm}\right){c}_{s,2}\\ \left({C}_{1}{\lambda}_{1}^{\pm}+{C}_{2}{\lambda}_{2}^{\pm}+{C}_{2}{\lambda}_{3}^{\pm}\right){c}_{s,4}\\ \left({C}_{1}{\lambda}_{1}^{\pm}+{C}_{2}{\lambda}_{2}^{\pm}+{C}_{2}{\lambda}_{3}^{\pm}\right){c}_{s,5}\\ \left({C}_{1}{\lambda}_{1}^{\pm}+{C}_{2}{\lambda}_{2}^{\pm}+{C}_{2}{\lambda}_{3}^{\pm}\right){c}_{s,6}\\ \left({C}_{1}{\lambda}_{1}^{\pm}+{C}_{2}{\lambda}_{2}^{\pm}+{C}_{2}{\lambda}_{3}^{\pm}\right){c}_{s,7}\end{array}\right\}$ ; (52)

Finally, the vibrational component is defined by the following equation:

${F}_{Vib}^{\pm}=\left\{\left({C}_{1}{\lambda}_{1}^{\pm}+{C}_{6}{\lambda}_{2}^{\pm}+{C}_{6}{\lambda}_{3}^{\pm}\right){e}_{v}\right\}$, (53)

where *e _{v}* is the vibrational energy. For the positive components of the flux vectors, the cell properties are employed, whereas for the negative components the neighboring cell properties are used.

The viscous formulation follows that of [25], which adopts the Green theorem to calculate primitive variable gradients. The viscous gradients at the flux interface are obtained by arithmetical average between cell (*i*,*j*) and its neighbors. As was done with the convective terms, there is a need to separate the viscous flux in three parts: dynamic viscous flux, chemical viscous flux, and vibrational viscous flux. The dynamic part corresponds to the first four equations of the Navier-Stokes ones plus the electron equation, the chemical part corresponds to the six equations immediately below the electron energy equation, and the vibrational part corresponds to the one that follows the last chemical equation. The resultant ordinary differential equation system can be written as:

${V}_{i,j}\frac{\text{d}{Q}_{i,j}}{\text{d}t}=-\left({R}_{i,j-1/2}+{R}_{i+1/2,j}+{R}_{i,j+1/2}+{R}_{i-1/2,j}\right)=-{C}_{i,j}$, (54)

where ${R}_{i+1/2,j}={R}_{i+1/2,j}^{Dyn}+{R}_{i+1/2,j}^{Chem}+{R}_{i+1/2,j}^{Vib}$, and the cell volume is given by:

$\begin{array}{c}{V}_{i,j}=0.5\left|\left({x}_{i,j}-{x}_{i+1,j}\right){y}_{i+1,j+1}+\left({x}_{i+1,j}-{x}_{i+1,j+1}\right){y}_{i,j}+\left({x}_{i+1,j+1}-{x}_{i,j}\right){y}_{i+1,j}\right|\\ \text{\hspace{0.17em}}+0.5\left|\left({x}_{i,j}-{x}_{i+1,j+1}\right){y}_{i,j+1}+\left({x}_{i+1,j+1}-{x}_{i,j+1}\right){y}_{i,j}+\left({x}_{i,j+1}-{x}_{i,j}\right){y}_{i+1,j+1}\right|\end{array}$ (55)

For the [7] and [8] schemes, Equation (54) can be used directly. For the [9] scheme, however, the *R* components are calculated by the following equation, considering the dynamic contribution:

${R}_{i,j-1/2}^{Dyn}={F}_{Dyn}^{+}\left({Q}^{-}\right)+{F}_{Dyn}^{-}\left({Q}^{+}\right)$, (56)

where *Q*^{−} and *Q*^{+} are computed considering *Q*^{−} the cell properties and *Q*^{+} the neighboring cell properties.

In the present study, the Euler backward method was employed to march the scheme in time. This method is first-order accurate in time, to the three types of complete flux. To the convective dynamic component, this method can be represented in general form by:

${Q}_{i,j}^{\left(n+1\right)}={Q}_{i,j}^{\left(n\right)}-\Delta {t}_{i,j}\times \left\{\frac{C\left[{Q}_{i,j}^{\left(n\right)}\right]}{{V}_{i,j}}-{S}_{E}\left[{Q}_{i,j}^{\left(n\right)}\right]\right\}$, (57)

with *S _{E}* being the modified electron source term, composed by the last line in Equation (14). To the convective chemical part, it can be represented in general form by:

${Q}_{i,j}^{\left(n+1\right)}={Q}_{i,j}^{\left(n\right)}-\Delta {t}_{i,j}\times \left\{\frac{C\left[{Q}_{i,j}^{\left(n\right)}\right]}{{V}_{i,j}}-{S}_{C}\left[{Q}_{i,j}^{\left(n\right)}\right]\right\}$, (58)

where the chemical source term *S _{C}* is defined by Equation (14) and is calculated with the temperature Trrc (reaction rate control temperature, see [18] [19] ). Finally, to the convective vibrational component, where

${Q}_{i,j}^{\left(n+1\right)}={Q}_{i,j}^{\left(n\right)}-\Delta {t}_{i,j}\times \left\{\frac{C\left[{Q}_{i,j}^{\left(n\right)}\right]}{{V}_{i,j}}-{S}_{V}\left[{Q}_{i,j}^{\left(n\right)}\right]\right\}$. (59)

6. Spatially Variable Time Step

The spatially variable time step has proved efficient gains in terms of convergence acceleration, as verified by [10] [11]. Initially, the parameter *σ* is determined, where:

${\sigma}_{s}=\frac{{c}_{s}}{{M}_{s}}$ and $\sigma ={\displaystyle {\sum}_{s=1}^{ns}{\sigma}_{s}}$, (60)

with *c _{s}* being the mass fraction and

${c}_{v,t}={\displaystyle {\sum}_{s=1}^{ns}{\sigma}_{s}{c}_{v,t,s}}$, (61)

where, for each gas constituent of the seven (7) species chemical model, the specific heat at constant volume, based on the kinetic theory of gases ( [26] ), is defined by

${c}_{v,t,N}=\frac{3}{2}{R}_{N}$, ${c}_{v,t,O}=\frac{3}{2}{R}_{O}$, ${c}_{v,t,{N}_{2}}=\frac{5}{2}{R}_{{N}_{2}}$, ${c}_{v,t,{O}_{2}}=\frac{5}{2}{R}_{{O}_{2}}$ ; (62)

${c}_{v,t,NO}=\frac{5}{2}{R}_{NO}$, ${c}_{v,t,N{O}^{+}}=\frac{5}{2}{R}_{N{O}^{+}}$ and ${c}_{v,t,{e}^{-}}=\frac{3}{2}{R}_{{e}^{-}}$. (63)

being *R _{s}* the specific gas constant. The total pressure of the gaseous mixture is determined by Dalton law, which indicates that the total pressure of the gas is the sum of the partial pressure of each constituent gas, resulting in:

${p}_{s}={c}_{s}\rho {R}_{s}T$ and $p={\displaystyle {\sum}_{s=1}^{ns-1}{p}_{s}}+{p}_{e}$, (64)

where *s *= 7 is related to the electron. The speed of sound to a reactive mixture considering thermochemical non-equilibrium is given by Equation (43), where
$\beta =\frac{{R}_{univ}\sigma}{{c}_{v,t}}$, with *R _{univ}* = 1.987 cal/(g-mol.K). Finally, the spatially variable time step is defined from the CFL definition:

$\Delta {t}_{i,j}=\frac{\text{CFL}\Delta {s}_{i,j}}{\sqrt{{u}_{i,j}^{2}+{v}_{i,j}^{2}}+{a}_{i,j}}$, (65)

where ∆*s _{i}*

7. Results

The initial conditions to the circumference, to the reentry capsule, to the blunt body, and to the double ellipse problems, for a seven species chemical model, are presented in Table 5. The Reynolds number is obtained from data of [27].

Tests were performed in a Core i7 processor of 2.8 GHz and 6.0 Gbytes of RAM microcomputer, in a Windows 10.0 environment. Three (3) orders of reduction of the maximum residual in the field were considered to obtain a converged solution. The residual was defined as the value of the discretized conservation equation.

Table 5. Initial conditions to the four problems studied in this work.

The attack angle was adopted equal to zero. For a matter of simplicity, the following abbreviations were used: [7] scheme = VL, [8] scheme = LS, [9] scheme = SW, Chebyshev-Gauss-Radau = CGR. Only some results obtained with these configurations are presented in this section. Converged results involve second to fifth orders accuracy and eighth and sixteenth orders of accuracy. It is presented the pressure contours, the translational/rotational temperature contours, electron temperature contours and temperature distributions at the body wall (translational/rotational, vibrational and electron temperatures).

7.1. Circumference

The first studied problem is the circumference configuration. Only the SW scheme has presented converged results for the inviscid case. Their solutions are presented in relation to the four types of results described above.

Figures 2-5 present the 2^{nd} order SW solutions in the CGR variant of the

Figure 2. Pressure contours (SW-CGR-2^{nd}).

Figure 3. Temperature contours (SW-CGR-2^{nd}).

Figure 4. Electron temperature contours (SW-CGR-2^{nd}).

Figure 5.Wall temperature distributions (SW-CGR-2^{nd}).

proposed spectral method. The stagnation pressure was estimated in 1585.09 units. The maximum translational/rotational temperature was 11,378.80 K. The maximum electron temperature was 817.79 K. Good symmetrical properties are observed in all contours.

The temperature distributions show that the vibrational temperature remains practically constant along the body wall. The electron temperature has an increasing behavior close to the circumference’s trailing edge. The translational/rotational temperature has a decreasing behavior along the body, increasing close to the trailing edge.

7.2. Reentry Capsule

For this problem, the inviscid solutions are described. The VL and SW schemes have presented converged results and only the CGR variant solutions of the spectral method are exhibited for the VL scheme.

Figures 6-9 present the 3^{rd} order VL solutions for the four types of results described above. The maximum pressure, stagnation pressure, was calculated in 1494.54 unities, at the body’s nose, and the maximum translational/rotational temperature, the stagnation temperature, was observed in 9364.49 K. The temperature distributions present the quasi-constant vibrational behavior and the increase of the electron temperature along the body.

Figure 6. Pressure contours (VL-CGR-3^{rd}).

Figure 7. Temperature contours (VL-CGR-3^{rd}).

Figure 8. Electron temperature contours (VL-CGR-3^{rd}).

Figure 9. Wall temperature distributions (VL-CGR-3^{rd}).

7.3. Blunt Body Problem

Solutions of second, third, fourth, fifth, eighth and sixteenth orders of accuracy are obtained with the Spectral method for the inviscid and viscous cases.

*Inviscid**Case.* Figures 10-13 exhibit the pressure and temperature contours as well the temperature distribution profiles along the body’s wall obtained with the SW scheme as using the Spectral method of fourth order of accuracy, with CGR weighting function. The maximum pressure is obtained at the body’s nose with the value of 141.94 unities, whereas the maximum temperature is obtained

Figure 10. Pressure contours (SW-CGR-4^{th}).

Figure 11. Temperature contours (SW-CGR-4^{th}).

Figure 12. Electron temperature contours (SW-CGR-4^{th}).

with the value of 9366.70 K. The maximum electron temperature is obtained with the value of 783.00 K. The wall temperature profiles are shown at Figure 13. In this figure the vibrational temperature is not constant and has a decreasing behavior. The electron profile has a typical constant profile along the wall. The translational/rotational temperature profile reaches a minimum at the body’s end, close to 5000 K. The solutions are symmetrical in relation to the body’s symmetry line.

*Viscous Case.* Figures 14-17 present the pressure and the temperature contours

Figure 13. Wall temperature distributions (SW-CGR-4^{th}).

Figure 14. Pressure contours (LS-CGR-8^{th}).

Figure 15. Temperature contours (LS-CGR-8^{th}).

Figure 16. Electron temperature contours (LS-CGR-8^{th}).

as well the temperature distributions along the body’s wall obtained with the LS scheme using CGR Spectral method of 8^{th} order of accuracy. The stagnation pressure is about 179.97 unites, whereas the stagnation translational/rotational temperature is estimated as 10,321.00 K. In Figure 16 the maximum electron temperature is obtained by the LS scheme with value 824.72 K. Some oscillations are perceptible in the electron temperature contours.

7.4. Double Ellipse

For the double ellipse problem, converged results were obtained with the CGR weighting function. It yields converged results for second, third, fourth, fifth, eighth and sixteenth orders of accuracy.

Figures 18-21 exhibit the pressure and temperature contours as well the

Figure 17. Wall temperature distributions (LS-CGR-8^{th}).

Figure 18. Pressure contours (LS-CGR-16^{th}).

temperature distribution profiles along the body’s wall obtained with the LS scheme as using the Spectral method of 16^{th} order of accuracy, with CGR weighting function.

The maximum pressure is obtained at the body’s nose with the value of 1822.17 unities, whereas the maximum temperature is obtained with the value of 11,025.30 K. The maximum electron temperature is obtained with the value of 302.40 K. The wall temperature profiles are shown at Figure 21. In this figure the vibrational temperature is quasi-constant and is unsymmetrical due to the asymmetry of the configuration. The electron profile has also a quasi-constant profile along the wall. The translational/rotational temperature profile reaches a minimum at

Figure 19. Temperature contours (LS-CGR-16^{th}).

Figure 20. Electron temperature contours (LS-CGR-16^{th}).

the body’s end, close to 5000 K.

7.5. Error Analysis

In order to perform an error analysis, the present reactive results are compared to the perfect gas solutions. The stagnation pressure at the blunt body nose, circumference nose and reentry capsule nose were evaluated assuming the perfect gas formulation. Such parameter calculated at this way is not the best comparison, but in the absence of practical reactive results, this constitutes the best available solution. To calculate the stagnation pressure ahead of the blunt body, for example, [28] presents in its B Appendix values of the normal shock wave properties ahead of the configuration. The ratio *pr*_{0}/*pr*_{∞} is estimated as function of the normal Mach number and the stagnation pressure *pr*_{0} can be determined

Figure 21. Wall temperature distributions (LS-CGR-16^{th}).

from this parameter. The value of *pr*_{∞} is determined by the following expression:

$p{r}_{\infty}=\frac{p{r}_{initial}}{{\rho}_{\infty}{a}_{\infty}^{2}}$, (66)

where, for example, for the blunt body problem, *pr _{initial}* = 687 N/m

Tables 7-11 compare the values of the stagnation pressure obtained from the simulations with these theoretical values and show the percentage errors. Table 7 has the circumference data, Table 8 has the reentry capsule data, and Table 9 has the inviscid blunt body data. Finally, Table 10 has the viscous blunt body data and Table 11 has the double ellipse data. These results are composed by all converged results obtained in this study, with the majority not exhibited in this article due to paper size.

As can be observed, the best solution to the circumference problem was obtained by the SW scheme using the CGR spectral method of second order, with the estimation error of the stagnation pressure of only 5.69%. To the reentry capsule problem, the best solution was again obtained with the SW scheme, addressing the minimum relative error of 0.11% in the estimative of the stagnation pressure. It was also obtained with the CGR spectral method of [6], in its fifth order of accuracy. To the third problem, the blunt body configuration, two error analysis were done: one for the inviscid case and the other for the viscous case.

Table 6. Values of theoretical stagnation pressure.

Table 7. Values of stagnation pressure and respective errors (Circumference Problem).

Table 8. Values of stagnation pressure and respective errors (Reentry Capsule Problem).

For the inviscid case, the best value of the stagnation pressure was obtained by the VL scheme using the CGR spectral method in its second order of accuracy, with an error of 0.33%. On the other hand, for the viscous case, the best estimative of the stagnation pressure was obtained by the LS scheme using CGR and sixteenth order of accuracy, presenting an error of 2.69%. Finally, for the fourth problem, the double ellipse problem, the best solution in terms of stagnation pressure estimation was obtained by the LS scheme using the CGR spectral method

Table 9. Values of stagnation pressure and respective errors (Blunt Body Inviscid Problem).

Table 10. Values of stagnation pressure and respective errors (Blunt Body Viscous Problem).

Table 11. Values of stagnation pressure and respective errors (Double Ellipse Problem).

in its second order accuracy, with an error of 1.77%.

7.6. Grid Convergence Index (GCI)

The discretization error ( [29] ) will be estimated using a method popular in the Computational Fluid Dynamics (CFD) community called the Grid Convergence Index (GCI), due to [30]. Two analyses are possible: constant grid refinement ratio and non-uniform grid refinement ratio. In the present study, the constant grid refinement ratio was considered and will be presented in this section. The property under analysis is the maximum pressure, stagnation pressure, of the blunt body problem, obtained by the [9] flux vector splitting scheme using the CGR spectral method of 2^{nd} order.

*Constant Grid Refinement Ratio.* For the special case where the meshes are constructed with a constant grid refinement ratio, *r*, e.g. *r* = *h*_{2}/*h*_{1} = *h*_{3}/*h*_{2} = constant where *h*_{1} < *h*_{2} < *h*_{3}, for a finite volume formulation, the convergence rate, *p*, can be estimated as

$p=\frac{\left|{\mathrm{log}}_{10}\left(\frac{{f}_{3}-{f}_{2}}{{f}_{2}-{f}_{1}}\right)\right|}{{\mathrm{log}}_{10}\left(r\right)}$. (67)

where *f*_{1}, *f*_{2} and *f*_{3} are numerical values of the property under investigation for each mesh. Once the observed order-of-convergence is known, an estimate of the error between the fine grid solution and the unknown exact solution can be done. The relative error between the two finest grids is given by

${e}_{21}=\left|\frac{{f}_{2}-{f}_{1}}{{f}_{1}}\right|$. (68)

The Grid Convergence Index (GCI) provides an estimate of the amount of discretization error in the finest grid solution relative to the converged numerical solution. Note the discretization error is with respect to the convergent numerical solution, as the exact solution is generally unknown. The GCI is given by

${\text{GCI}}_{21}={F}_{s}\frac{{e}_{21}}{{r}_{21}^{p}-1}$. (69)

where *F _{s}* is a “safety factor” multiplying the relative error term. The safety factor,

*The converged numerical solution lies in the interval *[4.90, 5.10]*, i.e. *[*f*_{1}(1 − GCI), *f*_{1} (1 + GCI)]*, with a *95%*confidence level. *

Finally, based upon the two finest mesh solutions, and the estimate of the observed convergence rate, an extrapolation of the numerical solution is possible:

${f}_{21}^{*}=\frac{{r}_{21}^{p}{f}_{1}-{f}_{2}}{{r}_{21}^{p}-1}$. (70)

The extrapolation solution, ${f}_{21}^{*}$, provides a useful estimate of the converged solution, [29].

For the present work, it was chosen three meshes to observe the behavior of the maximum pressure in the simulations. The meshes are composed of 14,518 cells (123 × 120), 3658 cells (63 × 60), and 986 cells (35 × 30). The numerical values of the maximum pressure observed in the simulations were: *f*_{1} = 188.89, *f*_{2} = 160.88, and *f*_{3} = 118.94. The ratio *r*_{21} = *h*_{2}/*h*_{1} = 2.0 and the ratio *r*_{32} = *h*_{3}/*h*_{2} = 2.0, where *r* is constant. Hence, the value of the convergence rate is calculated by

$p=\frac{\left|{\mathrm{log}}_{10}\left(\frac{{f}_{3}-{f}_{2}}{{f}_{2}-{f}_{1}}\right)\right|}{{\mathrm{log}}_{10}\left(r\right)}=\frac{\left|{\mathrm{log}}_{10}\left(\frac{118.94-160.88}{160.88-188.89}\right)\right|}{{\mathrm{log}}_{10}\left(2.0\right)}=\frac{0.17532}{0.30103}=0.58240$. (71)

The relative error *e*_{21} is calculated as

${e}_{21}=\left|\frac{{f}_{2}-{f}_{1}}{{f}_{1}}\right|=\left|\frac{160.88-188.89}{188.89}\right|=0.14829$. (72)

and GCI_{21} can be calculated as

${\text{GCI}}_{21}={F}_{s}\frac{{e}_{21}}{{r}_{21}^{p}-1}=1.25\times \frac{0.14829}{{2.0}^{0.58240}-1}$ $\therefore {\text{GCI}}_{21}=0.37231$. (73)

The interval of solution for this parameter is

$S=\left[{f}_{1}\left(1-{\text{GCI}}_{21}\right),{f}_{1}\left(1+{\text{GCI}}_{21}\right)\right]=\left[118.56436,259.21564\right]$, (74)

with a 95% confidence level. Finally, the extrapolation solution is calculated by:

${f}_{21}^{*}=\frac{{r}_{21}^{p}{f}_{1}-{f}_{2}}{{r}_{21}^{p}-1}=\frac{{2.0}^{0.58240}\times 188.89-160.88}{{2.0}^{0.58240}-1}=245.21$. (75)

This value corresponds to the expected numerical solution to be obtained. This also suggests that the more realistic value to the stagnation pressure, taking into account the chemical effects, should be higher than 200.00 unities in place of 170.87 unities. The solution presented in this manuscript assumed the value *f*_{2} = 160.88 and the error in relation to the extrapolation solution is 34.39%. However, it pertains to the 95% confidence level, which suggests that a more refined mesh is optional.

7.7. Three-Temperature Effect

The three-temperature effect can be seen in Table 12 and Table 13, where the values of the minimum errors committed by each scheme are compared with the present paper and the reference [31], which was an old work of the present author. Comparing the two tables, it is possible to highlight the better performance of the three-temperature model in relation to the two-temperature model of [31]. Of course, this comparison is not the best because such study should be accomplished with 4^{th} order of accuracy by both works, for example. But as the 4^{th} order study of [31] incorporated the best results for the inviscid case for each problem, so, the present results compare better than the old work. The unique loss of accuracy in relation to the old work (0.23%) was for the blunt body problem in its inviscid case (0.33%). All other comparisons exhibit the best performance of the three-temperature model, with maximum error of the order of

Table 12. Values of stagnation pressure and respective errors (Inviscid case).

Table 13. Values of stagnation pressure and respective errors (Viscous case).

2.69%. Finally, it is important to highlight better performance of the three-temperature model in the determination of the stagnation pressure for the studied configurations, than the two-temperature model.

In summary, the paper has presented results for the four studied cases and has compared them with literature theoretical results of [28]. The advantage of the present implementation is evident by the minimum errors found in each case in the present study and for the comparison of results with another reference [31], also of this author, that presented worse values for the stagnation pressure estimation in relation to this paper. The three-temperature model performance was highlighted in this section, emphasizing this one as more accurate them other two-temperature models. Particularly, the results of [31] were obtained with a computational code written in FORTRAN language and the present results with a code written in OBJECT PASCAL language. So, the performances are different and the realistic comparison between codes are evident, with better behavior for the later implementation with the proposed three-temperature model.

The advantages of the present formulation with a three-temperature model is highlighted by the good numerical results obtained in terms of percentage errors, with values close to 0.11% in the best case, and in comparison with other references as [31]. About the effectiveness of the proposed model, this can be observed that the implemented spectral method, used to high order analysis, has presented converged results for the majority of the studied cases, evidencing the appropriated formulation presented here. Moreover, the three algorithms have presented converged results, emphasizing the correct implementation and, as they are upwind and flux vector splitting schemes, don’t requiring artificial dissipation to guarantee convergence, they are also cheaper than the flux difference splitting counterparts. The fact of not requiring artificial dissipation is a critical advantage of these schemes, eliminating the phase of user trial and error approaches to convergence and providing a solution with the appropriated level of dissipation.

8. Conclusions

In this work, a study involving the fully coupled Euler and Navier-Stokes reactive equations is performed. These equations, in conservative and finite volume contexts, employing structured spatial discretization, on a condition of thermochemical non-equilibrium, are analyzed. High-order studies are accomplished using the spectral method of Street *et al*. The high enthalpy hypersonic flows around a circumference, around a reentry capsule, along a blunt body, and along a double ellipse in two-dimensions are simulated. The Van Leer, Liou and Steffen Jr., and the Steger and Warming flux vector splitting algorithms are applied to execute the numerical experiments. The Euler backward integration method is employed to march the schemes in time. The convergence process is accelerated to steady state condition through a spatially variable time step procedure, which has proved effective gains in terms of computational acceleration (see Maciel). The reactive simulations involve Earth atmosphere chemical model of seven species and eighteen reactions, based on the Blottner model. Three temperatures are used to accomplish the numerical comparisons.

The results have demonstrated that the Liou and Steffen Jr. and Steger and Warming schemes were the best in the estimative of the stagnation pressure ahead of the hypersonic geometries studied in this work. For the circumference problem, the best solution was obtained by the Steger and Warming scheme using the CGR spectral method of second order, with the estimation error of the stagnation pressure of only 5.69%. To the reentry capsule problem, the best solution was again obtained with the Steger and Warming scheme, addressing the minimum relative error of 0.11% in the estimative of the stagnation pressure. It was also obtained with the CGR spectral method of Street *et al*., in its fifth order of accuracy. To the third problem, the blunt body configuration, two error analysis were done: one for the inviscid case and the other for the viscous case. For the inviscid case, the best value of the stagnation pressure was obtained by the Van Leer scheme using the CGR spectral method in its second order of accuracy, with an error of 0.33%. On the other hand, for the viscous case, the best estimative of the stagnation pressure was obtained by the Liou and Steffen Jr. scheme using CGR weighting function and sixteenth order of accuracy, presenting an error of 2.69%. Finally, for the fourth problem, the double ellipse problem, the best solution in terms of stagnation pressure estimation was obtained by the Liou and Steffen Jr. scheme using the CGR spectral method in its second order accuracy, with an error of 1.77%.

The main contribution of this work to the CFD community was the implementation of three-temperature model coupled with three flux vector splitting algorithms to the solution of reentry flow problems in 2D and the use of a spectral method to error analysis. With it, numerical tools are available to study thermochemical non-equilibrium flows in two-dimensions more realistically. The detailed source terms implemented to consider electron temperature and more realistic simulations were remarkable contributions from this work. Finally, this in house code was implemented in Object PASCAL language, with the use of Delphi software produced by Borland.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

[1] |
Howe, J.T. (1985) Introductory Aerothermodynamics of Advanced Space Transportation Systems. Journal of Spacecraft and Rockets, 22, 19-26. https://doi.org/10.2514/3.25705 |

[2] |
Madjarevic, D., Ruggeri, T. and Simic, S. (2014) Shock Structure and Temperature Overshoot in Macroscopic Multi-Temperature Model of Mixtures. Physics of Fluids, 26, Article ID: 106102. https://doi.org/10.1063/1.4900517 |

[3] | Cai, C. and Huang, X. (2011) A Kinetic BGK Scheme for Non-Equilibrium Gaseous Flow Computation. 20th AIAA Computational Fluid Dynamics Conference, Honolulu, 27-30 June 2011, AIAA Paper 2011-3543. |

[4] |
Zhong, X., MacCormack, R.W. and Chapman, D.R. (1993) Stabilization of the Burnett Equations and Application to Hypersonic Flows. AIAA Journal, 31, 1036-1043. https://doi.org/10.2514/3.11726 |

[5] |
Gottlieb, D. and Orszag, S. (1987) Numerical Analysis of Spectral Methods: Theory and Applications. Society for Industrial and Applied Mathematics, Philadelphia. https://doi.org/10.1137/1.9781611970425 |

[6] |
Streett, C.L., Zang, T.A. and Hussaini, M.Y. (1984) Spectral Methods for Solution of the Boundary-Layer Equations. AIAA Paper 84-0170. https://doi.org/10.2514/6.1984-170 |

[7] |
Van Leer, B. (1982) Flux-Vector Splitting for the Euler Equations. 8th International Conference on Numerical Methods in Fluid Dynamics, Vol. 170, Aachen, 28 June-2 July 1982, 507-512. https://doi.org/10.1007/3-540-11948-5_66 |

[8] |
Liou, M. and Steffen Jr., C.J. (1993) A New Flux Splitting Scheme. Journal of Computational Physics, 107, 23-39. https://doi.org/10.1006/jcph.1993.1122 |

[9] |
Steger, J.L. and Warming, R.F. (1981) Flux Vector Splitting of the Inviscid Gasdynamic Equations with Application to Finite-Difference Methods. Journal of Computational Physics, 40, 263-293. https://doi.org/10.1016/0021-9991(81)90210-2 |

[10] | Maciel, E.S.G. (2015) Simulations in 2D and 3D Applying Unstructured Algorithms, Euler and Navier-Stokes Equations—Perfect Gas Formulation. Lambert Academic Publishing (LAP), Saarbrücken, Ch. 1, 26-47. |

[11] | Maciel, E.S.G. (2015) Simulations in 2D and 3D Applying Unstructured Algorithms, Euler and Navier-Stokes Equations—Perfect Gas Formulation. Lambert Academic Publishing (LAP), Saarbrücken, Ch. 6, 160-181. |

[12] |
Blottner, F.G. (1969) Viscous Shock Layer at the Stagnation Point with Nonequilibrium Air Chemistry. AIAA Journal, 7, 2281-2288. https://doi.org/10.2514/3.5528 |

[13] | Maciel, E.S.G. (2021) A Spectral Method Applied to Re-entry Flow Problems. Volume 2—Chemical and Thermochemical Non-Equilibrium Applications to Reentry Flows. Lambert Academic Publishing (LAP), Republic of Moldova. |

[14] | Davis, P.A. and Rabinowitz, P. (1967) Numerical Integration. Blaisdell Publishing Co., Waltham. |

[15] |
Canuto, C., Hussaini, M.Y., Quarteroni, A. and Zang, T.A. (2007) Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics. Scientific Computation Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-540-30728-0 |

[16] | Prabhu, R.K. (1994) An Implementation of a Chemical and Thermal Nonequilibrium Flow Solver on Unstructured Meshes and Application to Blunt Bodies. National Aeronautics and Space Administration, Washington DC, NASA CR-194967. |

[17] |
Saxena, S.K. and Nair, M.T. (2005) An Improved Roe Scheme for Real Gas Flow. 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, 10-13 January 2005, AIAA Paper 2005-0587. https://doi.org/10.2514/6.2005-587 |

[18] | Maciel, E.S. G. (2015) Hypersonic Reactive Flow Simulations in Two-Dimensions, Chemical and Thermochemical Non-Equilibrium Conditions. Lambert Academic Publishing (LAP), Saarbrücken, Ch. 5, 333-406. |

[19] | Maciel, E.S.G. (2015) Hypersonic Reactive Flow Simulations in Two-Dimensions, Chemical and Thermochemical Non-Equilibrium Conditions. Lambert Academic Publishing (LAP), Saarbrücken, Ch. 6, 407-490. |

[20] |
Candler, G.V. and MacCormack, R.W. (1991) Computation of Weakly Ionized Hypersonic Flows in Thermochemical Nonequilibrium. Journal of Thermophysics and Heat Transfer, 5, 266-273. https://doi.org/10.2514/3.260 |

[21] |
Lee, J.H. (1985) Basic Governing Equations for the Flight Regimes of Aeroassisted Orbital Transfer Vehicles. Thermal Design of Aeroassisted Transfer Vehicles. Progress in Astronautics and Aeronautics, 96, 3-53. https://doi.org/10.2514/5.9781600865718.0003.0053 |

[22] |
Kim, M., Gülhan, A. and Boyd, I.D. (2011) Modeling of Electron Temperature in Hypersonic Flows. 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Orlando, 4-7 January 2011, AIAA Paper 2011-1028. https://doi.org/10.2514/6.2011-1028 |

[23] |
Ait-Ali-Yahia, D. and Habashi, W.G. (1997) Finite Element Adaptive Method for Hypersonic Thermochemical Nonequilibrium Flows. AIAA Journal, 35, 1294-1302. https://doi.org/10.2514/2.260 |

[24] |
Radespiel, R. and Kroll, N. (1995) Accurate Flux Vector Splitting for Shocks and Shear Layers. Journal of Computational Physics, 121, 66-78. https://doi.org/10.1006/jcph.1995.1179 |

[25] |
Long, L.N., Khan, M.M.S and Sharp, H.T. (1991) Massively Parallel Three-Dimensional Euler/Navier-Stokes Method. AIAA Journal, 29, 657-666. https://doi.org/10.2514/3.10635 |

[26] | Vincenti, W.G. and Kruger Jr., C.H. (1965) Introduction to Physical Gas Dynamics. John Wiley & Sons, New York, London, Sidney. |

[27] | Fox, R.W. and McDonald, A.T. (1988) Introducao à Mecanica dos Fluidos. Guanabara Editor, Rio de Janeiro. |

[28] | Anderson Jr., J.D. (2010) Fundamentals of Aerodynamics. 5th Edition, McGraw-Hill, Inc., New York. |

[29] | Schwer, L.E. (2008) Is Your Mesh Refined Enough? Estimating Discretization Error using GCI. LS-DYNA Anwenderforum, Bamberg. |

[30] |
Roache, P.J. (1994) Perspective: A Method for Uniform Reporting of Grid Refinement Studies. ASME Journal of Fluids Engineering, 116, 405-413. https://doi.org/10.1115/1.2910291 |

[31] | Maciel, E.S.G. (2017) Spectral Method Applied to Thermochemical Non-Equilibrium Reentry Flows in 2D: Seven Species—High Order Analysis. Journal of Scientific and Engineering Research, 4, 302-339. |

Journals Menu

Contact us

+1 323-425-8868 | |

customer@scirp.org | |

+86 18163351462(WhatsApp) | |

1655362766 | |

Paper Publishing WeChat |

Copyright © 2024 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.