Construction of a coarse-grain quasi-classical trajectory method. II. Comparison against the direct molecular simulation method

An accurate description of non-equilibrium phenomena in high speed flows is crucial to the accurate prediction of heat loads on vehicles traveling in this regime. The extreme temperature gradients produced by these vehicles as they move through a planetary atmosphere result in strong thermal and chemical non-equilibrium. Previous work to model this type of flowfield has relied on the assumption of an equilibrium or Maxwell-Boltzmann distribution of internal energy modes (e.g., rotational, vibrational, electronic). Moreover, these models have to be calibrated to match experimental data.^1–5 However, in strong non-equilibrium, these assumptions are known to break down.^6,7 This motivated the application of state-to-state (StS) models to this type of flow which treat internal states as pseudo-species and solve directly for their population using the microscopic kinetic data.⁷ Despite their physical accuracy, the computational cost associated with this type of model rapidly becomes infeasible, limiting their application to zero- or one-dimensional simulations.^8–25 A nitrogen molecule in the ground electronic state contains nearly 10 000 rovibrational energy levels. Therefore, modeling all the possible interactions between two colliding molecules would require the kinetic data for on the order of 10¹⁵ transitions, an impossible and unnecessary task.

In order to overcome the computational cost associated with StS models, the coarse grain method (CGM) was developed. In this approach, energy states are grouped together, and instead of solving for the population of each state, properties of the groups, such as population and energy, are solved for. The actual distribution of states within the group can be reconstructed using a function derived using the maximum entropy principle. This approach has been applied to study both electronic atomic states and rovibrational molecular states.^26–36 The coarse-grain quasi-classical trajectory (CG-QCT) method described in Paper I extended the CGM to study the interaction between two diatomic molecules.³⁷ The CG-QCT method couples the construction of the reduced order model with the determination of the grouped kinetic properties, using the quasi-classical trajectory method. The CG-QCT method is designed to bridge the gap between computational quantum mechanics and computational fluid dynamics (CFD) by incorporating information about the microscopic kinetics into a computationally tractable model. However, in Paper I,³⁷ it was found that the manner in which states are grouped together has a profound impact on the resulting predictions (i.e., composition, temperature, …). Moreover, because the microscopic StS kinetics are not available for the $N2(XΣg+1)−N2(XΣg+1)$ system due to its intractable size, the model cannot be evaluated by comparing to the StS solution, as was done for the $N2(XΣg+1)−N(Su4)$ system.³³ 

The direct molecular simulation (DMS) method is a technique for directly determining the transient non-equilibrium behavior of a gas.³⁸ The DMS method is similar to the direct simulation Monte-Carlo (DSMC) approach for modeling flows, but instead of relying on pre-computed cross-sectional data, trajectory calculations are performed within the simulation, making use of the potential energy surface (PES).^38–44 As a result, the only model input to the DMS method is the PES, making the method extremely accurate due to the limited number of assumptions. Therefore, the DMS method can be used as a benchmark solution on which to evaluate other models. The DMS method has been used previously to study the N₂–N₂ PES by Jaffe et al.⁴⁵ 

In this work, CG-QCT method is applied to study the non-equilibrium state of a nitrogen mixture considering the interaction between two nitrogen molecules in the ground electronic state, the $N2(XΣg+1)−N2(XΣg+1)$ system. For the CG-QCT method, two grouping strategies will be considered: vibrational- and energy-based grouping. Vibrational-based grouping lumps states with the same vibrational quantum number together, and the rotational distribution is then reconstructed under the assumption of equilibrium between rotation and translation (T_rot = T). Energy-based grouping lumps states near in energy together, regardless of quantum structure. Using each grouping strategy, the CG-QCT model is compared with the DMS method to study the non-equilibrium energy transfer and dissociation process in a zero-dimensional reactor simulation. The objective of this work is to determine which grouping strategy better captures the non-equilibrium phenomena for $N2(XΣg+1)−N2(XΣg+1)$ interactions. The paper is structured as follows: first, a description of the physical problem to be studied will be presented in Sec. II, followed by a brief description of the DMS and CG-QCT methods in Sec. III, followed by the results of the comparison in Sec. IV, and followed by an analysis of the dissociation and energy transfer processes in Sec. V.

In this work, two approaches to studying non-equilibrium chemistry are applied to a simple isothermal, isochoric reactor simulation. This benchmark problem considers a constant volume box containing initially only nitrogen molecules at a cold temperature, 2000 K. Instantaneously the translational temperature of the box is raised to between 10 000 K and 25 000 K and held fixed for the duration of the simulation, with a Maxwellian velocity distribution. The molecules then become rovibrationally excited and eventually dissociate into nitrogen atoms in some finite time. This process is representative of the thermochemical non-equilibrium behind a strong shock wave. Accurately predicting both the non-equilibrium excitation and dissociation behavior of the molecules requires an accurate description of the non-equilibrium kinetics. Two methods for studying this are presented in this work: the Coarse-Grain Quasi-Classical Trajectory (CG-QCT) model and Direct Molecular Simulation (DMS). Both methods make use of the PES provided by NASA Ames.^46,47

In Secs. III A and III B, the CG-QCT and DMS methods are briefly discussed. Further details on the CG-QCT method can be found in Paper I.³⁷ A more extensive description of the DMS method can be found in Refs. 43 and 44.

The coarse grain method has been detailed previously in Paper I.³⁷ The coarse grain model is constructed by grouping energy states together based on some property. In this work, two grouping strategies are considered: energy-based grouping and vibrational-based grouping. Energy-based grouping groups states together according to their internal energy and whether the states are bound or quasi-bound (or pre-dissociated). The energy width of the energy bins is uniform for the bound groups and for the quasi-bound groups. Vibrational-based grouping groups states according to their vibrational quantum numbers: all rotational states which share a vibrational state are lumped together, regardless of energy. In this approach, groups of states can span the entire energy spectrum. Within each group, the distribution is assumed to be in equilibrium with the translational mode. It is important to note that the vibrational grouping model and a conventional vibrational state-to-state model assuming equilibrium between rotation and translation are identical.

The macroscopic governing equations for the coarse grain model are obtained by taking moments of the microscopic master equations for states within a group. In this work, the governing equations are the conservation of mass equation for the groups,

where the set $I$ denotes the full set of states or groups, (n_p, n_q, n_r, n_s) denote the number density of groups (p, q, r, s), n_N denotes the number density of atomic nitrogen, $Kpq−rsE$ and $Krs−pqE$ denote the excitation and de-excitation rate coefficients, respectively, $Kpq−rED$ and $Kr−pqER$ denote, respectively, the combined excitation-dissociation and combined excitation-recombination rate coefficients, and $KpqD$ and $KpqR$ denote, respectively, the dissociation and recombination rate coefficients.

The DMS method also known as classical trajectory calculation DSMC (CTC-DSMC) was first introduced by Koura^39–41 and was implemented with modern DSMC algorithms by Norman et al.⁴² Recently, the DMS method has been applied to nitrogen dissociation using ab initio PESs.^43,44 Similar to the DSMC method,^48,49 the DMS method simulates molecular collisions in a time accurate flowfield. However, instead of using stochastic collision models like DSMC, the DMS method simulates particle collisions by performing trajectory calculations using a PES. Therefore, in a DMS simulation, the PES is the only model input. Following this approach, particle trajectories are carried out at every time step and the post-trajectory state of the particle becomes its initial state for the next interaction. Hence, the DMS method constantly updates the state and composition of the simulated gas and can capture transient behaviors and non-equilibrium physics without any decoupling of rotational and vibrational energy.^38,43,44

Just like DSMC, only a fraction of the gas particles are simulated in DMS. Hence, every simulated particle represents a larger population of near-identical physical particles in the simulation volume. The ratio of actual particles in the simulation volume to the simulated particles is called particle weight (W_p). The simulations are carried out with a time step of the order of the mean collision time (τ_c) and cell volumes are of the order of the mean collision path (λ_c). For the zero-dimensional isothermal and isochoric cases presented in this paper, the particle weight was set to unity (W_p = 1), and the DMS time step was chosen to be a hundredth of the mean collision time (Δt_DMS = τ_c/100). The isothermal and isochoric case at T = 25 000 K was carried out with 10⁶ nitrogen molecules at t = 0, and for the case at T = 10 000 K, there were 6 × 10⁴ nitrogen molecules at t = 0 in the simulation volume. Since the simulation was carried out in zero-dimensions, particle movement is not needed and the simulation volume was chosen such that for the number of particles present, the particle weight would be unity and the density ρ = 1.28 kg/m³. A detailed description of the DMS method applied to zero-dimensional problems can be found in Refs. 38, 43, and 44.

In Secs. IV A and IV B, the comparison between the macroscopic and microscopic properties predicted by the DMS and CG-QCT methods will be shown. Results will be shown for two temperatures: 10 000 K and 25 000 K. The density is 1.28 kg/m³, and the molecules are initially in equilibrium at 2000 K, corresponding to a pressure of 760.137 kPa. The CG-QCT model has been applied with both grouping strategies detailed earlier: 60 energy-based groups and 61 vibrational-based groups. In both approaches, the groups are assumed to be in equilibrium with the translational temperature. First, the properties and analysis relevant to the dissociation process will be presented, followed by the properties and analysis of the energy transfer process.

Proper characterization of the dissociation process is crucial for CFD calculations because this highly endothermic process absorbs a significant amount of energy from the flowfield. In addition, the atomic species produced are highly reactive and can degrade the thermal protection systems of hypersonic vehicles. Figure 1 shows the mole fraction of atomic nitrogen as a function of time at 10 000 K. The DMS prediction is well matched by the CG-QCT energy-based grouping model. The vibrational specific CG-QCT model predicts significantly slower dissociation which stems from the inaccurate treatment of high energy states in this model: the vibrational grouping enforces mode separation across the entire energy spectrum. However, this effect is expected to breakdown for the high energy states due to the rovibrational coupling of high energy states. Moreover, because all rotational states within a vibrational group are in equilibrium at a common rotational temperature (T_rot = T), the quasi-bound states are lumped with bound states, despite the distinct differences in the kinetics that characterize the dissociation process from these internal states (e.g., the state (⁠$v$⁠, J) = (0, 0) is very unlikely to dissociate, while (⁠$v$⁠, J) = (0, 273) is extremely likely to dissociate). The discrepancy between the vibrational-based grouping and the DMS data highlights the importance of considering a rotational state for predicting dissociation. This effect is amplified at 10 000 K because molecules are more likely to climb to quasi-bound states before dissociating at lower temperatures.⁵⁰ In addition, because the DMS method cannot simulate recombination reactions, it can reach the quasi-steady-state (QSS) distribution, but not the equilibrium composition. The energy-based CG-QCT method is included without recombination reactions. Until 10⁻⁶ s, the effect of recombination is negligible; however, after this, the CG-QCT results demonstrate the effect of recombination in forcing the system to equilibrium.

The distribution function predicted by both methods can be compared both during the energy transfer and dissociation process. It was found that at the conditions studied, the dissociation process occurs when the distribution is in QSS, meaning that the relative distribution of states (or groups) is not changing in time. The resulting QSS distribution predicted by the energy-based CG-QCT method is shown in Fig. 2(a) at 10 000 K. The DMS data are grouped using the same grouping strategy for comparison with the energy binning strategy. The QSS distribution of groups predicted by the energy-based CG-QCT method is in excellent agreement with that predicted by the DMS method. At the dissociation energy, 9.75 eV, the distribution of groups turns down, indicating that the quasi-bound states are significantly depleted in QSS. At 10 000 K using the energy-based CG-QCT method, only 0.006% of the molecules are in the quasi-bound states [blue squares in Fig. 2(a)], while at 25 000 K, the quasi-bound states are over 100 times more populated, with 0.7% of the molecules in quasi-bound states [blue squares in Fig. 3(b)]. Figure 2(b) shows the vibrational distribution predicted by the CG-QCT method. Despite the differences in the dissociation rate observed by the mole fraction of atomic nitrogen (Fig. 1), the vibrational energy distribution deviates only slightly from the DMS data. This suggests that the overall dissociation rate is quite sensitive to the energy state populations and that the rotational energy distribution (not shown explicitly in Fig. 2) may play a role.

In order to compare the two grouping strategies directly, the state specific distribution was reconstructed with the vibrational-based CG-QCT method and regrouped according to the energy bins. This maps the vibrational-based CG-QCT data to the energy-based CG-QCT data for comparison and is shown in Fig. 3 at both temperatures. In both cases, the quasi-bound groups are overpopulated in the vibrational binned CG-QCT method compared to the DMS data. This stems from the assumption of equilibrium of all rotational states within a vibrational state made in this model. Despite the agreement between the vibrational distributions, capturing the rovibrational distribution is necessary for an accurate prediction of the dissociation process.

The QSS dissociation rate predicted from the CG-QCT and DMS methods is shown in Fig. 4 at various temperatures. At 10 000 K, the QSS dissociation rate predicted by the vibrational specific CG-QCT model is significantly slower than that predicted by both the energy-based CG-QCT and DMS methods. This is due to the improper lumping of quasi-bound states with bound states in this approach, which hinders dissociation, particularly at lower temperatures. Across the entire range of temperatures, the energy binned CG-QCT method is in excellent agreement with the DMS data. By contrast, at 25 000 K, the vibrational specific model overshoots the dissociation rate predicted by the DMS data. Due to the excellent agreement across the entire range of temperatures, the energy bins can accurately capture the dissociation process in QSS regardless of temperature. Table I contains the coefficients for the modified Arrhenius fits to the data computed by the CG-QCT and DMS methods. The Arrhenius form is1t given by

where A, η, and E_A are the coefficients given in Table I, and the units of $KQSSD$ are cm³/s.

In order to understand the relative importance of energy states for dissociation, the distribution of dissociating molecules in QSS was computed. For the DMS method, this is computed by counting the molecules which dissociate at each time step and taking an average over several time steps in the QSS region to obtain the distribution of molecules dissociating from various energies. For the CG-QCT model, this was done by weighting the dissociation rate from each group with the distribution of groups in QSS. This comparison is shown in Fig. 5 for the energy-based CG-QCT and in Fig. 6 for the vibrational-based CG-QCT. In the energy-based CG-QCT method, the trend matches with the DMS results, showing that molecules climb to high energy states before dissociating, with the highest probability of dissociation occurring at the dissociation energy. However, the actual peak and width of the distribution are not well matched. In Fig. 6, the vibrational CG-QCT method distribution of dissociating molecules matches well at 10 000 K despite the overall dissociation rate being significantly slower than that predicted by the DMS method. At higher temperatures, the vibrational distribution of dissociating molecules deviates significantly from the DMS results, with significant dissociation occurring from the low vibrational states predicted by the CG-QCT method. The high proportion of dissociation from low vibrational states in the CG-QCT method is most likely caused by the increased weight of the high rotational states at higher temperatures. Therefore, the states which are highly probable to dissociate now hold higher weight from the Boltzmann factor, contributing to an increased total dissociation rate from the low vibrational states.

To quantify the contribution to dissociation from each mode, the fraction of energy lost from each group from the rotational and vibrational modes are shown in Fig. 7 at T = 10 000 K. The qualitative agreement between the CG-QCT and DMS results is good, with both approaches losing a significant amount of energy from the rotational mode. The DMS results predict that 34% of the energy for dissociation comes from the rotational mode, while CG-QCT predicts 40%. The DMS results show a distinct bump in the rotational energy contribution at the dissociation energy. This occurs because quasi-bound states are primarily low-$v$/high-J states with a significant amount of rotational energy. Therefore, the contribution of the rotational energy to dissociation from the quasi-bound states overshoots the vibrational contribution for high energy states.

The energy transfer process serves to excite the internal states of the molecule until they have enough energy for dissociation to take over. Therefore, understanding the time scale in which energy transfer occurs is necessary to accurately predict the onset of dissociation. Moreover, the $N2(XΣg+1)−N2(XΣg+1)$ energy transfer process is very important for air chemistry because in many situations of interest there are initially very few atoms present. The internal energy at both temperatures predicted by the DMS method, and energy-based and vibrational-based CG-QCT models is shown in Fig. 8. At 10 000 K, the internal energy relaxation predicted by the energy-based CG-QCT model is significantly faster than both the DMS method and vibrational specific CG-QCT model. The vibrational specific CG-QCT method starts with higher internal energy because the rotational mode is already excited; however, the energy of the molecules does not increase until 10⁻⁹ s, indicating that the molecules are not gaining energy from the translational mode until that point. When the vibrational excitation starts, corresponding to the second increase in internal energy in the DMS method, the vibrational CG-QCT method matches well with the DMS data. Similarly, at 25 000 K, the internal energy relaxation predicted by the energy-based CG-QCT model is significantly faster than the DMS method. Again, the vibrational-based CG-QCT model starts with higher internal energy. At a higher temperature, the distinction between rotational and vibrational excitation is not present, indicating that the two excitation processes are more closely coupled at this temperature. In both cases, the internal energy at the final time, corresponding to the QSS energy, differs by less than 5% due to differences in the rotational and vibrational energy in the QSS region.

In addition to the total internal energy, the rotational and vibrational temperatures predicted from the DMS and vibrational-based CG-QCT model were computed. For the vibrational-based CG-QCT model, the rotational temperature is assumed to be frozen at the translational temperature. The temperatures for both cases are shown in Fig. 9. The vibrational specific CG-QCT model matches very well with the vibrational temperature predicted by the DMS method. However, this is helped by the assumption of equilibrium between the rotational and translational temperatures. If this assumption were relaxed, the vibrational relaxation time predicted by the vibrational specific model would most likely become much slower because the rotational mode would need to become excited first. Moreover, the assumption of equilibrium between rotation and translation is seen to breakdown particularly at higher temperatures, where the rotational temperature predicted by the DMS method only reaches 21 200 K.

Figure 10(a) shows the energy-based CG-QCT group distribution during the relaxation process at 25 000 K. The group distribution is significantly different from the DMS data as it appears to be relaxing faster than the DMS data, as observed in the internal energy shown in Fig. 8(b). In particular, the low energy groups which contain most of the internal energy show significantly different behavior between the DMS and energy-based binning data. A similar comparison for the vibrational binned data is shown in Fig. 10(b), showing the vibrational distribution during the relaxation process at 25 000 K. In this case, the low energy vibrational states are in excellent agreement between the CG-QCT method and the DMS method. Although the agreement for the higher energy levels is not as good, the trend predicted by the two models is similar, showing a bimodal distribution. Moreover, statistical noise is present in the high energy states from the DMS method. The agreement between the DMS method and the vibrational-based CG-QCT method in predicting energy transfer is a result of the mode separation known to be present for low energy states. For low energy states (⁠$v$ = 0, …, 3), which generally dictate the internal energy of the molecules, the separation of rotational and vibrational energy prevails, resulting in a strand structure in these states.

The vibrational relaxation time predicted by both the DMS method and vibrational specific CG-QCT model is shown in Fig. 11. At low temperatures, as seen in Fig. 9(a), the vibrational relaxation time is well matched between the DMS method and the vibrational specific CG-QCT model. At higher temperature, the vibrational relaxation time predicted by the vibrational specific CG-QCT model is approximately 50% slower than the DMS method.

This work presents an analysis of two grouping strategies used in the CG-QCT model compared to the DMS method for analyzing energy transfer and dissociation in an isothermal isochoric reactor. At the conditions studied, the energy transfer and dissociation processes are decoupled, with the gas relaxing to the QSS distribution and then dissociating from this state. Accurately describing these two different processes requires very different considerations. For the energy transfer process, the low energy states are very important, as most of the internal energy is contained in approximately the first 0-2 eV. However, for dissociation, the high energy states (E_i ≳ 6 eV) are crucial because they are most likely to dissociate and contribute significantly to the global dissociation rate.

The vibrational specific CG-QCT model predicted the energy transfer process across the range of temperatures. This is due to the fact that it captures well the vibrational strand structure caused by mode separation previously observed when studying the $N2(XΣg+1)−N(Su4)$ system.⁵³ For low energy states, the states which share vibrational quantum numbers tend to equilibrium with each other; however, these strands are slower to equilibrate with each other. Initially, the energy transfer proceeds through a series of vibrational-translational (VT) energy exchange reactions, which are significantly more efficient for vibrationally excited states. Therefore, the population of the high vibrational states is pumped up, creating a bi-modal distribution observed by Sharma et al.⁵⁴ It was observed in Paper I that initially the rate of excitation from the first few vibrational states is very slow, resulting in a significantly slower relaxation process than that predicted by the energy-based CG-QCT method.³⁷ By contrast, the energy-based CG-QCT model excitation reactions neglect all this information and excitation proceeds through a series of small energy jumps which occur very quickly. This is due to the lumping of different vibrational states together: the vibrational strand structure is lost, and all the states are assumed to equilibrate to some average temperature.

By contrast, the dissociation process in QSS is well captured across the range of temperatures by the energy-based CG-QCT model. The dissociation process depends heavily on the accurate prediction of the high energy (including quasi-bound) states because molecules tend to climb to high energy states before dissociating. Therefore, grouping the high energy states considering only the energy (not the quantum configuration) results in accurate predictions of the dissociation process. By contrast, the vibrational specific CG-QCT model lumps states together across a large range of energies. As a result, the high energy states are assumed to be in equilibrium with low energy states. Because the dissociation behavior from the low and high energy states is very different, this averaging results in significant underestimation of the dissociation from the low-$v$/high-J states. At high temperatures, because of the increased weight from these states, the Boltzmann factor seems to take over, artificially enhancing the dissociation rate from the low-$v$ states. Therefore, the vibrational grouped CG-QCT model cannot accurately account for dissociation as the energy grouped CG-QCT model can.

This work presents a detailed comparison between two methods for describing the non-equilibrium energy transfer and dissociation in an isothermal and isochoric reactor considering $N2(XΣg+1)−N2(XΣg+1)$ reactions: the CG-QCT model and the DMS method. Due to the minimal number of assumptions made in the DMS method, it was used as a benchmark solution on which to compare the two grouping strategies in the CG-QCT method. The two grouping strategies, energy-based grouping and vibrational-based grouping, were used in the CG-QCT method to study the isothermal and isochoric reactor. The comparison was made at 10 000 K and 25 000 K using the same potential energy surface for a fair comparison.

It was found that the dissociation process was well captured by the energy based CG-QCT method. The molecules tend to climb to high energy states before dissociating, and the high energy states are properly characterized by the energy-based groups. By contrast, the vibrational-based grouping failed to reproduce the dissociation process at low temperatures due to the assumption of equilibrium of rotational states, a strong assumption because of the wide range of energies spanned by the rotational states. This work highlights the importance of rotational states for dissociation processes, with 34%-40% of the energy to dissociate a molecule coming from the rotational mode confirmed by the DMS results. Moreover, the dissociation process was found to proceed under QSS conditions, and the QSS distribution predicted by the energy-based CG-QCT method was in excellent agreement with the DMS data.

By contrast, the internal energy transfer process was well captured by the vibrational-based CG-QCT method. Although the vibrational CG-QCT method assumes equilibrium between the rotational and translational modes, the vibrational relaxation was accurately predicted at both temperatures. Moreover, during the relaxation process, the vibrational distribution of low energy states, found from DMS, was well matched by the vibrational CG-QCT method. However, the relaxation time predicted by the energy-based CG-QCT method was significantly faster than what was observed in the DMS calculations. This disagreement stems from the importance of mode separation for the low energy states during the relaxation process.

Due to the highly endothermic nature of dissociation processes, they carry increased importance for CFD predictions. Therefore, the authors recommend prioritizing this process over excitation and using an energy-based binning strategy for the prediction of non-equilibrium energy-transfer and dissociation. However, in order to remedy the shortcomings of both binning models, future work will focus on the development of a hybrid grouping model in which low energy states are lumped according to vibrational quantum numbers and high energy states are lumped according to energy. This approach was previously proposed by Munafò et al. and applied for the $N2(XΣg+1)−N(Su4)$ system.³³ Additionally, computing the energy transfer coefficients for the CG-QCT model will enable a dramatic reduction in the number of groups required to be possibly only 3-5 groups.³³ 

The authors would like to thank David Schwenke and Richard Jaffe from NASA Ames Research Center for their helpful discussions. Dr. M. Panesi was supported by the Air Force Office of Scientific Research Young Investigators Program No. FA9550-15-1-0132 with Program Officer Ivett Leyva. Ms. R. Macdonald was supported by the Department of Defense (DoD) through the National Defense Science and Engineering Graduate (NDSEG) Fellowship. Dr. T. Schwartzentruber is supported by Air Force Office of Scientific Research (AFOSR) under Grant No. FA9550-16-1-0161. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the AFOSR or the U.S. government. Mr. M. Grover was supported by the Doctoral Dissertation Fellowship at the University of Minnesota.