A vibrational state-specific model for dissociation and recombination reactions within the direct simulation Monte Carlo method is introduced to study the energy level dynamics of the O2 + O system. The state-resolved cross sections for vibrational relaxation and dissociation reactions are obtained from a rotationally averaged quasi-classical trajectory database based on the Varandas and Pais O2(  3Σg)+O(  3P) potential energy surface. A two-step binary collision framework is outlined to characterize the vibrational state-resolved recombination probabilities, which are constrained by detailed balance for orbiting pair formation, and microscopic reversibility applied to the dissociation cross sections for orbiting pair stabilization. The vibrational state-to-state (STS) model is compared to the phenomenological total collision energy (TCE) and quantum kinetic (QK) models through a series of 0-d non-equilibrium relaxation calculations. A quasi-steady state (QSS) region is established in the vibrational temperature profiles of the TCE, QK, and STS models under non-equilibrium heating. This QSS region is a result of the competition between vibrational relaxation by vibrational-translational (VT) transitions and O2 dissociation. The duration of QSS predicted by the STS model is approximately ten and four times that of the TCE and QK model predictions, respectively, and the total time to reach equilibrium is approximately 3.5 times that of the TCE model and 1.5 times that of the QK model. A distinct QSS region is not observed in the non-equilibrium cooling case. This is attributed to the relatively rapid VT transitions that work to equilibrate the vibrational energy distribution upon recombination, which is comparatively slow. The total time to reach equilibrium by the STS model in the non-equilibrium cooling case is five times and three times greater than those of the QK and TCE models, respectively.

Modeling the flowfield surrounding hypersonic vehicles poses a significant challenge owing to the complex thermochemical processes in the shock-heated region and in the wake. These regions are characterized by both thermal and chemical non-equilibrium owing to the disparity in characteristic relaxation times. Models based on state-to-state (STS) kinetics have received considerable attention and have been employed to study relaxation behind strong shock waves1,2 in nozzle flows3,4 and in boundary layer flows.5 In particular, the N2(  1Σg+)+N2(  1Σg+) system,6–10 the N2(  1Σg+)+N(  4Su) system,11–14 and the O2(X3Σg)+O(  3P) system15–17 have been studied extensively owing to their importance in Earth atmospheric reentry. These efforts have demonstrated the inadequacy of legacy multi-temperature models, but most importantly, these studies have provided a wealth of information regarding the complex evolution of the internal energy distributions under thermal and chemical non-equilibrium. A recent study by Macdonald et al.8 examined the non-equilibrium processes during dissociation of the N2(  1Σg+)+N2(  1Σg+) system. It was established that dissociation processes are driven by the high-energy tail of the rovibrational distribution. Furthermore, nearly all dissociation occurs under quasi-steady state (QSS) conditions at high temperatures, while approximately half of the dissociation occurs after QSS at low (<13 000 K) temperatures. This is in contrast to dissociation processes for the N2(  1Σg+)+N(  4Su) system, where more than half of dissociation occurs before QSS conditions are established.13 Panesi et al.13 also found that rotational-translational (RT) and vibrational-translational (VT) relaxation times converge at high temperatures in nitrogen, again contrary to the conventional idea that vibrational relaxation occurs over comparatively longer characteristic time scales.

Recent work by Andrienko and Boyd15 outlines a comprehensive master equation study of the O2 + O system based on the O2(  3Σg)+O(  3P) potential energy surface by Varandas and Pais.18 The STS rovibrationally resolved cross sections for inelastic (non-reactive), exchange, and dissociation collisional processes employed in the master equation calculation were obtained from quasi-classical trajectory (QCT) calculations.15 The cross sections were computed on a single surface corresponding to the ground electronic state (accounting for 1/27 degeneracy). VT processes resulting from this surface were thus assumed to represent processes occurring on the other eight surfaces, while a constant scaling factor of 16/3 proposed by Nikitin was used to account for dissociation from electronic excited states.19 A number of important results emerged from the analysis of this system. Exchange reactions were found to play a significant role in the thermalization of the internal energies during rovibrational heating, particularly for highly excited quasi-bound states. Furthermore, equilibration of rotational and vibrational degrees of freedom was found to occur over similar time scales for a broad range of temperatures, which is in contrast to recent studies of N2–N and O2–Ar molecular systems.13,20 Thermochemical relaxation under dissociation, recombination, and exchange reactions was also studied in a 0-d heat bath calculation for a range of translational temperatures. A QSS region is established owing to the competition between internal energy relaxation (via RVT non-reactive transitions and exchange reactions) and dissociation. Remarkably, however, the characteristic QSS rotational and vibrational temperatures exhibit a maximum at 10 000 K and decrease for higher heat bath temperatures. This is in stark contrast to both the N2–N and N2–N2 systems, which show a monotonic increase in the internal temperatures under QSS.

State-specific models have also been developed for use in the direct simulation Monte Carlo (DSMC) method for N2(  1Σg+)+N(  4Su), N2(  1Σg+)+N2(  1Σg+), and O2(  3Σg)+O(  3P)21–25 in place of the widely used phenomenological models, such as the Larsen-Borgnakke (LB) model for internal energy relaxation26 and the total-collision-energy (TCE) and quantum-kinetic (QK) models for chemical reactions.27–29 The state-specific cross sections for internal energy exchange favor monoquantum transitions, in contrast to the Larsen-Borgnakke model in which transitions to the ground state are most favored.24,25 A recent study by Borges et al.23 compared the state-specific dissociation (SSD) rate coefficients obtained from O2(  3Σg)+O(  3P) QCT calculations with rate coefficients from QCT and the TCE model in which the TCE model parameters were calibrated with the QCT data.23 It was found that in non-equilibrium conditions, the TCE model overestimates the QCT dissociation rate coefficient by more than an order of magnitude in vibrationally cold conditions and underestimates the dissociation rate coefficient in vibrationally hot conditions.

Although it is necessary to include recombination reactions to capture the correct mixture composition and temperature in the shock layer and in the wake of hypersonic flows,27,30 a formal treatment of the recombination reaction in state-resolved hypersonic DSMC calculations has been lacking. Recombination modeling in DSMC is included in both TCE and QK models,31 in which the probability of recombination is computed from the dissociation rate coefficient and an equilibrium constant. The internal state of the recombining molecule is then selected based on an assumed Boltzmann distribution at a given temperature, which cannot capture the up-pumping effect observed in detailed state-resolved calculations.4,32,33 Recently, however, Gimelshein and Wysong proposed a new chemistry model in DSMC, inclusive of both dissociation and recombination reactions, that utilizes the vibrational bias model of dissociation and an empirical microscopic recombination model.34–36 The dissociation probability, Pd, is computed for each dissociation collision based on the translational and internal energies of the colliding pair and includes two parameters, λ and Ad. Ad is calibrated to match a known rate constant of the specific reaction, and λ is used to specify the degree of vibrational favoring. The recombination probability, Pr, is computed from the equilibrium recombination rate coefficient kr of Arrhenius form for a rovibrational state (v, j), which is obtained as the ratio of the corresponding dissociation rate coefficient kd to the equilibrium constant (further details may be found in Ref. 34). An important feature of this model is that by its construction, both the principle of detailed balance and microscopic reversibility are satisfied inherently. Detailed balance is achieved through the determination of kr as mentioned above, while microscopic reversibility is satisfied indirectly through the selection of the internal energy states of recombining pairs; the internal energy state of dissociated molecules is stored in a list in each computational cell, and the internal state of each recombining molecule is selected from this list. Thus, the inelastic recombination (backward) cross section for each rovibrational state (v, j) is inherently reconstructed through the sampling of this list, which is populated according to the dissociation (forward) processes. A similar approach was also employed for the TCE model.1 After assigning the internal state to the recombining molecule, the relative translational energy between the recombined product and the third body is determined based on energy conservation. While the chemistry model by Gimelshein and Wysong addresses both dissociation and recombination reactions in an elegant way that satisfies detailed balance and microscopic reversibility, the computed dissociation (and, therefore, recombination) probability relies on the calibration of Ad and λ to match the known equilibrium or non-equilibrium reaction rate constants.

This work introduces a new vibrational state-to-state DSMC chemistry model that has been developed for dissociation/recombination and internal energy relaxation of the O2 + O system. This model is based on vibrational state-specific VT and dissociation cross sections obtained from QCT calculations16 involving the Varandas and Pais potential energy surface (PES),18 which has been previously used for studying the O2 + O system.15,23,37,38 The database utilized in this work was developed by Esposito and Capitelli which includes rotationally averaged cross sections.38 As will be detailed later in Sec. III, this database is a suitable first step for the extension of the two-step binary collision (TSBC) model to state-resolved cross sections,35 which is the aim of the present work. It is noted that a rovibrationally resolved cross-sectional database has recently been developed by Andrienko and Boyd15 based on the Varandas and Pais PES and new accurate PES by Varga et al.39 and Paukku et al.40,41 for the O3 and O4 systems is also recently available. This paper aims to first assess the vibrational state-resolved DSMC chemistry model based on the O3 database from Esposito and Capitelli. This will enable important comparisons to the fully resolved rovibrational state-to-state model based on the O3 database of Andrienko and Boyd and to the new O3 PES by Varga et al. Recent comparisons of the Varandas and Pais PES to the 11A′ surface by Varga show significant differences in the structure of the potential as well as the vibrational relaxation rates. A full discussion of these comparisons may be found in Grover et al.42 

The present model employs the two-step binary collision (TSBC) approach for DSMC recombination, which was first introduced by Koura35 for a modified forced harmonic oscillator. The TSBC approach treats recombination through two binary interaction processes, including orbiting pair formation (constrained by detailed balance) followed by subsequent stabilization involving a third body collision. As will be presented later, the internal energy of the recombined molecule formed through the stabilization collision is determined through the recombination cross section via microscopic reversibility, but in such a way that this principle is satisfied directly from the known vibrational state-resolved QCT dissociation cross sections. The internal energy level dynamics of the O2 + O system are investigated for 0-d vibrational heating and cooling cases in the presence of both dissociation and recombination reactions, and these results are compared against existing phenomenological models. To the authors’ best knowledge, this is the first time that both dissociation and recombination reactions have been included in the study of energy level dynamics of the O2 + O system in DSMC. The results from the relaxation simulations capture a distinct QSS region in the vibrational heating cases examined, but a QSS region is not established during vibrational cooling.

The rest of the paper is organized as follows: a brief overview of the Larsen-Borgnakke, TCE and QK phenomenological reaction models is presented in Sec. II, followed by a detailed outline of the TSBC approach employed in the current chemistry model in Sec. III. In Sec. IV A, a comparison of the dissociation and recombination rate coefficients computed from QCT with those of the TCE, QK, and STS models in equilibrium (T = Tv) and non-equilibrium (TTv) conditions is presented. Results of the 0-d vibrational heating and cooling cases are presented in Sec. IV B, and the observed differences in the QSS regions predicted by each of the models are discussed, followed by conclusions.

In this section, a brief overview of the phenomenological DSMC chemistry models for internal energy exchange, dissociation, and recombination reactions is presented. Although many chemistry models have been developed for DSMC, we limit our discussion to the Larsen-Borgnakke (LB), TCE, and QK models for later comparison to the STS chemistry model.

The Larsen-Borgnakke model has been widely used for describing energy transfer in non-reactive inelastic collisions26 owing to its simplicity and its ability to reproduce relaxation with the Landau-Teller form. A comparison of the LB cross sections with the QCT cross sections for VT transitions, however, shows significant differences, as seen in Fig. 1. The LB cross sections were evaluated using the variable hard sphere (VHS) model and a vibrational collision number based on the Millikan-White form with Park’s high-temperature correction43 and the Gimelshein DSMC collision number correction factor.44 The graphs in the figure show a comparison of the LB and QCT cross sections for VT transitions (including inelastic and exchange collisions; excluding elastic collision) of the O2 + O system from initial vibrational states v = 0, v = 20, and v = 40 to all possible final states. A rotational temperature of Trot = 300 K and a relative collision energy of Etr = 2 eV are specified. It can be seen that the LB cross section (and hence the probability for transition into the final vibrational state, v′, indicated by the x-axis) is the largest for transitions into the lowest state, regardless of the initial vibrational state. In contrast, QCT cross sections favor mono-quantum transitions over multi-quantum transitions. Thus, the most probable post-collision vibrational state is in general not the ground state and depends upon the initial vibrational state. This QCT cross-sectional behavior is well-established in the literature and has been observed by others for the O2 + O PES employed in this work.17 

FIG. 1.

VT transition cross sections of O2 + O system with pre-collision vibrational states (a) v = 0, (b) v = 20, and (c) v = 40, Trot = 300 K and Etr = 2 eV from QCT and LB models.

FIG. 1.

VT transition cross sections of O2 + O system with pre-collision vibrational states (a) v = 0, (b) v = 20, and (c) v = 40, Trot = 300 K and Etr = 2 eV from QCT and LB models.

Close modal

One of the main advantages of the QK model is that it does not require information about the rates of reaction.45,46 Starting from the LB model, the post-collision vibrational states are chosen from a Boltzmann distribution corresponding to a collision energy equivalent temperature. The collision energy is the sum of the relative translational energy and the vibrational energy. For instance, in a collision of O2 and O, the collision energy used to redistribute the vibrational energy of O2 is given by Ecoll = Etr + Ev, where Ev is the pre-collision vibrational energy of O2. The post-collision vibrational state v′ is chosen using a uniform random number from 0 to vmax = ⌊Ecoll/(kBθ)⌋ and accepted with a probability of

PvibQK(v)=1vkBθvibEcoll3/2ω.
(1)

In the phenomenological models considered here, a simple harmonic oscillator (SHO) assumption is used, meaning the vibrational levels are spaced evenly by a characteristic vibrational temperature θvib. A dissociation event occurs if the maximum possible vibrational state vmax satisfies vmax > θd/θvib, where θd is the characteristic dissociation temperature. The QK dissociation rate can be expressed explicitly for VHS gas at equilibrium temperature T,

kdQK(T)=2σrefπ1/2ϵTTref1ω2kBTrefμ1/2v=0vmax1Γ(5/2ω,(vmaxv)θvib/T)exp(vθvib/T)(1exp(θvib/T))Γ(5/2ω),
(2)

where σref is the cross section at temperature Tref and ϵ is the symmetry factor. Recombination reactions occur when there is a third body particle within a ternary collision volume Vcoll. The corresponding recombination probability can be expressed as

PrQK=nTVcoll,
(3)

where nT is the number density of the third body particles and Vcoll is suggested by Bird28 in the form Vcoll=a(T/θvib)bVref, where the reference volume Vref is conveniently set to the volume of the sphere with the radius equal to the sum of the reference radii of the three atoms.28 For an equilibrium VHS gas, the recombination rate coefficient can be expressed as

krQK=2σrefπ1/2ϵTTref1ω2kBTrefμ1/2nTa(T/θvib)bVref.
(4)

The ratio of the dissociation and recombination rate coefficients from Eqs. (2) and (4) is equated to the equilibrium constant to determine the parameters a and b. The QK model was also recently improved to satisfy detailed balance.47 

The TCE model27 expresses the probability of dissociation in terms of an Arrhenius dissociation rate constant of the form

kdTCE=ΛdTηdexp(Ed/kBT),
(5)

with Arrhenius parameters Λ and η. Ed denotes the dissociation energy. Unlike the QK model, the collision energy of the TCE model takes the rotational energy into account, which is given as a sum of the internal and relative translational energies; i.e., for the O2 + O collision, Ecoll = Etr + Ev + Ej, where Ev and Ej are the vibrational and the rotational energy of pre-collision O2, respectively. The probability of dissociation is given by

PdTCE=π1/2ϵΛdTrefηd2σref(kBTref)ηd1+ωΓ(ζ¯+5/2ω)Γ(ζ¯+ηd+3/2)×μ2kBTref1/2(EcollEd)ηd+ζ¯+1/2ζ¯+3/2ω,
(6)

where ζ¯ is the average number of internal degrees of freedom contributing to the collision energy. Again, an equilibrium constant Keq can be determined to cast the recombination rate in terms of the dissociation rate. Since it is assumed that the recombination reactions have zero activation energy, the recombination rate can be expressed as

krTCE=ΛrTηr.
(7)

The recombination probability which yields the rate in Eq. (7) is

PrTCE=ϵπ1/2nTΛrTrefηr2σrefΓ(5/2ω)Γ(ηr+3/2)μ2kBTref1/2EcollkBTrefηr1+ω.
(8)

Similar to the QK model, the TCE model was recently improved to satisfy detailed balance.1 

To provide a more detailed description of dissociation and recombination in DSMC, the phenomenological QK and TCE models are replaced by a high-fidelity vibrationally resolved state-to-state (STS) description using QCT collision cross sections. The cross sections are computed by performing scattering simulations of individual atoms or molecules for a prescribed pre-collision state to determine the resulting post-collision state. This provides a set of state-resolved dissociation cross sections for direct use in DSMC, while the recombination cross sections are inferred according to the TSBC model. In this section, we first introduce the STS cross-sectional database employed in this work. We then provide an overview of the TSBC model for the O3 system and assess the appropriateness of an indirect (two step) recombination mechanism for O + O interactions, including O + O orbiting pairs and their associated lifetime. We then develop in detail the mathematical descriptions of dissociation and recombination kinetics for the TSBC model, followed by their implementation in DSMC.

This work employs the available cross-sectional database from Esposito and Capitelli37,38 for the O2(  3Σg)+O(  3P) system.18 The potential energy surface was constructed using the double many-body expansion (DMBE) method for the ground electronic state of ozone. Rovibrational levels were computed using the WKB approximation for the O2 potential.18 The cross sections are available at pre-collision vibrational states: vinit = 0, 1, 3, 5, 7, 10, 15, 20, 25, 30, 35, 39, 40–46; rotational temperatures: Trot = 300, 3000, 20 000 K; and relative translational energies Etr = 0.006–9 eV in increments of 0.01 eV. Details of this database may be found in Refs. 16, 18, 37, and 38.

It should be noted that the reactions in these QCT calculations proceed on a single surface (accounting for 1/27 of the degenerate states) corresponding to the ground electronic state. VT processes involving the other eight surfaces, representing the remaining degeneracies, were not considered in the generation of this database. Dissociation processes, however, are unlikely to proceed adiabatically because the required dissociation energy (>5 eV) is much higher than the 1Δ first excited molecular state, which is about 1 eV above the molecular ground state. Accordingly, Esposito and Capitelli introduced a variable multisurface factor38 in the current database to account for dissociation from electronically excited states. This variable multisurface factor is an extension of the constant 16/3 factor introduced by Nikitin,19 in which the dissociation cross section at an electronically excited state is estimated as its degeneracy divided by the degeneracy of the electronic ground state. The variable multisurface factor employed in this work is computed as the summation of the degeneracy of each excited electronic state which lies at or below the energy of the rovibrational level under consideration, divided by the degeneracy of the ground state. Further details of this variable multisurface factor and of the QCT calculations are provided in Ref. 38. In summary, there is only one PES considered in the current work, which is the representative of 1/27 of O2 + O collisions, and the variable multisurface factor proposed by Esposito and Capitelli is employed to account for dissociation from electronically excited states.

As mentioned earlier, the database employed in this work provides rotationally averaged cross sections. The QCT trajectories are carried out for a specified relative translational energy, Etr, and an initial vibrational state, v. The initial rotational state, j, is drawn from a Boltzmann distribution at temperature Trot, based on the WKB approximation of rovibrational levels for the O2 potential. The resulting cross sections are obtained by summing over all final rotational states, j, corresponding to a given final vibrational state, v (full details may be found in Ref. 38),

σdv(Etr,Trot)=1Qrjgj2j+1expEv,jkTrotfMSσdv(Etr,j).
(9)

In this model, the influence of rotational energy on the reaction probability is limited to cell-averaged Trot. Because of this, only Etr and v have a collision-specific influence on the collision probability, and the reaction of interest is expressed in terms of v only,

O2(v)+O3O.
(10)

Since the cross sections contained in this QCT database have significant scatter, a polynomial with the form below is used to fit the cross sections as a function of relative translational energy,

σ(x)=a2x2+a1x1+a0+a1x+a2x2+a3x3+a4x4+a5x5.
(11)

A bilinear interpolation technique is applied to determine cross sections for vibrational states (or rotational temperatures) not explicitly available in the database. Although the cross sections employed in this work are rotationally averaged, the study by Andrienko and Boyd suggests that O2 rotational nonequilibrium is an important factor at low temperatures, where recombination is most significant. Further studies should investigate the role of rotational nonequilibrium in recombination using a set of rotationally resolved cross sections.

Although the current database does not provide rotationally resolved cross sections, it is well suited for Koura’s two-step binary collision (TSBC) model,35 which considers only vibrational non-equilibrium. Our objective is to thus extend Koura’s approach for the current QCT database. The TSBC model describes dissociation and recombination reactions in the following two steps:

O2(v)+O(3)kfvkvfOP(f)+O(3),
(12)
OP(f)kOPksplitO(1)+O(2).
(13)

The first step of the dissociation reaction [the forward direction of Eq. (12)] is the collision of O2 with a third-body oxygen atom, O(3), which forms an orbiting pair, OP, in the vibrational free state f. The OP is formed from the atoms of the O2 molecule and is assumed as a unique state in the chemistry model. The second step of dissociation [the forward direction of Eq. (13)] is the splitting of the OP into its two oxygen atoms, referred to as O(1) and O(2).

Recombination reactions are modeled as the reverse of the dissociation reaction detailed above. The first step of recombination [the backward direction of Eq. (13)] is the formation of an OP from oxygen atoms O(1) and O(2), and the third body O(3) remains unchanged. The second step [the backward direction of Eq. (12)] is the stabilization of the OP through collision with the third-body oxygen atom O(3). This produces the recombined O2, while the third-body O(3) remains unchanged. A diagram outlining the TSBC process is shown in Fig. 2.

FIG. 2.

Schematic of TSBC model35 for dissociation and recombination reactions.

FIG. 2.

Schematic of TSBC model35 for dissociation and recombination reactions.

Close modal

Before discussing the details of the TSBC model as applied to the QCT database, we first assess the appropriateness of an indirect (two step) mechanism for the recombination of O2. We consider specifically the evidence for O + O orbiting and the associated lifetime of these orbiting pairs. It is well-known that orbiting occurs if there exists an attractive portion of the interaction potential,48–50 and such orbiting has been observed by Weaver and Alexeenko51 and Subramaniam and Stephani52 for the O3 system based on the Varandas and Pais PES considered in this work. A recent study by Elchinger and co-workers53 also demonstrates long-lived O + O orbiting interactions based on the Hulbert-Hirschfelder pairwise potential. Although this potential is different from the pairwise potential employed in the Varandas and Pais PES, both feature an attractive portion of the interaction potential, and thus, long-lived orbiting interactions are observed.

The time spent in orbit prolongs the interaction between the O atoms, making stabilization of the orbiting pair (i.e., O2 recombination) by a third body more likely as the orbiting time increases. Obviously, as the orbiting time increases, the range of impact parameters (i.e., the cross section) that ultimately results in appreciable orbiting decreases. If appreciable orbiting is defined as greater than three revolutions, then the corresponding orbiting pair cross section is approximately σOP = O(10−32. In what follows, by employing Eq. (29) with Etr = 0.1 eV, the corresponding orbiting lifetime is approximately τOP = O(10−12)s. The orbiting time is thus very short, which is important for the development of the TSBC model below.

In order to compute the recombination probability in DSMC, several assumptions are introduced in the development of the TSBC model: (1) the orbiting pair is assumed as a unique state, and the mechanisms in Eqs. (12) and (13) are assumed as the only possible recombination channels. Therefore, only the two-step (indirect) recombination pathway is considered; other possible channels, such as the direct three-body recombination mechanism, are not considered. (2) It is assumed that the OP splitting and its formation [forward and backward reactions of Eq. (13)] are sufficiently rapid such that the formation and splitting processes may be assumed to be in equilibrium. This assumption implies that the excitation cross section from the vibrational state v to the OP free state may be approximated by the dissociation cross section. (3) The velocity distribution function of oxygen atoms follows the Maxwellian distribution, which allows us to derive the analytical form of expressions used in the TSBC model. (4) The orbiting lifetime is assumed to be constant in this work as suggested by Koura.35 As will be discussed below, the value of orbiting lifetime itself is never evaluated in the current TSBC model; rather, only the product of the orbiting pair cross section and its lifetime is evaluated.

In Koura’s TSBC model, the dissociation reaction can be thought of as the following reaction:

O2+O(3)O(1)+O(2)+O(3).
(14)

Note that the reaction in Eq. (14) is simply the result of the two consecutive reactions in Eqs. (12) and (13). Energy conservation is applied to dissociation by

Etr+Ev=Etr+Ed,
(15)

where Etr is the relative translational energy of O2 with O(3) [Eq. (12)]. Although the rotational energy does not directly influence the dissociation probability in this model, Erot goes into the post-reaction scattering of O(1) and O(2) (i.e., the splitting of the OP) through energy and momentum conservation.

The description of dissociation through Eqs. (12) and (13) is necessary in order to establish microscopic reversibility for the reverse recombination process described in Sec. III B. The approach outlined above, in which the rotational energy has no influence on the dissociation probability but is used to split the dissociating molecule, is identical to the treatment of the QK model and similar to the TCE model. Elastic collisions are not included in the current QCT database and are accounted for using the phenomenological VHS model. Finally, the probability of a DSMC dissociation for each collision is simply the ratio of the dissociation cross section to the total cross section,

Pd(Etr,v,Trot)=σdv(Etr,Trot)σtot.
(16)

The total cross section σtot is the sum of all possible interaction cross sections for O + O2 collisions, including elastic, inelastic (non-reactive), exchange, and dissociation reactions.

Finally, the state-specific dissociation rate coefficient can be computed by integrating the state-specific dissociation cross section over the velocity space,

kdv(T,Trot)=fM(c)fA(c)gσdv(Etr,Trot)dcdc,
(17)

where fM and fA are the velocity distribution functions of the colliding molecule and atom, respectively, and are assumed to be Maxwellian for a given temperature T. It should be noted that Eq. (17) is only computed in post-processing for the purpose of analysis; the dissociation rate coefficient is not required as input in DSMC simulations.

Because the cross sections for the reverse recombination reaction 3O → O2 + O are not included in the QCT database, these cross sections must be backed out from the forward dissociation reaction through detailed balance and microscopic reversibility. Again following the TSBC method by Koura,35 recombination reactions occur through two steps: (i) OP formation [the backward direction of Eq. (13)] and (ii) OP stabilization [the backward direction of Eq. (12)], as shown in Fig. 2. The stabilization cross sections are determined from the state-specific dissociation cross sections through microscopic reversibility, while the orbiting pair cross sections are derived from detailed balance using equilibrium rates. The derivations of OP formation and stabilization cross sections are introduced next, followed by the state-specific recombination probability for DSMC.

The forward process of OP splitting and the backward process of OP formation are assumed to be rapid and in equilibrium,35 where

ksplitnOP=kOPnO2.
(18)

Therefore, the sum of the stabilization rate over all vibrational states can be referred to as the overall recombination rate according to the above assumption, which gives a recombination rate

dnO2dtkrnO3=vkfvnOPnO.
(19)

Substituting for nOP from Eq. (18) into Eq. (19) gives the recombination rate coefficient at a given equilibrium temperature T,

kr(T)=kOP(T)τOPvkfv(T).
(20)

The lifetime of the orbiting pair (τOP) is introduced in Eq. (20) where τOP=ksplit1, which is assumed to be constant.

In order to determine the recombination probability, the recombination rate coefficient must be expressed in terms of the dissociation rate coefficient. The equilibrium dissociation rate coefficient is expressed by the sum of the vibrational state-specific dissociation rate coefficient weighted by the Boltzmann vibrational distribution at temperature T,

kd(T)=vkdv(T)exp(Ev/kT)vexp(Ev/kT).
(21)

Recalling the second assumption outlined in the beginning of this section, the rate coefficient from state vf is assumed to be approximately equal to the dissociation rate coefficient (kdvkvf). The law of mass action can then be used to relate the dissociation and recombination rate constants as54 

kr(T)kd(T)=QMtrQMrQMvQMel(QAtrQAel)2exp(Ed/kT),
(22)

where QM and QA are partition functions of molecular and atomic oxygen, respectively. The superscripts tr, r, v, and el denote the translational, rotational, vibrational, and electronic mode, respectively. Combining Eqs. (22), (20), and (21) and employing detailed balance between kvf(T) and kfv(T),

kvf(T)exp(Ev/kT)=kfv(T)exp(Ed/kT),
(23)

we obtain

kOP(T)τOP=QMtrQMrQMel(QAtrQAel)2.
(24)

It is important to note that since the forward and backward reactions of Eq. (13) are assumed to be in equilibrium, the TSBC model does not capture the pressure dependence for the recombination rate coefficient and displays the standard fall-off region with respect to pressure. Such a pressure dependence can however be captured if the steady-state approximation is used, which will be addressed in future work. It is worth noting that the recombination rate coefficient in the current work is consistent with the low-pressure limit, for which the equilibrium and steady-state approximations yield the same rate coefficient. To compute the probability of recombination in DSMC, consider the expression for kOP from kinetic theory,

kOP=12fA(cA1)fA(cA2)gσOP(Etr)dcA1dcA2,
(25)

where the OP formation rate coefficient is obtained from the OP formation cross section when integrated over the velocity distribution functions for the reactant atoms. Assuming that the oxygen atoms follow a Maxwellian distribution at the temperature T, we have the following:

kOP(T)=128kTπμOO1/2(kT)2EtrσOPexp(Etr/kT)dEtr,
(26)

where μO–O is the reduced mass of two oxygen atoms. Here, we introduce the average value of the quantity EtrσOP, which can be expressed as

EtrσOP¯=EtrσOP(T)exp(Etr/kT)dEtrexp(Etr/kT)dEtr.
(27)

Combining Eqs. (24) and (26) with the ratio of the electronic partition function ηel=(QAel)2/QMel, the translational partition function for species s, Qstr=(2πkTms/h2)3/2, and the rotational partition function QMr=T/2θr and applying Eq. (27), we obtain

EtrσOP¯τOP=h34πmOkθrηel,
(28)

where mO is the molecular weight of an oxygen atom and h is Planck’s constant. To implement in DSMC, the product σOPτOP for each O + O collision pair is computed based on its relative collision energy Etr,

σOPτOP=h34πmOkBθrηel1Etr.
(29)

The total cross section of OP stabilization [the backward direction of Eq. (12)], σstab, is the sum of the cross section of all possible fv transitions,

σstab=vσfv(Etr).
(30)

The cross section of stabilization into each vibrational state v (σfv) can be determined in terms of the corresponding excitation cross section σvf (from QCT database) via microscopic reversibility and energy conservation,55 

EtrOPOσfv(EtrOPO)=EtrO2(v)Oσvf(EtrO2(v)O),
(31)
EtrOPO+Ed=EtrO2(v)O+Ev,
(32)

where EtrOPO and EtrO2(v)O are the relative translational energy between OP and O and between O2(v) and O, respectively. After computing EtrO2(v)O by Eq. (32) for the specific vibrational state v and substituting this value into Eq. (31), σfv can then be solved. Note that the rotational energy of the stabilized O2 is assumed equal to the orbiting energy (the relative translational energy between two atoms with respect to the center-of-mass frame), consistent with the dissociation process,35 while the bond energy goes into vibration.

Previous studies have proposed various approaches to three-body collisions in DSMC,56 including NTC-based methods.57 The collision frequencies outlined in Koura’s TSBC method use a time counter approach, which may include null collisions (candidates which are not selected for collision).35 The no-time-counter (NTC) scheme for the TSBC model is introduced instead for use in this model. For a binary collision between species p and q, the number of collision candidates selected is

Ncollpq=1/ϵNpNqF(σpqg)maxΔt/V,
(33)

where g denotes the relative velocity of the collision pair, F is the number of real molecule represented by a simulation particle, and ϵ is the symmetric factor [(1) if p and q are the same species; (2) otherwise]. The candidates are selected for collision with probability

Pcollpq=σpqg(σpqg)max.
(34)

To extend the NTC scheme for a three-body collision based on the TSBC approach, the number of collision candidates of OP + O is computed first, which is expressed as

Ncollrec=NOPNOF(σstabgOPO)maxΔt/V
(35)

with NOP defined35 as

NOP=(NOnO/2)(σOPgOOτOP)max,
(36)

and nO = NOF, we substitute Eq. (36) into Eq. (35) to obtain the number of recombination collision candidates,

Ncollrec=1/2NO3F2(σOPgOOτOP)max(σstabgOPO)maxΔtV.
(37)

This form is similar to the NTC scheme for binary collisions but involves a third N term to account for the three oxygen atoms involved in the recombination reaction. These recombination collision candidates are then selected for recombination with probability

Prec=(σOPgOOτOP)(σstabgOPO)(σOPgOOτOP)max(σstabgOPO)max.
(38)

The vibrational and rotational energies of the recombined molecule are determined according to

Prec(v)=σfvvσfv,
(39)
Erot=1/2μOO(gOO)2.
(40)

The relative translational energy of the recombined species and the third body is then determined through energy conservation.

The no-time-counter scheme of the TSBC algorithm is shown schematically in Fig. 3, and the steps for implementation are as follows:

  1. Three particles (atoms) are randomly selected from the DSMC particle list (first two particles are OP candidates; third particle is a stabilization candidate).

  2. Compute σOPτOP [Eq. (29)] from the first two particles.

  3. Compute σfv for each possible recombined vibrational state v [Eqs. (31) and (32)]. For example, for v ∈ [0–46], take v = 0, solve Eq. (32) for EtrO2(v)O (all other quantities known), and then solve Eq. (31) for σfv (all other quantities known); repeat for other v.

  4. Compute σstab [Eq. (30)].

  5. Compute Prec [Eq. (38)]; use acceptance-rejection to determine if recombination occurs.

  6. Select vibrational state of the recombined molecule [Eq. (39) via acceptance-rejection].

  7. Assign rotational energy of the recombined molecule equal to orbiting energy [Eq. (40)].

  8. Assign relative translational energy of the recombined molecule and third body through energy conservation.

FIG. 3.

No-time-counter scheme of TSBC algorithm in DSMC.

FIG. 3.

No-time-counter scheme of TSBC algorithm in DSMC.

Close modal

The results presented in this section are focused on the relaxation of the state-specific vibrational mode through both thermal and chemical relaxation processes. Although the rotationally averaged cross sections are obtained by assuming that the pre-collision rotational energy is Boltzmann distributed at a temperature Trot, the rotational mode is not constrained to be in equilibrium with either translational or vibrational modes, both in the cross-sectional calculation and in the DSMC simulations.

The TCE, QK, and STS models are first tested in equilibrium conditions, where T = Tv = 4000–20 000 K. In Fig. 4, the QCT rate coefficients are computed by integrating the QCT cross-sectional data directly. All of the TCE parameters for dissociation and recombination reactions are calibrated from the QCT data for a consistent comparison. The experimental data from Ibraguimova58 are also plotted in Fig. 4, in which the low-temperature (2000–6000 K) and high-temperature (6000–11 000 K) ranges of the dissociation rate coefficient are provided for comparison in Fig. 4(a). We note that the rates reported by Ibraguimova et al. include contributions from O2–O2 collisions, and contributions from O2–O interactions are considered negligible. The corresponding recombination coefficients [Fig. 4(b)] are computed based on the experimental dissociation rate coefficients and the equilibrium constant via detailed balance. The dashed line of the experimental data denotes the extrapolation from the valid temperature region.

FIG. 4.

Equilibrium total (a) dissociation and (b) recombination coefficients compared among the STS, TCE, and QK models and experimental data.58 Solid lines correspond to experimental conditions; dashed lines correspond to extrapolation beyond experimental conditions. “QCT + detailed balance” denotes that the value is computed via detailed balance and dissociation rate coefficient based on QCT database.

FIG. 4.

Equilibrium total (a) dissociation and (b) recombination coefficients compared among the STS, TCE, and QK models and experimental data.58 Solid lines correspond to experimental conditions; dashed lines correspond to extrapolation beyond experimental conditions. “QCT + detailed balance” denotes that the value is computed via detailed balance and dissociation rate coefficient based on QCT database.

Close modal

The agreement between the STS implementation in DSMC and the QCT data is excellent. The phenomenological TCE model also shows good agreement with the QCT results in the equilibrium conditions since the TCE parameters are calibrated from the QCT data. The QK model shows reasonable agreement with the dissociation coefficients from the QCT calculations, but comparison with the recombination coefficients is relatively poor, particularly at high temperatures [Fig. 4(b)]. Recall that the QK parameters are not explicitly calibrated to match the reaction rate coefficients; rather, the parameters a and b are calibrated to match the equilibrium constant (as discussed earlier in Sec. II). Therefore, it is not surprising to have some discrepancy in the equilibrium dissociation/recombination rate coefficients between QK and STS models.

In comparison with the experimental data, the dissociation rate coefficients of the STS model show good agreement for T < 7000 K. The STS recombination rate coefficients show reasonable agreement with the low-temperatureexperimental data. At high temperatures (T > 7000 K), the STS model shows a higher dissociation rate coefficient than the experiment. This is also observed by Kulakhmetov et al.17 Similar discrepancies are observed for the recombination rate coefficient at high temperatures. This overprediction may be a result of the use of rotationally averaged cross sections from the Varandas and Pais PES as well as the nonadiabatic correction factor that is used to account for dissociation through electronically excited reaction pathways. It is also noted that the experiments do not directly measure dissociation rates. Rather, the rates are reconstructed from the absorption data for O2–O2 collisions, using gas dynamic relations and determination of a vibrational temperature from the low-lying vibrational states using the Landau-Teller (harmonic oscillator) model.58 To obtain the final dissociation rate coefficient of O2–O, a factor is employed based on the O2–O2 rate coefficient, kdO2O=3.5kdO2O2.58 

Figure 5 shows the vibrational state-specific dissociation and recombination coefficients at T = 7000 K for the STS model compared with the TCE and QK models. This comparison provides a more detailed picture of how the total rate coefficients (shown previously in Fig. 4) compare by examining the contribution to the equilibrium total rate coefficient by each vibrational state. At low vibrational states, the STS model is characterized by a relatively low dissociation and recombination rate coefficient compared to the TCE and QK models, while at higher vibrational states, the STS rate coefficients exceed those of the TCE/QK models. This indicates a preference for dissociation/recombination from high-lying internal (vibrational) states.

FIG. 5.

Equilibrium state-specific (a) dissociation and (b) recombination coefficients at T = 7000 K comparing STS, TCE, and QK models.

FIG. 5.

Equilibrium state-specific (a) dissociation and (b) recombination coefficients at T = 7000 K comparing STS, TCE, and QK models.

Close modal

The rate coefficients of the TCE, QK, and STS models are compared under non-equilibrium conditions in Fig. 6. The coefficients are computed for T = 10 000 K, with the vibrational temperature ranging from Tv = 4000–20 000 K. Results from the QCT data are also shown for comparison. It is observed in general that the rate coefficients of the TCE and QK phenomenological models have weak vibrational temperature dependence for the dissociation rate coefficient, and these models perform poorly in strong non-equilibrium conditions. The recombination rate coefficients show no vibrational dependence and are approximately constant across the vibrational temperature range based on the current model.

FIG. 6.

Comparison of (a) dissociation and (b) recombination coefficients at T = 10 000 K with variable Tv using STS, TCE, and QK models. “QCT + detailed balance” denotes that the value is computed via detailed balance and dissociation rate coefficient based on QCT database.

FIG. 6.

Comparison of (a) dissociation and (b) recombination coefficients at T = 10 000 K with variable Tv using STS, TCE, and QK models. “QCT + detailed balance” denotes that the value is computed via detailed balance and dissociation rate coefficient based on QCT database.

Close modal

In this section, the dynamics of the vibrational distribution function and the state-resolved reaction rates arestudied through a series of 0-d non-equilibrium heat bath calculations. The simulations track the chemical relaxation of the gas in an isochoric reactor, in which the translational temperature is held constant. The vibrational temperature is initialized at the beginning of the simulation below or above the translational temperature such that either vibrational heating or cooling occurs as a result of the relaxation process. Dissociation, recombination, and VT transitions are modeled in the STS DSMC simulations using the QCT cross-sectional data and the TSBC recombination model, while RT transitions and elastic collisions are modeled using the LB and VHS models, respectively. It should be noted that the vibrational collision number, Zv, which controls the rate of vibrational relaxation in the TCE and QK model implementations was fit to a constant value to match the relaxation rate of the state-based model for a consistent comparison. This step was taken in efforts to minimize the influence of differences in non-reactive inelastic collisions when interpreting the current results. The chemical reaction parameters for the TCE and QK models are shown in Table I. The parameters for the TCE model are computed by fitting the reaction rate coefficients from the QCT database with the Arrhenius form, and the correction suggested by Gimelshein and co-workers59 is applied for employing the discrete vibrational energy in DSMC. The parameters a and b in the QK model28 are used to compute the collision volume for the third-body particle using the same equilibrium constant as the STS model.

TABLE I.

TCE and QK chemistry model parameters for the O/O2 mixture.

ModelΛdηdΛrηrab
TCE 3.7784 × 10−11 −0.713 373 7.892 03 × 10−46 0.161 092 … … 
QK … … … … 0.05 −1.1 
ModelΛdηdΛrηrab
TCE 3.7784 × 10−11 −0.713 373 7.892 03 × 10−46 0.161 092 … … 
QK … … … … 0.05 −1.1 

1. Vibrational heating case (dissociation regime)

The first case presented is the 0-d vibrational heating case, in which the system is held at a constant translational temperature, T = 7000 K, and the initial vibrational and rotational temperatures are set to Tv = Trot = 2000 K. The molecular/atomic oxygen gas mixture is initialized with number density 6.79 × 1024 and 1.48 × 1024, respectively, such that the system undergoes net dissociation during relaxation toward chemical equilibrium. Approximately 40 000 simulation particles were employed to obtain good statistics.

The results of the simulation are shown in a series of figures to highlight the important stages of relaxation and the energy level dynamics corresponding to the observed macroscopic behavior. Figures 7(a) and 7(b) show the evolution of vibrational temperature and mole fraction of atomic and molecular oxygen as a function of time. The system eventually reaches thermal and chemical equilibrium at T = 7000 K, where the equilibrium mole fraction of the molecular oxygen is about 0.6%. Figures 8(a) and 8(b) indicate the simulation times [(a)–(e)] at which vibrational distribution functions (Fig. 9) and state-specific dissociation/recombination rates (Fig. 10) are sampled. From Fig. 7(a), the evolution of the system can be divided into 4 distinct stages: (i) vibrational heating, (ii) quasi-steady state (QSS), (iii) final relaxation, and (iv) equilibrium. In stage (i), the energy added to the system by the isothermal reservoir goes into heating the vibrational mode. This is seen as a rapid initial increase in Tv at the beginning of the simulation. This increase in vibrational temperature is also seen in Figs. 9(a)–9(c) as distribution broadening from t = 0 s to t = 1.2 × 10−7 s.

FIG. 7.

Comparison of (a) vibrational temperature and (b) mole fraction as a function of time during vibrational heating with T = 7000 K.

FIG. 7.

Comparison of (a) vibrational temperature and (b) mole fraction as a function of time during vibrational heating with T = 7000 K.

Close modal
FIG. 8.

Schematic indicating the simulation times [(a)–(e)] during (a) thermal and (b) chemical relaxation at which vibrational distribution functions (Fig. 9), and state-specific dissociation/recombination rates (Fig. 10) are sampled.

FIG. 8.

Schematic indicating the simulation times [(a)–(e)] during (a) thermal and (b) chemical relaxation at which vibrational distribution functions (Fig. 9), and state-specific dissociation/recombination rates (Fig. 10) are sampled.

Close modal
FIG. 9.

Evolution of the vibrational distribution function for the heating case using the STS model at time (a) 4.0 × 10−9 s, (b) t = 1.2 × 10−8 s, (c) t = 1.2 × 10−7 s, (d) t = 2.2 × 10−7 s, and (e) t = 8.8 × 10−7 s, and corresponding Boltzmann distributions fitted to the data in (a)–(e) are shown in (f). Initial vibrational temperature is indicated by the black dashed line; vibrational temperature(s) of each distribution are indicated by solid lines.

FIG. 9.

Evolution of the vibrational distribution function for the heating case using the STS model at time (a) 4.0 × 10−9 s, (b) t = 1.2 × 10−8 s, (c) t = 1.2 × 10−7 s, (d) t = 2.2 × 10−7 s, and (e) t = 8.8 × 10−7 s, and corresponding Boltzmann distributions fitted to the data in (a)–(e) are shown in (f). Initial vibrational temperature is indicated by the black dashed line; vibrational temperature(s) of each distribution are indicated by solid lines.

Close modal
FIG. 10.

Evolution of state-specific dissociation and recombination rates for vibrational heating using the STS model at time (a) 4.0 × 10−9 s, (b) t = 1.2 × 10−8 s, (c) t = 1.2 × 10−7 s, (d) t = 2.2 × 10−7 s, and (e) t = 8.8 × 10−7 s.

FIG. 10.

Evolution of state-specific dissociation and recombination rates for vibrational heating using the STS model at time (a) 4.0 × 10−9 s, (b) t = 1.2 × 10−8 s, (c) t = 1.2 × 10−7 s, (d) t = 2.2 × 10−7 s, and (e) t = 8.8 × 10−7 s.

Close modal

Figure 7(a) shows that all models spend approximately the same duration of time in the initial vibrational heating stage, which is reasonable since Zv in both the TCE and QK models is determined from the STS vibrational relaxation rate. The system then enters into the quasi-steady state (QSS) region in stage (ii), where the dissociation rate becomes comparable to the vibrational relaxation rate. This results in a competition between dissociation and vibrational heating, which work to deplete and replenish the vibrational states, respectively, and a nearly stationary vibrational distribution from t = 1.2 × 10−7 s to t = 2.2 × 10−7 s is observed in Fig. 9(e). It is noted that the TCE and QK phenomenological models exhibit strongernet dissociation, resulting in comparatively short QSS regions compared to the STS model, whose QSS region is longer than TCE and QK by a factor of ten and four, respectively. These results are consistent with the fact that the phenomenological models overestimate the dissociation coefficient for Tv < T, as shown in Fig. 6(a). The total time to reach equilibrium for the STS model is approximately 3.5 times that of the TCE models and 1.5 times that of the QK model.

The evolution of the vibrational distribution function and the state-resolved dissociation and recombination rates of STS model are plotted in Figs. 9 and 10. The local vibrational temperature is computed by the energy-equivalent approach.13 Again, the instances in time indicated in the figure captions correspond to times indicated by labels (a)–(e) in Fig. 8. At t = 4.0 × 10−9 s, the VT transition pumps the molecules from low to high vibrational states and leads to theoverpopulation at the high-lying vibrational states [Fig. 9(a)]. At t = 1.2 × 10−8 s, the vibrational distribution continues to broaden and becomes closer to the Boltzmann distribution. The start and end of the QSS region correspond with t = 1.2 × 10−7 s and t = 2.2 × 10−7 s, respectively. At this stage, dissociation becomes more important and starts to deplete the tail population of the distribution. At the same time, VT transitions are working to populate the tail, resulting in a competition between the two processes. As seen in Figs. 9(c) and 9(d), the populations of the vibrational states at these times are nearly stationary, and the vibrational temperature [indicated by the slope of the green and blue lines fit to the data in (c) and (d)] is nearly constant [Fig. 9(f)]. It is also noted that the vibrational distribution function can be represented by two distinct temperatures, where the vibrational temperature at high-lying states is about 2000 K lower than the low-lying states during QSS. The corresponding state-specific dissociation and recombination rates during QSS from the STS model are shown in Figs. 10(c) and 10(d). During QSS, the state-specific dissociation and recombination rates are nearly stationary. It is also noted that at these times, the dissociation/recombination rates of the high-lying vibrational states are already balanced during QSS (indicated by the overlapping of the red and blue symbols). As expected, these rates are balanced across all vibrational states at equilibrium as required by detailed balance [Fig. 10(e)].

2. Vibrational cooling case (recombination regime)

The final case presented is the 0-d vibrational cooling case, in which the system is held at a constant translational temperature, T = 2000 K, and the initial vibrational and rotational temperatures are set to Tv = Trot = 7000 K. The molecular/atomic oxygen gas mixture is initialized with mole fractions of 0.2 and 0.8, respectively, such that the system undergoes net recombination during relaxation toward chemical equilibrium. The results of the simulation are again shown in a series of figures to highlight the energy level dynamics corresponding to the observed relaxation behavior. Figures 11(a) and 11(b) show the temperature and mole fraction as a function of time during the vibrational cooling process. The evolution of the vibrational distribution function and the state-resolved dissociation/recombination rates are plotted in Figs. 12 and 13.

FIG. 11.

Comparison of (a) vibrational temperature and (b) mole fraction as a function of time during vibrational cooling with T = 2000 K. Simulation times [(a)–(d)] are indicated at which vibrational distribution functions (Fig. 12) and state-specific dissociation/recombination rates (Fig. 13) are sampled.

FIG. 11.

Comparison of (a) vibrational temperature and (b) mole fraction as a function of time during vibrational cooling with T = 2000 K. Simulation times [(a)–(d)] are indicated at which vibrational distribution functions (Fig. 12) and state-specific dissociation/recombination rates (Fig. 13) are sampled.

Close modal
FIG. 12.

Evolution of the vibrational distribution function for the cooling case using the STS model at time (a) 2.5 × 10−8 s, (b) t = 3.5 × 10−8 s, (c) t = 9.5 × 10−7 s, and (d) t = 8.5 × 10−4 s. The evolution of the corresponding Boltzmann distributions fitted to the data in (a)–(d) are shown in (e).

FIG. 12.

Evolution of the vibrational distribution function for the cooling case using the STS model at time (a) 2.5 × 10−8 s, (b) t = 3.5 × 10−8 s, (c) t = 9.5 × 10−7 s, and (d) t = 8.5 × 10−4 s. The evolution of the corresponding Boltzmann distributions fitted to the data in (a)–(d) are shown in (e).

Close modal
FIG. 13.

Evolution of state-specific dissociation and recombination rates for vibrational cooling using the STS model at time (a) 2.5 × 10−8 s, (b) t = 3.5 × 10−8 s, (c) t = 9.5 × 10−7 s, and (d) t = 8.5 × 10−4 s.

FIG. 13.

Evolution of state-specific dissociation and recombination rates for vibrational cooling using the STS model at time (a) 2.5 × 10−8 s, (b) t = 3.5 × 10−8 s, (c) t = 9.5 × 10−7 s, and (d) t = 8.5 × 10−4 s.

Close modal

Unlike the heating case, there is no apparent plateau in the vibrational temperature profile during cooling [Fig. 11(a)], and a QSS region is not observed for any of the models. The competing processes which would contribute to QSS in the cooling case are VT transitions, which work to deplete the vibrational distribution tail, and recombination, which works to replenish the tail population. The recombination rate, however, is not large enough to balance the rate of depletion by VT transitions, which is comparatively faster. Recall that for the vibrational heating case, QSS was observed in which the dissociation rates shown in Fig. 10 were on the order of 1030 s−1, while the recombination rates shown in Fig. 13 during vibrational cooling are more than an order of magnitude lower. This suggests that the VT transitions dominate the thermal relaxation process, and the temperature profile observed in Fig. 11(a) resembles a relatively rapid thermal relaxation, with a comparatively slow chemical relaxation. This is also reflected in Fig. 11(b), which shows the mole fractions of molecular and atomic oxygen as a function of time. Relaxation toward chemical equilibrium requires more time in the cooling case, owing to slow recombination rates.

It is observed that the STS model reaches chemical equilibrium more slowly than the phenomenological TCE and QK models. This is due to the vibrational favoring of dissociation/recombination in the STS model, as shown in Fig. 5. Themolecules at high-lying vibrational states tend to be dissociated, while the recombining molecules have a preference to enter into the high-lying vibrational states for the STS model. Therefore, the molecules produced by recombination are very readily re-dissociated. For the phenomenological models, the trend of the state-resolved dissociation coefficients is similar to the STS model; however, the recombining molecules tend to enter into low-lying vibrational states. Therefore, the net recombination rate of STS model is lower than those of the phenomenological models and requires more time to reach chemical equilibrium.

The evolution of the vibrational distribution function and the state-resolved dissociation and recombination rates are plotted in Figs. 12 and 13, respectively. The labels (a)–(d) correspond to the results at time 2.5 × 10−8 s, t = 3.5 × 10−8 s, t = 9.5 × 10−7 s, and (d) t = 8.5 × 10−4 s, respectively. During the vibrational cooling stage, the vibrational temperature drops rapidly due to VT transitions which drive the tail population to lower-lying vibrational energy states. In the first stages of rapid vibrational cooling [Figs. 12(a) and 12(b)], the low-lying vibrational states relax to lower temperatures first. At these early times, both dissociation/recombination reactions are important due to the relatively high vibrational temperature and the high mole fraction of atomic oxygen, respectively, as shown in Fig. 13(a).

As the vibrational temperature is further decreased and approaches the translational temperature, the vibrational distribution becomes closer to Boltzmann; however, the very high-lying vibrational states are in strong non-equilibrium, as shown in Figs. 12(b) and 12(c). For instance, in Fig. 12(c), the vibrational distribution function can be represented by three distinct temperatures. The vibrational temperature of the intermediate energy states is much hotter (≈14 000 K) due to recombination, while the vibrational temperature of the tail (≈2000 K) is close to the translational temperature due to the relatively strong dissociation in the tail. The overpopulated high-lying states assist the dissociation; therefore, even though the mole fraction of molecular oxygen is low (χO ≈ 0.2) and the vibrational temperature is low (≈2000 K), there are still appreciable amounts of dissociation. There is an obvious overpopulation at high-lying vibrational states in Figs. 13(b) and 13(c), which is mainly due to the up-pumping by the recombination. In Fig. 13(d), the gas has nearly reached chemical equilibrium.

In this work, a vibrational state-specific model for dissociation and recombination reactions in DSMC is presented based on QCT cross-sectional calculations. The two-step binary collision framework is applied for recombination reactions in which both microscopic reversibility and state-specific detailed balance are used to constrain the orbiting pair and stabilization cross sections. The dissociation and recombination rate coefficients are discussed in both equilibrium and non-equilibrium conditions. In comparing the STS model with the phenomenological models, in the equilibrium condition, TCE and STS show good agreement in the overall reaction rate coefficients.

The TCE, QK, and STS models were compared in 0-d vibrational heating and cooling cases. The vibrational heating case shows a distinct QSS region during relaxation, while QSS is not established during vibrational cooling. In all cases, the vibrational distribution function exhibited strong non-equilibrium behavior during the vibrational heating and cooling processes. In the heating case, the phenomenological models start to dissociate ahead of the STS model because of the overestimation of the dissociation coefficient at low-lying vibrational states, which also leads to a shorter time to reach thermal and chemical equilibrium. In the cooling case, the phenomenological models again reach chemical equilibrium faster than STS model. A stronger dissociation rate is observed in the STS model which is due to the assistance of the recombination, resulting in a longer time to reach equilibrium. Finally, the balance between dissociation and recombination in the STS model is first achieved in the high vibrational states, which indicates that the net dissociation/recombination during the heating/cooling processes is driven by the low-lying vibrational states. Because of its importance at low temperatures for O2 relaxation,15 rotational nonequilibrium should be addressed by extending the current work to a rotationally resolved cross-sectional database.

This work was supported by an Early Career Faculty grant from NASA’s Space Technology Research Grants Program and the NASA Space Technology Research Fellowship, and by the Air Force Office of Scientific Research under award number FA9550-17-1-0127. The authors gratefully acknowledge Dr. F. Esposito for providing the QCT database and for many useful discussions.

1.
I.
Borges Sebastião
and
A.
Alexeenko
, “
Consistent post-reaction vibrational energy redistribution in DSMC simulations using TCE model
,”
Phys. Fluids
28
,
107103
(
2016
).
2.
M.
Panesi
,
A.
Munafò
,
T.
Magin
, and
R.
Jaffe
, “
Nonequilibrium shock-heated nitrogen flows using a rovibrational state-to-state method
,”
Phys. Rev. E
90
,
013009
(
2014
).
3.
A.
Munafo
,
M.
Panesi
,
R.
Jaffe
,
G.
Colonna
,
A.
Bourdon
, and
T.
Magin
, “
QCT-based vibrational collisional models applied to nonequilibrium nozzle flows
,”
Eur. Phys. J. D
66
,
188
(
2012
).
4.
A.
Guy
,
A.
Bourdon
, and
M.-Y.
Perrin
, “
Consistent multi-internal-temperatures models for nonequilibrium nozzle flows
,”
Chem. Phys.
420
,
15
24
(
2013
).
5.
I.
Armenise
and
F.
Esposito
, “
Dissociation–recombination models in hypersonic boundary layer O2/O flows
,”
Chem. Phys.
398
,
104
110
(
2012
).
6.
R. L.
Jaffe
,
D. W.
Schwenke
,
G.
Chaban
, and
W.
Huo
, “
Vibrational and rotational excitation and relaxation of nitrogen from accurate theoretical calculations
,” AIAA Paper 2008-1208,
2008
.
7.
G.
Chaban
,
R. L.
Jaffe
,
D. W.
Schwenke
, and
W.
Huo
, “
Dissociation cross-sections and rate coefficients for nitrogen from accurate theoretical calculations
,” AIAA Paper 2008-1209,
2008
.
8.
R. L.
Macdonald
,
R. L.
Jaffe
,
D. W.
Schwenke
, and
M.
Panesi
, “
Construction of a coarse-grain quasi-classical trajectory method. I. Theory and application to N2–N2 system
,”
J. Chem. Phys.
148
,
054309
(
2018
).
9.
J. D.
Bender
,
P.
Valentini
,
I.
Nompelis
,
Y.
Paukku
,
Z.
Varga
,
D. G.
Truhlar
,
T.
Schwartzentruber
, and
G. V.
Candler
, “
An improved potential energy surface and multi-temperature quasiclassical trajectory calculations of N2 + N2 dissociation reactions
,”
J. Chem. Phys.
143
,
054304
(
2015
).
10.
N.
Parsons
,
D. A.
Levin
,
A. C. T.
van Duin
, and
T.
Zhu
, “
Modeling of molecular nitrogen collisions and dissociation processes for direct simulation Monte Carlo
,”
J. Chem. Phys.
141
,
234307
(
2014
).
11.
A.
Munafò
,
S.
Venturi
,
R. L.
Macdonald
, and
M.
Panesi
, “
State-to-state and reduced-order models for recombination and energy transfer in aerothermal environments
,” AIAA Paper 2016-0505,
2016
.
12.
O.
Kunova
,
E.
Kustova
, and
A.
Savelev
, “
Generalized Treanor–Marrone model for state-specific dissociation rate coefficients
,”
Chem. Phys. Lett.
659
,
80
87
(
2016
).
13.
M.
Panesi
,
R. L.
Jaffe
,
D. W.
Schwenke
, and
T. E.
Magin
, “
Rovibrational internal energy transfer and dissociation of N2(1Σg+)–N(4Su) system in hypersonic flows
,”
J. Chem. Phys.
138
,
044312
(
2013
).
14.
M. S.
Grover
,
N.
Singh
,
T. E.
Schwartzentruber
, and
R. L.
Jaffe
, “
Dissociation and internal excitation of molecular nitrogen due to N2-N collisions using direct molecular simulation
,” AIAA Paper 2017-0660,
2017
.
15.
D. A.
Andrienko
and
I. D.
Boyd
, “
Rovibrational energy transfer and dissociation in O2–O collisions
,”
J. Chem. Phys.
144
,
104301
(
2016
).
16.
F.
Esposito
,
I.
Armenise
,
G.
Capitta
, and
M.
Capitelli
, “
O–O2 state-to-state vibrational relaxation and dissociation rates based on quasiclassical calculations
,”
Chem. Phys.
351
,
91
98
(
2008
).
17.
M.
Kulakhmetov
,
M.
Gallis
, and
A.
Alexeenko
, “
Ab initio-informed maximum entropy modeling of rovibrational relaxation and state-specific dissociation with application to the O2 + O system
,”
J. Chem. Phys.
144
,
174302
(
2016
).
18.
A.
Varandas
and
A.
Pais
, “
A realistic double many-body expansion (DMBE) potential energy surface for ground-state O3 from a multiproperty fit to ab initio calculations, and to experimental spectroscopic, inelastic scattering, and kinetic isotope thermal rate data
,”
Mol. Phys.
65
,
843
860
(
1988
).
19.
E.
Nikitin
, in
Theory of Elementary Atomic and Molecular Processes in Gases
(
Clarendon Press
,
Oxford
,
1974
), Chap. 7.
20.
J. G.
Kim
and
I. D.
Boyd
, “
Thermochemical nonequilibrium modeling of electronically excited molecular oxygen
,” AIAA Paper 2014-2963,
2014
.
21.
J. G.
Kim
and
I. D.
Boyd
, “
Monte Carlo simulation of nitrogen dissociation based on state-resolved cross sections
,”
Phys. Fluids
26
,
012006
(
2014
).
22.
C.
Zhang
and
T. E.
Schwartzentruber
, “
Consistent implementation of state-to-state collision models for direct simulation Monte Carlo
,” AIAA Paper 2014-0866,
2014
.
23.
I.
Borges Sebastião
,
M.
Kulakhmetov
, and
A.
Alexeenko
, “
DSMC study of oxygen shockwaves based on high-fidelity vibrational relaxation and dissociation models
,”
Phys. Fluids
29
,
017102
(
2017
).
24.
T. J.
Wilson
and
K. A.
Stephani
, “
State-to-state vibrational energy modeling in DSMC using quasiclassical trajectory calculations for O + O2
,” AIAA Paper 2016-3839,
2016
.
25.
T. J.
Wilson
,
T.-J.
Pan
, and
K. A.
Stephani
, “
State-to-state dissociation and recombination modeling in DSMC using quasi-classical trajectory calculations for O + O2
,” AIAA Paper 2018-1487,
2018
.
26.
C.
Borgnakke
and
P. S.
Larsen
, “
Statistical collision model for Monte Carlo simulation of polyatomic gas mixture
,”
J. Comput. Phys.
18
,
405
420
(
1975
).
27.
G.
Bird
,
Molecular Gas Dynamics and the Direct Simulation of Gas Flows
(
Clarendon
,
1994
).
28.
G.
Bird
,
The DSMC Method
(
Create Space Independent Publishing Platform
,
2013
).
29.
I.
Wysong
,
S.
Gimelshein
,
Y.
Bondar
, and
M.
Ivanov
, “
Comparison of direct simulation Monte Carlo chemistry and vibrational models applied to oxygen shock measurements
,”
Phys. Fluids
26
,
043101
(
2014
).
30.
I. D.
Boyd
, “
Analysis of vibration-dissociation-recombination processes behind strong shock waves of nitrogen
,”
Phys. Fluids A
4
,
178
185
(
1992
).
31.
M. A.
Gallis
,
J. R.
Torczynski
,
S. J.
Plimpton
,
D. J.
Rader
, and
T.
Koehler
, “
Direct simulation Monte Carlo: The quest for speed
,”
AIP Conf. Proc.
1628
,
27
36
(
2014
).
32.
G.
Colonna
,
L. D.
Pietanza
, and
M.
Capitelli
, “
Recombination-assisted nitrogen dissociation rates under nonequilibrium conditions
,”
J. Thermophys. Heat Transfer
22
,
399
406
(
2008
).
33.
D. A.
Andrienko
,
K.
Neitzel
, and
I. D.
Boyd
,
“Vibrational relaxation and dissociation in O2–O mixtures
,” AIAA Paper 2016-4021,
2016
.
34.
S.
Gimelshein
and
I.
Wysong
, “
DSMC modeling of flows with recombination reactions
,”
Phys. Fluids
29
,
067106
(
2017
).
35.
K.
Koura
, “
A set of model cross sections for the Monte Carlo simulation of rarefied real gases: Atom–diatom collisions
,”
Phys. Fluids
6
,
3473
3486
(
1994
).
36.
D. C.
Wadsworth
and
I. J.
Wysong
, “
Vibrational favoring effect in DSMC dissociation models
,”
Phys. Fluids
9
,
3873
3884
(
1997
).
37.
F.
Esposito
and
M.
Capitelli
, “
The relaxation of vibrationally excited O2 molecules by atomic oxygen
,”
Chem. Phys. Lett.
443
,
222
226
(
2007
).
38.
F.
Esposito
and
M.
Capitelli
, “
Quasiclassical trajectory calculations of vibrationally specific dissociation cross-sections and rate constants for the reaction O + O2(v) → 3O
,”
Chem. Phys. Lett.
364
,
180
187
(
2002
).
39.
Z.
Varga
,
Y.
Paukku
, and
D. G.
Truhlar
, “
Potential energy surfaces for O + O2 collisions
,”
J. Chem. Phys.
147
,
154312
(
2017
).
40.
Y.
Paukku
,
K. R.
Yang
,
Z.
Varga
,
G.
Song
,
J. D.
Bender
, and
D. G.
Truhlar
, “
Potential energy surfaces of quintet and singlet O4
,”
J. Chem. Phys.
147
,
034301
(
2017
).
41.
Y.
Paukku
,
Z.
Varga
, and
D. G.
Truhlar
, “
Potential energy surface of triplet O4
,”
J. Chem. Phys.
148
,
124314
(
2018
).
42.
M. S.
Grover
,
T. E.
Schwartzentruber
,
Z.
Varga
, and
D.
Truhlar
, “
Dynamics of vibrational energy excitation and dissociation in oxygen from direct molecular simulation
,” AIAA Paper 2018-0238,
2018
.
43.
C.
Park
, “
Review of chemical-kinetic problems of future nasa missions. I-Earth entries
,”
J. Thermophys. Heat Transfer
7
,
385
398
(
1993
).
44.
N.
Gimelshein
,
S.
Gimelshein
, and
D.
Levin
, “
Vibrational relaxation rates in the direct simulation Monte Carlo method
,”
Phys. Fluids
14
,
4452
4455
(
2002
).
45.
G.
Bird
, “
A comparison of collision energy-based and temperature-based procedures in DSMC
,”
AIP Conf. Proc.
1084
,
245
250
(
2008
).
46.
G.
Bird
, “
The QK model for gas-phase chemical reaction rates
,”
Phys. Fluids
23
,
106101
(
2011
).
47.
G.
Bird
, “
Setting the post-reaction internal energies in direct simulation Monte Carlo chemistry simulations
,”
Phys. Fluids
24
,
127104
(
2012
).
48.
J. P.
Toennies
,
W.
Welz
, and
G.
Wolf
, “
Molecular beam scattering studies of orbiting resonances and the determination of van der Waals potentials for H–Ne, Ar, Kr, and Xe and for H2–Ar, Kr, and Xe
,”
J. Chem. Phys.
71
,
614
642
(
1979
).
49.
J. H.
Ferziger
and
H. G.
Kaper
,
Mathematical Theory of Transport Processes in Gases
(
North-Holland
,
1972
).
50.
M.
Capitelli
,
D.
Bruno
, and
A.
Laricchiuta
,
Fundamental Aspects of Plasma Chemical Physics: Transport
(
Springer
,
2013
), Vol. 74.
51.
A. B.
Weaver
, “
Assessment of high-fidelity collision models in the direct simulation Monte Carlo method
,” Ph.D. thesis,
Purdue University
,
2015
.
52.
S.
Subramaniam
and
K. A.
Stephani
, “
Characterization of state-resolved transport for O + O2 collisions
,” AIAA Paper 2018-1484,
2018
.
53.
A.
Mahfouf
,
P.
André
,
G.
Faure
, and
M.-F.
Elchinger
, “
Theoretical and numerical study of transport collision integrals: Application to O (3P)–O(3P) interaction
,”
Chem. Phys.
491
,
1
10
(
2017
).
54.
W. G.
Vincenti
and
C. H.
Kruger
,
Introduction to Physical Gas Dynamics
(
Wiley
,
New York
,
1965
), Vol. 246.
55.
B. H.
Mahan
, “
Microscopic reversibility and detailed balance. An analysis
,”
J. Chem. Educ.
52
,
299
(
1975
).
56.
M.
Gallis
,
R.
Bond
, and
J.
Torczynski
, “
Assessment of reaction-rate predictions of a collision-energy approach for chemical reactions in atmospheric flows
,”
AIAA Paper
2010
-
4499
,
2010
.
57.
P.
Norman
and
T.
Schwartzentruber
, “
Classical trajectory calculation direct simulation Monte Carlo: GPU acceleration and three body collisions
,”
AIAA Paper AIAA-2013-1200
,
2013
.
58.
L.
Ibraguimova
,
A.
Sergievskaya
,
V. Y.
Levashov
,
O.
Shatalov
,
Y. V.
Tunik
, and
I.
Zabelinskii
, “
Investigation of oxygen dissociation and vibrational relaxation at temperatures 4000–10 800 K
,”
J. Chem. Phys.
139
,
034317
(
2013
).
59.
S.
Gimelshein
,
N.
Gimelshein
,
D.
Levin
,
M.
Ivanov
, and
I.
Wysong
, “
On the use of chemical reaction rates with discrete internal energies in the direct simulation Monte Carlo method
,”
Phys. Fluids
16
,
2442
2451
(
2004
).