Advanced models for vibrational–chemical coupling in multi-temperature flows

Chemical Physics 464 (2016) 1–13

Contents lists available at ScienceDirect

Chemical Physics journal homepage: www.elsevier.com/locate/chemphys

Advanced models for vibrational–chemical coupling in multi-temperature flows E. Kustova, E. Nagnibeda, G. Oblapenko ⇑, A. Savelev, I. Sharafutdinov Saint Petersburg State University, 198504, Universitetskiy pr., 28, Saint Petersburg, Russia

a r t i c l e

i n f o

Article history: Received 2 September 2015 In final form 28 October 2015

Keywords: Kinetic theory Vibrational–chemical coupling Non-equilibrium flows Relaxation rates

a b s t r a c t In this paper, self-consistent models for coupled vibrational relaxation and dissociation in multitemperature gas mixture flows are proposed on the basis of the kinetic theory methods. The Treanor–Marrone dissociation model is generalized taking into account the dependence of model parameter on the vibrational state. Multi-temperature dissociation rate coefficients are calculated on the basis of the improved Treanor–Marrone model using the Boltzmann and Treanor vibrational distributions; comparison with traditional models shows a significant difference in the dissociation rate coefficients for high temperatures. Generalization of the well known Landau–Teller model overcoming limitations of the original model is proposed. The developed models are validated against experimental data and assessed in non-equilibrium shock heated flows of O2/O and N2/N mixtures. It is shown that the proposed models provide a good accuracy in the wide range of flow conditions while being simple and computationally efficient. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Modeling of thermochemically non-equilibrium flows remains a challenging problem of hypersonic aerothermodynamics during the last decades. Although a wide variety of models for vibrational relaxation and chemical reactions have been proposed, the accuracy of fluid dynamic variables and heat transfer predictions still requires assessment and improvement. The most detailed description of strongly non-equilibrium reacting gas flows in the frame of continuum media approach is given by the state-to-state model [1–3]. This model treats each vibrational (or rovibrational) state of a molecule as a separate chemical species. Solving the master equations coupled to the fluid dynamic equations one can find the flow-field and the populations of internal states. The main advantage of the state-to-state approach is that it does not adopt any quasi-stationary distribution over internal energy states and therefore can be applied for arbitrary deviations from thermochemical equilibrium. This model is widely used now for simulations of stationary 1D [4–12] and simple 2D [13–15] flows. However its implementation for modeling complex 2D and 3D problems is hardly feasible due to very high computational costs, especially for viscous flow simulations requiring calculations of the state-to-state transport properties [16,17]. To the authors’ knowledge, 3D simulations using the

⇑ Corresponding author. http://dx.doi.org/10.1016/j.chemphys.2015.10.017 0301-0104/Ó 2015 Elsevier B.V. All rights reserved.

state-to-state approach have been done only in [18]. As a consequence, development of reduced multi-temperature models is of importance for engineering applications. Multi-temperature flow description is based on the assumption that a part of vibrational energy transitions proceeds faster than vibration–translation (VT) exchange and chemical reactions. Fast vibration–vibration (VV) transitions result in establishing quasistationary distributions over vibrational energy. Commonly used multi-temperature models are based on the Boltzmann vibrational distributions for harmonic oscillators [19,20]; their application is limited to flow conditions with weak impact of anharmonicity. Advanced multi-temperature models have been developed for the Treanor distribution taking into-account near-resonant VV exchange [21–23] and for more sophisticated combined Treanor– Gordiets–Boltzmann distribution accounting for different rates of vibrational transitions at various vibrational levels [24–26]. Reduced models based on the analysis of state-to-state distributions are discussed in [27,28]. An interesting model combining several Boltzmann distributions for different parts of the vibrational ladder is proposed recently in [29]. Development of a multi-temperature flow description requires modeling of vibrational relaxation, chemical reactions and coupling of these processes. For the state-dependent rate coefficients of vibrational transitions many models have been proposed, both analytical [30,4,31] and those based on molecular dynamic simulations [32–41]. Nonetheless in the multi-temperature flow simulations the Landau–Teller relaxation equation [42] with the

2

E. Kustova et al. / Chemical Physics 464 (2016) 1–13

Millikan–White relaxation time [43] is still commonly used, sometimes with empirical corrections [20]. It is worth noting that although the original Landau–Teller equation has been derived for weak deviations from thermal equilibrium under assumption of harmonic oscillations, at the present time it is widely implemented into the computational fluid dynamics (CFDs) for strongly non-equilibrium conditions and arbitrary model of vibrational spectra. Recently, a self-consistent model for the rate of vibrational relaxation based on the Chapman–Enskog method has been proposed in the frame of the state-to-state [44] and multi-temperature [45] flow descriptions. It is shown that in viscous flows, the rate of vibrational relaxation is specified by the rates of all other slow processes; the cross-coupling effect between the stress tensor and reaction rates is demonstrated. In the present study we use the model developed in [45] to assess the mixture thermochemical properties and establish the limits of applicability for the Landau–Teller equation. Rates of chemical reactions are of crucial importance for the accurate prediction of the non-equilibrium flow dynamics and heat transfer. In vibrationally non-equilibrium flows, it is necessary to account for strong coupling of vibrational relaxation and chemistry, and one-temperature Arrhenius rate coefficients are not applicable. There are numerous models of cross-coupling between vibrational and chemical non-equilibrium such as CVDV [46–50] (coupled vibration–dissociation–vibration), CVCV [51,50] (coupled vibration–chemistry–vibration), but they are basically applicable only to inviscid flows (except [49]) and rely on quasi-empirical expressions for the coupling terms derived on the basis of the Boltzmann vibrational distributions. On the other hand, as is demonstrated in [6,8,52,53], the averaged reaction rate coefficients depend strongly on the shape of non-equilibrium distributions. The model developed in [45] overcomes these limitations, however the coupling terms derived in this paper have a rather complicated form and it is hard to implement them directly into 3D viscous flow solvers. Therefore it is interesting to assess simpler models. For this purpose we have chosen the Treanor–Marrone CVDV model. It is quite simple, efficient and is widely used in many hypersonic applications. The model contains an adjustable parameter which is basically set to constant. A thorough analysis performed in [54] showed that setting this parameter to the constant value or linear function of temperature cannot provide satisfactory agreement with dissociation rate coefficients obtained from the molecular dynamic calculations. Accurate definition of the parameter in the Treanor–Marrone model can significantly improve the predicted dissociation rates in non-equilibrium flows. The objectives of the present paper are:
To generalize the Treanor–Marrone model taking into account the dependence of its parameter on the vibrational state and temperature.
To estimate the impact of the model and different vibrational distributions on the multi-temperature dissociation rate coefficients.
To derive a rigorous model of VT relaxation and establish the limits of applicability for the Landau–Teller model.
To validate the proposed models.
To study the shock heated flows of binary mixtures O2/O and N2/N.

calculated under various free stream conditions on the basis of modified Treanor–Marrone and Landau–Teller expressions are compared with those found experimentally [55] as well as with the solutions obtained using the traditional models such as original Landau–Teller and Treanor–Marrone models, the forced harmonic oscillator model (FHO) [31], and analytical expressions of Shwartz, Slawsky and Herzfeld (SSH) [30]. The results allow estimating the limits of applicability for the proposed simple and numerically efficient models. 2. Generalized Treanor–Marrone model Modeling of vibration–dissociation coupling in reacting flows still represents a challenging task. Experimental measurement of dissociation rate coefficients in vibrationally non-equilibrium flows is fairly complicated, and consequently, experimental data on the state-to-state or multi-temperature dissociation rate coefficients are scarce (some recent experimental results for oxygen can be found in [55]). In this situation, molecular dynamic methods are suitable to provide missing information on the reaction cross sections and state-specific rate coefficients [35–40]. However direct implementation of trajectory calculations into hypersonic flow simulations is not feasible at the present time due to extremely high computational efforts. Moreover, the state-to-state dissociation rate coefficients provided in the database [40] are given only for a few selected vibrational states, and the proposed interpolation formulas are valid for the temperatures lower than 10,000 K. Therefore accurate analytical models for vibration–dissociation coupling are of importance for computational fluid dynamics. Theoretical models for multi-temperature dissociation rate coefficients have been proposed by many authors. The simplest empirical Park’s model widely used in CFD is described in [20]; analytical Treanor–Marrone, Kuznetsov, Macheret–Fridman models [46–48] are better justified but sometimes contain cumbersome expressions hardly implementable to CFD solvers. Theoretical models generally include adjustable parameters, and correct choice of these parameters can significantly improve the model accuracy. For the present study. we have chosen the Treanor–Marrone model of preferential dissociation first proposed for two-temperature flows of harmonic oscillators in [46] and later generalized for multi-temperature flows of anharmonic oscillators [23] and for the state-to-state model [54]. This simple model is based on clear physical reasoning and can be easily applied for non-equilibrium flow simulations. The Treanor–Marrone model involves one adjustable parameter U having the dimension of temperature and characterizing the decrease in the dissociation rate with decreasing in the vibrational level. In CFD, this parameter is basically set to constant, U ¼ 1 (which corresponds to equiprobable dissociation from each vibrational level), U ¼ D=ð6kÞ (D is the dissociation energy, k is the Boltzmann constant), or it is defined as a linear function of temperature, for instance, U ¼ 3T. The latter two cases describe preferential dissociation from high vibrational levels. Analysis performed in [54,56] shows that constant values of U (or linear function of T) cannot provide satisfactory agreement of the state-specific dissociation rate coefficients with those obtained by quasi-classical trajectory calculations in the whole range of temperature and vibrational states. Let us consider first the state-to-state modification of the TreM

The paper is organized as follows. First we develop the generalized Treanor–Marrone model and study multi-temperature dissociation rate coefficients. Then we introduce the accurate model for the vibrational–chemical coupling and discuss its relation to the Landau–Teller equation. Then the models are applied to the high-temperature binary mixture flows behind the shock waves. Distribution of fluid dynamic variables in O2/O and N2/N mixtures

anor–Marrone model. The rate coefficient ki;diss of dissociation from the vibrational state i after a collision with a partner M can be calculated in the form: M

M

ki;diss ¼ Z M i kdiss; eq ðTÞ M

ð1Þ

where kdiss; eq ðTÞ is the thermal equilibrium dissociation rate coefficient,

ZM i

is the non-equilibrium factor determined as:

3

E. Kustova et al. / Chemical Physics 464 (2016) 1–13





ZM i ¼ Z i ðT; UÞ ¼

Z vibr ðTÞ ei 1 1 þ exp Z vibr ðUÞ k T U

 ;

ð2Þ

ei is the vibrational energy of the ith state, Z vibr ðTÞ is the equilibrium vibrational partition function Z vibr ðTÞ ¼

X

 exp 

i

ei  kT

:

ð3Þ

The thermal equilibrium dissociation rate coefficient can be calculated on the basis of the Arrhenius law:


 E M kdiss; eq ¼ AT b exp  ; T

ð4Þ

A; b are the parameters, E ¼ D=k is the dissociation energy in temperature units. The accuracy of the Treanor–Marrone model depends on two factors, and the first of them is the choice of parameters in the Arrhenius law. As is shown in [56], using the parameters A; b from [20] prevents to determine the parameter U providing satisfactory agreement for the dissociation rate coefficients with the molecular dynamic calculations at high temperatures and high vibrational levels. After assessment of various parameter sets, we can recommend those from [38] for O2/O and [57] for N2/N mixtures. These parameters are given in Table 1:

Table 1 Recommended parameters in the Arrhenius law. A (m3/s)

b

E (K) 60496

O þ O2

6:1953  10

0:55724

O2 þ O2

5:33  1011

1:0

59370

N þ N2

1:0  1011

0:68

113500

N2 þ N2

4:1  1012

0:62

113500

12

Table 2 Coefficients in Eq. (5) for O2 . n; k

an

bk ~ei < 0.67 eV

2:9177  10

1:0034  10

1

5:7075  104

1:3047  105

4

5

~ ei > 0.67 eV

4

0

1:2267  104 3:0141  105 –

5

2

2:7321  10

3

6:6291  105

2:4290  10 –

4

8:8254  105

–

–

5

6:8447  105

–

–

6

3:2435  105

–

–

7

9:5511  104

–

–

8

1:7081  104

–

–

9

1:7002  103

–

–

10

7:2293  101

–

–

–

Once the parameters A; b are specified, we can try to look for the optimal choice of the parameter U depending on the temperature and vibrational state. For this purpose we use the results of quasi-classical trajectory calculations [36,38–40]. In these references, the rate coefficients of dissociation reactions O2 ðiÞ þ O O þ O þ O and N2 ðiÞ þ N N þ N þ N are calculated for selected vibrational levels and interpolation formulas are proposed. However for our study, a few vibrational levels considered in [36,38–40] are not enough. Moreover, the interpolation formulas are derived only for the temperature range limited by 10,000 K. To solve our problem, first of all, we have found optimal parameter value for all energy levels and temperatures given in [40] using bisection method with low border 1 and upper border 2  106 and accuracy 103 . Then we have carried out exponential interpolation of obtained values in the whole temperature interval for fixed energy state with successive interpolation of the results for all energy states. As a result, we get piece-wise continuous functions for the parameter U depending on the vibrational energy and temperature T. The final formula takes the form:

! N K X X n k Uði; TÞ ¼ an ~ei exp T bk ~ei ; n¼0

ð5Þ

k¼0

where ~ei is the vibrational energy in eV, N and K depend on the reaction and vibrational energy range. The coefficients an ; bk for oxygen and nitrogen are presented in Tables 2 and 3 respectively. Note that this formula can be used for both harmonic and anharmonic oscillators and for any vibrational ladder, since it depends on the vibrational energy and not on the vibrational quantum number. To assess the proposed model, we have calculated the state-to-state dissociation rate coefficients using Eqs. (1), (2), and (5) and compared them to those given in [40]. Rate coefficients calculated using our model show a very good agreement with ki;diss obtained in [40] by means of trajectory calculations in the whole range of temperatures and vibrational states except for the highest vibrational levels. In the worst cases (for the last vibrational state), calculated ki;diss and those given in [40] do not differ more than twice whereas if we use traditional values for the parameter U such as U ¼ 1; U ¼ D=ð6kÞ; U ¼ 3T, the discrepancy can reach two orders of magnitude. In Figs. 1 and 2, the parameter U calculated in oxygen and nitrogen using Eq. (5) for different vibrational states and temperatures is presented. The vibrational energy levels from [40] are used for these plots. For comparison, the lines corresponding to U ¼ 3T and U ¼ D=ð6kÞ are also plotted. The parameter U increases monotonically with temperature for N2 as well as for low and intermediate states of O2 while for the states with i P 35 in oxygen it decreases. Strongly non-monotonic variation of U with vibrational level can be seen from Fig. 2. It is worth noting that in all cases, neither constant nor linear approximation of U does not describe correctly its behavior in the whole range of temperature and energy states. Let us discuss now the two-temperature dissociation rate coefficients. Once the state-to-state reaction rate coefficients are

Table 3 Coefficients in Eq. (5) for N2 . n; k

~ei < 2.85 eV

2.85 eV < ~ei < 9.27 eV

~ei > 9.27 eV

an

bk

an

bk

0

1:8586  104

5:4119  105

1:5970  104

1:7421  104

1

1:2897  103 –

7:3290  106 –

8:4580  102

2

6:9535  10

2:7691  105

3

–

–

1:1342  102

3:3879  106

4

–

–

8:9603  105 2

0

6:9833  10

1:3575  107

D U ¼ 3k

4

E. Kustova et al. / Chemical Physics 464 (2016) 1–13 200000 U, K

(a)

O2

180000 160000

i=3 i=15 i=35 U=3T U=D/6k

140000 120000 100000

U, K

60000

(b)

N2 i=3 i=20 i=40 U=3T U=D/6k

40000

80000 20000

60000 40000 20000

T, K

0 0

5000

10000

15000

0

20000

0

5000

10000

15000

20000

T, K

Fig. 1. Parameter U as a function of T for selected fixed vibrational states in O2 (a) and N2 (b).

300000

U, K

(a)

O2

70000 U, K

T=2000 K T=10000 K T=20000 K U=D/6k

60000

250000

T=2000 K T=10000 K T=20000 K U=D/6k

200000

(b)

N2

50000 40000

150000

30000 100000

20000

50000

10000

0 0

10

20

30

0

i

40

0

10

20

30

40

50

60

70

i

Fig. 2. Parameter U as a function of i for fixed temperature values in O2 (a) and N2 (b).

specified, one can easily calculate the averaged rate coefficients using any quasi-stationary vibrational distribution, M

kdiss ðT; T v Þ ¼

1 X M ni ki;diss ðTÞ nm i

ð6Þ

where ni is the number density of molecules on the vibrational level P i; nm ¼ i ni is the number density of molecules, T v is the vibrational temperature of the molecular species. For the present study, we use the Treanor and Boltzmann distributions over vibrational energy. The Treanor distribution is introduced by the expression [21,22]


 ni 1 ei  ie1 ie1 ; ¼ exp   nm Z vibr ðT; T 1 Þ kT kT 1

ð7Þ

e1 is the vibrational energy of the first level, T 1 is the temperature of

The non-equilibrium Boltzmann distribution for harmonic oscillators can be derived if we put ei ¼ ie1 and T 1 ¼ T v in Eqs. (7) and (8). Substituting distribution (7) into Eq. (6) and taking into account Eqs. (1) and (2) one can easily get the two temperature dissociation rate coefficient based on the Treanor distribution M

kdiss ðT; T 1 Þ ¼

ð9Þ Similarly, if we use the Boltzmann distribution with the vibrational temperature T v , we obtain M

kdiss ðT; T v Þ ¼

the first level, and the non-equilibrium partition function reads

Z vibr ðT; T 1 Þ ¼

X i


 ei  ie1 ie1 : exp   kT kT 1

ð8Þ

(a)

3

-18

10

-20

10

kdiss, m /s

,

N2 T1=4000 K

Treanor distr.



 M X kdiss; eq Z vibr ðTÞ ei ie1 1 1 :  exp þ kU k T T1 Z vibr ðUÞZ vibr ðT; T 1 Þ i



 M kdiss; eq Z vibr ðTÞ X ei 1 1 1 exp : þ  k U T Tv Z vibr ðUÞZ vibr ðT v Þ i

Comparison of two-temperature dissociation rate coefficients calculated on the basis of Eqs. (9) and (10) in O2 and N2 using the Treanor–Marrone model with the constant U ¼ D=ð6kÞ and -17

10 10

-20

10

-23

10

-26

10

-29

10

-32

10

-35

(b)

3

kdiss, m /s N2, T1=1421 K Treanor distr.

-22

10

-24

10

U=D/6k U(i,T), Eq.(5)

-26

10

-28

10

T, K 5000

10000

15000

20000

ð10Þ

U=D/6k U(i,T), Eq.(5)

5000

10000

15000

20000 T, K

Fig. 3. Dissociation rate coefficients based on the Treanor distribution as functions of T for N2 at T 1 ¼ 4000 K (a) and 1421 K (b).

5

E. Kustova et al. / Chemical Physics 464 (2016) 1–13 3

3

kdiss, m /s

kdiss, m /s

-17

10

-28

10

-33

10

-19

10

N2

-38

10

O2

-21

10

T1=4000K, Boltzmann T1=4000K, Treanor T1=1421K, Boltzmann T1=1421K, Treanor

-43

10

-48

10

-23

10

-53

10

1000

T1=4000K, Boltzmann T1=4000K, Treanor T1=1421K, Boltzmann T1=1421K, Treanor

-25

2000

3000

4000

10

6000 T, K

5000

5000

10000

20000 T, K

15000

Fig. 4. Dissociation rate coefficients for N2 and U ¼ Uði; TÞ as functions of T. Comparison for the Boltzmann and Treanor distributions.

Fig. 6. Dissociation rate coefficients for O2 as functions of T. Comparison for the Boltzmann and Treanor distributions. U ¼ Uði; TÞ.

the generalized model (Eq. (5)) is presented in Figs. 3 and 5. One can notice considerable difference, especially at high temperature; the discrepancy reaches two orders of magnitude in N2, and more than one order of magnitude in O2. While for nitrogen using the constant U yields overpredicted kdiss ðT; T v Þ values, for oxygen the opposite trend is found. It can be explained analyzing Eqs. (9) and (10). As we can see, the parameter U enters there twice, in the partition function and in the exponentials as one of the terms. For low vibrational levels making the major contribution to the distributions, the value of 1=U can be neglected with respect to 1=T. Therefore the behavior of the dissociation rate coefficients is mainly determined by the term Z vibr ðUÞ. For oxygen, Uði; TÞ > D=ð6kÞ in the whole temperature range (see Fig. 2(a)), which leads to lower values of Z vibr ðUÞ and consequently higher kdiss ðT; T v Þ. For nitrogen (Fig. 2(b)), the dependence is more complicated, and some compensation effects result in the decrease in kdiss ðT; T v Þ calculated using Uði; TÞ compared to those obtained with constant U. Comparison of two-temperature dissociation rate coefficients calculated on the basis of the Treanor and Boltzmann distributions (Figs. 4 and 6) shows that for oxygen, anharmonic effects are of importance. At low temperatures dissociation rate coefficients for anharmonic oscillators demonstrate essentially non-monotonic behavior contrarily to the case of harmonic oscillators. For nitrogen, the impact of anharmonicity is weak; it becomes noticeable only for low temperature T and high vibrational temperature.

specific internal energy. In mixtures of anharmonic oscillators, relaxation equations are written for the numbers of vibrational quanta per unit mass [22,45]. In a binary mixture of molecules and atoms, the set of equations describing the flow includes the following relaxation equation:

qm

dW þ r  qW ¼ RW  Wmm Rreact þ W r  ðqm Vm Þ; dt

ð11Þ

where W is the number of vibrational quanta in the molecular species per unit mass, qm is the density of the molecular species, qW is the flux of vibrational quanta, mm is the mass of the molecular species, Vm is the diffusion velocity of the molecules, RW characterizes the change in number of vibrational quanta due to slow vibrational energy transitions and chemical reactions, Rreact characterizes the change in the number of molecules due to chemical reactions. The number of vibrational quanta per unit mass W is defined as follows:

WðT; T 1 Þ ¼


 X 1 ei  ie1 ie1 : i exp   mm Z vibr ðT; T 1 Þ i kT kT 1

ð12Þ

When the vibrational spectrum of molecular species in a flow is considered to be harmonic, Eq. (11) can be multiplied by the energy of the first vibrational level e1 and this yields the equation describing the relaxation of vibrational energy:

qm

dEvibr þ r  qvibr ¼ Rvibr  Evibr mm Rreact þ Evibr r  ðqm Vm Þ: dt

ð13Þ

Here Evibr is the specific vibrational energy, qvibr is the vibrational 3. Revision of Landau–Teller relaxation equation

energy flux, and Rvibr characterizes the change in vibrational energy due to slow vibrational energy transitions and chemical reactions.

In multi-temperature flows, conservation equations and equations of chemical kinetics are coupled to relaxation equations for

The terms RW and Rreact can be expressed as functions of reaction rates n_ c;r [45]:

10

-16

10

-18

10

-20

10

-22

10

-24

(a)

3

kdiss, m /s

-26

10

-28

(b)

3

kdiss, m /s

-17

10

-18

,

10

O2 T1=4000 K

,

O2 T1=4000 K

-19

10

Treanor distr.

10

-16

10

Boltzmann distr.

-20

10

U=D/6k U(i,T), Eq.(5) T, K 5000

10000

15000

U=D/6k U(i,T), Eq.(5)

-21

10

20000

-22

10

T, K 5000

10000

15000

20000

Fig. 5. Dissociation rate coefficients based on the Treanor (a) and Boltzmann (b) distributions as functions of T for O2 at T 1 ¼ 4000 K.

6

E. Kustova et al. / Chemical Physics 464 (2016) 1–13

Rreact ¼ Na

2¡3; XNX

mr;ci n_ 2¡3;r ¼ Na

r¼1

i

X

N 2¡3;

n_ 2¡3;r ;

ð14Þ

r¼1

Nc Nc X X X X X i mr;ci n_ c;r ¼ Na Din_ c;r ;

RW ¼ N a

c¼VT;2¡3 r¼1

i

ð15Þ

c¼VT;2¡3 r¼1

where N VT and N 2¡3 are the amount of VT transitions and dissociation–recombination reactions occuring in the flow (for each specific reaction r, the chemical species of the collision partner and the vibrational levels before and after the transition are fixed), N a is the Avogadro constant, n_ c;r is the reaction rate of a reaction r of type

c; mr;ci is the global stoichiometric coefficient for the rth transition, 0

where Eeq vibr ðTÞ is the equilibrium vibrational energy per unit of mass and sVT is the VT relaxation time. In [45] it was shown that the Landau–Teller expression can be obtained from a strict kinetic-theory definition of the relaxation terms (by expanding the expressions in a power series and keeping only linear terms) under the following assumptions: (1) the vibrational spectrum of molecular species is harmonic; (2) deviations from vibrational equilibrium are weak (T=T v  1); (3) vibrational specific heats are constant. An intermediate version of the Landau–Teller expression was also obtained:

RVT ¼

0

Di ¼ i  i, and i denotes the vibrational level after the inelastic collision. For a harmonic oscillator model the term Rvibr ¼ e1 RW may be written as

Rvibr ¼ Na

Nc X X Den_ c;r ;

ð16Þ

c¼VT;2¡3 r¼1

where De ¼ ei0  ei . In the zero-order approximation of the Chapman–Enskog method, the reaction rates n_ c;r are linear functions of corresponding rate coefficients: r ð0Þ n_ ð0Þ r ¼ C kf ;r

mðrÞ lc
YY nci r;ci : Na c¼m;a i¼1

ð17Þ

ð0Þ

Here kf ;r is the rate coefficient of the process r; lc is the number of vibrational levels in chemical species c; mr;ci is the reactant stoiðrÞ

chiometric coefficient, and Cr is a function of the corresponding generalized affinity Ar [45]:




 A Cr ¼ 1  exp r : kT

ð18Þ

X nM T ðT  T v Þqm cvibr : Tv nsVT M M

ð24Þ

where nM is the number density of the partner species M; n is the mixture number density, cvibr is the specific heat of vibrational degrees of freedom of the molecular species, sVT M is the VT relaxation time for collisions with chemical species M. It was shown that such an expression for the relaxation terms provides significantly better agreement with strict kinetic theory results for a harmonic oscillator case. Thus, it is of interest to compare results obtained using different versions of the Landau–Teller formula, since the modified expression can be considered as an easy-to-implement and more accurate replacement of the classic Landau–Teller formula. A model for calculation of vibrational relaxation times is also needed. The Millikan–White formula is commonly used, but as shown in [55], it does not provide a satisfactory agreement with experimental data at high temperatures (T > 6000 K) even with the correction proposed by Park [20]. VT relaxation times can be calculated using a rigorous kinetic-theory definition, if the crosssections of various VT transitions are known:

1

¼

sVT d

 2 VT 4kn DE vibr ; mm cvibr md

ð25Þ

DE vibr ¼ ðei0  ei Þ=kT; i and i are the vibrational levels before and after the VT transition and averaging operators for molecule–molecule and molecule–atom collisions are defined for any function F ij of molecular velocity by: 0

In the first-order approximation of the Chapman–Enskog method, the rates of slow process are found in the form

_ ð1Þ n_ r;c ¼ n_ ð0Þ r;c þ nr;c :

ð19Þ

The first-order corrections to the reaction rates n_ r;c are not independent quantities, and cross-coupling effects have to be taken into account [45]: ð1Þ

Na n_ r;ð1Þc ¼ kT

X lc;r;b;s Cb;s  lc;r;v r  v :

ð20Þ

ð26Þ

dEvibr ¼ Rvibr  Evibr mm Rreact : dt

ð21Þ

The relaxation term Rvibr can be presented as follows:

R

vibr

¼ RVT þ Rchem ; VT

where R

ð22Þ

is the term describing the change in vibrational energy

due to VT transitions and Rchem is the term describing the change in vibrational energy due to chemical reactions. In computational fluid dynamics, the relaxation terms RVT for various chemical species are commonly modeled using the Landau–Teller formula [42]:

RVT ¼ qm

1=2 X Z kT sj sl F ij g 30 2 2pmmm Z ini0 jlj0 l0 int

 0 eij enl ði þ nÞe1 1 1 2 j0nl0  exp g 20     rimm;ijnl d Xdg 0 ; T1 T kT kT k

b;s

Here lc;r;b;s and lc;r;v are integral brackets of scalar terms occurring in the first-order correction to the distribution function. In the present study, however, we limit our consideration to inviscid flows, and therefore, neglect the influence of cross-coupling and velocity divergence on the reaction rates. In the case of an inviscid flow, the Eq. (13) is reduced to

qm


hFiVT mm ¼

Eeq vibr ðTÞ  Evibr ðT v Þ

s

VT

;

ð23Þ


hFiVT ma ¼

1=2 X Z kT sj F ij g 30 2pmma Z int ii0jj0

 0 0 eij ie1 1 1 2 j  exp g 20   rima;ij d Xdg 0 ;  kT k T1 T

ð27Þ

where mcd ¼ mc md =ðmc þ md Þ is the collision reduced mass, sj is the degeneracy of the molecular state with rotational energy ej (since we only consider diatomic molecules, the degeneracy of vibrational states is equal to 1), Z int is the internal partition function, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi g 0 ¼ mcd =2kT g is the dimensionless relative velocity, g is the relative velocity of the colliding particles, i and n are the vibrational levels of the colliding molecules,

0 0

0

j nl rimm;ijnl is the cross section of a

VT transition in which the vibrational level i and rotational level j 0 0 of a molecule change to i and j after colliding with a molecule with vibrational level n and rotational level l, the latter of which changes 0 0

j to l ; rima;ij is the cross section of a VT transition in which the vibra0

0

0

tional level i and rotational level j of a molecule change to i and j after colliding with an atom.

7

E. Kustova et al. / Chemical Physics 464 (2016) 1–13

In Fig. 7 the ratio of relaxation times for molecule–atom collisions obtained using the strict definition (25) and (27) to those calculated using the Millikan–White formula is plotted for various values of T and T 1 . Note that contrarily to Eq. (25), the Millikan– White formula does not include any dependence of sVT on the vibrational temperature. It can be seen that the choice of vibrational spectrum (harmonic and anharmonic) does not play a significant role except for oxygen at high temperatures. At high values of T, the relaxation times calculated using formula (27) are noticeably larger than those given by the Millikan–White formula, especially for large ratios T=T 1 . The calculation of integrals (26) and (27) requires much computational efforts which is not favorable for their implementing to CFD solvers. On the other hand, pre-computed values of sVT obtained by using the strict definition of the VT relaxation time can be approximated numerically in a wide range of temperatures and these numerical approximations can be used as more accurate substitutes for the currently used Millikan–White formula.

In Fig. 8, the ratio of relaxation terms obtained using expressions (16) (considering only VT transitions) and (24) to those obtained by the Landau–Teller formula (23) are shown under various conditions for pure nitrogen and oxygen flows. For calculations, the VT process cross-sections rVT;r were assumed to be in the form

rVT;r ¼ PVT;r rtot;el ;

ð28Þ

where PVT;r is the probability of the VT transition and rtot;el is an elastic total cross-section. The forced harmonic oscillator (FHO) model [4] was used to calculate the probabilities PVT;r , while the elastic cross-sections were calculated with the Variable Soft Sphere (VSS) model [58] with data from [59]. VT relaxation times were calculated using the strict definition (25). As can be seen from Fig. 8, the modified Landau–Teller formula (24) provides better agreement with strict kinetic theory results, especially in cases of moderately strong vibrational non-equilibrium. For flows with T v < T (conditions typical for flows behind shock waves), the original Landau–Teller expression

9

4.5

(a)

4.0

(b)

8

7

3.5

3.0

6

2.5

5

2.0

4

1.5

3

1.0

2

0.5

1

0.0 2000

3000

4000

5000

6000

7000

8000

9000

0 2000

10000

3000

4000

5000

6000

7000

8000

9000

10000

Fig. 7. The ratio of VT relaxation times for N2 + N (a) and O2 + O (b) collisions obtained using formula (27) to relaxation times calculated using the Millikan–White formula as a function of temperature T for various values of T 1 .

5

12

Kinetic theory, N2

Kinetic theory, N2

Kinetic theory, O2

Kinetic theory, O2 10

Modified Landau-Teller, N2

Modified Landau-Teller, N2

Modified Landau-Teller, O2

4

Modified Landau-Teller, O2

8

3

6

2 4

(b)

1 2

(a) 0 0

0 2000

4000

6000

8000

10000

12000

14000

16000

0

2000

4000

6000

8000

10000

12000

14000

16000

Fig. 8. The ratio of relaxation terms in O2 and N2 for various models to the relaxation term given by the Landau–Teller expression (23) as a function of temperature T with T v ¼ 1000 K (a) and as a function of vibrational temperature T v with T ¼ 10; 000 K (b).

8

E. Kustova et al. / Chemical Physics 464 (2016) 1–13

underestimates vibrational relaxation rates, while for flows with T v > T, (typical conditions in expanding flows), the original Landau–Teller expression overestimates vibrational relaxation rates, and in general, gives accurate results only near the point of vibrational equilibrium (T v ¼ T). 4. Application to a binary mixture flow behind shock wave 4.1. Governing equations In this section we apply the models discussed above for simulations of non-equilibrium flows of binary mixtures O2 /O and N2 /N behind shock waves. Experimental data (see [60,50,61,55]) show that in shock heated gases translational and rotational energies equilibrate much faster than vibrational energy relaxation and chemical reactions proceed. Therefore, in the present study equilibrium distributions over translational and rotational energies are supposed to be established within a thin shock front and then they keep this form whereas vibrational energy distributions are strongly non-equilibrium in the relaxation zone behind a shock. It is also known that collisions with VV vibrational energy exchanges between molecules of the same species occur more often than those with vibrational and translational energy exchanges. Simultaneous translational, rotational and vibrational energy exchanges are not considered in our paper as less probable. Rapid VV exchanges of vibrational quanta between colliding molecules make it possible to consider level populations of N2 and O2 molecules along the relaxation zone in the form of Treanor distributions (7) depending on the gas temperature and vibration temperature of the first level T 1 . The total numbers of vibrational levels are set to lN2 ¼ 46; lO2 ¼ 35 for anharmonic oscillators and lN2 ¼ 33; lO2 ¼ 26 for harmonic oscillators. Since in this paper only binary mixtures of molecules and atoms are considered, subscripts ‘‘m” and ‘‘a” indicate molecules and atoms correspondingly. Treanor distributions are reduced to non-equilibrium Boltzmann distributions with vibrational temperatures T v for harmonic oscillators (ei ¼ ie1 ) and come to Boltzmann distributions with the gas temperature T in thermally equilibrium flows when T 1 ¼ T. In the frame of multi-temperature kinetic model [3] based on Treanor distributions, the set of macroscopic flow parameters for each binary mixture considered in the present paper includes five r; tÞ and unknown functions: number densities of molecules nm ð~ atoms na ð~ r; tÞ, temperature Tð~ r; tÞ, temperature of the first vibrational level T 1 ð~ r; tÞ, and gas velocity v ð~ r; tÞ. In the case of a onedimensional steady-state non-viscous flow behind a plane shock wave the governing equations take the form:

dðnm v Þ ¼ Rreact ; dx

ð29Þ

dðna v Þ ¼ 2Rreact ; dx

ð30Þ

dðqm W v Þ ¼ RW ; dx

ð31Þ ð32Þ

q0 v 20 þ p0 ¼ qv 2 þ p;

ð33Þ

v 20 2

¼hþ

v2 2

ha ¼

5 R ea Tþ ; 2 la ma

hm ¼

7 R 1 X Tþ ne; 2 lm qm i i i

ð35Þ

Y m ¼ qqm ; Y a ¼ qqa are the mass fractions of molecules and atoms, ea is

the atom formation energy, ma is the atom mass, lm and la are molar masses, ni are given by distributions (7), R is the universal gas constant. Right-hand sides in Eqs. (29) and (30) characterize variation of number densities of molecules and atoms due to dissociation and recombination reactions behind a shock wave and are written as follows:

Rreact ¼

 X  M M nM n2a krec ðTÞ  nm kdiss ðT; T 1 Þ :

ð36Þ

M M

Two-temperature dissociation rate coefficients kdiss ðT; T 1 Þ are given by Eq. (9), whereas the recombination rate coefficients M

krec ðTÞ are calculated using the detailed balance principle [3]. Eq. (31) describes variation of the number of vibrational quanta along the relaxation zone due to VT(TV) vibrational energy transitions and dissociation–recombination reactions. Following [3] (see also Section 3) we consider source terms in Eqs. (31) in the form:

RW ¼ RW;VT þ RW;react where the terms RW;VT and RW;react are responsible for VT transitions and reactions. It is worth noting that the most natural (but computationally expensive) way to calculate these terms avoiding uncertainty introduced by relaxation times is to use corresponding state-to-state production terms multiplied by the quantum number i. In this case, the production terms take the form

RW;VT ¼

 XX  M M M M inM ni1 ki1 i þ niþ1 kiþ1 i  ni ðki i1 þ ki iþ1 Þ ; M

RW;react ¼

ð37Þ

i

 XX  M M inM n2a ki;rec  ni ki;diss : M

ð38Þ

i

These equations contain state-dependent rate coefficients M

M

ki i1 ; ki iþ1 for VT(TV) vibrational energy transitions for N2 and O2 M

M

molecules and ki;diss ; ki;rec for dissociation and recombination processes. Coefficients of forward and backward processes are found using the detailed balance principle [3]. In calculations, we take into account only single-quantum vibrational energy transitions neglecting multi-quantum jumps as less probable. More computationally efficient but less accurate method for the calculation of VT production terms is to use the Landau–Teller Eq. (23) or its modification (24). In the next section we assess the limits of applicability for the simplified models. Eqs. (30)–(34) are solved for various conditions in a free stream. Level populations ahead a shock front are supposed to be equilibrium and described by Boltzmann distributions with the gas temperature. To avoid confusing, below we denote the temperatures before ð0Þ

q0 v 0 ¼ qv ;

h0 þ

h ¼ hm Y m þ ha Y a ;

:

ð34Þ

Here p; q are the gas pressure and density, x is the distance from the shock front, index ”0” specifies flow parameters in a free stream, h is the total enthalpy per unit mass:

the shock front as T ð0Þ ; T 1 , and immediately behind the shock ð1Þ

ð0Þ

ð0Þ

as T ð1Þ ; T 1 . For all considered cases, nN2 ¼ nO2 ¼ nð0Þ ¼ pð0Þ =kT ð0Þ nN

ð0Þ nO

ð0Þ

;

¼ ¼ 0. Flow parameters just behind a shock front in the considered two-temperature approach are found from the set of equations for mass, momentum and total energy conservation in the form (32)–(34) with frozen chemical reactions as well as conservation of vibrational quanta numbers:

  ð1Þ nð0Þ v ð0Þ WðT ð0Þ Þ ¼ nð1Þ v ð1Þ W T ð1Þ ; T 1 :

With the use of Eqs. (12) and (32) we can write:

9

E. Kustova et al. / Chemical Physics 464 (2016) 1–13

Z

ð0Þ

ðT Þ

i

ei kT

ð0Þ

 ¼

Z

vibr

X 1   i ð1Þ ð1Þ T ; T1 i

 exp 

ei  ie1 kT

ð1Þ



ie1 ð1Þ

kT 1

! :

ð39Þ

For harmonic oscillators T 1 ¼ T v and level populations in the relaxation zone have the form of Boltzmann distributions with vibrational temperature T v . Due to resonant VV exchanges between harmonic oscillators, level populations do not change within the shock front and just behind a shock are defined by the free stream gas temperature T ð0Þ . For calculation of flow parameters immediately behind a shock, in the conservation Eqs. (32)–(34) chemical reactions and vibrational distributions are assumed to be frozen. Below we present the comparison of shock heated gas flow parameters for considered mixtures calculated with the use of different models for rate coefficients of TV(VT) transitions and dissociation as well as on the basis of different vibrational distributions. 4.2. O2/O mixture.Validation against experimental data In this section, we study the post-shock flows of O2/O mixtures for the conditions corresponding to experiments carried out in [55]. Two test cases are considered: (1) Mð0Þ ¼ 9:3; v ð0Þ ¼ 3070 m/s, T ð0Þ ¼ 295 K, pð0Þ ¼ 2 Torr;

(2) Mð0Þ ¼ 13:4; v ð0Þ ¼ 4440 m/s, T ð0Þ ¼ 295 K, pð0Þ ¼ 0:8 Torr.

To validate our models, we compare the calculated vibrational temperature distributions with those measured in [55], see Fig. 9. The following five models for vibrational energy production terms have been used in calculations: (1) classical Landau–Teller model (LT), (2) modified LT model given by Eq. (24) with Millikan–White expression for the VT relaxation time (LT-m (MW)), (3) modification of LT model on the basis of the accurate kinetic theory Eq. (24) with Eq. (25) for the relaxation time (LT-m), (4) Eq. (37) with forced harmonic oscillator model proposed by Adamovich et al. [31] (FHO), (5) Eq. (37) with generalized Shwartz, Slawsky, Herzfeld analytical expressions [30] (SSH). For the dissociation rate, the two-temperature Treanor–Marrone (TM) model with constant U ¼ D=6k as well as its modification with variable Uði; TÞ given in Eq. (5) have been used. One can see that for M ¼ 9:3, the influence of dissociation model is rather weak whereas the model of VT relaxation is of importance. Implementation of the original and modified Landau–Teller models with the Millikan–White relaxation time as well as the SSH model result in very fast relaxation and higher value for the maximum of vibrational temperature. The best agreement for

4500

the location of T v maximum is given by the FHO and modified LT model with vibrational relaxation time given by Eq. (25). For M ¼ 9:3, the vibrational temperatures obtained using the FHO and modified LT model with correct relaxation time practically coincide regardless the dissociation model. For M ¼ 13:4, the influence of dissociation model is much more important. Again, SSH and LT models with the Millikan–White relaxation time significantly overpredict the rate of vibrational relaxation compared to experiments. The best agreement for the location of T v maximum is provided by the FHO + generalized Treanor–Marrone model with variable parameter U. Modified LT model with vibrational relaxation time given by Eq. (25) + generalized TM model also provides a good agreement to measured vibrational temperature. Thus we can conclude that the proposed modifications of the LT and TM models give a simple, accurate and numerically efficient alternative to the rigorous way of calculation of the production terms using the state-to-state rate coefficients of vibrational transitions and dissociation (Eqs. (37) and (38)). In Fig. 10 we plot the maximum vibrational temperature attained in the relaxation zone as a function of the temperature T ð1Þ just behind the shock front. Surprisingly, the best prediction of the T v maximum is given by the SSH model. Other models yield underpredicted values of the T v maximum. This conclusion is in contradiction with the previous one about the poor accuracy of the SSH model in the prediction of the vibrational relaxation rate. However it is worth noting that there is an uncertainty in the experimental detection of T v maximum value, and its location seems to be a better indicator for the comparison of measured and calculated vibrational temperatures. 7500

6500 6000

3500 5000

LT-m (M-W) LT-m, U(i,T), Eq.(5) LT-m, U=D/6k experiment

500 0.000

0.002

0.004

0.006

x, m

0.008

0.010

8000

9000

10000

11000

Fig. 10. Maximum of vibrational temperature in O2/O mixture as a function of temperature T ð1Þ . Comparison with experiments.

T1, K

T1, K

1000

FHO, U(i,T), Eq.(5) FHO, U=D/6k SSH LT

7000

(1)

3000

1500

6000

T ,K

(a)

M=9.3

5000

4000

3500

2000

5500

4500

4000

2500

FHO, U=D/6k FHO, U(i,T), Eq.(5) SSH LT LT-m experiment

7000

max Tv, K

1 vibr


X i exp 

7500 7000 6500 6000 5500 5000 4500 4000 3500 3000 2500 2000 1500 1000 500 0.000

(b) M=13.4

FHO, U(i,T), Eq.(5) FHO, U=D/6k SSH LT

0.002

0.004

LT-m (M-W) LT-m, U(i,T), Eq.(5) LT-m, U=D/6k experiment

0.006

0.008

0.010

x, m

Fig. 9. Vibrational temperature in O2/O mixture as a function of x. Comparison with experiments. (a) M = 9.3 and (b) M = 13.4.

10

E. Kustova et al. / Chemical Physics 464 (2016) 1–13

5500

4500

T, K

5000

T, K

(a)

M=9.3 FHO, U(i,T), Eq.(5) FHO, U=D/6k SSH LT LT-m (M-W) LT-m, U(i,T), Eq.(5) LT-m, U=D/6k

4000

3500 0.000

0.002

0.004

0.006

0.008

0.010

11000 10500 10000 9500 9000 8500 8000 7500 7000 6500 6000 5500 5000 4500 4000 0.000

(b)

M=13.4 FHO, U(i,T), Eq.(5) FHO, U=D/6k SSH LT LT-m (M-W) LT-m, U(i,T), Eq.(5) LT-m, U=D/6k

0.002

0.004

x, m

0.006

0.008

0.010

x, m

Fig. 11. Temperature for different models in O2/O mixture as a function of x. (a) M = 9.3 and (b) M = 13.4.

(a)

0.10 n /n O 0.08

0.4

M=13.4

0.06

0.3

M=9.3 FHO, U(i,T), Eq.(5) FHO, U=D/6k SSH LT LT-m (M-W) LT-m, U(i,T), Eq.(5) LT-m, U=D/6k

0.04

0.02

0.00 0.000

(b)

nO/n 0.5

0.002

0.004

0.006

0.008

0.010

FHO, U(i,T), Eq.(5) FHO, U=D/6k SSH LT LT-m (M-W) LT-m, U(i,T), Eq.(5) LT-m, U=D/6k

0.2 0.1 0.0 0.000

0.002

0.004

0.006

0.008

0.010

x, m

x, m

Fig. 12. Atom molar fractions in O2/O mixture as functions of x. (a) M = 9.3 and (b) M = 13.4.

Now let us discuss the variation of temperature and oxygen atom molar fraction along the relaxation zone shown in Figs. 11 and 12. We compare different models for energy transitions and dissociation considered above. It can be noticed that modified Treanor–Marrone dissociation model, FHO model for energy transitions or modification of the Landau–Teller model based on the kinetic theory (with Eq. (25) for relaxation time) provide close results for T and nO =n for M ¼ 9:3, however the difference increases for M ¼ 13:4. In the frame of these models we obtain higher values for the gas temperature and lower values for nO =n behind a shock compared to other models for both M numbers whereas SSH and Landau–Teller expressions in both classical and modified form with the Millikan–White formula for the relaxation time provide noticeably lower values for the gas temperature and higher values for atom molar fractions thus overpredicting the rate of dissociation. In the case M ¼ 9:3 dissociation is rather slight and vibrational relaxation plays a dominating role. Therefore a dissociation model occurs not so important and results obtained in the frame of different dissociation models appear to be close to each other whereas quite different values are obtained with the use of various models for energy transitions. For M ¼ 13:4, the choice of dissociation model becomes important: noticeable distinctions between the values obtained using the original TM model and its modification with U ¼ Uði; TÞ are found. The difference between the temperature and molar fraction values found with the use of various models decreases with x.

4.3. N2/N mixture In this section, we study the shock heated flow of N2/N mixture using the models discussed above. The free stream conditions are the following: M ð0Þ ¼ 13; T ð0Þ ¼ 271 K, pð0Þ ¼ 100 Pa. First we compare the results obtained using the anharmonic oscillator model based on the Treanor vibrational distributions (7) and those for harmonic oscillators and non-equilibrium Boltzmann distribution. The use of harmonic oscillator vibrational energies makes it possible to simplify considerably kinetic models and governing equations. Estimates for anharmonic effects under considered conditions are shown in Figs. 13–15(a). Figs. 13(a), 14(a) and 15(a) demonstrate the influence of anharmonicity of molecular vibrations on the gas temperature T, vibrational temperature T 1 and nitrogen atom molar fraction in the relaxation zone. Neglecting anharmonic effects leads to slightly underpredicted T 1 values as well as overpredicted T and, as a consequence, higher atom molar fractions. The difference between values found in the considered test case for anharmonic oscillators and neglecting anharmonicity does not exceed 1.5% for T, 5% for T 1 but can reach higher values for species molar fractions. Since anharmonic effects are weak, in the following discussion we limit our consideration to harmonic oscillators. Figs. 13(b) and 14(b) show the impact of models for dissociation and vibrational energy transitions proposed in Sections 2 and 3 on temperatures T and T v for harmonic oscillators. In this case very close temperature values are obtained using TM dissociation model

11

E. Kustova et al. / Chemical Physics 464 (2016) 1–13

and its modification TM-m with U ¼ Uði; TÞ with FHO model for vibrational transitions. At the same time, LT and LTm (M-W) models with TM expressions provide a little bit lower T values. For vibrational temperature T v in the beginning of the relaxation zone the LT model or its modification LTm (M-W) with TM dissociation model give overestimated values compared to those found with the use of FHO model; different dissociation models do not affect the vibrational temperature calculated using the FHO model, corresponding curves cannot be distinguished in Fig. 14(b). For

x > 0:01 m all four considered models provide rather close T v values. Fig. 15(b) depicts the variation of N atom molar fraction for harmonic oscillators found in the frame of the same models as in Figs. 13(b) and 14(b). Again, FHO model with TM and TMm Uði; TÞ expressions give rather close values for molar fractions. LT and LTm (M-W) models lead to higher nN =n values. We would like to emphasize that for nitrogen, the impact of dissociation model is found to be rather weak for M ¼ 13. However we expect more important effects at higher Mach numbers.

9000

harm. osc. anh. osc. 8500

FHO, U(i,T), Eq.(5) FHO, U=D/6k LT, U=D/6k LT-m (M-W), U=D/6k

8500

T, K

T, K

(b)

9500

(a)

9000

8000

8000 7500

7500

7000

0.00

0.01

0.02

0.00

0.01

0.02

0.03

x, m

0.04

0.05

x, m

Fig. 13. Temperature distribution in N2/N mixture. (a) Comparison of harmonic (Boltzmann distribution) and anharmonic (Treanor distribution) oscillator models. (b) Harmonic oscillator, different VT relaxation and dissociation models.

(a)

7500 7000 6500 6000 5000

T1, K

T1, K

5500

harm. osc. anh. osc.

4500 4000 3500 3000 2500 2000 0.00

0.01

0.02

7500 7000 6500 6000 5500 5000 4500 4000 3500 3000 2500 2000 1500 1000 500 0.00

(b)

FHO, U(i,T), Eq.(5) FHO, U=D/6k LT, U=D/6k LT-m (M-W), U=D/6k

0.01

0.02

x, m

0.03

0.04

0.05

x, m

Fig. 14. Vibrational temperature distribution in N2/N mixture. (a) Comparison of harmonic (Boltzmann distribution) and anharmonic (Treanor distribution) oscillator models. (b) Harmonic oscillator, different VT relaxation and dissociation models.

(a)

0.05

(b)

0.06

0.04 0.04

nN/n

nN/n

0.03

0.02

0.01

0.00 0.00

0.01

0.02

0.03

x, m

0.04

FHO, U(i,T), Eq.(5) FHO, U=D/6k LT, U=D/6k LT-m (M-W), U=D/6k

0.02

harm. osc. anh. osc.

0.05

0.00 0.00

0.01

0.02

0.03

0.04

0.05

x, m

Fig. 15. Atom molar fraction distribution in N2/N mixture. (a) Comparison of harmonic (Boltzmann distribution) and anharmonic (Treanor distribution) oscillator models. (b) Harmonic oscillator, different VT relaxation and dissociation models.

12

E. Kustova et al. / Chemical Physics 464 (2016) 1–13

5. Conclusions Simple self-consistent, efficient and accurate models of coupled vibrational relaxation and dissociation for multi-temperature flow simulations are proposed on the basis of the kinetic theory methods. The improved Treanor–Marrone model taking into account statespecific parameter Uði; TÞ provides very good agreement with the state-to-state dissociation rate coefficients obtained in the frame of molecular dynamics methods. Two-temperature dissociation rate coefficients are calculated using the new parameter on the basis of the Treanor and Boltzmann distributions. For oxygen, anharmonic effects are significant in the whole temperature range, however the discrepancy decreases with temperature; at low temperatures dissociation rate coefficients for anharmonic oscillators show essentially non-monotonic behavior contrarily to the case of harmonic oscillators. For nitrogen, the impact of anharmonicity is weak; it become noticeable only for low temperature T and high vibrational temperature. Comparison of the two-temperature dissociation rate coefficients calculated using the constant and varying parameter U demonstrates considerable difference, especially at high temperature; the discrepancy reaches two orders of magnitude in N2. While for nitrogen using the constant parameter yields overpredicted kdiss ðT; T v Þ values, for oxygen the opposite trend is found. The generalization of the Landau–Teller model for the rate of VT relaxation is proposed on the basis of the kinetic theory. The improved model includes the ratio T=T v missing in the original model developed for the case T  T v and the relaxation time introduced rigorously as integral of the vibrational energy variation depending on both T and T v . Comparison of the relaxation rates calculated using the original and new models demonstrates significant difference at high temperature and big ratio T=T v . The developed models are validated using the experimental measurements of vibrational temperature in shock heated oxygen flows. It is shown that combination of generalized Treanor–Marrone and Landau–Teller models (with correct expression for the relaxation time) provides a very good agreement with experimental T v distributions. Classical LT and SSH models overestimate considerably the rate of VT relaxation. At M ¼ 9:3, the impact of dissociation model is weak whereas at M ¼ 13:4 it plays an important role in the vibrational kinetics. For nitrogen at M ¼ 13 the effect of dissociation model is weak but we expect stronger influence for higher Mach numbers. Finally we conclude that the developed models represent an accurate and numerically efficient alternative to the rigorous (but time and resource consuming) method of the production terms calculation based on the state-to-state rate coefficients of vibrational transitions and dissociation. Since the models are derived on the basis of the kinetic theory and molecular dynamics, their implementation is not restricted to shock heated flows only, and the models can be applied for arbitrary non-equilibrium flows such as nozzle expansions and boundary layers. Conflict of interest There is no conflict of interest. Acknowledgments This study is supported by the Russian Science Foundation (Project 15-19-30016). References [1] E. Montroll, K. Shuler, J. Chem. Phys. 26 (3) (1957) 454. [2] M. Capitelli, C. Ferreira, B. Gordiets, A. Osipov, Plasma kinetics in atmospheric gases, Springer Series on Atomic, Optical and Plasma Physics, vol. 31, SpringerVerlag, Berlin, 2000.

[3] E. Nagnibeda, E. Kustova, Nonequilibrium Reacting Gas Flows, Kinetic Theory of Transport and Relaxation Processes, Springer-Verlag, Berlin, Heidelberg, 2009. [4] I. Adamovich, S. Macheret, J. Rich, C. Treanor, J. Thermophys. Heat Transfer 12 (1) (1998) 57. [5] S. Ruffin, C. Park, Vibrational relaxation of anharmonic oscillators in expanding flows, AIAA Paper 92-0806. [6] G. Colonna, M. Capitelli, M. Tuttafesta, D. Giordano, J. Thermophys. Heat Transfer 13 (3) (1999) 372. [7] I. Armenise, M. Capitelli, G. Colonna, C. Gorse, J. Thermophys. Heat Transfer 10 (3) (1996) 397. [8] E. Kustova, E. Nagnibeda, I. Armenise, M. Capitelli, J. Thermophys. Heat Transfer 16 (2) (2002) 238. [9] A. Orsini, P. Rini, V. Taviani, D. Fletcher, E. Kustova, E. Nagnibeda, J. Thermophys. Heat Transfer 22 (3) (2008) 390. [10] J. Kim, I. Boyd, Chem. Phys. 415 (2013) 237. [11] O. Kunova, E. Nagnibeda, Chem. Phys. 441 (2014) 66. [12] M. Panesi, A. Munafò, T.E. Magin, R.L. Jaffe, Phys. Rev. E 90 (2014) 013009. [13] D. Giordano, V. Bellucci, G. Colonna, M. Capitelli, I. Armenise, C. Bruno, J. Thermophys. Heat Transfer 11 (1) (1997) 27. [14] E. Josyula, J. Burt, E. Kustova, P. Vedula, Influence of state-to-state transport coefficients on surface heat transfer in hypersonic flows, AIAA Paper 20140864, doi:10.2514/6.2014-0864. [15] E. Josyula, J. Burt, E. Kustova, P. Vedula, M. Mekhonoshina, State-to-state kinetic modeling of dissociating and radiating hypersonic flows, AIAA Paper 2015-0475, doi:10.2514/6.2015-0475. [16] E. Kustova, E. Nagnibeda, Chem. Phys. 233 (1998) 57. [17] I. Armenise, E. Kustova, Chem. Phys. 428 (2014) 90. [18] L. Cutrone, M. Tuttafesta, M. Capitelli, A. Schettino, G. Pascazio, G. Colonna, 3D nozzle flow simulations including state-to-state kinetics calculation, in: J. Fan (Ed.), Rarefied Gas Dynamics, vol. 1628 of AIP Conference Proceedings, 2014, pp. 1154–1161. [19] R. Brun, Transport et relaxation dans les écoulementsgazeux, Masson, Paris, New York, Barselone, Milan, Mexico, Saõ Paulo, 1986. [20] C. Park, Nonequilibrium Hypersonic Aerothermodynamics, J. Wiley and Sons, New York, Chichester, Brisbane, Toronto, Singapore, 1990. [21] C. Treanor, I. Rich, R. Rehm, J. Chem. Phys. 48 (1968) 1798. [22] A. Chikhaoui, J. Dudon, E. Kustova, E. Nagnibeda, Physica A 247 (1–4) (1997) 526. [23] A. Chikhaoui, J. Dudon, S. Genieys, E. Kustova, E. Nagnibeda, Phys. Fluids 12 (1) (2000) 220. [24] E. Kustova, E. Nagnibeda, Chem. Phys. 208 (3) (1996) 313. [25] A. Chikhaoui, E. Nagnibeda, E. Kustova, T. Alexandrova, Chem. Phys. 263 (2001) 111. [26] E. Kustova, E. Nagnibeda, T. Alexandrova, A. Chikhaoui, Chem. Phys. 276 (2) (2002) 139. [27] G. Colonna, I. Armenise, D. Bruno, M. Capitelli, J. Thermophys. Heat Transfer 20 (3) (2006) 477. [28] G. Colonna, L.D. Pietanza, M. Capitelli, J. Thermophys. Heat Transfer 22 (3) (2008) 399. [29] A. Guy, A. Bourdon, M. Perrin, Chem. Phys. 420 (2013) 15. [30] R. Schwartz, Z. Slawsky, K. Herzfeld, J. Chem. Phys. 20 (1952) 1591. [31] I. Adamovich, AIAA J. 39 (10) (2001) 1916. [32] G. Billing, E. Fisher, Chem. Phys. 43 (1979) 395. [33] G. Billing, R. Kolesnick, Chem. Phys. Lett. 200 (4) (1992) 382. [34] F. Esposito, M. Capitelli, C. Gorse, Chem. Phys. 257 (2000) 193. [35] F. Esposito, I. Armenise, M. Capitelli, Chem. Phys. 331 (1) (2006) 1. [36] I. Armenise, F. Esposito, M. Capitelli, Chem. Phys. 336 (1) (2007) 83. [37] F. Esposito, M. Capitelli, Chem. Phys. Lett. 443 (2007) 222. [38] F. Esposito, I. Armenise, G. Capitta, M. Capitelli, Chem. Phys. 351 (1–3) (2008) 91. [39] I. Armenise, F. Esposito, Chem. Phys. 398 (2012) 104. [40] Planetary entry integrated models, . [41] R. Jaffe, D. Schwenke, G. Chaban, Vibration–rotation excitation and dissociation in N2–N2 collisions from accurate theoretical calculations, in: Proc. 10th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, June 2010, 2010, AIAA Paper 2010-4517. [42] L. Landau, E. Teller, Phys. Z. Sowjetunion 10 (1936) 34. [43] R. Millikan, D. White, J. Chem. Phys. 39 (1963) 3209. [44] E. Kustova, G.M. Kremer, Chem. Phys. 445 (2014) 82. [45] E. Kustova, G. Oblapenko, Phys. Fluids 27 (2015) 016102. [46] P. Marrone, C. Treanor, Phys. Fluids 6 (9) (1963) 1215. [47] N.M. Kuznetsov, Kinetics of Monomolecular Reactions, Nauka, Moscow, 1982. in Russian. [48] S.O. Macheret, A.A. Fridman, I.V. Adamovich, J.W. Rich, C.E. Treanor, Mechanisms of nonequilibrium dissociation of diatomic molecules, AIAA Paper 94-1984. [49] N. Belouaggadia, R. Brun, J. Thermophys. Heat Transfer 12 (4) (1998) 482. [50] G. Chernyi, S. Losev, S. Macheret, B. Potapkin, Physical and Chemical Processes in Gas Dynamics, vol. 1, 2, American Institute of Aeronautics and Astronautics, 2004. [51] O. Knab, H. Frühauf, E. Messerschmid, J. Thermophys. Heat Transfer 9 (2) (1995) 219. [52] E. Kustova, E. Nagnibeda, T. Alexandrova, A. Chikhaoui, Chem. Phys. Lett. 377 (5–6) (2003) 663.

E. Kustova et al. / Chemical Physics 464 (2016) 1–13 [53] O. Kunova, E. Nagnibeda, Chem. Phys. Lett. 625 (2015) 121. [54] F. Esposito, M. Capitelli, E. Kustova, E. Nagnibeda, Chem. Phys. Lett. 330 (2000) 207. [55] L.B. Ibraguimova, A.L. Sergievskaya, V.Y. Levashov, O.P. Shatalov, Y.V. Tunik, I.E. Zabelinskii, J. Chem. Phys. 139 (2013) 034317. [56] A.S. Savel’ev, E. Kustova, Limits of applicability of the Treanor–Marrone model for N2 and O2 state-to-state dissociation rate coefficients, Vestnik St. Petersburg. Universiteta, Math., Mech., Astr. 2 60 (2) (2015) 268–279.

13

[57] T.J. Scanlon, C. White, M.K. Borg, R.C. Palharini, E. Farbar, I.D. Boyd, J.M. Reese, R.E. Brown, AIAA J. 53 (6) (2015) 1670–1680. [58] K. Koura, H. Matsumoto, Phys. Fluids A 3 (10) (1991) 2459. [59] K. Koura, H. Matsumoto, Phys. Fluids A 4 (5) (1992) 1083–1085. [60] Y. Stupochenko, S. Losev, A. Osipov, Relaxation in Shock Waves, SpringerVerlag, Berlin, Heidelberg, New York, 1967. [61] R. Brun (Ed.), High Temperature Phenomena in Shock Waves, Shock Wave Science and Technology Reference Library, vol. 7, Springer, Berlin, 2012.