A two-temperature open-source CFD model for hypersonic reacting :

: A two-temperature CFD (computational ﬂuid dynamics) solver is a prerequisite to any spacecraft re-entry numerical study that aims at producing results with a satisfactory level of accuracy within realistic timescales. In this respect, a new two-temperature CFD solver, hy2Foam , has been developed within the framework of the open-source CFD platform OpenFOAM for the prediction of hypersonic reacting ﬂows. This solver makes the distinct juncture between the trans-rotational and multiple vibrational-electronic temperatures. hy2Foam has the capability to model vibrational-translational and vibrational-vibrational energy exchanges in an eleven-species air mixture. It makes use of either the Park TTv model or the coupled vibration-dissociation-vibration (CVDV) model to handle chemistry-vibration coupling and it can simulate ﬂows with or without electronic energy. Veriﬁcation of the code for various zero-dimensional adiabatic heat baths of progressive complexity has been carried out. hy2Foam has been shown to produce results in good agreement with those given by the CFD code LeMANS (The Michigan Aerothermodynamic Navier-Stokes solver) and previously published data. A comparison is also performed with the open-source DSMC (direct simulation Monte Carlo) code dsmcFoam . It has been demonstrated that the use of the CVDV model and rates derived from Quantum-Kinetic theory promote a satisfactory consistency between the CFD and DSMC chemistry modules.


Introduction
The study of high-speed vehicles re-entering the Earth's atmosphere is of current interest as witnessed by the ongoing tests on the Orion capsule (see Figure 1a) [1,2].Access to space continues to be a challenging area with significant economic and scientific implications in the near-future for the leading countries.Mastering the art of the high-speed regime is not solely limited to space missions though and new prospects may emerge such as those related to hypersonic civilian transportation.This future vision of air-space transportation is embodied by vehicles, such as the cFASTT-1 [3] shown in Figure 1b designed in the Centre for Future Air-Space Transportation Technology (cFASTT) at the University of Strathclyde, flying at cruise-altitudes around ten times higher than today's aircraft.However, many scientific and technological hurdles remain to be cleared for turning this future into a reality.The technological challenges remain considerable from an engineering, environmental, and societal standpoint [4].In particular, previous missions have highlighted the harsh environment experienced by a spacecraft during the descent phase and the essential role played by thermal protection systems to preserve the integrity of the vehicle.The accurate prediction of aerodynamic and thermal loads are thus a vital prerequisite to any computational code.Despite the significant amount of studies focusing on the subject area since the 1960s [5] a thorough understanding and characterisation of the aerothermodynamic flow conditions during re-entry is yet to be resolved as no definitive models exist to describe the wide range of physical phenomena encountered by the re-entry craft.The highly complex flow field surrounding a re-entry vehicle needs to be captured by means of numerical simulations with a satisfactory level of accuracy and within realistic timescales.Between 60 and 80 km altitude, the flow field consists of a mixture of continuum and rarefied regions and that influences the numerical approach to be chosen.Non-continuum regions cannot be modelled using conventional computational fluid dynamics (CFD) as the Navier-Stokes-Fourier (NSF) equations are inadequate for accurately predicting regions with such flow behaviour [6].Conversely, the direct simulation Monte Carlo (DSMC) [7] particle-based methodology is the predominant approach for solving rarefied flows, however, significant computational resources are required when computing in the continuum region.The chemically reacting environment over the craft has also to be accounted for as it is of importance for the correct prediction of aerodynamic loads and surface heating.
The requirements in terms of CFD translate into the need for a solver that makes the distinction between the trans-rotational and the vibrational-electronic temperatures, namely the two-temperature model of Park [8,9].These temperatures are derived from the decomposition of the internal energy into its elementary energy modes.Translational-vibrational and vibrational-vibrational energy exchanges are then determined and coupled with an appropriate chemistry module.Among such CFD codes dedicated to the hypersonic regime are NASA's DPLR (Data-Parallel Line Relaxation from the National Aeronautics and Space Administration) [10], LAURA (Langley Aerothermodynamic Upwind Relaxation Algorithm) [11], VULCAN (Viscous Upwind aLgorithm for Complex flow ANalysis) [12], LeMANS (The Michigan Aerothermodynamic Navier-Stokes solver) [13,14], and US3D (UnStructured 3D) [15,16] from the University of Minnesota.For rarefied flow conditions, a DSMC solver called dsmcFoam has been developed and validated within the framework of the open-source CFD platform OpenFOAM [17] at the University of Strathclyde [18,19].The dsmcFoam code contains capabilities to effectively model the vibrational internal energy mode and chemical reactions, adopting the recent Quantum Kinetic (QK) approach introduced by Bird [20].Post-collision energy redistribution between the internal energy modes is executed using the quantum Larsen-Borgnakke procedure [18,21].
Within the same software suite, OpenFOAM, the open-source two-temperature CFD solver, named hy2Foam, to compute hypersonic reacting flows is under development.The ultimate objective being the coupling of these two distinct numerical approaches to form a hybrid CFD-DSMC solver.The work presented in this paper describes the initial stages of the verification of hy2Foam for various zero-dimensional adiabatic heat baths of increasing complexity, highlighting the contributions of competing physical mechanisms.The verification process is supported by comparisons with data from the literature, results from the CFD solver LeMANS and the DSMC code dsmcFoam.This paper concludes with the prescription of an appropriate combination of chemistry-vibration/chemical rate models for use in a hybrid CFD-DSMC (QK) solver.

Methodology
The two-temperature formulation is considered in this section with particular focus on the models that have been implemented in hy2Foam.

Two-Temperature Model
The total energy, E, can be written as the sum of the kinetic, translational, internal, electron, and chemical energies as follows [6] where indices t, r, v, el, and e refer to translational, rotational, vibrational, electronic and electron energy modes, respectively.ρ is the density and u i for i to vary between 1 and 3 represents one of three components of the velocity vector.s is an index to vary between 1 and the number of species in the mixture, except electrons, such that ρ s is the partial density of species s and h • s denotes the standard enthalpy of formation of species s.The relation between the total energy of a specific mode and the energy per unit mass of species s is expressed in Equation (2) for the vibrational energy mode The relaxation towards thermal equilibrium of the translational and rotational energy modes is usually achieved within a small number of particle collision events [20].In a two-temperature description, the translational and rotational temperatures are thus considered to be in equilibrium with each other at all times and are equal to the trans-rotational temperature, denoted by T tr .Similarly, electron and electronic energy modes are assumed to be in equilibrium with the vibrational energy mode [20].Hence, all three temperatures are forced to follow the vibrational-electron-electronic temperature designated by T ve in a single vibrational temperature formulation.Considering a multiple vibrational temperature model, which is precisely the focus of this work, each vibrationally-excited species has its own vibrational-electronic temperature, T ve,s , while the electronic and electron temperatures are set to the vibrational-electronic temperature of a reference molecule, T ve,ref .
The different energy modes per unit mass of species are detailed below for each mode.
e r = R s T tr (4) e e = R e T ve,ref R s is the specific gas constant, θ v,s is the characteristic vibrational temperature of the particle of interest.θ el,i,s and g i,s are the characteristic electronic temperature and the degeneracy degree of a given electronic state i for species s, respectively.The values employed in this paper are to be found in Tables A1 and A2.The common model of the simple harmonic oscillator which assumes equally-spaced vibrational levels is adopted [20].
Two source terms are introduced into the NSF equations to quantify energy redistributions between the trans-rotational and the vibrational-electronic energy modes (V-T) and in between the vibrational-electronic energy modes of the different vibrationally-excited species (V-V).
The source term to model V-T energy exchange, and denoted by Q m, V−T , is governed by the Landau-Teller equation [22].This equation prescribes the V-T energy exchange rate and thus the time rate of change in vibrational energy as follows where N m is the list of the molecules in the mixture, τ m, V−T is the molar-averaged V-T relaxation time and e ve,m is the vibrational-electronic energy per unit mass of species m, for m ∈ N m , such that e ve,m = e v,m + e el,m .Hence, it also follows that E ve,m = E v,m + E el,m for m ∈ N m .Millikan and White [23] introduced a semi-empirical correlation to evaluate the V-T relaxation time for individual or mixture of species.The temperature range of validity of the proposed formula is from 300 K up to 8000 K. Within this range, the estimation doesn't deviate from experimental measurements by more than a factor of 5. Park added a correction factor to the Millikan-White formula [9] to take into account the inaccurate estimation of the collision cross-section at high temperatures.For a mixture, this produces to where X is the species molar fraction and the interspecies relaxation time, τ m−s, V−T , is expressed as the sum of the Millikan and White (MW) and Park (P) contributions where firstly Coefficients A m,s and B m,s can either be calculated following Equations ( 12) and ( 13) or be read from a table.
M m,s is the reduced mass of species m and s, and is given by where M m is the molecular weight of species m in g/mol.The tabulated values of A m,s and B m,s used in this paper originates from the work of Park [5].Secondly, Park's counterpart may be written as where n m,s is the number density of the colliding pair (m, s), cm represents the average molecular speed, and σ v,m is the limited collision cross-section defined as follows and In the original formulation, σ m is taken as 10 −21 m 2 .In this work, σ m is set to 3.0 × 10 −21 m 2 for N 2 and O 2 , and 3.0 × 10 −22 m 2 for NO [24].
The vibrational energy pools of the different molecules are not supposed to be tightly coupled.Although the vibrational-vibrational (V-V) energy redistribution is known to play only a secondary role in vibrational energy exchange, the multi-vibrational temperature option has been retained to allow a consistent coupling between dsmcFoam and hy2Foam.The formulation that has been implemented to account for V-V energy transfer is the one of Knab et al. [25,26].The source term designated by Q m, V−V may be written as where N A is the Avogadro number, R is the universal gas constant, and P m,l represents the assessed exchange probability whose recommended value is a constant equal to 10 −2 [27].σ m,l is the collision cross-section of the pair (m, l) based on the respective molecular diameters.
For validation purposes in particular, it is useful to evaluate the temperature that the gas particles actually feel.This temperature is considered as the overall temperature, T ov , and can be determined considering the different molecular temperatures and their respective degrees of freedom.This may be written as Y is the mass-fraction and ζ is the number of degrees of freedom in relation to one specific energy mode.Both translational and rotational energy modes are supposed to be fully excited; therefore, three degrees of freedom are associated with the translational mode and ζ t,s = 3. Neglecting the degree of freedom of rotation along the nuclear axis for diatomic particles, two degrees of freedom are allocated to the rotational energy mode and ζ r,s = 2.For the vibrational and electronic modes, ζ v,s and ζ el,s are derived from the following equations and and finally for electrons, ζ e = 3.

Non-Equilibrium Navier-Stokes-Fourier Equations
The Navier-Stokes-Fourier (NSF) equations are typically solved to compute transient compressible reacting flows where the continuum hypothesis holds.Let N s be the number of species and N m the number of molecules composing the mixture.Both total and all individual species continuity equations are solved for robustness as explained in [28].In a spatially-zero flux-divergence form, the continuity equation, the N s species transport and reaction equations, the momentum equations, the N m vibrational energy equations, and the total energy equation can be written as [6] ∂U ∂t = Ẇ (22) where the vector of conserved quantities, U , is given by where u, v, and w are the three components of the velocity vector.The total pressure, p, is computed as the sum of the partial pressures as follows For ions, no additional vibrational energy equations are solved but instead, the temperature of an ion is forced to follow the temperature of the non-ionised molecule.Moreover, temperatures of atoms and electrons are set to the vibrational-electronic temperature of a reference molecule.
Finally, the source term vector, Ẇ, is written as Ẇ = 0, ωs , 0, 0, 0, ωv,m , 0 where ωs is the net mass production of species s and ωv,m , for m ∈ N m , is given by Q m, C−V being the vibrational-electronic energy added or removed by reactions to the species m.
In Equation ( 26), the last three source terms that deal with electrons are omitted and they will be introduced in future work.

Generalities
The law of mass action stipulates that the net source of chemical species s due to chemical reactions, previously denoted by ωs , is computed as the sum of the reaction sources over the N R reactions in which the species participates in.It is written as [29] where ν s,r and ν s,r are the forward and backward stoichiometric coefficients of species s in the reaction r, and N r is the number of species in reaction r.
The forward rate constant, k f , is assumed to follow the Arrhenius law, and is thus given by where A is a pre-exponential factor, β is the temperature exponent, and T a is the temperature of activation derived from the activation energy.T c, f is the controlling temperature of the forward reaction.

Chemistry-Vibration Coupling: The Park TTv Model
For a dissociation reaction, the controlling temperature is Park's temperature, defined as An exponent of 0.7 in favour of the trans-rotational temperature is commonly adopted [6].Other controlling temperatures corresponding to the various types of reaction that may occur in hypersonic flows are also implemented in hy2Foam and are shown in Table 1.Throughout this paper chemical reactions are considered to be irreversible to be in line with the dsmcFoam study.The relation giving the backward rate constant, k b , is thus omitted here but one could refer to [30] to compute its value.
The vibrational-electronic energy added or removed by chemical reactions to species m is translated into an additional source term to appear in the NSF equations as seen in Equation (25).It takes the following form [30,31] The coefficient D m can either be determined from a non-preferential or a preferential model.The non-preferential model postulates that molecular creation and destruction take place at the average vibrational energy such that [30] whereas in the preferential model molecules at the higher vibrational energy states are more likely to undergo dissociation [32] D m is the dissociation potential of species m and it can be found in Table A1.The fraction of the dissociation energy, α m , is a constant usually taken as 0.3 [6].However, heat bath DSMC studies from Holman and Boyd have shown that α m follows a linear trend as a function of the translational temperature [33].In hy2Foam, both non-preferential and preferential chemistry-vibration models are available for use and for the latter model, the choice is left between a constant and a linear fit with translational temperature for α m .

Chemistry-Vibration Coupling: The Coupled Vibration-Dissociation-Vibration Model
The coupled vibration-dissociation-vibration (CVDV) model of Marrone and Treanor [34] is another popular model to describe the preferential energy removal from the upper vibrational energy states due to dissociation [32].The forward rate constant introduced in Equation ( 28) is modified as follows where U m is a factor defined as and T F,m is a modified temperature given by In Equation (33), Z(T) is the vibrational partition function for the particle m at temperature T computed as in which α is a positive integer and N is a cut-off value for the number of vibrational levels considered for m.Using the simple harmonic oscillator model, the energy of the α-th vibrational level of the molecule m is given by Introducing now and and splitting the net source of chemical species into its forward and backward components, the CVDV chemistry-vibration source term eventually becomes where m m is the mass of one particle of species m.

Implementation in OpenFOAM 2.3.0
The version 2.3.0 [35] of the OpenFOAM CFD platform possesses one single-temperature solver that is specifically dedicated to resolve high-speed compressible flows [36].This solver named rhoCentralFoam uses the central-upwind interpolation schemes of Kurganov and Tadmor [37].It has been tested and assessed, and has been shown to produce similar results to those given by the MISTRAL flow solver [38] and the TAU code from the German Aerospace Center (DLR) [39] for the study of a non-reacting hypersonic flow past a hollow cylinder in the continuum regime [40].Basic chemistry features have first been added to rhoCentralFoam by taking advantage of another OpenFOAM solver, reactingFoam, that covers subsonic combustion.This allows the modeling of chemical reactions at low temperatures and the introduction of quantities either related to the mixture itself or to each of the individual species composing the mixture.
The newly blended single-temperature solver has been named hyFoam and validated considering the reversible reaction of hydrogen and iodine to produce hydrogen-iodide in the reaction I 2 + H 2 − − − − 2HI.This reaction takes place at a constant temperature of 700 K and at an initial pressure of 0.528 atm.The species concentration versus time can be easily derived [41] to produce the analytical solution shown in Figure 2. hyFoam results are in very satisfactory agreement thus validating the first stage of the new solver implementation.The flexibility of the C++ programming language then enables the user to incorporate all of the new two-temperature aspects outlined in the previous paragraphs.Based on the foundations of hyFoam this new solver has been given the name hy2Foam.The next section will be dedicated to the step-by-step validation of hy2Foam for a five-species air mixture.

Results and Discussion
The newly coded two-temperature solver will now be validated considering the contributions of competing mechanisms in isolation and a zero-dimensional adiabatic heat bath will be used for this purpose.From an initially disturbed system, the V-T relaxation of a single-species gas to recover equilibrium conditions will first be investigated.This will be followed by the V-T and V-V thermal relaxations of several mixtures.The relaxation towards equilibrium of a chemically-reacting mixture will finally be carried out.The verification is supported by comparisons against results from DSMC codes dsmcFoam and MONACO [42], and the CFD two-temperature solver LeMANS.Throughout this entire section, the heat bath is composed of a single cubic cell of length 1 × 10 −5 m.The time-step for CFD and DSMC computations has been set at 1 × 10 −9 s.At least 50,000 DSMC particles were used for each dsmcFoam run [43] and data shown are the ensemble average resulting from three statistically-independent computations.It should be noted that the dsmcFoam code used here has three temperatures: translational, rotational, and vibrational.However, the rotational collision number, Z r , can be tuned and set to unity to allow a direct comparison between dsmcFoam and two-temperature CFD solvers in the absence of electronic energy.

Vibrational-Translational Relaxation of a Single-Species Gas
The hy2Foam solver has already been validated by considering the thermal relaxation of a single-species gas [44].Since new features have been added such as corrected species-dependent coefficients in the interspecies relaxation time (Park 1993, Ref. [5]) and the inclusion of the electronic energy mode, it has been decided to further validate the solver.

Case without Electronic Energy
Following the work of Boyd and Josyula [45], vibrational heating (T 0 v < T 0 tr ) and vibrational cooling (T 0 v > T 0 tr ) heat bath simulations involving nitrogen have been performed.Results from hy2Foam and dsmcFoam, and published data using LeMANS and MONACO are shown in Figure 3a,b.In the heating case, the initial trans-rotational temperature is set at 10,000 K while the vibrational temperature is lowered down to 1000 K.The initial pressure is set to 1 atm.From Equations ( 5) and ( 20), the number of vibrational of degrees freedom associated with the simple harmonic model for nitrogen at T 0 v is equal to ζ 0 tot = 0.240.This corresponds to an initial overall temperature of 9588 K according to Equation (19).In Figure 3a the hy2Foam solver is run twice, first in a default configuration with the A N 2 ,N 2 coefficient taken from Park's 1993 table, and secondly in a configuration that represents a best-fit to dsmcFoam results by adjusting the value of A N 2 ,N 2 , as shown with the red solid lines.
An excellent agreement is observed between the default configuration of hy2Foam, the CFD code LeMANS and the DSMC code MONACO.dsmcFoam predicts a noticeably longer thermal relaxation that is 48% slower (defined as the time to recover 99% of the equilibrium temperature).
It can be seen that all simulations converge towards an equilibrium temperature of T eq = 7623.3K.The number of vibrational of degrees freedom for nitrogen at T eq is 1.590.Hence, the total number of degrees of freedom at equilibrium is eq tot = 6.590.A straightforward calculation that ensures energy equipartition enables the determination of the initial overall temperature The calculation produces an initial overall temperature of 9588.0K which is therefore consistent with the aforementioned value and in agreement with both DSMC codes.
In the cooling case, the initial trans-rotational temperature is set to 3000 K while the vibrational temperature is increased to 10,000 K.The hy2Foam results match exactly with the prediction given by MONACO in Figure 3b.Conversely, there are significant differences between dsmcFoam and hy2Foam.The post-collisional energy redistribution in dsmcFoam is handled by the quantum Larsen-Borgnakke procedure [18].MONACO adopts the standard Larsen-Borgnakke procedure but post-collisional energy redistribution is ensured to be consistent with the CFD approach at a macroscopic scale [45,46], which explains why the results are comparable to hy2Foam and LeMANS.

Case with Electronic Energy
In the preceding heating case, the electronic energy has no significant effect on thermal relaxation if the initial trans-rotational temperature is maintained at 10,000 K and it has thus been increased up to T 0 tr = 30,000 K.The effect of the electronic energy mode on the thermal relaxation of the nitrogen molecules is displayed in Figure 4.The equilibrium temperatures reached in both cases agree with energy equipartition and are specified on the right-hand side.
Without electronic energy, the hy2Foam and dsmcFoam results are in very satisfactory agreement considering the degree of empiricism that the two-temperature solution carries.With the inclusion of the electronic energy the relaxation to equilibrium is observed to be around 26% faster.These runs highlight the importance of the electronic energy mode that is often overlooked in re-entry case studies [20].Indeed, the difference in equilibrium temperature is not of the order of a few percent but is in the present case as large as 23.7%.Accounting for the electronic mode of nitrogen thus becomes necessary above 10,000 K to accurately predict thermal quantities for example.
The results shown in Figures 3a,b and 4 therefore validate the implementation of hy2Foam for a single-species gas.

Vibrational-Translational Relaxation of a Non-Reacting Multi-Species Gas
The non-reacting heat bath is now composed of gas mixture made of N 2 and N in equal proportions.The initial temperatures remain unchanged with regard to the previous paragraph and number densities are equal to 5.0 × 10 22 m −3 for both species.In Figure 5, the equilibrium temperatures specified on the right-hand side are once more consistent with energy equipartition and the same values are recovered using dsmcFoam and LeMANS.The disaccord between the hy2Foam and dsmcFoam solutions appear to be slightly greater that in the case considered without atomic nitrogen, nonetheless, they do remain satisfactory.The results using LeMANS are not in agreement with hy2Foam; however, a different convention is adopted for Equation (15) in LeMANS where the number density n m,s represents the mixture number density and σ m equals 1 × 10 −20 m 2 .For clarity, these modifications have been temporarily implemented in hy2Foam and are shown by the red solid line in Figure 5.The temperature profiles are now shown to be nearly superimposed, thus verifying the hy2Foam implementation for a mixture.Compared to the full N 2 configuration shown in Figure 4, the increase in equilibrium temperature due to the loss of half of the mixture vibrational energy is less pronounced with the inclusion of the electronic mode.This can be explained by the relative importance of the electronic mode of N that brings 1.39 degrees of freedom to the mixture at T eq , and thus compensates part of the vibrational energy loss.

Vibrational-Translational and Vibrational-Vibrational Relaxations
A vibrational cooling scenario with a high initial vibrational temperature is chosen to better illustrate the role of vibrational-vibrational energy exchange.The mixture is composed of oxygen and nitrogen molecules in equal proportions.The initial trans-rotational temperature is set at 5000 K, the vibrational temperature at 30,000 K, and the pressure at 1 atm.In Figure 6, the thermal relaxation process in hy2Foam is presented with and without V-V energy transfer using dashed-dotted and solid lines, respectively.As it is supposed to, V-V energy transfer promotes vibrational energy redistribution between the vibrational energy pools and thus a quicker relaxation towards equilibrium is achieved.In the present case, the N 2 and O 2 vibrational energy pools remain distinct throughout the calculation.The solver LeMANS has only one vibrational temperature.Logically, this unique vibrational temperature profile should be somewhere between the T v,N 2 (t) and T v,O 2 (t) profiles when the V-V transfer is enabled and this is precisely what is observed in Figure 6.
Figures 5 and 6 provide sufficient evidence of the verification of hy2Foam for a non-reacting multi-species gas.

Relaxation of a Chemically-Reacting Mixture
The chemical reaction considered is the irreversible molecule-molecule dissociation of nitrogen The rate coefficients are shown in Table 2 in which the units of A and T a are given in m 3 •mol −1 •s −1 and Kelvin, respectively.They come from the Quantum Kinetic (QK) theory [18], and Park's rates for the use of a two-temperature model [9].The next simulations are run from an initial trans-rotational temperature of 30,000 K and vibrational temperature of 1000 K.The initial number density is equal to 5.0 × 10 22 m −3 for both species.In the following, the hy2Foam configuration named Park shown with solid lines in Figure 7 makes the use of the Park TTv model and Park's 1993 chemical rate constants, which is the most commonly adopted set-up in the hypersonic community.The second configuration shown with dashed-dotted lines and called CVDV-QK results in the combination of the CVDV chemistry-vibration model and the QK rates.can first be the conventionally-adopted Park set-up, the decrease in trans-rotational temperature is not well-captured.The chemical reactions take about ten times longer to have a significant effect on the mixture composition when compared to the DSMC solution.Using the CVDV-QK configuration in hy2Foam, a much improved agreement is achieved with regard to the DSMC simulations.The trans-rotational temperature and species concentrations are now correctly estimated over time, while the early increase in vibrational temperature can be imputed to a departure from a Boltzmann distribution and certain aspects of current multi-temperature CFD modelling approaches used to describe N 2 −N 2 and N 2 −N interactions [47].
Both CFD and DSMC chemistry modules must deliver consistent results the development of any hybrid CFD-DSMC code.The appropriateness of the CVDV-QK configuration is thus further evaluated for a mixture initially set in a state of thermal equilibrium.Figure 8 confirms the trend observed in Figure 7.The Park configuration first initiates a thermal relaxation before any chemistry takes place which results in a lag in the vibrational temperature decrease.Conversely, the profiles predicted by dsmcFoam are once again much better approximated by the CVDV-QK configuration.It can be concluded that is a more suitable choice for use in a CFD solver than the Park configuration when compared to a DSMC code that uses Quantum-Kinetic chemistry.

Chemically-Reacting Air
To conclude on the validation of the hy2Foam solver using adiabatic heat baths, the complete set of 19 reactions to occur a pre-heated 5-species is considered.15 dissociation reactions and four are listed in Table 1 in Ref. [18].The initial pressure is set at 0.063 atm and the initial temperature is equal to 10,000 K.It is supposed in this case that all molecular temperatures are in thermal equilibrium at all times.In hy2Foam, the possibility is offered to the user to downgrade the two-temperature solver to a single-temperature solver.By doing so, there is again a very good agreement between hy2Foam and LeMANS temperature and number density fields as shown in Figure 9.

Conclusions
The open-source CFD OpenFOAM has been enhanced with a new solver to compute hypersonic reacting flows in near-equilibrium conditions.The code, called hy2Foam, introduces a trans-rotational temperature and multiple vibrational-electronic temperatures.It makes the use of state-of-the-art vibrational-translational and vibrational-vibrational energy transfers and offers the choice of whether to include the electronic energy into the calculations.Chemistry-vibration coupling is operated using either the Park TTv model or the CVDV model.The demonstration of the successful implementation using a zero-dimensional benchmark case has been carried out by considering a single-species gas, a non-reacting mixture, and a reacting mixture initially set in a state of thermal non-equilibrium.In comparison with previously heat bath data, hy2Foam has shown perform equally well.Further testing against provided by the DSMC code dsmcFoam has highlighted again the discrepancies within CFD for flows that depart significantly from a Boltzmann distribution.
The assessment of hy2Foam in an 11-species air mixture configuration will be considered for future work together with the further exploration of multi-dimensional case scenarios.Finally, the authors will aim at the development of an open-source hybrid CFD-DSMC solver within the OpenFOAM framework.As a key feature of the hybrid code, hy2Foam will be able to correctly model high-speed rarefied gas flows using state-of-the-art numerical techniques.In this respect, it has been demonstrated that the use of the CVDV model and rates derived from Quantum-Kinetic theory promote excellent agreement between the CFD and DSMC chemistry modules.

Figure 2 .
Figure 2. Species concentration versus time for a chemically reacting H 2 −I 2 reservoir.

Figure 4 .
Figure 4. Thermal relaxation of a N 2 heat bath with and without electronic energy.

Figure 5 .
Figure 5. V-T relaxation of a N 2 −N heat bath.

Figure 6 .
Figure 6.V-T and V-V relaxations a N 2 −O 2 heat bath.

Figure 7 .
Figure 7. Influence of the chemistry-vibration model and chemical rate constants on a chemically-reacting N 2 −N heat bath.(a)Temperature versus time; (b)Normalised number density versus time.

Figure 8 .
Figure 8. Chemically-reacting N 2 −N heat bath in an initial state of thermal equilibrium.(a) Temperature versus time; (b) Normalised number density versus time.

Table 1 .
Controlling temperature depending on the type of the chemical reaction.

2 .
Parameters for the evalution of k f .