Investigations of vibrational kinetics relaxation within air shock wave plasma

A vibrationally detailed kinetic model is used to study the relaxation behind shock waves in air. The role of recently published data for the rate coefficients of the Zeldovich reactions of NO formation is studied in detail. Results allow to study the radiation emitted from the shock-heated gas. Comparison with some emission spectroscopy study performed in a shock tube facility shows only qualitative agreement with the model predictions but it allows to identify directions for model improvement.


Introduction
During the re-entry of a spacecraft into a planetary atmosphere, the relative stream velocity is so high that a bow shock is created in front of the vehicle. After the strong shock wave, in the relaxing flow, a massive amount of the free stream kinetic energy is converted through inelastic collisions, into internal energy of the surrounding gases, which results in intense convective and radiative heating over the spacecraft structure. Inter-particle collisions promote energy exchanges among the internal modes of atoms and molecules, and produce chemical reactions including dissociation and ionization. Since the time scale of energy exchange processes for internal levels and the chemical time scale are comparable with the fluid dynamics one, the flow behind the shock is in thermal and chemical nonequilibrium. Understanding such nonequilibrium processes in shock-heated flows is a key issue for the design of spacecrafts. Although progress has been achieved during the last decades, modeling of thermo-chemically nonequilibrium flows is still a challenging problem in both academia and industry.
The state of the gas, i.e. the concentration of the gas species and their internal energy (rotational, vibrational and electronic) level populations can be estimated by either multitemperature (MT) models [1] or state-to-state (StS) approaches [2]. The MT models make the assumption that each energy mode is characterized by its own temperature, which is determined by solving an individual energy conservation equation. These models rely on the common assumption that the relaxation of the internal energy modes occurs through a series of Boltzmann distributions at their respective temperatures. However, in many flow conditions, the populations deviate from Boltzmann distributions [1]. Conversely, the StS approaches treat each molecule with different internal state as a separate pseudo-species, for which a continuity equation is solved. Without adopting any quasi-stationary distribution over the internal energy modes these approaches allow representing arbitrary nonequilibrium situations. The StS approaches, however, require a large number of reaction rate constants corresponding to the different physico-chemical processes. In this study, the nonequilibrium relaxation in shockheated air is simulated with a StS approach, since the vibrational state-specific rate coefficients from ab initio methods are becoming available to describe vibrational relaxation and reactive processes in N 2 -O 2 mixtures. Particularly sensitive is the population of NO whose formation is governed by a strong coupling between chemical and vibrational kinetic and consequently tightly depending on the concentration of other molecular species.
Making use of the recent quasi-classical trajectory (QCT) results [3,4,5] for the Zeldovich exchange reaction of NO formation on theoretically calculated potential energy surfaces (PES), the thermal and chemical kinetics of NO is discussed in detail.
As an assessment of the modeling, the post shock radiative signature has been calculated under the adopted assumptions and critically compared to intensity spectra measured in shock tube experiments [6] in conditions representative of Earth re-entries (u ∼ 5-8 km/s).

Modelling
In this work, we consider an air gas mixture consisting of the five chemical species N 2 , O 2 , NO, N and O. For the shock conditions of interest, the ionization degree is expected to be small and the focus is on the coupling of translational, vibrational and chemical kinetics. Thereby, ionization is neglected and excitation of vibrational and electronic levels occurs mainly through heavy -heavy particle collisional impact.
A vibrational state-specific description is adopted for the diatomic molecules in their ground states. Vibrational levels are obtained from a RKR analysis on the low energy part of spectroscopically determined potential energy curves and the extrapolation to the dissociation limit [7]. This gives a total of 61 bound levels for N 2 , 46 for O 2 and 48 for NO. The rotational mode is assumed in equilibrium at the gas temperature, T . In the following paragraphs, the physical model used for the simulation of the relaxation behind shock waves in air is presented.

Governing Equations
In the frame of reference moving with the shock wave, the flow is described by the stationary 1D Euler equations: with H = {N 2 , O 2 , NO, N, O}. ρ, u, p and h are the density, velocity, total pressure and specific enthalpy of the gas mixture. y v k andω v k are the mass fraction and source term of the molecular species k with vibrational level v. v k max is the maximum vibrational quantum number for species k and v k max = 0 for atoms. The species mass fractions and related chemical source terms are obtained by summing over all vibrational levels, i.e.
The conservation equations are completed by closure relations obtained from ideal gas assumption and then solved by the ODEPACK library [8].

Vibration-Chemical Kinetic Mechanism and the Reaction Rates
The most important part of the modelling is the kinetic mechanism that strongly affects the thermal and chemical relaxation of the flow. The kinetic mechanism drives the evolution of the composition before equilibrium is achieved. Three different categories of physico-chemical processes are included in this work. By denoting AB and CD the diatomic molecules, and X any of the 5 species, as well as v and w the specific vibrational levels, the processes include: 1) the vibrational-translational (VT) energy exchange that drives the vibrational level populations: 2) the dissociation-recombination (DR) reactions, responsible for the production of atomic species: and 3) the two Zeldovich exchange (ZE) reactions [9] which are the essential mechanism for NO formation: The generic chemical reaction i can be expressed in the general form: This allows to express the chemical source termω i k,v using the law of mass action: where, s k,i r , resp. s k,i p are the reactant, resp. product stoichiometric coefficients of species k in reaction i. k i f , resp. k i b are the forward, resp. backward reaction rate constants, which satisfy the detailed balance condition k i b = k i f /k i eq , k i eq being the equilibrium constant. M k is the molar mass of species k.
[·] denotes molar density. N s is the total number of species in mixtures. N r is the total number of chemical reactions.
Values for the state-specific rates are obtained, where available, from quasi-classical trajectory (QCT) calculations and the semi-analytical forced harmonic oscillator (FHO) model. For the remaining processes, where none of the above is available, the classical models of Landau-Teller (LT) relaxation and Park's two temperature model are used. Table 1 summarizes the state-specific reactions used in this work and the methods used to calculate the reaction rate constants.
2.2.1. State-specific rates from FHO and QCT calculations. In the present work, the statespecific rates for vibrational excitation and dissociation in diatom-diatom collisions are taken from the STELLAR Database [10]; it collects results from the work of Lino da Silva et al., who applied the FHO model [11] to the production of a multi-quantum datasets of VT and DR rates [12] in a wide temperature range. The excitation and dissociation rates for N 2 − N and O 2 − O systems are from the work of Esposito et al., based on the QCT approach modelling the evolution of the vibrational and rotational states during a collision by means of classical mechanics [13,14]. For the two Zeldovich exchange reactions, we have two alternative datasets using QCT calculations: one from Bose and Candler [15,16] and one from the very recent work by Esposito [3,4,5]. The dissociation rates in N 2 − O and O 2 − N collisions are also from this recent work.
State-specific rates from classical models. The classical models are used to define the state-specific reaction rates in a way to recover the macroscopic chemical and energy source terms. The LT model assumes vibrational energy relaxes as a consequence of VT energy transfer collisions AB − X as: where, e AB vib,eq (T ) is the specific vibrational energy in equilibrium at gas temperature T ; τ AB−X is the relaxation time, modelled from the Millikan-White's formula [17]. In a StS framework, this behaviour is recovered by choosing the state-specific rates for reactions AB(v)+X AB(w)+X as: where, N A is the Avogadro's number; is the reduced equilibrium populations. Backward rates are calculated from the detailed balance condition, with k v→w Next, for reactions as eq. (5), according to Park's model [1], the chemical reaction rates are functions of different mode temperatures, depending on the type of reactions. The controlling temperatures T f , T b in forward and backward reactions is given by where T vib is the energy-representative vibrational temperature for which a Boltzmann distribution would give the same vibrational energy as the non-equilibrium distribution. Then the macroscopic rate coefficients can be expressed by modified Arrhenius functions: where, the constants α, β, γ are adjusted to reproduce experimental determinations of K f and K eq . The total chemical and vibrational energy source terms are obtained by: where, E k f , E k b are the average vibrational energy lost/gained in reactions. In practice, it is by assuming the reaction has a non-preferential character. In the current StS approach, the above behaviour is reproduced by choosing: Values for the macroscopic rate constants K f and K b are taken from Ref. [6].

Post-shock radiative signature modelling
In the following we present the approach adopted to rebuild the radiative signature measured in shock tube experiments (see [6,18]). As in the experiments, the radiative signature of the postshock plasma is considered in the range 200-400 nm and is dominated by transitions between bound energy levels of diatomic molecules. The spectral emission coefficient of the plasma is expressed as the sum of bound-bound transitions: where n u is the population of the upper level, A ul is the Einstein coefficient for spontaneous emission, σ ul is the wavenumber of the transition, h is Planck constant, c is the light velocity and f (σ − σ ul ) is the line profile shape approximated by a Doppler profile to account for the thermal broadening (in our conditions the collisional broadening remains marginal).
The emission coefficient has been calculated accounting for the electronic systems of N 2 , O 2 and NO molecules listed in Table 2. The spectroscopic data have been provided by the high resolution spectroscopic database HTGR considering an exhaustive set of vibrational bands, taken into account up to the last vibrational levels below the dissociation limit for each electronic state (for details see [19]). Also, although in practice the Second positive system of N + 2 is a (0:4 -0:22) significant contributor above 300 nm, it is neglected in the present work which is focused on neutral molecules. The population of molecular energy levels has been calculated in the frame of a two temperature model. For a well defined energy level of a diatomic molecule (n, v, J), characterized by the electronic state n and the vibrational, resp. rotational, quantum number v, resp. J, the population can be expressed as: where Q(T ve , T rot ) is the two temperatures partition function [20], N is the total population of the considered molecule and g nvJ is the degeneracy of the level (n, v, J); E el , E vib (n, v) and E rot (n, J) are respectively the pure electronic, vibrational and rotational term and E inter (n, v, J) is the vibrational-rotational interaction term (in practice much smaller that the other terms).
The total population density of the species, N and the temperatures T ve and T rot are output of the flow simulation discussed above. The measurements dealt with in this study concern the intensity emitted by the flow in direction perpendicular to the its axis (for details on the data acquisition system and the methodology see [18]). Hence, the intensity collected by Optical Emission Spectroscopy corresponds to the local emission integrated along the flow diameter, where uniform properties are assumed in the transverse direction. Since the plasma medium is in practice optically thin, the transverse intensity I σ (x), at a position x behind the shock, is merely calculated as : where the shock diameter D is assumed to be equal to the shock-tube test section diameter 50 mm and η λ (x) is the emission at the location x (along the flow axis). Also, the acquisition time was set higher than the radiating flow lasting time, meaning that the collected radiation profile is not spatially resolved, but instead integrated along the flow axis. The intensity collected by the system is then calculated as: where ∆x = v shock × t rad is the length of the plasma axial extent, and t rad refers to the effective duration of the plasma radiation when passing through the measurement test section of 6 and 17 µs respectively for the 0.25 and 1 Torr cases. Finally, for comparisons with measured spectra, the so-called collected intensity is convoluted with the apparatus function approximated by a normalized triangle function of 2 nm full width half maximum.

Flow properties: results and discussion
We consider shock waves propagating at u 1 = 5.0 to 8.06 km/s in mixtures of N 2 and O 2 (80%-20% molar fraction composition) at T 0 = 298 K, p 0 = 0.25 or 1.0 Torr. The conditions just behind the shock are obtained assuming frozen chemistry and vibrational kinetics. The detailed State-to-State approach is used to solve for the evolution of the flow behind the shock waves. We employed two different models of state-specific rates for the Zeldovich exchange reactions, i.e. Esposito's and Bose's rates, which are labeled as 'Espo' and 'Bose' in the following figures respectively. Calculations show that the relaxation behaviors are similar in all cases. Therefore, in the following paragraphs, we firstly analyze the flow properties, vibration excitation, DR and exchange reactions, as well as non-Boltzmann vibrational distributions, and compare the Esposito's and Bose's rates based in the case with u = 8.06 km/s, p 0 = 0.25 Torr. Then we discuss the formation of NO as a function of the shock conditions.
3.1. Relaxation flow for u 1 = 8.06 km/s Figure 1 shows the profiles of the gas temperature and representative vibrational temperatures of diatomic species (Figure 1(a)), and species mole fractions (Figure 1(b)) up to 1 m behind the shock. The solid lines are the results obtained by employing Esposito's ZE rates, and dashed lines are those by Bose's ZE rates. Differences are only noticeable for the vibrational temperature and abundance of NO. At the shock, the gas temperature increases rapidly and reaches about 31500 K. After the shock, the vibrational temperatures of N 2 and O 2 increase up to a maximum due to VT transitions. Then, once sufficient dissociation occurs, the loss of vibrational energy by dissociation takes over the energy gain due to excitation, the two temperatures fall and gradually equilibrate with the gas temperature. The two models for the ZE reactions result in large differences in the predicted evolution of the NO vibrational temperature. The composition profiles show the dissociation of N 2 and O 2 and the formation of NO, N and O. At the early phase of the relaxation, the low vibrational temperatures of N 2 and O 2 prevent the dissociation and exchange processes to occur. Therefore, their mole fractions remain at the initial levels before 1 × 10 −5 m. Then, due to VT excitation, the molecules gradually climb the vibrational ladder, which initializes the dissociation of upper levels. As the vibrational temperatures increase, the dissociation progressively goes deep into the lower energy levels. As a result, the molecules N 2 and O 2 are dramatically reduced, while the atoms N and O are abundantly produced. At about 1 × 10 −1 m, the mole fraction of N 2 is reduced from 80% to 10%. Dissociation of O 2 is much faster so that already at 1 × 10 −3 m, the mole fraction of O 2 falls below 0.1%. NO, instead, is produced in the Zeldovich exchange reactions. The onset of these processes requires enough atomic nitrogen and oxygen from dissociation of N 2 and O 2 . Therefore, a long incubation distance is observed for the appearance of NO. Esposito's rates predict a relatively larger formation than Bose's rates. At about 2 × 10 −5 m, Esposito's rates predict more than 0.1% mole fraction of NO, which reaches its maximum value of 0.65% at 8 × 10 −5 m. Bose's rates obtained more than 0.1% at 4 × 10 −5 m, and then get the maximum value of 0.24% at 1.7 × 10 −4 m. After the peak, dissociation prevails over formation, and the mole faction of NO begins to curb. At 1.0 × 10 −3 m, the mole fraction of NO falls back below 0.1%. Moreover, for the reaction O 2 (v) + N NO(w) + O, once NO is formed in sufficient amount, the reverse reaction takes over and leads to a production of O 2 . This slightly increases the vibration temperature of O 2 and slows down the consumption of O 2 by dissociation, as observed between 3 × 10 −4 m to 5 × 10 −4 in Figure 1.
For the purpose of comparison, we also apply the MT method to simulate the relaxation flow; in this case, vibrational excitation is modeled by the Landau-Teller relation with Millikan-White's formula for the relaxation time, and the vibrational-chemical coupling is from Park's two-temperature model. The comparisons of resulting temperature and concentration profiles are illustrated in Figure 2. The MT method predicts different behaviors. All the vibrational  temperatures for diatomic molecules start from low temperatures. Then they increase up to the maximum values and undergo overshoot to the gas temperature during their decrease when dissociation is predominant. Dissociation of N 2 and O 2 , as well as formation of NO are slower in the MT simulation, and longer incubation distances are observed. Finally, the MT model also predicts less amount of N O.
It is very interesting to plot the vibrational distribution functions (VDF) for diatomic species in the flow. Figure 3 shows the locations at which VDF for N 2 are plotted; x 0.79 , x 0.7 , x 0.4 and x 0.1 denote the locations at which the N 2 mole fraction is reduced to 79%, 70%, 40% and 10%, respectively. The calculated VDF together with the Boltzmann distribution functions at the equivalent vibrational temperatures are illustrated in Figure 4. At x 0.79 , a location close to the shock, the VDF highly departs from the Boltzmann distribution with large overpopulation of the high-lying vibrational levels. Then, as excitation and dissociation proceeds, the VDF gradually tends towards a Boltzmann distribution at x 0.7 . However, at x 0.4 when half of the N 2 dissociation has occurred, the VDF departs from the Boltzmann distribution again due to the depletion of the high-lying levels caused by dissociation. Finally, when equilibrium has been reached at x 0.1 , the VDF is close to the Boltzmann distribution. Figures 5 and 6 show the locations and VDF for O 2 ; x 0.19 , x 0.1 , and x 0.001 are the locations at which the O 2 molar fraction is 19%, 10% and 1%. We also plot the VDF at x max O where the O has reached its maximum. Similar behaviors are observed here as the ones of N 2 . Interestingly, differences in the high-lying level populations predicted by Esposito's and Bose's rates are observed at x max O and x 0.001 . The Esposito's rates predict larger underpopulation of the high-lying vibrational levels. The discrepancy appears after the exchange reaction O 2 (v) + N NO(w) + O starts proceeding predominatly in the backward direction. Figure 7 and 8 show the locations and VDF for NO, and x F 0.001 is where the NO mole fraction has reached 0.1% during formation; x max is where the NO mole fraction has reached its maximum; x D 0.001 is where the NO mole fraction has reduced to 0.1% during dissociation. It shows that both the Esposito's and Bose's results are close to the Boltzmann distribution functions at all three locations. However, VDF predicted by Esposito's rates shows underpopulation of the high-lying vibrational levels.
A detailed investigation shows that Bose's rates to form high-lying levels of NO are significantly larger than Esposito's. This explains the vibrational temperature overshoot predicted with the former model. On the other hand, Esposito's rates are larger for the formation of NO in intermediate-lying states. This leads to a larger production of NO. Note, however, that the original work [15,16], only reported calculated rates for low and intermediate vibrational levels. Data for high-lying vibrational levels are obtained by extrapolation [7].

Formation of NO as a function of shock speed
We calculate shock waves propagating at velocities from 5.0 km/s to 8.06 km/s. All the upstream gas mixtures contain 80% N 2 and 20% O 2 in mole fraction at T 1 = 298 K. The upstream pressure and shock speed conditions considered are summarized in Table 3. Figure 9 illustrates the maximum NO mole fraction and the its position behind the shock waves. The position has been normalized by the upstream mean free path λ 1 , which is estimated as 4.71 × 10 −5 m and 1.88 × 10 −4 m for 1 Torr and 0.25 Torr respectively. Esposito's rates produce the largest    amounts of NO, while MT method has the smallest amounts. Besides, the predicted position of the maximum mole fraction by Esposito's rates is closest to the front of the shock wave, while the one by MT method is located farthest from the shock wave. Generally speaking, as the velocity of the shock wave increases, the formation of NO is reduced and the position of the maximum mole fraction shifts towards the shock front. Both effects can be explained by the strong temperature dependence of the dissociation reactions: faster shocks will result in larger temperatures, so the initial dissociation of N 2 and O 2 will occur faster, which in turn triggers NO formation closer to the shock front. At the same time, larger temperatures will promote NO dissociation, so that the maximum mole fraction decreases.

Assessment with shock tube experiments
The emission in the UV range, 200-400 nm, has been calculated at a dozen positions behind the shock using the three models investigated here referred to as MT, ES and BO. The resulting emission profiles for 5.56 and 6.94 km/s are plotted in Figure 10     the differences between MT and StS calculations become important, reflecting differences in the predicted flow properties and the emission intensity is dominated by the thermal equilibrium region. However, whatever the model, the relative contribution of N 2 FN has increased to the level of O 2 SR in the emission peak region. In all cases, the contribution of NO emission remains much smaller than the rest and has consequently been reported on a distinct plot for the sake of readability (see Figure 10, right plots), detailing the contribution of the NO systems listed in Table 2. With increasing velocity the contribution of NO radiation slightly increases, however more in the case of MT than in StS models.
Comparisons between rebuilt spectral intensities in the range 200-300 nm and measured spectra are gathered in Figure 11 for the 1 Torr case. Although the experimental and simulated Below 250 nm, the discrepancies between spectra increases with increasing velocity, whereas MT calculations tends toward a better agreement for the bands ∆v = 1 of N 2 FN around 315 nm. For all the cases in Table 3, the total emission has been calculated in the range 200-300 nm and compared in Figure 11 to the experimental data to focus on the contributions to post shock emission due to NO systems and O 2 SR. As already appreciated from the spectral analysis in the 1 torr pressure case, all calculations under-estimate the measured total emission and the discrepancy increases with velocity. MT calculations seem to reproduce the main trend of the measurements, whereas StS calculations exhibit an opposite evolution. This is connected to the decay of the O 2 SR band which contributes up to 90% of the emission in the examined spectral range.

Conclusions
A detailed kinetic mechanism has been studied for describing the relaxation behind shock waves in air. The mechanism is vibrationally specific for the electronic ground state of the neutral molecules involved. In particular, two sets of reaction rate coefficients for the Zeldovich reactions of NO formation, obtained from accurate QCT calculations, have been critically analyzed and the results compared to a standard MT formulation. Results show that the NO formation kinetics is sensitive to the details of the adopted rate coefficients. Results from the calculations have been then used to estimate the radiation emitted assuming the electronic excited molecular states are populated according to a Boltzmann distribution at the N 2 vibrational temperature. This allows some comparison with recent optical emission spectroscopy measurements performed in a shock tube facility documenting the spectral intensity of the overall shock (without any spatial resolution). A poor agreement is found in all cases, and the contribution of O 2 SR in the range 200-300 nm was found to be dominant, calling for additional studies. However the lack of information on the spatial distribution of the intensity behind the shock makes it difficult to assess specific routes of improvement. Also, indications on the merits of different kinetic models of NO formation are shown to be masked by the predominant contribution of O 2 emission in the same spectral range. Calculations with a 11 species air model, that take into account the ionization kinetics and the energy balance of the free electron component have also been carried out (not reported here), but results are 14

APhM2016
IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 815 (2017) 012026 doi:10.1088/1742-6596/815/1/012026 seen not to differ appreciably. This is to be expected since, at these relatively low shock speeds, the free electron concentration is expected to be small and excitation processes are driven mainly by heavy particle collisions. A better picture is expected from a detailed kinetic mechanism describing the population of the electronically excited states responsible for the emitted radiation. Future work will therefore be focused on the development of such a Collisional-Radiative model.