SU2-NEMO: An Open-Source Framework for High-Mach Nonequilibrium Multi-Species Flows

Abstract

SU2-NEMO, a recent extension of the open-source SU2 multiphysics suite’s set of physical models and code architecture, is presented with the aim of introducing its enhanced capabilities in addressing high-enthalpy and high-Mach number flows. This paper discusses the thermal nonequilibrium and finite-rate chemistry models adopted, including a link to the Mutation++ physio-chemical library. Further, the paper discusses how the software architecture has been designed to ensure modularity, incorporating the ability to introduce additional models in an efficient manner. A review of the numerical formulation and the discretization schemes utilized for the convective fluxes is also presented. Several test cases in two- and three-dimensions are examined for validation purposes and to illustrate the performance of the solver in addressing complex nonequilibrium flows.

1. Introduction

Despite continued research efforts, numerical simulation of high-Mach flows remains a significant challenge in computational fluid dynamics (CFD). From a design standpoint, accurate prediction of the nonequilibrium aerothermodynamic environment is necessary for the development and optimization of vehicle systems to sustain the large aerodynamic and thermal loads experienced during hypersonic flight. This prediction is made challenging by the presence of finite-rate processes that arise due to highly-energetic molecular collisions, affecting both fluid bulk and transport properties. The development of efficient and robust computational models that effectively capture the coupling between fluid bulk motion, chemical kinetics, and thermodynamic processes remains an active area of research within the CFD community.

The hypersonic regime is characterized by high-temperature effects associated with high-flow kinetic energy relative to the static thermal enthalpy. This can give rise to finite-rate kinetic processes when the flow is rapidly decelerated, such as across a strong shock. When the relaxation time of these processes is on the order of the fluid mechanical residence time, a state of thermochemical nonequilibrium is established. In this state, the thermodynamic properties and species concentrations must be spatially resolved to accurately model the flow. High fluid-specific energy densities make recreating the mission-relevant time and length scales difficult and costly, particularly in ground-based experimental facilities [1,2]. The inability to completely capture these physical phenomena can lead to catastrophic failure, such as the case of the NASA X-15 [3], where an unanticipated shock impingement and the associated localized heating caused structural failure. Numerical simulation provides an additional, and often necessary, method to predict the flow environment. It is commonly employed for test planning and data reduction in ground-based test facilities [4], as well as in analyses during the vehicle design process. As flow enthalpy increases, nonequilibrium effects become increasingly significant and complex, requiring additional models to capture the underlying physics. The addition of these models to CFD codes gives rise to new challenges in solver stability, robustness, and efficiency that have been documented in the literature.

Design of vehicles for atmospheric entry has provided a major historical impetus for the development of nonequilibrium modeling capabilities. An early effort to incorporate chemical nonequilibrium effects into numerical methods using upwinding schemes was carried out by Gnoffo and McCandless [5], who implemented loose coupling of chemical source terms that provided accurate predictions of the chemical species concentrations for high-speed flows around spheres. Increased computational power and improved numerical methods have allowed for higher-fidelity solutions, including the implementation of individual conservation equations for each species present in the flow and each molecular energy mode. The LAURA solver [6] developed at the NASA Langley Research Center represents a well-validated solver incorporating these reacting Navier–Stokes equations.

The DPLR solver [7], developed at the NASA Ames Research Center, takes advantage of the uniform geometry of blunt body entry vehicles for simulation on structured grids. Unlike LAURA, DPLR employs a modified Gauss–Seidel line relaxation method to efficiently solve along lines perpendicular to the surface, which shows rapid convergence even in the presence of large gradients [7]. The primary challenge in these cases is efficiently generating structured grids, particularly for complex geometries.

US3D [8,9] is a solver developed as an extension of the DPLR algorithm for unstructured grids. US3D uses a hybrid grid approach with line relaxation in the structured region of the grid, and point relaxation everywhere else. US3D represents the state-of-the-art in hypersonic simulation capabilities, including a full suite of physical models for nonequilibrium effects, turbulence, wall-modeled large eddy simulation, grid alignment, and material thermal response and shape change [10].

Several open-source codes for modeling flows in thermochemical nonequilibrium have also been developed. Within the widely used OpenFOAM CFD framework, a coupled DSMC-CFD solver, hy2Foam [11], has been deployed. The hy2Foam solver utilizes a two-temperature model and has been validated for a range of operating conditions employing a five-species air model. Another open-source code suite for simulating hypersonic flows is COOLFluiD [12], developed at the Von Karman Institute for Fluid Dynamics and at the KU Leuven Center for mathematical Plasma Astrophysics. COOLFluiD leverages state-of-the-art numerical methods to solve multi-physics problems on unstructured methods, and has been applied to problems in re-entry aerothermodynamics, aeroacoustics, turbulence modeling, and plasma dynamics.

Recent interest in sustainable spaceflight has motivated the design of reusable space systems capable of accommodating the entry environment many times over. These developments highlight the need for design and optimization software capable of modeling the complex physics of nonequilibrium flows across a wide range of mission profiles and vehicle geometries. One such design-optimization-focused code is SU2. SU2 [13,14] is an open-source software suite developed to solve partial differential equations (PDEs) and PDE-constrained optimization problems. While SU2 was developed primarily for computational fluid dynamics, it has also been applied to problems in heat conduction and radiation. In order to capture the complex physics of nonequilibrium flows [15], an effort has been undertaken to enhance the nonequilibrium modeling capabilities in SU2, culminating in the NonEquilibrium MOdeling (NEMO) code base within the SU2 framework. This effort has extended the SU2 base code to include additional conservation equations to capture the physics of multi-component reacting flows in thermodynamic and chemical nonequilibrium.

Past work in SU2 has included implementation of a continuous adjoint formulation for design sensitivity analysis in hypersonic flows [15] and study of shock interaction mechanisms in other-than-air gases [16,17]. Other developments include integration of anisotropic mesh adaptation using the discrete adjoint for nonequilibrium flows [18]. Because of its modular design, SU2 is an excellent framework for the test and development of simulation and physical modeling capabilities for hypersonic, reacting flows. Recent efforts have focused on coupling SU2-NEMO with the Mutation++ thermochemical library [19,20] to provide an alternative to the set of chemistry models included in the code. Moreover, SU2 is built with design optimization in mind, utilizing the discrete adjoint methodology to efficiently extract design sensitivities. This will allow SU2-NEMO to not only study the complex flow physics, but to also be used to aid the design of hypersonic systems. This paper discusses the implementation of SU2-NEMO, including the code structure, constitutive physical models, and the numerical implementation. Verification and validation (V&V) test cases are presented, demonstrating the code’s capability of modeling a variety of flow physics in comparison to numerical and experimental results. Numerical results are presented for complex and realistic vehicle geometry, demonstrating a direction for future application and continued development of the code base.

2. Code Architecture and Design

SU2-NEMO is built upon the existing open-source, multiphysics code suite, SU2. SU2 was designed for the analysis of PDE-based problems using state-of-the-art numerical methods [13]. The primary code base is written in C++ and includes Python scripts for analysis, design, and optimization tasks. SU2 utilizes the object-oriented nature of the C++ programming language to efficiently implement new capabilities within its class structure. The overarching themes within the software are polymorphism and modularity; functions can be adapted or added without impacting the rest of the code. SU2-NEMO relies on the higher-level functions provided by the SU2 framework (geometry, integration, and output class structures), but has added numerical methods (convective schemes, viscous fluxes, and source terms) unique to the modeling of the nonequilibrium hypersonic environment.

The modularity of the SU2-NEMO architecture facilitates the rapid implementation of new models, providing an efficient framework for the test and development of new schemes and thermochemical libraries without having to significantly alter existing codes [15]. Additional advantages of creating SU2-NEMO within the class structure of the SU2 framework include the ability to utilize and modify a variety of different implicit solvers and the ability to leverage the mature discrete adjoint infrastructure already available in SU2 for nonequilibrium flows. An overview of the general SU2 class architecture is given in Figure 1. The authors invite the readers to see the original SU2 paper [13] for a more in-depth discussion of the code architecture and philosophy.

The main extensions of the code for SU2-NEMO occur within the CSolver, CFluidModel, and CNumerics classes. Like the SU2 framework, SU2-NEMO includes solver classes that define the physical phenomena being modeled: Euler Equations, Navier–Stokes, RANS, and so forth. To account for the additional conservation equations and differing data structures required to model the flow physics, the solver classes are extended with NEMO-specific functions. Similarly, this logic is extended to the numerics classes to include NEMO-specific convective schemes, such as MSW and AUSM. These extensions to the solver and numerics classes can be seen in Figure 2a,b. The NEMO-specific classes have been generalized for any nonequilibrium flow. The additions to the framework allow for an arbitrary gas model to be simulated without further changes to the CNumerics, CVariable, or CSolver structures.

Gas-specific properties are computed using the CFluidModel class. A major benefit of the SU2-NEMO framework is the ability to alter or add different thermochemical models using the CNEMOGas child class within the CFluidModel, as seen in Figure 3. The CNEMOGas class contains all data and functions necessary to compute flow thermochemical and transport properties. Currently, CNEMOGas has two additional child classes: CSU2TCLib, the native thermochemical library, and CMutationTCLib, a class implementing a connection with Mutation++.

2.1. Mutation++

Due to the modular implementation of SU2, the code suite is amenable to interfacing with external libraries. SU2-NEMO takes advantage of this through an abstraction of thermochemical source terms using the CMutationTCLib class. SU2-NEMO is linked to the MUlticomponent Thermodynamic And Transport properties for IONized gases library in C++ (Mutation++) [19], a well-validated physio-chemical library. CMutationTCLib implements calls to Mutation++ for computation of mixture thermodynamic, chemical kinetic, and transport properties. Mutation++ is developed, updated, and configuration-managed by the von Karman Institute, utilizing high-fidelity models for nonequilibrium processes in reacting mixtures. Mutation++ provides the capability to simulate a wide array of different gas compositions.

2.2. SU2 Native Thermochemical Library

SU2-NEMO also contains a native thermochemical library, CSU2TCLib. The native SU2 library includes three standard gas models: Air-5, N₂, and Argon. CSU2TCLib demonstrates similar robustness and convergence behavior to CMutationTCLib, and permits easier integration with other SU2 capabilities, but has been shown to be less computationally efficient than CMutationTCLib on benchmark test cases.

3. Physical Modeling

Continuum hypersonic flows are modeled in SU2-NEMO with the Navier–Stokes equations extended for reacting flows in thermochemical nonequilibrium. These are a set of coupled nonlinear partial differential equations, given below in their conservation form in Equation (1).

$$\mathbf{R}(\mathbf{U},\nabla \mathbf{U})=\frac{\partial \mathbf{U}}{\partial t}+\nabla ·{\mathbf{F}}^c\left(\mathbf{U}\right)-\nabla ·{\mathbf{F}}^v(\mathbf{U},\nabla \mathbf{U})-\mathbf{Q}\left(\mathbf{U}\right)=\mathbf{0},$$

(1)

where the conservative variables, convective flux, viscous flux, and source terms are given using standard notation for a number of species $n_s$ by

$$\mathbf{U}=\left\{\begin{array}{c}{\rho }_1\\ ⋮\\ {\rho }_{n_s}\\ \rho \mathbf{u}\\ \rho e\\ \rho e^{ve}\end{array}\right\}, {\mathbf{F}}^c=\left\{\begin{array}{c}{\rho }_1\mathbf{u}\\ ⋮\\ {\rho }_{n_s}\mathbf{u}\\ \rho \mathbf{u}{\mathbf{u}}^T+P\mathcal{I}\\ \rho h\mathbf{u}\\ \rho e^{ve}\mathbf{u}\end{array}\right\}, {\mathbf{F}}^v=\left\{\begin{array}{c}-{\mathbf{J}}_1\\ ⋮\\ -{\mathbf{J}}_{n_s}\\ \sigma \\ {\mathbf{u}}^T\sigma -\sum _k{\mathbf{q}}^k-\sum _s{\mathbf{J}}_s{h}_s\\ -{\mathbf{q}}^{ve}-\sum _s{\mathbf{J}}_s{e}^{ve}\end{array}\right\}, \mathbf{Q}=\left\{\begin{array}{c}{\dot{\omega }}_1\\ ⋮\\ {\dot{\omega }}_{n_s}\\ \mathbf{0}\\ 0\\ {\dot{\mathsf{\Theta }}}^{tr:ve}+\sum _s{\dot{\omega }}_s{e}_s^{ve}\end{array}\right\},$$

(2)

where the details of the notation will be described in subsequent sections. The thermochemical nonequilibrium closure terms (Q), transport properties, and mixture energies introduced in the next section are provided by the CFluidModel class, either through the CMutationTCLib class linked to the Mutation++ library or CSU2TCLib class containing the native SU2 library.

3.1. Two-Temperature Model

Energy carried within polyatomic molecules in a gas-phase species is partitioned among translational, rotational, vibrational, and electronic degrees of freedom. The energy stored in each mode is quantized, and the allowable energy levels are given by eigenstates of the time-independent Schrödinger equation [21].

A complete description of the interaction among energy modes is too complex for implementation in a CFD code, so simplifying assumptions are used. In the models implemented in SU2-NEMO, rotational energy is determined assuming rigid body dynamics (fixed bond length) for polyatomic molecules, and vibrational energy is approximated as a simple harmonic oscillator. Inter-mode coupling between the energy modes, such as coupling between the rotational and vibrational modes due to bond length variation, is neglected and the energy modes are assumed to be independent.

To accommodate differences in the number of collisions required to reach equilibrium, separate temperatures are used to track the translational–rotational energy modes and vibrational–electronic energy modes. This two-temperature model assumes equilibrium between the translational and rotational energy modes, and between the vibrational and electronic modes, but also assumes that these two sets are not necessarily in equilibrium with each other. The rigid-rotor harmonic oscillator (RRHO) two-temperature model has been shown to accurately predict flow and transport properties from high supersonic to atmospheric entry conditions [15]. More sophisticated correction terms can be used where modes are excited sufficiently far from their ground state, such that non-ideal and anharmonic behavior becomes significant to solution accuracy.

Through the independence of the energy levels, the total energy and vibrational–electronic energy per unit volume can be expressed as

$$\rho e=\sum _s{\rho }_s\left(e_s^{tr}+e_s^{rot}+e_s^{vib}+e_s^{el}+e_s^{\circ }+\frac{1}{2}{\mathbf{u}}^T\mathbf{u}\right),$$

(3)

and

$$\rho e^{ve}=\sum _s{\rho }_s\left(e_s^{vib}+e_s^{el}\right).$$

(4)

Considering a general gas mixture consisting of polyatomic, monatomic, and free electron species, expressions for the energy stored in the translational, rotational, vibrational, and electronic modes are given as

$$e_s^{tr}=\left\{\begin{array}{cc}\frac{3}{2}\frac{R}{M_s}T\hfill & \mathrm{for} \mathrm{monatomic} \mathrm{and} \mathrm{polyatomic} \mathrm{species},\hfill \\ 0\hfill & \mathrm{for} \mathrm{electrons},\hfill \end{array}\right.$$

(5)

$$e_s^{rot}=\left\{\begin{array}{cc}\frac{\xi }{2}\frac{R}{M_s}T\hfill & \mathrm{for} \mathrm{polyatomic} \mathrm{species},\hfill \\ 0\hfill & \mathrm{for} \mathrm{monatomic} \mathrm{species} \mathrm{and} \mathrm{electrons},\hfill \end{array}\right.$$

(6)

where $\xi$ is an integer specifying the number of axes of rotation,

$$e_s^{vib}=\left\{\begin{array}{cc}\frac{R}{M_s}\frac{{\theta }_s^{vib}}{exp\left({\theta }_s^{vib}/T^{ve}\right)-1}\hfill & \mathrm{for} \mathrm{polyatomic} \mathrm{species},\hfill \\ 0\hfill & \mathrm{for} \mathrm{monatomic} \mathrm{species} \mathrm{and} \mathrm{electrons},\hfill \end{array}\right.$$

(7)

where ${\theta }_s^{vib}$ is the characteristic vibrational temperature of the species, and

$$e_s^{el}=\left\{\begin{array}{cc}\frac{R}{M_s}\frac{{\sum }_{i=1}^{\infty }{g}_{i,s}{\theta }_{i,s}^{el}exp(-{\theta }_{i,s}^{el}/T_{ve})}{{\sum }_{i=0}^{\infty }{g}_{i,s}exp(-{\theta }_{i,s}^{el}/T_{ve})}\hfill & \mathrm{for} \mathrm{polyatomic} \mathrm{and} \mathrm{monatomic} \mathrm{species},\hfill \\ \frac{3}{2}\frac{R}{M_s}{T}^{ve}\hfill & \mathrm{for} \mathrm{electrons},\hfill \end{array}\right.$$

(8)

where ${\theta }_s^{el}$ is the characteristic electronic temperature of the species and $g_i$ is the degeneracy of the $i\mathrm{th}$ state.

Vibrational–Electronic Relaxation

Vibrational relaxation is computed using a standard Landau–Teller [22] relaxation time with Park high-temperature correction

$${\dot{\mathsf{\Theta }}}^{tr:ve}=\sum _s{\rho }_s\frac{d{e}_s^{ve}}{dt}=\sum _s{\rho }_s\frac{e_s^{ve*}-e_s^{ve}}{{\tau }_s},$$

(9)

where ${\tau }_s$ is computed using a combination of the Landau–Teller relaxation time, ${\langle {\tau }_s\rangle }_{L-T}$, and a limiting relaxation time from Park, ${\tau }_{ps}$ using

$${\tau }_s={\langle {\tau }_s\rangle }_{L-T}+{\tau }_{ps},$$

(10)

and

$${\langle {\tau }_s\rangle }_{L-T}=\frac{{\sum }_r{X}_r}{{\sum }_r{X}_r/{\tau }_{sr}}.$$

(11)

The interspecies relaxation times are taken from experimental data from Millikan and White [23], expressed as

$${\tau }_{sr}=\frac{1}{P}exp\left[A_{s}r\left(T^{-1/3}-0.015{\mu }_{sr}^{1/4}\right)-18.42\right].$$

(12)

A limiting relaxation time, ${\tau }_{ps}$, is used to correct for under-prediction of the Millikan–White model at high temperatures [23]. ${\tau }_{ps}$ is defined as

$${\tau }_{ps}=\frac{1}{{\sigma }_s{c}_{s}n},$$

(13)

where ${\sigma }_s$ is the effective collision cross-section.

3.2. Finite-Rate Chemical Kinetics

Energetic collisions also result in chemical reactions taking place in the flow. The source terms in the species conservation equations described in Equation (2) are the species’ volumetric mass production rates which are governed by the forward and backward reaction rates, $R^f$ and $R^b$, for a given reaction r, and can be expressed as

$${\dot{\omega }}_s=M_s\sum _r({\beta }_{s,r}-{\alpha }_{s,r})(R_r^f-R_r^b).$$

(14)

From kinetic theory, the forward and backward reaction rates are dependent on the molar concentrations of the reactants and products, as well as the forward and backward reaction rate coefficients, $k^f$ and $k^b$, respectively [21], and can be expressed as

$$R_r^f=k_r^f\prod _s{\left(\frac{{\rho }_s}{M_s}\right)}^{{\alpha }_{s,r}},$$

(15)

and

$$R_r^b=k_r^b\prod _s{\left(\frac{{\rho }_s}{M_s}\right)}^{{\beta }_{s,r}}.$$

(16)

For an Arrhenius reaction, the forward reaction rate coefficient can be computed as

$$k_r^f=C_r{\left(T_r\right)}^{{\eta }_r}exp\left(-\frac{{ϵ}_r}{k_B{T}_r}\right),$$

(17)

where $C_r$ is the pre-factor, $T_r$ is the rate-controlling temperature for the reaction, ${\eta }_r$ is an empirical exponent, and ${ϵ}_r$ is the activation energy per molecule. A list of these values can be seen in

Table A1. Arrhenius parameters for forward rate reactions.

Reaction	Cr ((cm3)/(mol-s))	ηr	${\mathit{epsilon}}_{\mathit{r}}^{\mathit{A}}$ (K)
N2 + e− ⇌ N + N + e−	3.0 × 1024	−1.6	113,200
N2 + N2 ⇌ 2N + O2	7.0 × 1021	−1.60	113,200
N2 + O2 ⇌ 2N + O2	7.0 × 1021	−1.60	113,200
N2 + NO ⇌ 2N + NO	7.0 × 1021	−1.60	113,200
N2 + N ⇌ 2N + NO	3.0 × 1022	−1.60	113,200
N2 + O ⇌ 2N + NO	3.0 × 1022	−1.60	113,200
O2 + N2 ⇌ 2O + N2	2.0 × 1021	−1.50	59,500
O2 + O2 ⇌ 2O + O2	2.0 × 1021	−1.50	59,500
O2 + NO ⇌ 2O + NO	2.0 × 1021	−1.50	59,500
O2 + N ⇌ 2O + N	1.0 × 1022	−1.50	595,200
O2 + O ⇌ 2O + O	1.0 × 1022	−1.50	595,200
NO + N2 ⇌ N + O + N2	5.0 × 1015	0.00	75,500
NO + O2 ⇌ N + O + O2	5.0 × 1015	0.00	75,500
NO + NO ⇌ N + O + NO	5.0 × 1015	0.00	75,500
NO + N ⇌ N + O + N	1.1 × 1017	0.00	75,500
NO + O ⇌ N + O + O	1.1 × 1017	0.00	75,500
N2 + O ⇌ NO + N	6.4 × 1017	−1.00	38,400
NO + O ⇌ O2 + N	8.4 × 1017	0.00	19,450
N + O ⇌ NO+ + e−	5.3 × 1012	0.0	31,900

.

The rate-controlling temperature of the reaction is calculated as a geometric average of the translation–rotational and vibrational–electronic temperatures,

$$T_r={\left(T\right)}^{a_r}{\left(T^{ve}\right)}^{b_r},$$

(18)

where the constants $a_r$ and $b_r$, shown in

Table A2. Reaction rate coefficient controlling temperature factors.

Reaction	${\mathbf{a}}_{\mathit{r}}^{\mathit{f}}$	${\mathbf{b}}_{\mathit{r}}^{\mathit{f}}$	${\mathbf{a}}_{\mathit{r}}^{\mathit{b}}$	${\mathbf{b}}_{\mathit{r}}^{\mathit{b}}$
N2 + M ⇌ 2N + O2	0.5	0.5	1.0	0.0
N2 + M ⇌ 2O + O2	0.5	0.5	1.0	0.0
NO + M ⇌ N + O + NO	0.5	0.5	1.0	0.0
N2 + O ⇌ NO + N	1.0	0.0	1.0	0.0
NO + O ⇌ O2 + N	1.0	0.0	1.0	0.0
N + O ⇌ NO+ + e−	1.0	0.0	1.0	0.0

, are dependent on the nature of the reaction. The forward and backward rate coefficients are related by the reaction equilibrium constant which is determined using curve fits from the Park 1990 model [24], shown in

Table A3. Parameters for the Park 1990 equilibrium reaction constants.

Reaction	N (1/cm3)	A0	A1	A2	A3	A4
N2 + M ⇌ 2N + M	1 × 1014	3.4907	0.83133	4.0978	−12.728	0.7487
	1 × 1015	2.0723	1.38970	2.0617	−11.828	0.015105
	1 × 1016	1.6060	1.57320	1.3923	−11.533	−0.004543
	1 × 1017	1.5351	1.60610	1.2993	−11.494	−0.00698
	1 × 1018	1.4766	1.62910	1.2153	−11.457	−0.00944
	1 × 1019	1.4766	1.62910	1.2153	−11.457	−0.00944
O2 + M ⇌ 2O + M	1 × 1014	1.8103	1.9607	3.5716	−7.3623	0.083861
	1 × 1015	0.91354	2.3160	2.2885	−6.7969	0.046338
	1 × 1016	0.64183	2.4253	1.9026	−6.6277	0.035151
	1 × 1017	0.55388	2.4600	1.7763	−6.5720	0.031445
	1 × 1018	0.52455	2.4715	1.7342	−6.55534	0.030209
	1 × 1019	0.50989	2.4773	1.7132	−6.5441	0.29591
NO + M ⇌ N + O + M	1 × 1014	2.1649	0.078577	2.8508	−8.5422	0.053043
	1 × 1015	1.0072	0.53545	1.1911	−7.8098	0.004394
	1 × 1016	0.63817	0.68189	0.66336	−7.5773	−0.011025
	1 × 1017	0.55889	0.71558	0.55396	−7.5304	−0.014089
	1 × 1018	0.5150	0.73286	0.49096	−7.5025	−0.015938
	1 × 1019	0.50765	0.73575	0.48042	−7.4979	−0.016247
N2 + O ⇌ N + O + M	1 × 1014	1.3261	0.75268	1.2474	−4.1857	0.02184
	1 × 1015	1.0653	0.85417	0.87093	−4.0188	0.010721
	1 × 1016	0.96794	0.89131	0.7291	−3.9555	0.006488
	1 × 1017	0.97646	0.89043	0.74572	−3.9642	0.007123
	1 × 1018	0.96188	0.89617	0.72479	−3.955	0.006509
	1 × 1019	0.96921	0.89329	0.73531	−3.9596	0.006818
NO + O ⇌ O2 + N	1 × 1014	0.35438	−1.8821	−0.72111	−1.1797	−0.30831
	1 × 1015	0.093613	−1.7806	−1.0975	−1.0128	−0.41949
	1 × 1016	−0.003732	−1.7434	−1.2394	−0.94952	−0.046182
	1 × 1017	0.004815	−1.7443	−1.2227	−0.95824	−0.45545
	1 × 1018	−0.009758	−1.7386	−1.2436	−0.949	−0.046159
	1 × 1019	−0.002428	−1.7415	−1.2331	−0.95365	−0.04585
N + O ⇌ NO+ + e−	1 × 1014	−2.1852	−6.6709	−4.2968	−2.2175	−0.50748
	1 × 1015	−1.0276	−7.1278	−2.637	−2.95	−0.0021
	1 × 1016	−0.65871	−7.2742	−2.1096	−3.1823	0.01331
	1 × 1017	−0.57924	−7.3079	−1.9999	−3.2294	0.016382
	1 × 1018	−0.53538	−7.3252	−1.937	−3.2572	0.01823
	1 × 1019	−0.52801	−7.3281	−1.9264	−3.2618	0.01854

.

3.3. Transport Properties

Mass, momentum, and energy transport in fluids are all governed by molecular collisions, and expressions for these transport properties can be derived from the kinetic theory. The mass diffusion fluxes, ${\mathbf{J}}_s$, are computed using Fick’s Law of Diffusion:

$${\mathbf{J}}_s=\rho D_s\nabla \left(c_s\right),$$

(19)

where $c_s$ is the species mass fraction and $D_s$ is the species multi-component diffusion coefficient. The values of $D_s$ are computed as a weighted sum of binary diffusion coefficients between all species in the mixture. These are obtained by solving the Stefan–Maxwell equations under the Ramshaw approximations [25]. The viscous stress tensor is written as

$$\sigma =\mu \left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\mathsf{T}}-\frac{2}{3}\mathbf{I}(\nabla ·\mathbf{u})\right),$$

(20)

where $\mu$ is the mixture viscosity coefficient. The conduction heat flux for each thermal energy mode, ${\mathbf{q}}^k$, is assumed to be given by Fourier’s Law of heat conduction:

$${\mathbf{q}}^k={\lambda }^k\nabla \left(T^k\right),$$

(21)

where $T^k$ is the temperature and ${\lambda }^k$ is the thermal conductivity coefficient of the kth energy mode. Viscosity and thermal conductivity are computed using Wilke’s mixing rule [26]. The species viscosity model is calculated using either Blottner’s three-parameter curve fits for high-temperature air [27], or Gupta–Yos [28]. Thermal conductivity is calculated using Eucken’s formula [29].

3.4. Turbulence Modeling

One of the core features in the SU2 suite is a Reynolds-Averaged-Navier–Stokes (RANS) solver to solve compressible turbulent flow [30]. The RANS equations are solved in SU2 using two main models: the one-equation Spallart-Allmaras (SA) model [31] and the two-equation Menter k-omega SST model [32]. Within the SU2-NEMO framework, the SA model has been implemented. The turbulent eddy viscosity is computed as

$${\mu }_t=\rho \widehat{\nu }{f}_{v1},$$

(22)

where $\widehat{\nu }$ is solved for using

$$\frac{\partial \widehat{\nu }}{\partial t}+{\mathbf{u}}_j\frac{\partial \widehat{\nu }}{\partial {\mathbf{x}}_j}=c_{b1}(1-f_{t2})\widehat{S}\widehat{\nu }-\left[c_{w1}{f}_w-\frac{c_{b1}}{k^2}{f}_{t2}\right]{\left(\frac{\widehat{\nu }}{d}\right)}^2+\frac{1}{\sigma }\left[\frac{\partial }{\partial {\mathbf{x}}_j}\left((\nu +\widehat{\nu })\frac{\partial \widehat{\nu }}{\partial {\mathbf{x}}_j}\right)+c_{b2}\frac{\partial \widehat{\nu }}{\partial {\mathbf{x}}_i}\frac{\partial \widehat{\nu }}{\partial {\mathbf{x}}_i}\right].$$

(23)

The SA model was not designed for multiple species being present in a given flow, and as such, does not account for the multiple species’ densities. This work takes the total density of the flow as

$$\rho =\sum {\rho }_s.$$

(24)

Hypersonic turbulence modeling remains a challenge in the community, and work is being done to improve and validate models within the NEMO framework.

3.5. Modeling of the Slip Regime

One of the most recent implementations in SU2-NEMO is the capability of simulating flows with moderate rarefaction, for which the Navier–Stokes continuum equations hold for the whole domain except at the vicinity of solid surfaces. Near the wall, due to the lack of collisions, the gas flow will not have the same properties as the surface, that is, wall temperature and velocity. Thus, there is the need to account for the molecular slippage. SU2-NEMO uses the Maxwell velocity slip [33] and Smoluchowski temperature jump [34] equations to compute the velocity and temperature of the gas in contact with the surface. The equations are given as

$$U_s=\frac{2-\sigma }{\sigma }\lambda \frac{\partial U}{\partial n}+\frac{3}{4}\frac{\mu }{\rho T}\frac{\partial T}{\partial x},$$

(25)

and

$$T-T_w=\frac{2-\alpha }{\alpha }\lambda \frac{2\gamma }{(\gamma +1)Pr}\frac{\partial T}{\partial n},$$

(26)

respectively, where $\mu$ is the flow viscosity, $\rho$ is the mixture density, $Pr$ is the Prandtl number, $\gamma$ is the specific heat ratio, T is the temperature of the gas, $T_w$ is the temperature of the surface, and $\lambda$ is the mean free path, calculated as [35]

$$\lambda =\frac{\mu }{\rho }\frac{\pi }{\sqrt{2RT}}.$$

(27)

The coefficients $\sigma$ and $\alpha$ are referred to as the Tangential Momentum Accommodation Coefficient (TMAC) and the Thermal Accommodation Coefficient (TAC), respectively. The values of the accommodation coefficients depend on the physical characteristics of the surface, and are usually determined empirically [36].

4. Numerical Implementation

This section highlights the numerical implementation of models within SU2-NEMO. This includes both the discretization of the governing equations and time-integration strategies. The basic numerical procedures are consistent with those implemented in the base SU2 software, and more details can be found in the Reference [13]. However, SU2-NEMO employs some specific convective schemes well-suited to high-speed flows, which are described in more detail.

4.1. Spatial Integration

SU2 utilizes the Finite Volume Method to solve the discretized governing equations on an edge-based median dual-grid numerical mesh. The discretized conservation equations can be written for a control volume ${\mathsf{\Omega }}_i$ as

$${\int }_{{\mathsf{\Omega }}_i}\frac{\partial {\mathbf{U}}_i}{\partial t} d\mathsf{\Omega }+\sum _{j\in N\left(i\right)}({\widehat{\mathbf{F}}}_{ij}^c+{\widehat{\mathbf{F}}}_{ij}^v)\mathsf{\Delta }{S}_{ij}-\mathbf{Q}\left|{\mathsf{\Omega }}_i\right|={\int }_{{\mathsf{\Omega }}_i}\frac{\partial {\mathbf{U}}_i}{\partial t} \mathrm{d}\mathsf{\Omega }+\mathbf{R}\left({\mathbf{U}}_i\right)=0,$$

(28)

where ${\mathbf{U}}_i$ is the vector of conservative variables and $\mathbf{R}\left({\mathbf{U}}_i\right)$ is the residual representing the integration of all spatial terms at node i. ${\widehat{\mathbf{F}}}_{ij}^c$ and ${\widehat{\mathbf{F}}}_{ij}^v$ are the numerical approximations of the convective and viscous fluxes, and $\mathbf{Q}$ is a vector of source terms. $\mathsf{\Delta }{S}_{ij}$ is the area of the face associated with the edge $ij$, ${\mathsf{\Omega }}_i$ is the volume of the dual control volume, and $N\left(i\right)$ is the set of neighboring nodes to node i. The fluxes are computed at the midpoint of each edge and added/subtracted to the residual for each of the two nodes making up a particular edge.

4.1.1. Convective Flux

SU2-NEMO, like the base SU2, can be discretized using both upwind and central convective schemes. This section will not discuss each scheme at length, but will focus on two schemes implemented specifically for high-speed flow simulation: the Modified Steger-Warming (MSW) and Advection Upstream Splitting Method (AUSM) flux-vector splitting schemes.

• Modified Steger-Warming

The original Steger-Warming flux-vector splitting scheme is an upwinding numerical method for resolving flows in the presence of large shocks. The scheme involves splitting the flux vector based on the direction of information propagation in the flow [37]. MSW is widely used due to its highly dissipative nature, mitigating convergence issues due to the stiffness of nonequilibrium equations. The approximated flux between nodes i and j is given by

$${\widehat{\mathbf{F}}}_{ij}={\left(A^c·\mathbf{n}\right)}_i^{+}{\mathbf{U}}_i+{\left(A^c·\mathbf{n}\right)}_j^{-}{\mathbf{U}}_j.$$

(29)

The flux Jacobian can be diagonalized through a similarity transformation, as

$${\left(A^c·\mathbf{n}\right)}^{\pm }=\mathcal{P}{\Lambda }^{\pm }{\mathcal{P}}^{-1}.$$

(30)

To address accuracy concerns of the Steger-Warming scheme, Modified Steger-Warming adds a weighting factor to reduce the numerical dissipation in the presence of large pressure gradients, such as near-shocks and in boundary layers, to ensure greater accuracy without impacting flow stability. The weighting factor is calculated as

$$w=\frac{1}{2}\frac{1}{{\left(ϵ\nabla p\right)}^2+1},$$

(31)

where $ϵ$ is a tunable parameter, taken as 5 by default in SU2-NEMO. The weighting factor is used to compute the corrected primal state vectors as

$${\mathbf{U}}_{ij}=\left(1-w\right){\mathbf{U}}_i+w{\mathbf{U}}_j,$$

(32)

and

$${\mathbf{U}}_{ji}=\left(1-w\right){\mathbf{U}}_j+w{\mathbf{U}}_i.$$

(33)

• Advection Upstream Splitting Method

AUSM is another flux-vector splitting method commonly used for high-speed flows. The standard AUSM scheme [38] is described in detail below. The scheme separates the flux vector into convective terms propagated at the local flow velocity and pressure terms propagated at the local speed of sound, expressed as

$${\widehat{\mathbf{F}}}_{ij}·\mathbf{n}=\left[\begin{array}{c}{\rho }_s\\ \rho \mathbf{u}\\ \rho h\\ \rho e^{ve}\end{array}\right]\left(\mathbf{u}·\mathbf{n}\right)+\left[\begin{array}{c}0\\ p\mathcal{I}·\mathbf{n}\\ 0\\ 0\end{array}\right],$$

(34)

where the Mach number can be split as

$$M_{ij}=M_i^{+}+M_j^{-},$$

(35)

and $M^{\pm }$ are defined by using Van Leer splitting [39]. Pressure is similarly discretized as

$$P_{ij}=P_i^{+}+P_j^{-},$$

(36)

such that

$$P^{\pm }=\left\{\begin{array}{c}\frac{p}{4}{(M\pm 1)}^2(2\mp M) \mathrm{for} \left|M\right|\le 1.\hfill \\ \frac{p}{4}(M\pm |M\left|\right)/M \mathrm{otherwise}.\hfill \end{array}\right.$$

(37)

The flux approximation is given by

$${\widehat{\mathbf{F}}}_{ij}=\frac{(M_{ij}+|M_{ij}\left|\right)}{2}{\left[\begin{array}{c}{\rho }_{s}c\\ \rho \mathbf{u}c\\ \rho hc\\ \rho e^{ve}c\end{array}\right]}_i+\frac{(M_{ij}-|M_{ij}\left|\right)}{2}{\left[\begin{array}{c}{\rho }_{s}c\\ \rho \mathbf{u}c\\ \rho hc\\ \rho e^{ve}c\end{array}\right]}_j+P_{ij}\left[\begin{array}{c}0\\ \mathcal{I}·\mathbf{n}\\ 0\\ 0\end{array}\right].$$

(38)

The resulting scheme is accurate and robust, with lower numerical dissipation than MSW. In addition to the scheme described above, several extensions of AUSM have also been implemented in SU2-NEMO. These include the AUSM+M [40] and AUSM⁺ up2 [41] schemes, which contain improvements on the stability characteristics and accuracy of the original AUSM scheme. The AUSM family of schemes offer superior shock-capturing and avoid the presence of carbuncles observed in stagnation regions around blunt bodies. SU2-NEMO, like the base SU2 code, has slope-limiters for higher-order accuracy solutions.

4.1.2. Viscous Flux

Viscous flux values are computed at the median dual-grid interfaces. Gradient information required for the viscous fluxes can be evaluated by either the Gauss-Green Theorem or Weighted Least-Squares approach, with appropriate corrections in areas of high cell skewness.

4.2. Boundary Conditions

SU2-NEMO employs several wall boundary conditions similar to those in the base SU2 code. These include non-catalytic isothermal and heat flux wall boundaries for viscous flow simulation, as well as the inviscid Euler wall. The Smoluchowski–Maxwell boundary condition, as described in [42], is also implemented for modeling rarefied flows. A more detailed description of the slip boundary condition implementation is done in Section 3.5.

Surface-catalyzed chemical reactions can have a significant impact in high-speed flows, and current efforts focusing on the development and validation of more complex gas–surface models in SU2-NEMO are ongoing.

4.3. Time Integration

SU2-NEMO takes advantage of the time integration available in the general SU2 framework, utilizing both explicit time integration (explicit forward Euler and Runge–Kutta explicit), and implicit time integration (implicit backward Euler).

Explicit time integration schemes are easily parallelizable, since they involve a local update at each cell. The drawback is that explicit schemes place requirements on the allowable Courant–Friedrichs–Lewy (CFL) number of the local time-step based on wave speed. The CFL number must be sufficiently small, which can result in a large number of iterations to converge, especially in regions of large gradients.

Implicit methods incorporate information at the next time-step, so they do not suffer from the same stability concerns as explicit methods. However, they require the simultaneous solution of all residual equations by means of a linear system, resulting in a higher computational cost per time-step. Several linear solvers are implemented in SU2, as well as a suite of preconditioners to reduce the computational cost of the solution.

5. Results and Discussion

This section presents results from several verification and validation cases highlighting the capabilities of the SU2-NEMO framework’s numerical schemes, thermochemical libraries, and flow models. All cases were run using Euler explicit time-integration, unless stated otherwise.

5.1. Zero-Dimensional Thermal Bath

The zero-dimensional adiabatic thermal bath is a canonical verification case in the hypersonics community [43,44], which isolates the impacts of thermal nonequilibrium and finite-rate chemistry. A heat bath in thermochemical nonequilibrium is simulated using the N₂ and Air-5 gas models within both the SU2 native library and the Mutation++ library. In the cases run using Mutation++, both preferential and non-preferential dissociation models were employed. The Air-5 gas results are compared to results generated using the LeMANS code [45]. The thermal bath is simulated using a 5 × 5 structured mesh with symmetry planes on each side, and an unsteady explicit time-stepping method is used for these simulations.

5.1.1. N₂ Gas Model

In the first case, a N₂ gas mixture is simulated with an initial concentration of pure diatomic nitrogen. The initial conditions are shown below in

Table 1. N₂ thermal bath initial conditions.

${\mathit{T}}_{\mathbf{\infty }}^{\mathit{tr}}$ (K)	${\mathit{T}}_{\mathbf{\infty }}^{\mathit{ve}}$ (K)	P∞ (atm)	$\mathit{Y}{\left[{\mathit{N}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathit{N}\right]}_{\mathbf{\infty }}$
20,000	300	27.25	1.0	0.0

. The comparison of the native (CSU2TCLib) and Mutation++ (CMutationTCLib) CFluidModel classes for temperature relaxation and dissociation are shown in Figure 4a,b, respectively.

The SU2-Mutation++ simulation was significantly faster, running in 26.38 s compared to 73.1 s for SU2-Native in serial on a single CPU. Overall, there is close agreement between the thermal relaxation results of both the SU2-Native and SU2-Mutation++ simulations, as can be seen in Figure 4a. This agreement is to be expected as both codes use similar thermochemical models. All the simulations performed using the SU2-Native and Mutation++ library (preferential and non-preferential) reach an equilibrium temperature of 7500 K at 10⁻⁶ s. The SU2-Native data, which uses a non-preferential model, is bounded by the SU2-Mutation++ preferential and non-preferential data, with $T^{ve}$ reaching a peak value of 11,870 K compared to 11,690 K and 12,100 K, respectively. Comparing the non-preferential data, the peak vibrational–electronic temperature varies by 2.8%. This difference is attributed to minor differences in the characteristic values chosen within the respective models. The non-preferential data exhibit a larger overshoot in both temperature relaxation time-evolutions compared to the preferential data. As a result, the non-preferential simulations have marginally higher rates of dissociation, as seen in Figure 4b. This behavior is caused by the non-preferential model driving the creation and destruction of species at the average $T^{ve}$ in the cell, rather than the maximum $T^{ve}$; the preferential model will lead to higher dissociation at higher $T^{ve}$. In the context of the N₂ thermal bath, the impact of the preferential and non-preferential models on dissociation is minimal.

5.1.2. Air-5 Gas Model

The second case simulates a standard five-species air model (Air-5) with initial conditions seen below in

Table 2. Air-5 thermal bath initial conditions.

${\mathit{T}}_{\mathbf{\infty }}^{\mathit{tr}}$ (K)	${\mathit{T}}_{\mathbf{\infty }}^{\mathit{ve}}$ (K)	P∞ (atm)	$\mathit{Y}{\left[{\mathit{N}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[{\mathit{O}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathbf{NO}\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathit{N}\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathit{O}\right]}_{\mathbf{\infty }}$
15,000	300	20.42	0.767	0.233	0.0	0.0	0.0

. The results for temperature relaxation and dissociation can be seen in Figure 5a,b.

Both the SU2-Mutation++ and LeMANS simulation data [46] show significant agreement, reaching an equilibrium temperature of 6200 K at 3.5 × 10⁻⁷ s, as shown in Figure 5a. Unlike the N₂ data, the differences between the non-preferential and preferential models are more significant. The non-preferential models predict greater peak $T^{ve}$ values relative to the preferential models: 11.3% for SU2-Mutation++, and 11.8% for LeMANS. Furthermore, both the SU2-Mutation++ models predict higher peak temperatures than the LeMANS models. For example, an increase of 3.5% in the peak non-preferential model vibrational temperature is observed. The differences directly correlate to the dissociation in the Air-5 models; the impact of the preferential dissociation is evident in Figure 5b. In general, the preferential models begin dissociating at a slower rate, due to the preference for higher $T^{ve}$. This is most notable in atomic oxygen, around 10⁻⁸ s. The percentage of oxygen in the flow is tabulated in

Table 3. Atomic oxygen percentage at 10⁻⁸ s.

	Non-Preferential	Preferential
SU2-Mutation++	14.91%	13.27%
LeMANS	17.48%	15.42%

.

The differences in thermal relaxation between the SU2 and LeMANS codes are mirrored in the dissociation results. Notably, NEMO-Mutation++ dissociates at a slower rate than LEMANs. This is believed to be caused by differences in implementations of each code. Although not shown above, simulations were also run using SU2-Native to characterize speed differences within SU2-NEMO. SU2-Mutation++ was ∼4.3 times faster than SU2-Native, running in 35.6 s compared to 154.6 s on a single CPU. In summary, SU2-NEMO’s thermochemical libraries, CSU2TCLib and CMutationTCLib, are compared to each other and verified against the LeMANS code. As expected, all models perform similarly, with minor differences in relaxation times seen between the SU2 and Mutation++ models and LeMANS. The CMutationTCLib was also shown to be 3–4.3 faster than CSU2TCLib for this test case.

5.2. HEG Cylinder

Experiments from the High-Enthalpy shock tunnel in Göttingen (HEG) [47] are commonly used as a basis for validation of the physio-chemical models implemented in CFD codes, where nonequilibrium chemical relaxation processes occur within the shock layer, affecting the density distribution of the species. SU2-NEMO is compared with experimental results of flow past a cylinder with a radius of 45 mm using the free-stream conditions in

Table 4. HEG free-stream conditions.

U∞ (m/s)	p∞ (Pa)	ρ∞ (kg/m3)	T∞ (K)	M∞	$\mathit{Y}{\left[{\mathit{N}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[{\mathit{O}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathbf{NO}\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathit{N}\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathit{O}\right]}_{\mathbf{\infty }}$
5956	476	1.547 × 10−3	901	8.98	0.7543	0.00713	0.01026	6.5 × 10−7	0.2283

.

A quadrilateral grid, consisting of 68,541 nodes and 68,000 elements, is used to avoid the generation of entropy errors due to the misalignment of local elements with the shock, which could propagate and contaminate the solution in the boundary layer region, thus invalidating the results. This issue is more concerning in the stagnation region of blunt bodies, since there is little dissipation of entropy gradients, and as such, higher levels of attention in the generation of the grid is required. An isothermal boundary condition with no catalytic effects is employed at the cylinder surface, with a wall temperature of 300 K. The resulting temperature contours can be observed in Figure 6. As expected, both the translational–rotational and vibrational–electronic temperatures are greatest immediately behind the shock, then quickly begin to relax, reaching equilibrium at the sphere surface.

The surface pressure and heat flux quantities are plotted in Figure 7 and compared with experimental data. The pressure plot is a near-exact fit with the available experimental data. However, heat flux quantities are slightly under-predicted, particularly in the stagnation region, when compared to the experimental data. This under-prediction is attributed to the lack of catalytic effects; according to Knight et al. [48], catalytic effects must be considered to accurately predict surface heat flux quantities. In comparison with other non-catalytic results from Lani et al. and Nompelis [48], the SU2-NEMO solution has the best agreement with experimental data in the stagnation region, showing a 2.79% peak increase relative to Lani’s heat flux peak value.

5.3. RAM-C II Test Vehicle

The RAM-C II vehicle is a blunted cone geometry that was used to study the impacts of plasma formation on radio-attenuation during atmospheric re-entry [49]. The nose radius is 0.1524 m, with a cone angle of 9 degrees. The RAM-C II flight test case has become a common verification and validation test for weakly-ionized flows. SU2-NEMO is used to simulate the RAM-C vehicle using an Air-7 Mutation++ gas model, consisting of N₂, O₂, N, O, NO, NO⁺, and e⁻, with viscous effects absent. A mesh comprised of 316,910 nodes and 267,241 elements is used, with 160 elements normal to the vehicle body. The free-stream conditions for the 61 km flight can be seen in

Table 5. RAM-C free-stream conditions.

M∞	P∞ (Pa)	ρ∞ (kg/m3)	T∞ (K)	$\mathit{Y}{\left[{\mathit{N}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[{\mathit{O}}_2\right]}_{\mathbf{\infty }}$
23.9	19.7	2.7024 × 10−4	254	0.767	0.233

.

The presence of the free-electron species in the Air-7 model allows for prediction of plasma density in the vicinity of the vehicle. While more complex chemistry models, like the Air-11 model, can also be used to more accurately model the flow, Air-7 is sufficient for this case, as NO⁺ production, due to dissociative recombination of monatomic nitrogen and oxygen, is the primary mechanism for ion formation at these flight conditions. Because the plasma density is low and the primary focus is on prediction of flow properties, magnetohydrodynamic effects are neglected in this simulation. Results can be seen in Figure 8 and Figure 9.

In Figure 8a, a strong bow shock that develops at the nose of the vehicle is observed. High temperatures attained across this shock drive chemical reactions, resulting in significant variation in the flow chemistry along the stagnation line, as illustrated in Figure 9a. Diatomic nitrogen and oxygen are depleted, while concentrations of monatomic oxygen and nitrogen increase. The mass fraction of nitric oxide is relatively small in comparison; it peaks in the region immediately behind the shock and quickly diminishes. While temperatures are sufficient to result in ionization through dissociative recombination of monatomic nitrogen and oxygen, the resulting concentrations are orders of magnitude smaller than the neutral species.

Figure 8b shows contours of electron density around the vehicle. Electron density is highest in the vicinity of the vehicle forebody, where the shock strength is largest. As the flow is accelerated around the forebody, electron density continues to be highest in the vicinity of the shock. Figure 9b shows a plot of maximum electron number density in the wall-normal direction, compared to other numerical results as well as experimental measurements. The experimental measurements of electron density were collected along the axial length of the vehicle using reflectometers and an electrostatic rake [50]. SU2-NEMO predicts similar electron number density values to Scalabrin [45] and Candler and MacCormack [51], however these predictions begin to diverge past an axial point of approximately x/R ≈ 1. Candler and MacCormack predict higher values along the flank of the vehicle, whereas SU2-NEMO and Scalabrin more closely match the reflectometer measurements. SU2-NEMO predictions show strong agreement with the reflectometer measurements past x/R ≈1.5. Josyula and Bailey utilized a six-temperature model, which allows for superior agreement with experimental results over the forebody [52]. All simulation results are within the uncertainty of the electrostatic rake measurement at the aft of the vehicle. The under-prediction of electron density in the nose region is believed to be a consequence of vibrational and electronic modes being considered together in the two-temperature model. It is expected that improved prediction in this region can be achieved by further decomposing the energy modes into a higher-temperature model, as Josyula and Bailey did, with the trade-off of increased computation time.

5.4. Axisymmetric Shock-Wave Boundary Layer Interaction

Accurate prediction of flow separation and turbulent heating augmentation is crucial for the analysis and design of hypersonic vehicles, particularly for endo-atmospheric flight. Computational prediction of these values continues to be an active area of research. Shock-wave boundary-layer interaction is a phenomenon known to encompass both separation and turbulent heating effects, making it an excellent candidate for V&V of turbulence modeling. The NASA Turbulence Modeling Resource’s (TMR) high-speed axisymmetric shock-wave boundary-layer interaction (ASWBLI) test case models flow around a conical/ogive cylinder forebody with an axisymmetric flared compression corner, depicted in Figure 10. This test case is based on a set of experiments performed by Kussoy and Horstman at the NASA Ames Research Center, where experimental measurements of surface pressure and heat flux were obtained for a range of compression angles between 20 and 35 degrees [53].

In this study, a three-dimensional mesh constituting a 1 degree wedge of a body of revolution with periodic boundary conditions was used. The mesh consists of 514,082 nodes and 256,000 quadrilateral volume elements with a wall spacing of 2.5 μm. This represents the finest grid available from the NASA TMR website, and has been demonstrated to produce grid-independent solutions [54]. The flare angle chosen for this study was 20 degrees. The surface of the body was modeled using a non-catalytic isothermal boundary condition at 311 K. Free-stream values were applied at the inlet and top boundary, while a standard extrapolation was used in the supersonic portion of the outlet. The free-stream conditions were chosen to be consistent with Georgiadis et al. [54] for simulating the flow without the conical forebody, and are provided in

Table 6. ASWBLI free-stream conditions.

M∞	P∞ (Pa)	T∞ (K)	Twall (K)	$\mathit{Y}{\left[{\mathit{N}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[{\mathit{O}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathbf{NO}\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathit{N}\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[\mathit{O}\right]}_{\mathbf{\infty }}$
7.11	550.13	80.0	311.0	0.77	0.23	0.0	0.0	0.0

. The Native Air-5 gas model is used for this simulation.

The AUSM convective scheme is utilized for spatial discretization, and diffusive fluxes are evaluated using a Weighted Least-Squares approach. Numerical simulation is performed using the SA turbulence model implemented in SU2-NEMO in conjunction with a Sutherland viscosity model for air.

Contours of Mach number over the compression corner are given in Figure 11a. A curved oblique shock caused by the axisymmetric compression corner can be observed. At the base of the compression, there is a region of re-circulation, shown in Figure 11b, associated with the boundary layer separation due to the impinging shock. Capturing this region is necessary to predict the downstream flowfield, as the separation region drives the creation of a reattachment shock/compression wave. Subsequent results for this study are normalized by the surface values 6 cm upstream from the compression corner. As such, boundary layer profiles of x-velocity, Figure 12a, and total translational temperature, Figure 12b, are taken at this point for comparison.

As expected, there is little difference between the boundary layer profiles of SU2-NEMO and WIND-US. Distributions of surface pressure and heat flux, normalized by the values 6 cm upstream of the compression corner (subscript ∞), are given in Figure 13a,b, respectively. The numerical results show close agreement in heat flux and surface pressure distribution in the vicinity of the shock, with some differences in the post-shock flow features. In particular, we note both sets of numerical results under-predict the normalized pressure on the surface of the flare (x > 7 cm). SU2-NEMO predicts surface heat fluxes in a similar trend as WIND-US until 3 cm from the compression. From here, SU2-NEMO predicts a consistently higher heat flux, reaching a peak value 9.85 times the reference value and 3.46% higher than WIND-US. However, both SU2-NEMO and WIND-US over-predict the local minimum at 7 cm and fail to predict the second peak in heat flux at 12 cm. Georgiadis et al. note a strong dependence of heat flux on mesh orthogonality near the wall, even for very refined grids, and more accurate heat flux predictions can be obtained by modifying the baseline grids utilized in this study [54]. Georgiadis et al. also found the turbulent Prandtl number has a significant impact on the heat flux predictions. Using a variable Prandtl number in this case may lead to better agreement with experimental data. Measurement uncertainties of ±10% are reported in the Reference [53] for pressure and surface heat flux. Differences in peak and post-shock heat flux prediction between SU2-NEMO and WIND-US may be attributed to differences in convective scheme, as well as code-to-code differences. While initial results show good agreement with numerical and experimental results, refinement and validation of the turbulence models remains an area of ongoing research and development effort in SU2-NEMO.

5.5. Slip Flow over a Cylinder

To verify the implemented slip regime models (see Section 3.5), a hypersonic flow of nitrogen over a cylinder with a radius of 0.1524 m is simulated. Two points in the slip regime are considered; the crossover between the continuum and slip regime at a Knudsen number (Kn) of 0.01 and a more rarefied case at a Knudsen number of 0.05. The Mutation++ N₂ gas model is used to simulate this case with free-stream flow conditions illustrated in

Table 7. Cylinder with slip flow conditions.

Kn∞	M∞	U∞ (m/s)	ρ∞ (kg/m3)	T∞ (K)	Twall (K)
0.01	10	2883	1.974 × 10−5	200	500
0.05	10	2883	3.949 × 10−6	200	500

.

For both cases considered, the global Knudsen number is inside the range allowed for the use of slip-regime models. A grid-independence analysis was conducted for both cases to reach a mesh-independent solution using hybrid grids, with a total of 210,656 nodes for the Kn = 0.01 case and 207,526 nodes for the Kn = 0.05 case. The temperature contours are illustrated in Figure 14.

As the flow becomes more rarefied, the vibrational excitation tends to decrease as a result of fewer molecular collisions, which are required to equilibrate the translational and vibrational energy modes. It can be observed in Figure 14b, that for the Kn = 0.05 case, the vibrational temperature does not exceed the imposed wall temperature. The vibrational modes are primarily activated by the cylinder surface, as opposed to the bow shock. It can also be observed that the increase in rarefaction leads to an increase in the shock stand-off distance. The skin friction (Figure 15), heat flux (Figure 16), and pressure (Figure 17) coefficients obtained with SU2-NEMO are compared with the DSMC results obtained by Lofthouse [36] using the MONACO code [44].

Observing the surface plots and the correspondent peak values in

Table 8. Peak values at the surface of the cylinder.

	Kn = 0.01			Kn = 0.05
	SU2-NEMO	DSMC	Δ	SU2-NEMO	DSMC	Δ
Peak Pressure Coefficient	1.7764	1.7634	0.74 %	1.8370	1.8201	0.93 %
Peak Heat flux Coefficient	0.1480	0.1503	−1.53 %	0.3254	0.3099	4.97 %
Peak Skin friction Coefficient	0.1026	0.1045	−1.81 %	0.2173	0.2049	6.05 %

, SU2-NEMO shows good agreement with DSMC for the Kn = 0.01 case, but begins to overestimate the heat flux and skin friction coefficients as rarefaction increases. The skin friction coefficient is the most sensitive of the surface properties to the amount of continuum breakdown, which increases as the flow becomes more rarefied. These differences are most noticeable in the wake region of the cylinder for the Kn = 0.05 case, between 120–180 degrees. However, this over-prediction regarding the peak value is below 5% for heat flux coefficient and 6% for skin friction coefficient, thus the slip model implemented in SU2-NEMO can be used to reasonably simulate the slip flow regime.

5.6. NASA X-43a (Hyper-X)

The NASA X-43a was an experimental aircraft used to test hypersonic air-breathing propulsion systems. It was the first vehicle to fly freely using scramjet technology in 2004. The first test flight flew at Mach 6.8 after being air-launched from a B-52 Stratofortress and accelerated by a boost rocket. The X-43a went on to break the air-breathing propulsion aircraft airspeed record at Mach 9.6. SU2-NEMO is used to simulate a cruise portion of an X-43a-like vehicle’s flight, demonstrating the capability to simultaneously simulate complex geometries while also utilizing mesh adaptation. The Mutation++ Air-5 gas model is used with free-stream flight conditions given in

Table 9. X-43a free-stream conditions.

M∞	P∞ (Pa)	T∞ (K)	$\mathit{Y}{\left[{\mathit{N}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{\left[{\mathit{O}}_2\right]}_{\mathbf{\infty }}$	$\mathit{Y}{[-]}_{\mathbf{\infty }}$
7	190	239	0.77	0.23	0.0

. At these conditions, thermal relaxation and chemical reactions do not play a significant role in the physics predictions. This case will serve as a baseline for future work incorporating more realistic propulsion models.

To account for the internal flow of the scramjet, the engine inlet is modeled with an outlet boundary condition and the engine outlet is modeled with an inlet boundary condition. The scramjet outlet flow parameters were set to the free-stream values. This case is simulated without viscous contributions. Throughout the simulation process, adaptive mesh refinement is performed to improve the shock resolutions. Five levels of mesh adaptation were used, where the adaptation was based on gradients of the local Mach number. The original mesh was comprised of 777,365 nodes and 4,559,568 elements, while the final mesh includes 5,114,158 nodes and 28,962,633 elements. The resulting mesh is shown in Figure 18.

Due to the large separation region and resulting recompression slipline seen behind the vehicle, much of the adaptation algorithm focused on adding grid elements in that region, shown in Figure 18, often overlooking other weaker Mach gradients. A more precise grid adaptation method would allow for equally resolved flow features on a coarser mesh.

As expected, the concentrations of N, O, and NO⁺ are negligibly small, confirming that thermochemical effects had little impact on the flow features under the simulated conditions. Select Mach number contours are shown in Figure 19 and Figure 20. Attached oblique shocks are formed by the leading edge of the vehicle. A second ramp on the underside of the vehicle creates a secondary shock, compressing the air further at the inlet of the scramjet. Without the proper amount of compression, there is a significant risk of a scramjet unstart and vehicle failure. The flow behind the scramjet is accelerated, leading to a high-Mach region on the underside of the vehicle and resulting in an expansion jet visible emanating from the rear surface. This region is flanked by lower Mach regions generated by a system of shock-waves that emanate from the vertical and horizontal tails, shown in Figure 20.

The shock-waves generated by the leading edge of the vertical tails impact the horizontal tail’s leading edge, resulting in shock–shock and shock–surface interactions. These shocks interact with each other, creating a localized high Mach region on the inboard section of the horizontal wing, with a low-pressure region on the outboard section of the wing. The mesh adaptation algorithm captured the interactions, increasing the surface element density along the discontinuities, seen in Figure 20. This vast range of flow speed and resulting stagnation quantities on a control surface can lead to a loss of control authority if not properly understood. These interactions must be accounted for and accurately predicted, particularly around the engine inlet and control surfaces, in the design process, as they can lead to significant heating augmentation, engine flame-out, and/or structural failure.

As a last test case for this initial V&V study of SU2-NEMO, the X-43a-like vehicle results showcase the potential for this newly developed numerical approach to simulate the flow around complex geometries leveraging mesh adaptation. The use of unstructured meshes with mesh adaptation in SU2-NEMO paves the way for more complex analysis and design optimization studies in the future.

6. Conclusions

This paper presented the SU2-NEMO framework, detailing the multiphysics extensions to SU2 for applications in high-temperature and high-Mach flows. This paper outlined the SU2-NEMO framework design, which provides the flexibility to rapidly modify existing thermochemical models, or to couple external software suites like Mutation++ with the addition of the CNEMOGas CFluidModel class. This paper also provided details on the numerical models implemented to account for nonequilibrium and chemical phenomena required to accurately model the hypersonic flight regime.

Additionally, the SU2-NEMO thermochemical models were verified and validated using several canonical hypersonic test cases, including a zero-dimensional heat bath, the HEG hypersonic cylinder, and the weakly ionized RAM-C vehicle. A low Knudsen number cylinder has been presented, verifying the slip flow regime boundary conditions. An axisymmetric turbulent cone-flare geometry has also been simulated to demonstrate the implementation of the RANS SA turbulence model in SU2-NEMO. Finally, SU2-NEMO was used to simulate an X-43a-like vehicle, capturing the shock–surface interactions on a complex body using anisotropically adapted unstructured meshes.

SU2-NEMO has been developed to take advantage of the capabilities present in the rapidly evolving SU2 mulitphysics suite. While this paper discussed the fundamental capabilities of SU2-NEMO, future work will include expanding the framework to incorporate the use of the adjoint solver for sensitivity analysis and design optimization applications, more complex gas–surface interaction boundary conditions, and improved turbulence models.