2.1. Vibrationally partitioned equation of state

Consider the $i{\rm th}$ species within some mixture, and assume it may be treated as an ideal gas. This species is considered to be in thermal equilibrium within some volume, $V$ , if the probability distribution function (PDF) value of the $ \boldsymbol{n} {\rm th}$ quantised molecular energy level, denoted $ \widetilde {f} _i( \boldsymbol{n} )$ , is determined via a Boltzmann distribution with a single temperature (Vincenti & Kruger Reference Vincenti and Kruger1965), i.e.

(2.1) \begin{align} \widetilde {f} _i( \boldsymbol{n} ) = \frac {{\rm e}^{- \widetilde { \varepsilon } _i( \boldsymbol{n} )/ T_i}}{Q_i(T)}, \quad Q_i(T) = \sum _{ \boldsymbol{n} } {\rm e}^{- \widetilde { \varepsilon } _i( \boldsymbol{n} )/ T_i}, \end{align}

where $Q_i(T)$ is the species partition function of molecular energy states, $ \widetilde { \varepsilon } _i( \boldsymbol{n} )$ is the species’ $ \boldsymbol{n} {\rm th}$ energy state expressed in units of temperature and $T_i$ is the species temperature. The $\, \widetilde {\boldsymbol{\cdot }} \,$ notation is used to denote state-specific quantities that depend on the quantum state.

If the system is not in thermal equilibrium, then the definition of $ \widetilde {f} _i( \boldsymbol{n} )$ in (2.1) must be reconsidered. Following convention (Vincenti & Kruger Reference Vincenti and Kruger1965; McQuarrie Reference McQuarrie2000), we assume that $ \widetilde { \varepsilon } _i( \boldsymbol{n} )$ can be partitioned into translational ( $t$ ), rotational ( $r$ ), vibrational ( $v$ ) and electronic ( $el$ ) energy components so that

(2.2) \begin{align} \widetilde { \varepsilon } _i( \boldsymbol{n} ) = \widetilde { \varepsilon } _{t,i}(n_t) + \widetilde { \varepsilon } _{r,i}(n_r) + \widetilde { \varepsilon } _{v,i}(n_v) + \widetilde { \varepsilon } _{el,i}(n_{el}), \end{align}

where the subscript denotes the component, and the energy state $ \boldsymbol{n}$ is now represented as a vector of quantum numbers for the translational, rotational, vibrational and electronic state, i.e. $ \boldsymbol{n} = (n_t,n_r,n_v,n_{el})$ . We assume that coupled rotational–vibrational energy can either be neglected (via the rigid-rotor approximation) or be subsumed into either the vibrational or rotational energy component. Substituting (2.2) into (2.1), and assuming each component is characterised by a component temperature, $T_{\alpha ,i}$ , the PDF becomes

(2.3) \begin{align} \widetilde {f} _i( \boldsymbol{n} ) = \prod _{\alpha } \left ( \frac {{\rm e}^{- \widetilde { \varepsilon } _{\alpha ,i}(n_\alpha )/ T_{\alpha ,i}}} {Q_{\alpha ,i}(T_{\alpha ,i})} \right ), \quad Q_{\alpha ,i}(T_{\alpha ,i}) = \sum _{n_\alpha } {\rm e}^{- \widetilde { \varepsilon } _{\alpha ,i}(n_\alpha )/T_{\alpha ,i}}, \end{align}

for $\alpha =t,\,r,\,v\ \mbox{and}\ el$ . When $T_{\alpha ,i} = T_i$ for all $\alpha$ , the species is in thermal equilibrium and (2.1) is attained. Otherwise, the system is in a state of non-equilibrium.

With a PDF that accounts for thermal non-equilibrium, we are able to determine non-equilibrium relations for the internal energy and pressure. From the component partition functions in (2.3) and Boltzmann’s entropy relation (Vincenti & Kruger Reference Vincenti and Kruger1965), the state-averaged (continuum-scale) specific internal energy and pressure components are given by

(2.4) \begin{align} e_{\alpha ,i} = \frac {\mathcal{R}_u T_{\alpha ,i}^2}{W_i} \,\boldsymbol{\cdot }\, \frac {\partial \left ( \ln Q_{\alpha ,i} \right )}{ \partial T_{\alpha ,i}} , \quad p_{\alpha ,i} = \frac {M_i \mathcal{R}_u T_{\alpha ,i}}{W_i} \,\boldsymbol{\cdot }\, \frac {\partial \left ( \ln Q_{\alpha ,i} \right )}{ \partial V} , \end{align}

where $\mathcal{R}_u$ is the universal gas constant, $W_i$ is the species molar mass and $M_i$ is the total mass of the species within the volume, $V$ . The details associated with the component forms of energy and pressure are described in Vincenti & Kruger (Reference Vincenti and Kruger1965). Here, we include a brief summary of assumptions and associated results. For the translational and rotational energy, we assume the spacing between $ \widetilde { \varepsilon } _{t,i}$ and $ \widetilde { \varepsilon } _{r,i}$ states is small at the temperatures expected in high-speed combustion (above 300 K). It follows that $Q_{t,i}$ and $Q_{r,i}$ may be approximated by integrals rather than summations. For the vibrational energy, we will assume that some temperature-dependent oscillator model can be employed (e.g. the harmonic oscillator assumption). Finally, we note that the rotational, vibrational and electronic partition functions are independent of $V$ , so $\partial ( \ln Q_{\alpha ,i} ) / \partial V = 0$ and $p_{\alpha ,i} = 0$ for $\alpha = r,\,v,\ \mbox{and}\ el$ . In other words, the pressure only depends on the translational component of the internal energy.

Under these assumptions, the specific internal energy and pressure are defined as the summations of their respective components as defined in (2.4) and reduce to

(2.5) \begin{align} e_i = e_{t,i}(T_{t,i}) + e_{r,i}(T_{r,i}) + e_{v,i}(T_{v,i}) + e_{el,i}(T_{el,i}), \quad p_i = \rho _i R_i T_{t,i}, \end{align}

where $e_{t,i}$ , $e_{r,i}$ , $e_{v,i}$ and $e_{el,i}$ are the translational, rotational, vibrational and electronic components of the specific internal energy, respectively, each of which is characterised by a corresponding species- and mode-specific temperature. We assume $e_{el,i}$ includes the reference energy state associated with the enthalpy of formation. Additionally, $\rho _i$ is the mass density, and $R_i = \mathcal{R}_u/W_i$ is the species-specific gas constant. At this stage, we do not prescribe a specific model for $e_{\alpha ,i}$ , except to note that each component is related to its corresponding temperature via a functional relationship. For a mixture of $N$ species, the equation of state for ideal gas mixtures is typically based on Dalton’s law of partial pressures (Powers & Paolucci Reference Powers and Paolucci2005), with the mixture’s internal energy and pressure given by

(2.6) \begin{align} e = \sum _{i=1}^N Y_i e_i, \quad p = \sum _{i=1}^N p_i, \end{align}

where $Y_i$ is the species mass fraction.

This mode-partitioned equation of state requires additional transport equations to evolve each mode independently, one for each independent mode. Additionally, we must account for the exchange of energy between modes and species, and determine how each energy component affects and responds to chemical reactions. For a broad range of systems in high-speed combustion, we can reduce the complexity considerably by applying a series of well-established assumptions. First, we assume that the translational energy of each molecular species relaxes within a few intermolecular collisions (Vincenti & Kruger Reference Vincenti and Kruger1965). It follows that the translational PDF for all molecular species will relax to a common translational temperature. Additionally, we also assume that rotational and electronic relaxation is rapid compared with the vibrational relaxation and chemical reaction processes in the high-speed combustion regime (Vincenti & Kruger Reference Vincenti and Kruger1965; Reising et al. Reference Reising, Haller, Clemens, Varghese, Fiévet and Raman2016; Boyd & Schwartzentruber Reference Boyd and Schwartzentruber2017). Even hydrogen, which has a relatively large characteristic rotational temperature of $\approx 90$ K (Atkins et al. Reference Atkins, De Paula and Keeler2023), relaxes $2{-}4$ times faster than vibrational modes at temperatures below 5000 K (Furudate et al. Reference Furudate, Fujita and Abe2006). Thus, we assume that translational, rotational and electronic PDFs for all species can all be characterised by a single temperature associated with energy components that do not experience non-equilibrium, denoted the translational–rotational temperature, or $T_{tr}$ .

These assumptions applied to (2.6) give

(2.7) \begin{align} e = e_{tr}(T_{tr}) + \sum _{i=1}^N Y_i e_{v,i}, \quad p = \rho R T_{tr}, \end{align}

where $e_{tr} = \sum _i Y_i [e_{t,i}(T_{tr}) + e_{r,i}(T_{tr}) + e_{el,i}(T_{tr})]$ represents the mixture translational–rotational(–electronic) energy associated with $T_{tr}$ ; $\rho = \sum _i \rho _i$ is the mass density; and $R = \mathcal{R}_u/W$ is the specific gas constant, with $1/W = \sum _i Y_i/W_i$ denoting the average molar mass. Since the vibrational energy modes are decoupled from the internal energy associated with $T_{tr}$ , additional equations must be employed to model the transport of the species’ vibrational energy.

2.2. Vibrationally partitioned transport equations

The following transport equations are derived for a one-dimensional inviscid flow. Extending this derivation to viscous and higher-dimensional systems does not materially alter this underlying formulation. For a mixture with $N$ species, the one-dimensional reactive Euler equations are given by

(2.8) \begin{align} \frac {\partial }{ \partial t} \left [ \begin{array}{c} \rho \\ \rho u \\ \rho E \\ \rho Y_i \end{array} \right ] + \frac {\partial }{ \partial x} \left [ \begin{array}{c} \rho u \\ \rho u^2 + p \\\ u( \rho E + p ) \\ \rho u Y_i \end{array} \right ] = \left [ \begin{array}{c} 0 \\ 0 \\ 0 \\ W_i \dot {\omega }_i \end{array} \right ], \end{align}

for $i = 1, \ldots , N$ . Here, $x$ and $t$ denote the spatial and temporal components of the solution vector, respectively; $u$ is the velocity; $E = e + u^2/2$ is the total energy; and $\dot {\omega }_i$ is the species molar production rate, with other variables defined in § 2.1.

Substituting the mixture internal energy (2.7) into the total energy transport equation (2.8) and rearranging gives

(2.9) \begin{align} \sum _{i=1}^{N}\left [ \frac {\partial }{ \partial t} \Big ( \rho Y_i e_{v,i} \Big ) + \frac {\partial }{ \partial x} \Big ( \rho u Y_i e_{v,i} \Big ) \right ] + \left [ \frac {\partial }{ \partial t} \Big ( \rho E_{tr} \Big ) + \frac {\partial }{ \partial x} \Big ( u ( \rho E_{tr} + p ) \Big ) \right ] = 0, \end{align}

where $E_{tr} = e_{tr} + u^2/2$ . Now, (2.9) may be posed as $N+1$ separate transport equations for the vibrational energy and total translational–rotational–electronic energy given by

(2.10a) \begin{align} &\frac {\partial }{ \partial t} \Big ( \rho Y_i e_{v,i} \Big ) + \frac {\partial }{ \partial x} \Big ( \rho u Y_i e_{v,i} \Big ) = S_i, \quad i = 1,\,\ldots ,\,N, \end{align}

(2.10b) \begin{align} &{\kern 22pt}\frac {\partial }{ \partial t} \Big ( \rho E_{tr} \Big ) + \frac {\partial }{ \partial x} \Big ( u ( \rho E_{tr} + p ) \Big ) = -\sum _{i=1}^{N} S_i, \end{align}

where $S_i$ is the species’ vibrational energy exchange term, which governs the coupling of vibrational energy with other modes and its creation or destruction through chemical reactions. We choose a form given by

(2.11) \begin{align} S_i = S_{{\scriptsize T-V},i} + S_{{\scriptsize V-V},i} + S_{{\scriptsize C-V},i} , \end{align}

where $S_{{\scriptsize T-V},i}$ represents vibrational energy exchange of each species with translational–rotational–electronic energy of the mixture (denoted T–V exchange), $ S_{{\scriptsize V-V},i}$ represents vibrational energy exchange of each species with vibrational energy of all the other species in the mixture (denoted V–V exchange) and $ S_{{\scriptsize C-V},i}$ represents vibrational energy exchange via chemical reactions (denoted C–V exchange). As described in § 1, T–V exchange controls the relaxation process from vibrational non-equilibrium to equilibrium, V–V exchange controls the redistribution of vibrational energy among each species and C–V exchange controls the creation/destruction of vibrational energy from chemical reactions. Here, we do not include V–V–T exchange, which would account for non-equal spacing of quantised vibrational energy states between molecules in vibrational exchange collisions, by assuming that the corresponding contribution to the translational–rotational energy is negligible compared with T–V exchange (see § 3.2).

The set of equations in (2.10a ) models the vibrational energy transport of each species. This, combined with the reactive Euler equations (2.8), completes the system of vibrationally partitioned transport equations. The translational–rotational energy transport equation (2.10b ) is degenerate in this formulation, i.e. the translational–rotational energy can be calculated directly from the total and vibrational energies. As a final note, we recognise that the internal energy definition may be relaxed to include other non-equilibrium components, such as free-electron motion or rotational non-equilibrium. As such, additional transport equations are needed to close the system. A framework for extending the model to account for these considerations is discussed in Appendix A.

3.1. Translational–vibrational energy exchange

Translational–vibrational energy exchange describes the coupling between translational and vibrational energy within the flow. In general, T–V exchange drives each species towards thermal equilibrium. For a harmonic oscillator, the exchange rate is of the Landau–Teller form (Landau & Teller Reference Landau and Teller1936) given by

(3.1) \begin{align} S_{{\scriptsize T-V},i} = \sum _{j}^N S_{{\scriptsize T-V},ij} , \quad S_{{\scriptsize T-V},ij} = \frac {\rho Y_i e_{v,i}^{tr}}{ \tau _{{\scriptsize T-V},ij} /X_{\!j}} \boldsymbol{\cdot }\left ( 1 - \frac {e_{v,i}}{e_{v,i}^{tr}} \right ), \end{align}

where $e_{v,i}^{tr}$ is the vibrational energy as defined in (2.5) evaluated at $T_{tr}$ , $X_{\!j} = Y_i\boldsymbol{\cdot }W/W_i$ is the species mole fraction and $ \tau _{{\scriptsize T-V},ij}$ is the species pairwise T–V relaxation time scale. The mixture-averaged time scale is given by $1/ \tau _{{\scriptsize T-V},i} = \sum _j X_{\!j}/ \tau _{{\scriptsize T-V},ij}$ . In general, $ \tau _{{\scriptsize T-V},ij}$ depends on the molecular properties of the colliding pair, the translational–rotational temperature and the pressure, each of which in turn dictates the collision rate and probability of T–V exchange, i.e.

(3.2) \begin{align} \tau _{{\scriptsize T-V},ij} = \frac {\tau _{c,ij}}{ P_{{\scriptsize T-V},ij}}, \quad \tau _{c,ij} = \frac { W}{\rho \sigma _{ij} \mathcal{N}_A } \left ( \frac {\pi W_{ij}}{ 8 \mathcal{R}_u T_{tr} }\right )^{1/2}, \end{align}

where $\tau _{c,ij}$ denotes the collision time scale, $ P_{{\scriptsize T-V},ij}$ is the exchange probability, $\sigma _{ij}$ is the collision cross-section, $\mathcal{N}_A$ is Avogadro’s number and $W_{ij} = W_i W_j/(W_i + W_j)$ is the reduced molar mass. Though this model was originally derived for harmonic oscillators, we assume that the relaxation time scale can be modified for more complex vibrational energy models by modifying $ P_{{\scriptsize T-V},ij}$ rather than altering the form of (3.1).

3.2. Vibrational–vibrational energy exchange

Vibrational–vibrational energy exchange models the coupling of vibrational energy among species within a mixture. In general, V–V exchange pushes each species to share the same relative vibrational state (Candler & MacCormack Reference Candler and MacCormack1988, Reference Candler and MacCormack1991). As such, $ S_{{\scriptsize V-V},i}$ may be either positive or negative, depending on the degree of vibrational non-equilibrium of each species in the mixture. The V–V exchange term is derived from the collision rate and the probability of exchange (Candler & MacCormack Reference Candler and MacCormack1988, Reference Candler and MacCormack1991), given by

(3.3) \begin{align} S_{{\scriptsize V-V},i} = \sum _{j\ne i}^N S_{{\scriptsize V-V},ij} , \quad S_{{\scriptsize V-V},ij} = W_i Z_{ij} \left ( P_{{\scriptsize V-V},ji} e_{v,j} - P_{{\scriptsize V-V},ij} e_{v,i} \right ), \end{align}

where $Z_{ij}$ is the (symmetric) collision rate and $ P_{{\scriptsize V-V},ij}$ and $ P_{{\scriptsize V-V},ji}$ are the (not necessarily symmetric) probability of species $i$ increasing the local vibrational energy of species $j$ and vice versa. In this form, $\sum _i S_{{\scriptsize V-V},i} /W_i = 0$ , so V–V exchange does not contribute to the production/destruction of molar-specific translational–rotational energy. At thermal equilibrium, we expect $ S_{{\scriptsize V-V},ij} = 0$ , so, following Knab et al. (Reference Knab, Frühauf and Jonas1992), $ P_{{\scriptsize V-V},ij}$ and $ P_{{\scriptsize V-V},ji}$ are constrained such that

(3.4) \begin{align} P_{{\scriptsize V-V},ji} e_{v,j}^{tr} = P_{{\scriptsize V-V},ij} e_{v,i}^{tr}. \end{align}

With this constraint, and with the collision rate as defined by Candler & MacCormack (Reference Candler and MacCormack1988, Reference Candler and MacCormack1991), (3.3) reduces to the following form:

(3.5) \begin{align} S_{{\scriptsize V-V},ij} = \frac {\rho Y_i e_{v,i}^{tr} }{ \tau _{{\scriptsize V-V},ij} /X_{\!j} } \boldsymbol{\cdot }\left ( \frac {e_{v,j}}{e_{v,j}^{tr}} - \frac {e_{v,i}}{e_{v,i}^{tr}} \right ), \quad \tau _{{\scriptsize V-V},ij} = \frac {\tau _{c,ij}}{ P_{{\scriptsize V-V},ij}}, \end{align}

where $ \tau _{{\scriptsize V-V},ij}$ is the species pairwise V–V redistribution time scale and $\tau _{c,ij}$ is the collision time scale as defined in (3.2). Similar to T–V exchange, the V–V exchange process pushes each species’ vibrational energy to a specific value at a time scale based on the molecular collision rate and the associated exchange probability.

It is important to note that this model assumes that the energy difference between any two vibrational state transitions is negligible. In other words, we do not account for non-equal spacing of quantised vibrational energy states between molecules. If we did, then V–V exchange collisions would either generate or require a small amount of excess translational–rotational energy as part of the exchange process, denoted V–V–T exchange. Here, we assume that this effect is negligible. Otherwise, we would need to employ a more extensive state-specific implementation that we compute and then average (at every unique temperature combination) to determine the state-averaged rate, significantly increasing the modelling/computational burden. To our knowledge, a reduced, state-averaged model that accounts for V–V–T exchange has not yet been developed, so to maintain a tractable modelling framework, and motivated by this historical approach, we neglect V–V–T exchange effects. We believe such consideration is worthy of further investigation and will investigate the limitations of this assumption in a future work.

3.3. Vibrationally informed chemical reactions

The influence of vibrational non-equilibrium on chemical reactions is based on the concept that chemical reactions are sensitive to the vibrational state, specifically that high-lying vibrational states tend to enhance reaction rates (Hammerling et al. Reference Hammerling, Teare and Kivel1959). If such high-lying states are over- or under-populated, the subsequent rate will increase or decrease, respectively. In other words, when $T_{v,i} \gt T_{tr}$ (vibrationally hot), high-lying states are more populated compared with the corresponding distribution at $T_{tr}$ , so the reaction rate should increase compared with when $T_{v,i} = T_{tr}$ . Likewise, when $T_{v,i} \lt T_{tr}$ (vibrationally cold), the rate should be suppressed. To model this effect, we build upon several historical modelling paradigms to develop a chemically–vibrationally consistent model for general, reversible reactions.

The molar production rates are calculated for $J$ reversible reactions of the form

(3.6) \begin{align} \sum _{i=1}^N \nu _{ij}^\prime \chi _i \,\rightleftharpoons \, \sum _{i=1}^N \nu _{ij}^{\prime \prime } \chi _i, \end{align}

where $\nu _{ij}^\prime$ and $\nu _{ij}^{\prime \prime }$ are the species forward and reverse stoichiometric coefficients, respectively, and $\chi _i$ denotes the species (e.g. H $_2$ ). The species molar production rate is given by

(3.7) \begin{align} \dot {\omega }_i = \sum _{j=1}^J \nu _{ij} \left ( \varLambda _{f,j} - \varLambda _{r,j} \right ), \quad \varLambda _{f,j} = k_{f,j} \,\boldsymbol{\cdot }\, \prod _{i=1}^N c_i^{\nu _{ij}^\prime }, \quad \varLambda _{r,j} = k_{r,j} \,\boldsymbol{\cdot }\, \prod _{i=1}^N c_i^{\nu _{ij}^{\prime \prime }}, \end{align}

where $\nu _{ij} = \nu _{ij}^{\prime \prime } - \nu _{ij}^\prime$ is the net stoichiometric coefficient matrix and $\varLambda _{f,j}$ and $\varLambda _{r,j}$ are the forward and reverse reaction rates, respectively. Together, these represent the bulk (state-averaged) rate of change of the species mole fractions and account for all the quantised vibrational states of both the reactants and products. Also, $k_{f,j}$ and $k_{r,j}$ are the forward and reverse rate constants, respectively, and $c_i = \rho Y_i/W_i$ is the species concentration. We choose to model the rate constants using a vibrational efficiency function formulation (Hammerling et al. Reference Hammerling, Teare and Kivel1959; Treanor & Marrone Reference Treanor and Marrone1962; Knab et al. Reference Knab, Frühauf and Messerschmid1995), i.e.

(3.8) \begin{align} k_{f,j} = k_{f,j}^{eq}(T_{tr}) \, \varPhi _f(T_{tr}, \boldsymbol{T_{v}} ), \quad k_{r,j} = k_{r,j}^{eq}(T_{tr}) \, \varPhi _r(T_{tr}, \boldsymbol{T_{v}} ), \end{align}

where $k_{f,j}^{eq}$ and $k_{r,j}^{eq}$ are the forward and reverse rate constants associated with thermal equilibrium, respectively, and $\varPhi _f$ and $\varPhi _r$ are the forward and reverse vibrational efficiency functions, respectively, which modify the rate to reflect the role of vibrational non-equilibrium. Here, $ \boldsymbol{T_{v}}$ denotes the set of vibrational temperatures associated with the system.

Ideally, the efficiency function would be calculated directly via quantum calculations or QCT methods (Voelkel et al. Reference Voelkel, Masselot, Varghese and Raman2016a ). However, this approach remains computationally prohibitive for many relevant reactions when using quantum-level methods, and the generation of the realistic potential energy surfaces necessary for QCT calculations is often laborious and artisanal in nature; for example, a realistic potential energy surface for $\textrm{N}_{2}+\textrm{O}_{2}$ requires over 2000 independent coefficients (Varga et al. Reference Varga, Meana-Pañeda, Song, Paukku and Truhlar2016). With this in mind, we look to derive a generalised formulation for the vibrational efficiency function.

3.3.1. Vibrational efficiency function formulation

The vibrational efficiency function can be derived from the state-specific description of a reaction, given by

(3.9) \begin{align} \sum _{i \in A} \chi _i (v_i) \,\rightleftharpoons \, \sum _{i \in B} \chi _i(w_i), \end{align}

where $A$ and $B$ are the set of reactant and product species, respectively, and $v_i$ and $w_i$ are the corresponding vibrational states of the reactant and product species, respectively ( $v_i$ and $w_i$ may be scalars or vectors depending on the species’ vibrational degrees of freedom). Here, we utilise $A$ (and $B$ ) rather than $N$ to allow for repeat reactants (and products), with each having their own unique vibrational state. Also, we have dropped reaction-specific notation (the $j$ subscript) in this section for brevity. Finally, we assume that this formulation is only vibrationally state-specific, and translational–rotational states are still characterised by $T_{tr}$ .

The state-specific forward reaction rate is given by

(3.10) \begin{align} \widetilde {\varLambda } _f( \boldsymbol v , \boldsymbol w ) = c_{\!A} \, \widetilde {f} _{A}( \boldsymbol v ) \, \widetilde {k} _f( \boldsymbol v , \boldsymbol w ), \end{align}

where $c_{\!A} = \prod _{i \in A} c_i$ is the product of species concentrations, $ \widetilde {f} _{A}( \boldsymbol v ) = \prod _{i \in A} \widetilde {f} _{v,i}(v_i)$ is the product of reactant vibrational PDFs evaluated at $T_{v,i}$ (as described in § 2.1) and $ \widetilde {k} _f( \boldsymbol v , \boldsymbol w )$ is the state-specific rate constant. As before, the $\, \widetilde {\boldsymbol{\cdot }} \,$ notation denotes that the quantity is state-specific, i.e. dependent on the vibrational state, and $ \boldsymbol v$ and $ \boldsymbol w$ are the set of reactant and product vibrational states, respectively. In general, we assume all quantities depend on both $T_{tr}$ and $ \boldsymbol{T_{v}}$ , with different notation used to indicate non-standard temperature dependence. We define the state-specific rate constant as

(3.11) \begin{align} \widetilde {k} _f( \boldsymbol v , \boldsymbol w ) = k_f^{eq} \, \widetilde {\varPhi } _f( \boldsymbol v , \boldsymbol w ), \end{align}

where $k_f^{eq}$ is the rate constant assuming thermal equilibrium (evaluated at $T_{tr}$ ) and $ \widetilde {\varPhi } _f( \boldsymbol v , \boldsymbol w )$ is the state-specific vibrational efficiency function. The state-averaged forward reaction rate, rate constant and vibrational efficiency function are defined by summing over the vibrational states, giving

(3.12) \begin{align} \varLambda _f = c_{\!A} \, k_f, \quad k_f = k_f^{eq} \, \varPhi _f, \quad \varPhi _f = \sum _{ \boldsymbol v , \boldsymbol w } \widetilde {f} _{A}( \boldsymbol v ) \, \widetilde {\varPhi } _f( \boldsymbol v , \boldsymbol w ), \end{align}

which is consistent with (3.7) and (3.8). Here, the summation over $ \boldsymbol v$ includes all reactant vibrational states and corresponds to an average sum weighted by the vibrational partition function, and the summation over $ \boldsymbol w$ includes all product vibrational states. Similar relations hold for the reverse rate, $\varLambda _r$ , as derived from a corresponding state-specific rate constant, $ \widetilde {k} _r( \boldsymbol v , \boldsymbol w )$ , and vibrational efficiency function, $ \widetilde {\varPhi } _r( \boldsymbol v , \boldsymbol w )$ .

For reversible reactions, it is important for the vibrationally informed model to be consistent with the idea of an equilibrium constant, i.e. balance of the forward and reverse rates. Thus, the state-specific efficiency function is now constrained to ensure consistency under detailed balance via the following relation:

(3.13) \begin{align} \widetilde {f} _{A}^{*}( \boldsymbol v ) \, \widetilde {k} _f( \boldsymbol v , \boldsymbol w ) = \mathcal{T} \, \widetilde {f} _B^{*}( \boldsymbol w ) \, \widetilde {k} _r( \boldsymbol v , \boldsymbol w ), \end{align}

where $ \widetilde {f} _A^{*}( \boldsymbol v )$ and $ \widetilde {f} _B^{*}( \boldsymbol w )$ are the products of the reactant and product vibrational PDFs evaluated at some temperature, $T_{*}$ , respectively. Summing over all vibrational states defines $\mathcal{T} = K_c^{eq} \, \varPhi _f^*/ \varPhi _r^*$ , where $K_c^{eq} = k_f^{eq}/k_r^{eq}$ is the equilibrium constant associated with thermal equilibrium, $\varPhi _f^* = \sum _{ \boldsymbol v , \boldsymbol w } \widetilde {f} _A^{*}( \boldsymbol v ) \, \widetilde {\varPhi } _f( \boldsymbol v , \boldsymbol w )$ and $\varPhi _r^* = \sum _{ \boldsymbol v , \boldsymbol w } \widetilde {f} _B^{*}( \boldsymbol w ) \, \widetilde {\varPhi } _r( \boldsymbol v , \boldsymbol w )$ . Now, summing (3.13) over $ \boldsymbol v$ results in the following relation:

(3.14) \begin{align} \frac {\sum _{ \boldsymbol v } \widetilde {f} _A^{*}( \boldsymbol v ) \, \widetilde {\varPhi } _f( \boldsymbol v , \boldsymbol w )}{\varPhi _f^*} = \frac {\sum _{ \boldsymbol v } \widetilde {f} _B^{*}( \boldsymbol w ) \, \widetilde {\varPhi } _r( \boldsymbol v , \boldsymbol w )}{\varPhi _r^*}, \end{align}

which must hold for the vibrationally aware rate constant to be consistent.

One approach to guarantee that (3.14) is satisfied is to assume that $ \widetilde {\varPhi } _f$ and $ \widetilde {\varPhi } _r$ are separable as follows:

(3.15) \begin{align} \widetilde {\varPhi } _f( \boldsymbol v , \boldsymbol w ) = a_f \, \widetilde {\omega } _f( \boldsymbol v ) \, \widetilde {s} _f( \boldsymbol w ), \quad \widetilde {\varPhi } _r( \boldsymbol v , \boldsymbol w ) = a_r \, \widetilde {\omega } _r( \boldsymbol w ) \, \widetilde {s} _r( \boldsymbol v ), \end{align}

where $a_f$ and $a_r$ are functions of $T_{tr}$ , defined to ensure consistency at thermal equilibrium; $ \widetilde {\omega } _f( \boldsymbol v )$ and $ \widetilde {\omega } _r( \boldsymbol w )$ are the forward and reverse efficiency scaling functions that determine how the reacting molecule’s vibrational state enhances the rate; and $ \widetilde {s} _f( \boldsymbol w )$ and $ \widetilde {s} _r( \boldsymbol v )$ are functions to be defined based on the detailed balance constraint (3.14). Substituting this form into (3.14) and rearranging gives

(3.16) \begin{align} \frac { \widetilde {s} _f( \boldsymbol w )}{\sum _{ \boldsymbol w } \widetilde {s} _f( \boldsymbol w )} = \frac { \widetilde {f} _B^{*}( \boldsymbol w ) \, \widetilde {\omega } _r( \boldsymbol w )} {\sum _{ \boldsymbol w } \widetilde {f} _B^{*}( \boldsymbol w ) \, \widetilde {\omega } _r( \boldsymbol w )}. \end{align}

A sufficient constraint to ensure this relation holds is to set $ \widetilde {s} _f( \boldsymbol w ) = \widetilde {f} _B^{*}( \boldsymbol w ) \, \widetilde {\omega } _r( \boldsymbol w )$ . Following the same logic with (3.13), but summing over $ \boldsymbol w$ instead of $ \boldsymbol v$ , results in a similar constraint for $ \widetilde {s} _r( \boldsymbol w )$ .

The subsequent state-averaged vibrational efficiency functions, with $a_f$ and $a_r$ defined to ensure $\varPhi _f = \varPhi _r = 1$ at thermal equilibrium, are given by

(3.17) \begin{align} \varPhi _f = \frac {\varOmega _f \, \varOmega _r^*}{\varOmega _f^{tr} \, \varOmega _r^{tr}}, \quad \varPhi _r = \frac {\varOmega _r \, \varOmega _f^*}{\varOmega _r^{tr} \, \varOmega _f^{tr}}, \end{align}

where

(3.18) \begin{align} &{\kern 4pt}\varOmega _f = \sum _{ \boldsymbol v } \widetilde {f} _A( \boldsymbol v ) \, \widetilde {\omega } _f( \boldsymbol v ), \quad \varOmega _f^* = \sum _{ \boldsymbol v } \widetilde {f} _A^{*}( \boldsymbol v ) \, \widetilde {\omega } _f( \boldsymbol v ), \quad \varOmega _f^{tr} = \sum _{ \boldsymbol v } \widetilde {f} _A^{tr}( \boldsymbol v ) \, \widetilde {\omega } _f^{tr}( \boldsymbol v ), \nonumber \\ &\varOmega _r = \sum _{ \boldsymbol w } \widetilde {f} _B( \boldsymbol w ) \, \widetilde {\omega } _r( \boldsymbol w ), \quad \varOmega _r^* = \sum _{ \boldsymbol w } \widetilde {f} _B^{*}( \boldsymbol w ) \, \widetilde {\omega } _r( \boldsymbol w ), \quad \varOmega _r^{tr} = \sum _{ \boldsymbol w } \widetilde {f} _B^{tr}( \boldsymbol w ) \, \widetilde {\omega } _r^{tr}( \boldsymbol w ), \end{align}

and where the $tr$ superscript denotes that the quantity is to be evaluated with all temperatures set to $T_{tr}$ . Any choice for $T_{*}$ is valid, but we choose to set $T_{*} = T_{tr}$ , which is similar to the convention established by Knab et al. (Reference Knab, Frühauf and Messerschmid1995). It follows that if $ \widetilde {\omega } _r( \boldsymbol w )$ is independent of $T_{v,i}$ , then $\varOmega _f^* = \varOmega _f^{tr}$ , and $\varPhi _f = \varOmega _f / \varOmega _f^{tr}$ is only dependent on the reactant state (with a similar result for the reverse reaction). In general, we do not enforce any specific temperature dependencies on $ \widetilde {\omega } _r( \boldsymbol w )$ , and instead look to define a generalised scaling function for a broad range of systems.

3.3.2. A generalised vibrational scaling function

The state-specific vibrational scaling function, $ \widetilde {\omega } _{f}( \boldsymbol v )$ , controls how vibrational energy enhances the reaction rate. The only specific guideline for defining this function is that high-lying states tend to increase the reaction rate, the effect of which lessens as $T_{tr}$ increases (Voelkel et al. Reference Voelkel, Raman and Varghese2016b ). With this in mind, and motivated by historical precedent (as discussed below) (Marrone & Treanor Reference Marrone and Treanor1963; Knab et al. Reference Knab, Frühauf and Messerschmid1995; Arsentiev et al. Reference Arsentiev, Loukhovitski and Starik2012; Voelkel et al. Reference Voelkel, Raman and Varghese2016b ), we choose a form given by

(3.19) \begin{align} \widetilde {\omega } _{f}( \boldsymbol v ) = \prod _{i \in A} \exp \bigg ( \frac { \widetilde {\varepsilon } _{v,i}(v_i)}{U_{i}} \bigg ) \,\boldsymbol{\cdot }\, \exp \bigg ( \frac { \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v ) }{ T_{tr} } \bigg ), \end{align}

where $ \widetilde {\varepsilon } _{v,i}(v_i)$ is the state-specific vibrational energy, normalised such that its units are temperature; $U_i$ is a model parameter that determines how the collisional cross-section of the species increases with increasing vibrational energy (Marrone & Treanor Reference Marrone and Treanor1963); and $ \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v )$ is the effective vibrational component of the reaction barrier. Again, in this section, we forgo reaction-specific notation for brevity, but recognise that each of these parameters may depend on both species and reaction. The effective vibrational component $ \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v )$ is defined to be a weighted sum given by

(3.20a) \begin{align}&{\kern 88pt} \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v ) = \sum _{ \boldsymbol w } \widetilde {f} _B( \boldsymbol w ) \, \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v , \boldsymbol w ), \end{align}

(3.20b) \begin{align} &\widetilde {\mathcal{E}} _{v,f}( \boldsymbol v , \boldsymbol w ) = \sum _{i \in A} \alpha _{v,i} \, \widetilde {\varepsilon } _{v,i}(v_i) - \delta ^{\prime } \max \Big [ 0, \, \widetilde {\varDelta } _{f}( \boldsymbol v ) - \delta ^{\prime \prime } \max \big ( 0, \, \widetilde {\varDelta } _{r}( \boldsymbol w ) \big ) \Big ], \end{align}

(3.20c) \begin{align} &\widetilde {\varDelta } _{f}( \boldsymbol v ) = \sum _{i \in A} \alpha _{v,i} \, \widetilde {\varepsilon } _{v,i}(v_i) - \beta _f \, \mathcal{E}_{a,f}, \, \widetilde {\varDelta } _{r}( \boldsymbol w ) = \sum _{i \in B} \alpha _{v,i} \, \widetilde {\varepsilon } _{v,i}(w_i) - \beta _r \, \mathcal{E}_{a,r}, \end{align}

where $\alpha _{v,i}$ is the vibrational efficiency parameter; $\beta _f$ and $\beta _r$ are parameters that scale forward and reverse activation energies, respectively; $\mathcal{E}_{a,f}$ and $\mathcal{E}_{a,r}$ are the forward and reverse activation energies in temperature units, respectively; and $\delta ^\prime$ and $\delta ^{\prime \prime }$ are constants equal to either zero or one. The effective vibrational component $ \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v )$ is defined so that vibrational energies above the forward activation energy are limited to the forward activation energy, which itself may increase for sufficiently large product vibrational energies greater than the reverse activation energy. For this generalised model, we assume that each of these parameters may have some temperature dependency. A corresponding form for the reverse efficiency scaling function $ \widetilde {\omega } _r( \boldsymbol w )$ can also be defined.

The motivation for the effective vibrational component of the reaction barrier, $ \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v )$ , is based on a simplified reaction pathway illustrated in figure 1. Here, we show how a reference reaction path (black line), defined assuming the vibrational contributions of reactant and product species are zero, is modified by the reactant and product vibrational energies (red line). In short, the reactant vibrational energy reduces the translational–rotational energy needed for reaction; as such, higher vibrational energy in the reactants leads to increased reaction rates. Here, we define the magnitude of that increase by $\sum _{i\in A} \alpha _{v,i} \widetilde {\varepsilon } _{v,i}(v_i)$ , limited by the reference activation energy barrier as defined in (3.20). The vibrational energy efficiency parameter, $\alpha _{v,i}$ , controls how effective the vibrational energy is at contributing to the activation energy. Here, $\alpha _{v,i}$ is intended to be a representative value associated with each chemical species, independent of its complexity (e.g. number of vibrational modes). In our analysis, because we have assumed that intramolecular vibrational relaxation is instantaneous (resulting in one vibrational temperature per species), we observe that a set of vibrational mode-specific efficiency parameters does not materially alter the modelling flexibility granted by assigning a single, representative value. The activation energies $\mathcal{E}_{a,f}$ and $\mathcal{E}_{a,r}$ are scaled by $\beta _f$ and $\beta _r$ , respectively, to account for the notion that the reference path assumes no contributions from vibrational energy, though $\mathcal{E}_{a,f}$ and $\mathcal{E}_{a,r}$ are determined independently of this assumption.

Figure 1.

A simplified reaction path diagram. The black line illustrates a reference reaction pathway assuming negligible vibrational energy contributions, and the red line shows how the path shifts as the reactant and product vibrational energies increase.

Several historical models can be reconstructed using this general model. Setting parameters so that $ \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v ) = \widetilde {\mathcal{E}} _{v,r}( \boldsymbol w ) = 0$ recovers the model presented by Marrone & Treanor (Reference Marrone and Treanor1963). Under this model, the efficiency function only varies with respect to $U_i$ , so the increased cross-section of reactants with high-lying vibrational states is accounted for, but the vibrational energy has no effect on the reaction barrier. Setting $\delta ^{\prime \prime } = 0$ recovers the CVCV model (Knab et al. Reference Knab, Frühauf and Messerschmid1995), though the original model did not include $\alpha _{v,i}$ (the $\alpha$ parameter from the original paper is replaced by $\beta _f$ in this formulation). Under this model, the reverse activation energy does not contribute to the modified reaction barrier. This neglects products with high-lying vibrational states greater than the reverse reaction barrier, thus potentially cutting off vibrational states that could enhance the rate. Setting $U_i = \infty$ recovers a model motivated by Arsentiev et al. (Reference Arsentiev, Loukhovitski and Starik2012). However, the original model did not impose an assumption that the efficiency function is separable into reactant and product states, resulting in a non-separable state-specific efficiency function, with $ \widetilde {\varPhi } _f( \boldsymbol v , \boldsymbol w ) \sim \exp ( \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v , \boldsymbol w )/T_{tr})$ . For general reactions, it is unclear how the approach by Arsentiev et al. (Reference Arsentiev, Loukhovitski and Starik2012) affects the detailed balance of reversible reactions. Finally, setting $U_i = \infty$ and $\alpha _{v,i} = 0.6 - T_{tr}/10^4$ recovers a model similar to that presented by Voelkel et al. (Reference Voelkel, Raman and Varghese2016b ), though this model had a slightly different approach for handling the reference energy barrier.

3.4. Chemical–vibrational energy exchange

Chemical–vibrational energy exchange represents the creation and destruction of vibrational energy from chemical reactions. For a state-specific rate as defined in § 3.3.1, the destruction of energy at a particular vibrational state is determined by the forward rate multiplied by the energy associated with that state (Knab et al. Reference Knab, Frühauf and Messerschmid1995). A similar definition applies for the creation of energy as determined via the reverse rate. A weighted average over the state-specific creation and destruction of energy results in the bulk C–V exchange. Thus, for a given reaction, $j$ , the C–V exchange rate for a reactant species is given by

(3.21) \begin{align} \frac { S_{{\scriptsize C-V},ij} }{W_i} = c_{B,j} \sum _{ \boldsymbol v , \boldsymbol w } \widetilde {f} _B( \boldsymbol w ) \, \widetilde {k} _{r,j}( \boldsymbol v , \boldsymbol w ) \, \widetilde {e} _{v,i}(v_i) - c_{A,j} \sum _{ \boldsymbol v , \boldsymbol w } \widetilde {f} _A( \boldsymbol v ) \, \widetilde {k} _{f,j}( \boldsymbol v , \boldsymbol w ) \, \widetilde {e} _{v,i}(v_i), \end{align}

where $ \widetilde {e} _{v,i}(v_i) = R_i \, \widetilde {\varepsilon } _{v,i}(v_i)$ is the state-specific vibrational energy per unit mass. A similar relation holds for product species. Accounting for generalised reactant and product combinations, and multiple reactions, the resulting state-averaged C–V exchange rate becomes

(3.22) \begin{align} & {\kern 114pt}S_{{\scriptsize C-V},i} = \sum _{j=1}^J S_{{\scriptsize C-V},ij} , \nonumber \\ & S_{{\scriptsize C-V},ij} = W_i \varLambda _{r,j} \left ( \nu _{ij}^\prime \, G_{f,ij}^* - \nu _{ij}^{\prime \prime } \, G_{r,ij} \right ) - W_i \varLambda _{f,j} \left ( \nu _{ij}^\prime \, G_{f,ij} - \nu _{ij}^{\prime \prime } \, G_{r,ij}^* \right ), \end{align}

where

(3.23) \begin{align} &{\kern 2pt} G_{f,ij} = \sum _{ \boldsymbol v } \widetilde {f} _A( \boldsymbol v ) \, \widetilde {\omega } _{f,j}( \boldsymbol v ) \, \widetilde {e} _{v,i}(v_i) / \varOmega _{f,j}, \quad G_{f,ij}^* = \sum _{ \boldsymbol v } \widetilde {f} _A^{*}( \boldsymbol v ) \, \widetilde {\omega } _{f,j}( \boldsymbol v ) \, \widetilde {e} _{v,i}(v_i) / \varOmega _{f,j}^*, \nonumber \\ &G_{r,ij} = \sum _{ \boldsymbol w } \widetilde {f} _B( \boldsymbol w ) \, \widetilde {\omega } _{r,j}( \boldsymbol w ) \, \widetilde {e} _{v,i}(w_i) / \varOmega _{r,j}, \quad G_{r,ij}^* = \sum _{ \boldsymbol w } \widetilde {f} _B^{*}( \boldsymbol w ) \, \widetilde {\omega } _{r,j}( \boldsymbol w ) \, \widetilde {e} _{v,i}(w_i) / \varOmega _{r,j}^*. \end{align}

Physically, $G_{f,ij}^*$ and $G_{r,ij}^*$ quantify the production of vibrational energy from reactions, whereas $G_{f,ij}$ and $G_{r,ij}$ quantify the destruction of vibrational energy.

4.1. Vibrational partitioning of the internal energy

For high-speed combustion applications, complex chemical mechanisms are typically developed under the assumption of a specific equation of state. For instance, the hydrogen–air chemical mechanism presented by Burke et al. (Reference Burke, Chaos, Ju, Dryer and Klippenstein2012) supports both thermodynamic state and molar production rate calculations via Chemkin-II (Kee et al. Reference Kee, Rupley and Miller1989). The associated equation of state is the NASA 7-coefficient polynomial model. Modifying the equation of state to, say, a rigid-rotor harmonic oscillator model would necessitate recalibration of the kinetic mechanism. With this in mind, we will assume that the internal energy (at thermal equilibrium) is defined via some functional form, denoted by $e_{i}^{{\scriptsize fit}}$ . Because we will eventually need the vibrational partition function for each species to define the vibrational efficiency function (see (4.7) below), we choose to define vibrational energy explicitly via the harmonic oscillator approximation and calculate the translational–rotational energy implicitly. Thus, the species internal energy as defined in (2.5) has translational–rotational and vibrational components given by

(4.2) \begin{align} e_{tr,i}(T_{tr}) = e_{i}^{{\scriptsize fit}}(T_{tr}) - e_{v,i}(T_{tr}), \quad e_{v,i}(T_{v,i}) = \sum _{j=1}^{N_{v,i}} \frac {R_i \theta _{v,ij}}{\textrm{e}^{\theta _{v,ij}/T_{v,i}} - 1}, \end{align}

where $N_{v,i}$ is the vibrational degrees of freedom and $\theta _{v,ij}$ is the characteristic vibrational temperature corresponding to the species’ $j{\rm th}$ vibrational mode. For a species with $N_{a,i}$ atoms, $N_{v,i} = 0$ if the species is monatomic ( $N_{a,i}=1$ ), $N_{v,i} = 3 N_{a,i} - 5$ if the species is linear and $N_{v,i} = 3 N_{a,i} - 6$ if the species is nonlinear.

Several other relatively simple options could be used to define the vibrational energy. First, we could employ an anharmonic oscillator model or truncated harmonic oscillator model and replace $e_{v,i}$ in (4.2) accordingly. Section 5.1 explores both of these choices with regard to the vibrational efficiency function. Second, for complex molecules, where characteristic vibrational temperatures are not fully defined, it may be preferred to define the translational–rotational–electronic energy explicitly (via a rigid-rotor approximation, for instance) and calculate the vibrational energy implicitly from $e_{i}^{{\scriptsize fit}}$ . Under this formulation, the vibrational partition function is calculated by integrating (2.4). Because we are calculating the vibrational energy implicitly, extra care must be taken to ensure that the vibrational energy and partition function are physically consistent (e.g. $e_{v,i} \gt 0$ , ${\rm d}e_{v,i}/{\rm d}T_{v,i} \gt 0$ , ${\rm d}Q_{v,i}/{\rm d}T_{v,i} \gt 0$ , etc.).

4.2. Translational–vibrational exchange

To fully specify T–V exchange, the pairwise T–V relaxation time scale, $ \tau _{{\scriptsize T-V},ij}$ , must be defined. Ideally, the time scale for each species’ pairwise combination could be individually calculated or measured. The literature provides an abundance of such studies, though they are typically focused on a relatively limited set of species’ pairwise interactions. For a wide range of diatomic molecules, Millikan & White (Reference Millikan and White1963) showed that the following correlation generally agrees to within 50 % of the empirical measurements:

(4.3) \begin{align} p \boldsymbol{\cdot }\tau _{{\scriptsize T-V},ij} = \exp \left [ A W_{ij}^{1/2} \theta _{v,i}^{4/3} \left ( T_{tr}^{-1/3} - B W_{ij}^{1/4} \right ) - C\right ] , \end{align}

where $A = 1.16 \times 10^{-3}$ , $B = 1.50 \times 10^{-2}$ and $C = 18.42$ are constants determined via calibration; and $W_{ij}$ is the reduced mass of the colliding pair. The units of $A$ , $B$ and $C$ are defined assuming $W_{ij}$ is in g mol⁻¹ and $T_{tr}$ and $\theta _{v,ij}$ are in K. The corresponding T–V exchange probability is given by $ P_{{\scriptsize T-V},ij} = \tau _{c,ij}/ \tau _{{\scriptsize T-V},ij}$ , where $\tau _{c,ij}$ is defined in (3.2). Many other choices for the T–V relaxation time scale have been proposed, including corrections for high temperatures (Park Reference Park1989). For temperatures below 10 000 K, these high-temperature corrections have negligible impact on the relaxation time scale (see figure 7.3 in Boyd & Schwartzentruber (Reference Boyd and Schwartzentruber2017)). The Millikan & White model is included here and used in the subsequent analysis due to its historical significance and to provide a general understanding of the approximate time scales in detonations.

4.3. Vibrational–vibrational exchange

The key parameters in defining V–V exchange are the collision cross-section, $\sigma _{ij}$ , which controls the collision time scale via (3.2), and the exchange probability, $ P_{{\scriptsize V-V},ij}$ . The collision cross-section is a measure of area that describes the propensity of two molecules to affect one another’s trajectory via a collision. The value for $\sigma _{ij}$ can be calculated via QCT methods, which may employ a potential energy surface determined from quantum chemistry computations, or from some intermolecular potential, such as a hard-sphere potential, inverse power-law potential or Lennard-Jones potential (Boyd & Schwartzentruber Reference Boyd and Schwartzentruber2017). The hard-sphere model represents the simplest such model, with $\sigma _{ij} = \pi d_{ij}^2$ , where $d_{ij}$ is the average diameter of the species pair. The challenge remains to determine this ‘diameter’ for each species. In general, the hard-sphere model can be used to calculate transport properties such as viscosity, which is approximately $5 (2 \pi W_i \mathcal{R}_u T)^{1/2}/(16 \pi \mathcal{N}_A d_i)$ , where $d_i$ is the species diameter (Boyd & Schwartzentruber Reference Boyd and Schwartzentruber2017). With viscosity measurements from Lemmon et al. (Reference Lemmon, Bell, Huber and McLinden2024), such a relation allows us to approximate the diameter for a range of molecules. For hydrogen–air mixtures below 5000 K, the calculated hard-sphere diameters range between $2$ and $4$ Å. Due to the complex nature of calculating the collision cross-section and the limited accuracy of the hard-sphere approximation itself, we set $d_{ij}$ to an average value of $3$ Å for each collision pair.

With regard to the exchange probability, $ P_{{\scriptsize V-V},ij}$ is typically set as a constant value depending on the system (Candler & MacCormack Reference Candler and MacCormack1988, Reference Candler and MacCormack1991; Knab et al. Reference Knab, Frühauf and Jonas1992, Reference Knab, Frühauf and Messerschmid1995; Voelkel et al. Reference Voelkel, Masselot, Varghese and Raman2016a ). However, such an assumption is inconsistent with the thermal equilibrium constraint described by (3.4) because $ P_{{\scriptsize V-V},ij}$ is not necessarily symmetric. Consider that if we assume $ P_{{\scriptsize V-V},ij}$ is symmetric, then (3.4) is valid if and only if e $_{v,i}^{tr} = e_{v,j}^{tr}$ , which itself is only valid under rare conditions. By invalidating (3.4), then $\sum _i S_{{\scriptsize V-V},i} /W_i$ , with $ S_{{\scriptsize V-V},i}$ defined in (3.5), is not necessarily zero. However, $\sum _i S_{{\scriptsize V-V},i} /W_i$ should be zero based on the base V–V exchange model as defined in (3.3). Thus, assuming $ P_{{\scriptsize V-V},ij}$ is symmetric inadvertently modifies how V–V exchange affects the translational–rotational energy, a process that we refer to as ‘leakage’ and one that we wish to avoid.

To rectify this inconsistency, we set $ P_{{\scriptsize V-V},ij}$ as follows:

(4.4) \begin{align} P_{{\scriptsize V-V},ij} = P_{{\scriptsize V-V},ij}^\prime \left ( \frac {(e_{v,i}^{tr} + e_{v,j}^{tr})/2}{e_{v,i}^{tr}} \right ), \end{align}

where $ P_{{\scriptsize V-V},ij}^\prime$ is a modified (symmetric) probability factor associated with species pairwise vibrational transfer. Setting $ P_{{\scriptsize V-V},ij}^\prime$ to a constant value is now consistent with (3.4) and $ S_{{\scriptsize V-V},ij} /W_i = - S_{{\scriptsize V-V},ji} /W_j$ , thus ensuring that $\sum _i S_{{\scriptsize V-V},i} /W_i = 0$ . Our choice for $ P_{{\scriptsize V-V},ij}$ is motivated by (3.4), wherein $ P_{{\scriptsize V-V},ij} e_{v,i}^{tr}$ and $ P_{{\scriptsize V-V},ji} e_{v,j}^{tr}$ are expected to increase monotonically with $T_{tr}$ . With $ P_{{\scriptsize V-V},ij}$ set to be proportional to the average species pairwise vibrational energy at thermal equilibrium, the left- and right-hand sides of (3.4) still increase monotonically with $ P_{{\scriptsize V-V},ij}^\prime$ set to a constant value.

In general, we expect $ P_{{\scriptsize V-V},ij}^\prime$ to be similar in magnitude to $ P_{{\scriptsize V-V},ij}$ and $ P_{{\scriptsize V-V},ji}$ , so the previous historical reference values are still valid. Originally, $ P_{{\scriptsize V-V},ij}$ was set to a constant value of $10^{-2}$ (Candler & MacCormack Reference Candler and MacCormack1988, Reference Candler and MacCormack1991; Knab et al. Reference Knab, Frühauf and Jonas1992, Reference Knab, Frühauf and Messerschmid1995; Fiévet et al. Reference Fiévet, Voelkel, Koo, Raman and Varghese2017). This value seemed to be sufficient for high-speed reentry flows of interest in those studies. However, in measuring the vibrational–rotational relaxation process of $\textrm{N}_{2}$ and $\textrm{O}_{2}$ within high-speed turbulent diffusion flames, Reising et al. (Reference Reising, Haller, Clemens, Varghese, Fiévet and Raman2016) observed that a value closer to $10^{-5}$ provided the best match between simulations and experimental measurements. It is also well established that $ P_{{\scriptsize V-V},ij}$ is not only species-dependent, but also temperature-dependent. For instance, Taylor et al. (Reference Taylor, Camac and Feinberg1967) measured V–V exchange probabilities between $10^{-4}$ and $10^{-3}$ for $\textrm{NO}-\textrm{CO}$ , $\textrm{NO}-\textrm{N}_{2}$ and $\textrm{CO}_{2}-\textrm{N}_{2}$ mixtures between 500 and 2000 K. Schwartz et al. (Reference Schwartz, Slawsky and Herzfeld1952) predicted values between approximately $10^{-5}$ and $10^{-4}$ for $\textrm{N}_{2}-\textrm{O}_{2}$ collisions. Voelkel (Reference Voelkel2016) directly calculated $\textrm{N}_{2}-\textrm{O}_{2}$ scattering rates via QCT and determined that the V–V exchange probability varies between $10^{-6}$ and $10^{-4}$ for temperatures between 1000 and 3000 K. These collective observations were used to motivate a change from $ P_{{\scriptsize V-V},ij} = 10^{-2}$ to $10^{-4}$ for the numerical study by Voelkel et al. (Reference Voelkel, Masselot, Varghese and Raman2016a ).

This breadth of studies demonstrates that the V–V exchange represents a complex, species-specific and temperature-dependent phenomenon, and that a robust model for the exchange probability is still a matter of considerable research. With this in mind, we explore a range of constant values for $ P_{{\scriptsize V-V},ij}^\prime$ to provide a broad view of how V–V exchange works and couples with other vibrational exchange processes, and we look towards potential future studies for the development of a generalised, species- and temperature-dependent exchange probability model.

4.4. Vibrationally informed chemical reactions

For modelling vibrationally-informed chemical reactions, it remains to define the model parameters associated with our form of the vibrational scaling function (3.19), which in turn defines the vibrational efficiency function (3.17). Under high-speed combustion conditions, where temperatures are typically below 5000 K, we assume vibrational energies exceeding the reaction activation energy are negligible in the flow. This is equivalent to setting $\delta ^\prime = 0$ in (3.20). The resulting effective vibrational component of the reaction barrier reduces to $ \widetilde {\mathcal{E}} _{v,f}( \boldsymbol v ) = \sum _{i \in A} \alpha _{v,i} \, \widetilde {\varepsilon } _{v,i}(v_i)$ , and the state-specific vibrational scaling function as defined in (3.19) becomes

(4.5) \begin{align} \widetilde {\omega } _{f}( \boldsymbol v ) = \prod _{i \in A} \exp \bigg [ \widetilde {\varepsilon } _{v,i}(v_i) \left ( \frac { \alpha _{v,i} }{ T_{tr} } + \frac {1}{U_{i}} \right ) \bigg ]. \end{align}

Recognising that $ \widetilde {f} _{A}( \boldsymbol v ) = \prod _{i \in A} \widetilde {f} _{v,i}(v_i)$ and $ \widetilde {f} _{v,i}(v_i) = \prod _{i \in A} \exp [- \widetilde {\varepsilon } _{v,i}(v_i)/T_{v,i} ]/{}Q(T_{v,i})$ , where $Q_{v,i}$ is the vibrational partition function as described in § 2.1, then $\varOmega _f$ and $\varOmega _f^{tr}$ as defined in (3.18) reduce to

(4.6) \begin{align} \varOmega _f = \prod _{i \in A} \frac {Q_{v,i}(\varGamma _i)}{Q_{v,i}(T_{v,i})}, \quad \varOmega _f^{tr} = \prod _{i \in A} \frac {Q_{v,i}(\varGamma _i^{tr})}{Q_{v,i}(T_{tr})}, \end{align}

where $1/\varGamma _{i} = 1/T_{v,i} - \alpha _{v,i}/T_{tr} - 1/U_i$ and $1/\varGamma _{i}^{tr} = (1 - \alpha _{v,i})/T_{tr} - 1/U_i$ are effective ‘gamma’ temperatures. Similar reductions exist for $\varOmega _r$ and $\varOmega _r^{tr}$ . The resulting efficiency function is given by

(4.7) \begin{align} \varPhi _{f,j} = \prod _{i=1}^N \varphi _{i}^{\nu _{ij}^{\prime }}, \quad \varPhi _{r,j} = \prod _{i=1}^N \varphi _{i}^{\nu _{ij}^{\prime \prime }}, \quad \varphi _{i} = \frac { Q_{v,i}(T_{tr}) \, Q_{v,i}(\varGamma _i) }{ Q_{v,i}(T_{v,i}) \, Q_{v,i}(\varGamma _{i}^{tr})}, \end{align}

where $\varphi _{i}$ is the species efficiency function. Here, we have assumed $\alpha _{v,i}$ and $U_i$ do not vary between reactions, resulting in a reaction-independent species vibrational efficiency function. The effective temperatures $\varGamma _i$ and $\varGamma _i^{tr}$ modulate the magnitude of the vibrational efficiency function’s scaling behaviour depending on the value of $\alpha _{v,i}$ and $U_i$ . Note that this model does not explicitly require a specific vibrational energy description, thus allowing for the use of harmonic, truncated harmonic and anharmonic oscillator models, as well as energy levels measured empirically or calculated directly from a potential energy surface, as long as the vibrational partition function is well defined.

For this vibrational efficiency function, we seek a model with the following physically consistent behaviour. For a given species denoted by $i$ , with $T_{tr}$ and $T_{v,k\ne i}$ fixed, a vibrationally cold state ( $T_{v,i} \lt T_{tr}$ ) has under-excited high-lying vibrational modes and should reduce the reaction rate, whereas a vibrationally hot state ( $T_{v,i} \gt T_{tr}$ ) has over-excited high-lying vibrational modes and should increase the reaction rate. Furthermore, the farther $T_{v,i}$ is from $T_{tr}$ , the more pronounced this effect will become. To ensure this behaviour, we simply require that $\varPhi _{f,j}$ increases monotonically with $T_{v,i}$ ( $\partial \varPhi _{f,j}/\partial T_{v,i} \ge 0$ ). We assume that the $Q_{v,i}$ is positive for physical (positive) temperatures, so it is sufficient to show that $\partial \varphi _{i}/\partial T_{v,i} \ge 0$ to satisfy the monotonicity condition.

This condition constrains the vibrational efficiency parameter, $\alpha _{v,i}$ . For our implementation, which utilises the harmonic oscillator approximation, a sufficient constraint requires $\alpha _{v,i} \ge -T_{tr}/U_i$ , as well as $\partial \alpha _{v,i}/\partial T_{v,i} \ge 0$ when $T_{v,i} \lt T_{tr}$ and $\partial \alpha _{v,i}/\partial T_{v,i} \le 0$ when $T_{v,i} \gt T_{tr}$ . Note that this is a sufficient constraint, with more relaxed constraints requiring a balance of $\partial \alpha _{v,i}/\partial T_{v,i}$ against a series of partition functions. This, in conjunction with a constraint requiring non-negative but finite gamma temperatures to ensure $ 0 \le \varPhi _{f,j} \lt \infty$ , results in the following constraints on the vibrational efficiency parameter:

(4.8) \begin{align} -\frac {T_{tr}}{U_i} \le \alpha _{v,i} \lt \frac {1}{\max (\,1, \,T_{v,i}/T_{tr} \,)} - \frac {T_{tr}}{U_i}. \end{align}

These bounds represent the minimum and maximum extent to which vibrational non-equilibrium can influence reaction kinetics. When $\alpha _{v,i} = -T_{tr}/U_i$ , then $\varGamma _{i} = T_{v,i}$ and $\varPhi _{f,j} = 1$ , indicating no vibrational non-equilibrium effect. When $\alpha _{v,i} = 1/\max (1,T_{v,i}/T_{tr}) - T_{tr}/U_i - \delta$ , where $\delta$ is a small perturbation, $\varGamma _i^{tr} = T_{tr}/\delta$ when $T_{v,i} \lt T_{tr}$ and the efficiency function goes to zero, while $\varGamma _i = T_{tr}/\delta$ when $T_{v,i} \gt T_{tr}$ and the efficiency function goes to infinity. Given these bounds, it is convenient to define $\alpha _{v,i}$ as some fractional value of the minimal and maximal vibrational effect, i.e.

(4.9) \begin{align} \alpha _{v,i} = \frac {\alpha _{v,i}^\prime }{\max (\,1, \,T_{v,i}/T_{tr} \,)} - \frac {T_{tr}}{U_i}, \end{align}

where $\alpha _{v,i}^\prime \in [0,1)$ is the modified vibrational efficiency parameter. The resulting gamma temperatures are now independent of $U_i$ , and only $\alpha _{v,i}^\prime$ needs to be defined.

As a final note, we recognise that introducing a vibrational temperature into the definition for $\alpha _{v,i}$ should be addressed. Rather than evaluating $\alpha _{v,i}$ at $T_{tr}$ when setting $\varGamma _i^{tr}$ , we choose to maintain the vibrational temperature dependency. This choice does not alter the underlying formulation of $\varPhi _f$ and $\varPhi _r$ in (3.17) because the equilibrium limit is still recovered. Consider the alternative, whereby $\varGamma _i^{tr}$ is evaluated using $\alpha _{v,i}^{tr}$ ( $\alpha _{v,i}$ evaluated with $T_{v,i}$ set to $T_{tr}$ ). The subsequent constraints (assuming $U_i \rightarrow \infty$ for brevity) to ensure $\partial \varphi _{i}/\partial T_{v,i} \ge 0$ and positive gamma temperatures are $0 \le \alpha _{v,i} \lt T_{tr}/T_{v,i}$ , $0 \le \alpha _{v,i}^{tr} \lt 1$ and $\partial \alpha _{v,i}/\partial T_{v,i} \ge -T_{tr}/T_{v,i}^2\boldsymbol{\cdot }Q(T_{v,i})/Q(T_{tr}/\alpha _{v,i})$ (or simply that $\partial \alpha _{v,i}/\partial T_{v,i} \ge 0$ ). It is straightforward to show that for vibrationally hot flows, $\alpha _{v,i}$ may potentially have to decrease to ensure $\alpha _{v,i} \lt T_{tr}/T_{v,i}$ , otherwise $\varGamma _i$ will be negative. However, as $\alpha _{v,i}$ decreases, the constraint for $\partial \alpha _{v,i}/\partial T_{v,i}$ will eventually be invalidated, resulting in $\partial \varphi _{i}/\partial T_{v,i} \lt 0$ (a rate that decreases as $T_{v,i}$ increases). It is this restriction that motivates our choice to simply evaluate $\varGamma _i^{tr}$ at $\alpha _{v,i}$ (and maintain the vibrational temperature dependency).

4.5. Chemical–vibrational exchange

For the vibrational efficiency model described in § 4.4, the resulting species- and reaction-dependent C–V exchange rate as defined in (3.22) reduces to

(4.10) \begin{align} S_{{\scriptsize C-V},ij} = \rho Y_i e_{v,i}^{tr} \left ( \frac {e_{v,i}(\varGamma _i^{tr})/e_{v,i}^{tr}}{ \tau _{{\scriptsize C-V(c)},ij} } - \frac {e_{v,i}(\varGamma _i)/e_{v,i}^{tr}}{ \tau _{{\scriptsize C-V(d)},ij} } \right ), \end{align}

where $ \tau _{{\scriptsize C-V(c)},ij} = \rho Y_i/[W_i(\nu _{ij}^{\prime } \varLambda _{r,j} +\nu _{ij}^{\prime \prime } \varLambda _{f,j})]$ and $ \tau _{{\scriptsize C-V(d)},ij} = \rho Y_i/[W_i (\nu _{ij}^{\prime } \varLambda _{f,j} + \nu _{ij}^{\prime \prime } \varLambda _{r,j})]$ are the C–V time scales for the creation (c) and destruction (d) of vibrational energy from chemical reactions. Here, we include the $\rho Y_i e_{v,i}^{tr}$ term to have C–V exchange in a similar form compared with T–V exchange via (3.1) and V–V exchange via (3.5). The corresponding C–V exchange time scale is proportional to the reaction time scales associated with the flow.

4.6. Guidance on selecting parameter values

Our proposed vibrational non-equilibrium model for high-speed combustion requires two sets of model parameters to close the system: $ P_{{\scriptsize V-V},ij}^\prime \gt 0$ and $\alpha _{v,i}^\prime \in [0,1)$ . For a typical hydrogen–air mixture with seven non-monatomic species, six of which are part of the reaction mechanism, then 21 $ P_{{\scriptsize V-V},ij}^\prime$ and 6 $\alpha _{v,i}^\prime$ values (or closures) are needed. Consider also that $ P_{{\scriptsize V-V},ij}^\prime$ generally depends on $T_{tr}$ . For instance, in V–V exchange collisions of $\textrm{N}_{2}$ and $\textrm{O}_{2}$ , Voelkel (Reference Voelkel2016) found that for any fixed $T_{tr}$ (below 3000 K), $ P_{{\scriptsize V-V},ij}^\prime$ was approximately constant. However, for $T_{tr}$ varied between 500 and 3000 K, the value of $ P_{{\scriptsize V-V},ij}^\prime$ varied by approximately three orders of magnitude. These results suggest that a temperature dependence is likely warranted for predicting accurate relaxation phenomena. The same likely applies to $\alpha _{v,i}^\prime$ (though using a different model, consider the discussion by Voelkel et al. (Reference Voelkel, Masselot, Varghese and Raman2016a )). It is also possible that a global model or correlation, akin to the T–V exchange time scale presented by Millikan & White (Reference Millikan and White1963), could significantly reduce the number of new parameters needed to close the system by relying on basic molecular properties of each species and colliding pair.

All of this points to the need for a modern reassessment of the V–V exchange mechanism. Though beyond the scope of this work, we believe that collating and systematically examining historical V–V exchange studies could potentially reveal a global correlation valid within the high-speed combustion regime. However, we also recognise that legacy studies only offer a patchwork of rates, probabilities and behaviours associated with specific binary mixtures. Thus, part of this work would also necessarily entail identification of data gaps that need to be addressed, through either state-specific modelling efforts or new relaxation experiments.

Once the relaxation mechanism is determined, then $\alpha _{v,i}^\prime$ can be addressed. Here, we would recommend two paths: (i) calculate state-specific reaction rates from studies motivated by first principles analysis, such as QCT-based methods, and (ii) calibrate $\alpha _{v,i}^\prime$ to recover empirical measurements associated with combustion, including the induction zone length and unsteady behaviour (cellular structure) of detonation waves. The first approach will provide valuable insights into the nature of $\alpha _{v,i}^\prime$ , specifically regarding the assumptions of reaction independence and potential dependencies on $T_{tr}$ . The second approach is more practical and recognises that complex chemical mechanisms are highly coupled systems and stem from a wide range of sources. Not every reaction can leverage QCT-based methods to evaluate $\alpha _{v,i}^\prime$ . Instead, we can set $\alpha _{v,i}^\prime$ based on a smaller subset of studies to achieve specific behaviours, and then validate that choice by ensuring that the model accurately predicts specific empirical measurements and behaviours. Essentially, we treat $\alpha _{v,i}^\prime$ as a reduced parametrisation and calibrate its value to recover important application-specific features, similar to other reduced modelling efforts (e.g. the chemical–diffusive model) (Lu et al. Reference Lu, Kaplan and Oran2021).

Another approach, consistent with historical non-equilibrium studies (Candler & MacCormack Reference Candler and MacCormack1988, Reference Candler and MacCormack1991; Knab et al. Reference Knab, Frühauf and Messerschmid1995; Voelkel et al. Reference Voelkel, Masselot, Varghese and Raman2016a ), would be to further assume that $ P_{{\scriptsize V-V},ij}^\prime$ and $\alpha _{v,i}^\prime$ can be set to a constant value representing an average effect, thus necessitating only two parameter values to close the system. This is the choice we make in Part 2 (Voelkel et al. Reference Voelkel, Short and Chiquete2026) of this work to study the role of non-equilibrium within inert shocks and steady detonation waves of hydrogen–air mixtures. Therein, $ P_{{\scriptsize V-V},ij}^\prime$ is set to a constant value, varied between $10^{-6}$ and $10^{-2}$ . Over this range, the relaxation dynamics shifts completely from no V–V effect (entirely dominated by T–V exchange) to being completely determined by the fastest-relaxing species in the mixture (species are ‘pulled’ towards equilibrium by the species with the greatest T–V exchange). A value of $10^{-4}$ is shown to provide a good balance of expected physical observations, namely that $\textrm{N}_{2}$ relaxation is significantly enhanced by $\textrm{H}_{2}\textrm{O}$ , but this effect is lessened when $\textrm{N}_{2}$ is mixed with $\textrm{O}_{2}$ .

Likewise, $\alpha _{v,i}^\prime$ is varied between a value of 0.0 and 0.8 in Part 2 (Voelkel et al. Reference Voelkel, Short and Chiquete2026), beyond which the induction zone length of steady detonations increases by a factor of more than 40 % (and increases exponentially as $\alpha _{v,i}^\prime$ grows further). A value of 0.0 reduces the induction zone by nearly 20 %. In recommending a value for future studies, we first recognise that detailed chemical mechanisms are derived from shock tube and detonation experiments that can neither completely isolate nor completely remove vibrational non-equilibrium effects. Thus, non-equilibrium effects are inherently present in such measurements, which in turn are used to calibrate models that traditionally assume vibrational equilibrium. This introduces a physical inconsistency, in that an equilibrium model is calibrated to a flow field with non-equilibrium effects. Even so, the equilibrium model is calibrated to the measured data. It follows that the non-equilibrium results should still recover the empirical measurement, specifically the induction zone length. Ideally, new chemical mechanisms could be calibrated without imposing vibrational equilibrium, thus rectifying this inconsistency. Until then, we postulate that the non-equilibrium model should accurately reproduce the established empirical measurements, which implicitly account for vibrational non-equilibrium flow physics. As such, we show that $\alpha _{v,i}^\prime$ between 0.4 and 0.6 represents an optimal choice.

As a final note regarding this recommendation for $\alpha _{v,i}^\prime$ , we assert that recovering the same steady detonation behaviour between equilibrium and non-equilibrium models does not imply that they will recover the same unsteady detonation behaviour. In particular, it is well established that equilibrium models have trouble accurately predicting detonation cell sizes (Taylor et al. Reference Taylor, Kessler, Gamezo and Oran2013a ; Voelkel et al. Reference Voelkel, Masselot, Varghese and Raman2016a ; Lu et al. Reference Lu, Kaplan and Oran2021). However, by introducing non-equilibrium into the model, it is straightforward to show that the signal speed associated with acoustic perturbations changes form (Chu Reference Chu1957; Clarke Reference Clarke1978). Essentially, the frozen sound speed is now driven by the ratio of specific heats associated with only translational–rotational motion (the equilibrium sound speed still includes the vibrational motion). Depending on the relaxation time scales, and considering that the ratio of specific heats used to define the signal speed is critical for determining the unstable behaviour of the flow (Erpenbeck Reference Erpenbeck1964; Short & Stewart Reference Short and Stewart1998), we expect the non-equilibrium processes to dramatically alter the unsteady behaviour.

5.1. Vibrational efficiency model features

With the $\alpha _{v,i}$ defined via (4.9), the general features of the vibrational efficiency function are now examined. The modified vibrational efficiency parameter, $\alpha _{v,i}^\prime$ , does not have an established value, except that it is bound between 0 and 1. Thus, we explore its relative effect by examining a range of values. Figure 2 shows the species vibrational efficiency function for common molecules in hydrogen–air combustion with $\alpha _{v,i}^\prime$ varied between 0.3 and 0.7. For each species, we evaluate $Q_{v,i}$ via the harmonic oscillator approximation. As expected, the efficiency function monotonically increases with $T_{v,\textrm{i}}$ for fixed $T_{tr}$ . Increasing $\alpha _{v,i}^\prime$ pushes the efficiency function farther from unity, amplifying the vibrational non-equilibrium effect, consistent with the reaction diagram shown in figure 1 and our previous discussion (§ 3.3.2). Additionally, molecules with more vibrational degrees of freedom exhibit a greater effect on the vibrational efficiency. This is consistent with the notion that molecules with a greater potential to absorb or expel vibrational energy in molecular collisions will have a greater impact on the reaction rate. For molecules with the same number of vibrational modes, we observe that the vibrational efficiency function scales with the molecule’s characteristic vibrational temperatures.

Figure 2.

Species vibrational efficiency function versus $T_{v,i}/T_{tr}$ for various molecules in hydrogen combustion, with $T_{tr} = 1500$ K.

The disparity between molecules with the same number of vibrational modes diminishes as $T_{tr}$ increases, as shown in figure 3, which again shows the species vibrational efficiency function for common molecules in hydrogen–air combustion, now with $T_{tr}$ varied between 1000 and 5000 K. This observation arises from the nature of the vibrational partition function for harmonic oscillators. As $T$ increases, $Q_{v,i}(T)$ becomes proportional to $\prod _{l=1}^{N_{v,i}} T/\theta _{v,il}$ , and $\varphi _{i}$ becomes proportional to $[(T_{tr} \varGamma _i)/(T_{v,i} \varGamma _i^{tr})]^{N_{v,i}}$ . This limit is shown by the black dashed lines in figure 3 at $T_{tr} = 5000$ K, with the various lines corresponding to $N_{v,i} = 1$ , 3 and 6, which corresponds to species with 2, 3 and 4 atoms, respectively. By $5000$ K, it is clear that the efficiency function has nearly reached its limiting behaviour.

Figure 3.

Species vibrational efficiency function versus $T_{v,i}/T_{tr}$ for various molecules in hydrogen combustion, with $\alpha _{v,i}^\prime = 0.5$ . The black dashed lines correspond to $[(T_{tr} \varGamma _i)/(T_{v,i} \varGamma _i^{tr})]^{N_{v,i}}$ , with $N_{v,i}$ varied between 1, 3 and 6, which corresponds to species with 2, 3 and 4 atoms, respectively.

Figure 3. Long description

Three line graphs depict the species vibrational efficiency function for various molecules in hydrogen combustion at different translational temperatures. Panel A: The line graph shows the species vibrational efficiency function for translational temperature of 1000 K. The x-axis represents the ratio of vibrational to translational temperature (T_v,i / T_tr) ranging from 0.3 to 1.2. The y-axis represents the vibrational efficiency function (φ_i) ranging from 0 to 1.4. The graph includes lines for H2, OH, O2, H2O, HO2, and H2O2, each with distinct colors. Panel B: The line graph shows the species vibrational efficiency function for translational temperature of 3000 K. The x-axis represents the ratio of vibrational to translational temperature (T_v,i / T_tr) ranging from 0.3 to 1.2. The y-axis represents the vibrational efficiency function (φ_i) ranging from 0 to 1.4. The graph includes lines for H2, OH, O2, H2O, HO2, and H2O2, each with distinct colors. Panel C: The line graph shows the species vibrational efficiency function for translational temperature of 5000 K. The x-axis represents the ratio of vibrational to translational temperature (T_v,i / T_tr) ranging from 0.3 to 1.2. The y-axis represents the vibrational efficiency function (φ_i) ranging from 0 to 1.4. The graph includes lines for H2, OH, O2, H2O, HO2, and H2O2, each with distinct colors. The black dashed lines correspond to φ_i = T_v,i / T_tr, with varied between 1, 3, and 6, which corresponds to species with 2, 3, and 4 atoms, respectively.

Figure 4 compares how varied partition function definitions affect the efficiency function for chain-branching reactions in hydrogen–air combustion at fixed $T_{v,i}/T_{tr}$ ratios and $T_{tr}$ varied between 1400 and 5400 K. Here, we compare the proposed model, as defined via (4.7), for different energy descriptions: harmonic oscillator (HO), truncated harmonic oscillator (THO) and anharmonic oscillator (AHO). With regard to $\textrm{H}_{2}\textrm{O}$ , only the HO and THO models are examined. The maximum vibrational energy of each molecule is defined as the dissociation/bond energy of the molecule as collated by Dean (Reference Dean1999).

Figure 4.

Species vibrational efficiency function versus $T_{tr}$ for various chain-branching reactions in hydrogen combustion, with $\alpha _{v,i}^\prime = 0.5$ .

For vibrationally hot and cold $\textrm{O}_{2}$ and $\textrm{H}_{2}$ , as shown in figures 4(a) and 4(b), the HO model reproduces similar efficiency function curves to those of the THO and AHO models, suggesting that vibrational energies exceeding the dissociation energy are negligible, and that anharmonic corrections only have a marginal effect. For vibrationally hot and cold $\textrm{H}_{2}\textrm{O}$ , as shown in figure 4(c), the THO and HO models diverge at approximately 3000 K, resulting in a 20 % discrepancy in $\varPhi _{f,j}$ at 5000 K, indicating that vibrational energies above the bond energy are non-negligible. This alludes to a significant discrepancy in the limiting behaviour of these two models. As the temperature increases, the partition function associated with the HO model approaches its limit as described for figure 3, whereas the THO model limits to a value proportional to the bond energy. This results in a vibrational efficiency function that limits to one. Thus, at sufficiently high temperatures, the vibrational effect is negligible. A similar limiting behaviour is observed for the AHO model. It is unclear which limiting behaviour is more favourable.

Beyond the HO, THO and AHO models examined in figure 4, we also include the results for the generalised form as defined via (3.19) and (3.20), which limits the vibrational scaling function by the forward and reverse activation energy. Here, we set $\delta ^\prime = \delta ^{\prime \prime } = 1$ and $\beta _f = \beta _r = 0.8$ . The general model adopts the AHO description for $\textrm{O}_{2}$ and $\textrm{H}_{2}$ and the THO model for $\textrm{H}_{2}\textrm{O}$ . Compared with the proposed model, where $\delta ^\prime = \delta ^{\prime \prime } = 0$ , the generalised model reduces the magnitude of the vibrational efficiency function by limiting the impact of the effective vibrational component of the reaction barrier as defined via (3.20). Thus, high-lying vibrational states may be limited by the activation energy, diminishing the vibrational influence on the reaction rate. This is particularly noticeable for larger molecules, as seen with $\textrm{H}_{2}\textrm{O}$ in figure 4(c), which have more vibrational modes and thus a greater vibrational energy.

5.2. Comparisons with historical models

Figure 5 shows the vibrational efficiency function for several chain-branching reactions in hydrogen–air detonations, with the proposed model compared with historical models, including the modified Park model presented by Voelkel et al. (Reference Voelkel, Raman and Varghese2016b ), the CVCV model presented by Knab et al. (Reference Knab, Frühauf and Messerschmid1995) and the model presented by Arsentiev et al. (Reference Arsentiev, Loukhovitski and Starik2012). For the historical models, dissociation and activation energies were taken from Dean (Reference Dean1999) and Burke et al. (Reference Burke, Chaos, Ju, Dryer and Klippenstein2012), respectively. For each reaction, we set $T_{tr} = 1500$ K and $\alpha _{v,i} = 0.5$ . The modified Park model (Voelkel et al. Reference Voelkel, Raman and Varghese2016b ) was directly calibrated to QCT results. Thus, we consider this model to be the target model for these specific reactions. The $\alpha$ parameter associated with the CVCV model was reduced from the recommended value of 0.8 to 0.3 for these comparisons, which provided a much closer match to the modified Park model.

Figure 5.

Vibrational efficiency function versus reactant vibrational temperature scaled by $T_{tr}$ for key reactions in hydrogen combustion, with $T_{tr} = 1500$ K and $\alpha _{v,i} = 0.5$ .

Figure 5. Long description

The image contains four line graphs labeled (a), (b), (c), and (d), each depicting the vibrational efficiency function (phi) versus the reactant vibrational temperature scaled by the translational temperature. Each graph compares data from different models: Proposed model (red line), Knab et al. (1995) (blue line), Arsentev et al. (2012) (orange line), and Voelkel et al. (2016b) (black line). Panel A shows the reaction H * O2, Panel B shows H2 * O, Panel C shows H2 * OH, and Panel D shows H2 * OH. The x-axis represents the ratio of vibrational to translational temperature (Tv,i / Ttr), and the y-axis represents the vibrational efficiency function (phi). Each graph includes a dashed horizontal line at phi = 1.0 and a dashed vertical line at Tv,i / Ttr = 1.0. The trends show how the vibrational efficiency function varies with the reactant vibrational temperature for different models.

In figure 5(a), the vibrational efficiency function of $\textrm{O}_{2}$ has a significant effect on the reaction rate. For $ T_{v,\textrm{O}_{2}} /T_{tr} = 0.4$ ( $ T_{v,\textrm{O}_{2}} = 600$ K), the rate is reduced by approximately $50$ % based on the target model. The proposed model closely matches the target model for $T_{v,\rm O_2}/T_{tr} \gt 0.8$ , below which the proposed model begins to flatten out, resulting in less reaction suppression for significantly vibrationally cold states. Both the CVCV and Arsentiev models overestimate the role of vibrational non-equilibrium for both vibrationally hot and cold states, until $T_{v,\rm O_2}/T_{tr} \lt 0.3$ . In figure 5(b), vibrational non-equilibrium of $\textrm{H}_{2}$ has a lesser effect, with the reaction reduced by $15$ % at $ T_{v,\textrm{H}_{2}} /T_{tr} = 0.4$ . For vibrational cold flows, each model is within 10 % of the others, with the proposed model exhibiting a slightly lesser vibrational effect. For vibrationally hot flows, the proposed model is significantly closer to the target model.

In figure 5(c), the vibrational efficiency function for our proposed model is the same as in figure 5(b). However, all of the other models show a reduced vibrational non-equilibrium effect, accounting for the reduced activation energy barrier associated with this reaction. Now, the reaction rate is only reduced by $8$ % at $ T_{v,\textrm{H}_{2}} /T_{tr} = 0.4$ . The Knab and Arsentiev models do a better job of capturing this subtle effect, but each model is still within a few percent. In figure 5(d), the target model exhibits a marginal non-equilibrium effect, suggesting that $\textrm{OH}$ vibrational non-equilibrium is not relevant for this reaction. None of the other models are able to account for this result, with the Knab model producing the best result.

Two final observations are made regarding figure 5. First, we note that the slope of the efficiency function for the proposed, CVCV and Arsentiev models decreases as $T_{v,i}$ decreases, in contrast to the target model. Subsequently, as $T_{v,i}$ decreases, these models will eventually diverge from the target model. We expect such states to only be observed directly behind the shock (where $T_{v,i}/T_{tr} \approx 0.2$ for a stoichiometric Chapman–Jouguet detonation at standard temperature and pressure). However, we expect the species to relax somewhat quickly to a more moderate vibrationally cold state, where the proposed model is more well suited. Second, the point of this comparison was to validate the choice of our proposed model. We feel that, compared with the significantly more complex historical models (which themselves are not necessarily accurate representations of reality), the proposed model performs well and captures the underlying physics associated with vibrational non-equilibrium’s effect on reaction rates.

5.3. General features of chemical–vibrational exchange

Chemical–vibrational exchange creates and destroys vibrational energy from chemical reactions and scales with the vibrational energy. Based on the proposed model (4.10), consider (2.10a ) in non-conservative form, given by

(5.1) \begin{align} \frac {\partial e_{v,i}}{ \partial t} + u \frac {\partial e_{v,i}}{ \partial x} = \frac { S_{{\scriptsize T-V},i} + S_{{\scriptsize V-V},i} }{\rho Y_i} + \frac {e_{v,i}(\varGamma _i^{tr}) - e_{v,i}}{ \tau _{{\scriptsize C-V(c)},ij} } - \frac {e_{v,i}(\varGamma _i) - e_{v,i}}{ \tau _{{\scriptsize C-V(d)},ij} }. \end{align}

In this form, it is straightforward to show that, for a given species $i$ , positive reactions that create new molecules push $T_{v,i}$ towards $\varGamma _i^{tr}$ ( $e_{v,i} \,\rightarrow \, e_{v,i}(\varGamma _i^{tr})$ ). This occurs at a time scale proportional to the species molar production rate scaled by $e_{v,i}(\varGamma _i^{tr})/e_{v,i}^{tr}$ . In contrast, negative reactions that destroy molecules push $T_{v,i}$ towards zero ( $e_{v,i} \,\rightarrow \, 0$ ). This occurs at a time scale proportional to the species molar production rate scaled by $e_{v,i}(\varGamma _i)/e_{v,i}^{tr}$ .

Figure 6 shows $\varGamma _i^{tr}/T_{tr}$ and $\varGamma _i/T_{tr}$ versus $T_{v,i}/T_{tr}$ for varied values of $\alpha _{v,i}^\prime$ . For all vibrational efficiency parameters, $\varGamma _{i}^{tr}/T_{tr} \ge 1$ as shown in figure 6(a), which indicates that vibrational energy is created in a vibrationally hot state. As $T_{v,i}/T_{tr}$ increases beyond 1, this effect becomes diminished. When $\alpha _{v,i}^\prime = 0$ (vibrational energy does not influence reactions), then $\varGamma _{i}^{tr}/T_{tr} = 1$ , so positive reactions push $T_{v,i}$ towards thermal equilibrium, consistent with (5.1).

Figure 6.

Gamma temperatures versus vibrational temperature, both relative to $T_{tr}$ .

Figure 6(b) shows how the destruction of vibrational energy scales with $\varGamma _i$ and $T_{v,i}$ . When the species is vibrationally hot ( $T_{v,i}/T_{tr} \gt 1$ ), then $\varGamma _i/T_{tr} \gg 1$ and $e_{v,i}(\varGamma _i)/e_{v,i}^{tr} \gg 1$ , leading to rapid depletion of vibrational energy. As the species becomes vibrationally cold, then the destruction rate of vibrational energy goes to zero. When $\alpha _{v,i}^\prime = 0$ , then $\varGamma _i = T_{v,i}$ and, similarly, the destruction rate goes to zero as shown in (5.1). Thus, negative reactions do not alter the local vibrational energy.