Chemical Physics ELSEVIER
Chemical Physics 195 (1995) 83 91
Vibrational relaxation models for dilute shock heated gases David A. Gonzales a,*, Philip L. Varghese b Department of Aeronautics and Astronautics, Massachusetts Institute of Technoh)gy, Cambridge, MA 02139, USA b Department of Aerospace Engineering and Engineering Mechanics. The University of Texas at Austin, Austin, TX 78712, USA
Received 6 July 1994; in final form 6 February 1995
Abstract This work presents a computational examination of two models of vibrational relaxation for coupled vibration-dissociation-recombination processes. These models are commonly used in state-of-the-art non-equilibrium flow codes. A shock heated gaseous mixture of a single diatomic species highly dilute in a rare gas heat bath is considered. The master equation is linearized about final equilibrium and provides remarkable accuracy to early simulation times. The computations show a two-step relaxation of the vibrational energy. This behavior is not predicted by a linear relaxation rate equation. An exact expression for the rate of change of vibrational energy is presented and it is observed that both the linear relaxation rate equation and the Park diffusion model do not accurately reproduce the computed results.
1. Introduction
In computational fluid dynamic simulations of vibrationally relaxing, chemically reacting flows of diatomic gases, two models are often used to calculate rates of vibrational relaxation. The first is the linear relaxation rate (LRR) equation. (The ERR equation is often referred to as the Landau-Teller equation in the literature.) The other model, given by Park [1], attempts to simulate the assumed diffusive nature of vibrational relaxation [2]. In this article we examine the accuracy of these relations for a gas consisting of anharmonic oscillators highly dilute in a rare gas. Large amounts of kinetic data have been obtained for gases at high temperatures from shock tube experiments. Vibrational relaxation times and rate coef-
* Corresponding author.
ficients for dissociation have been determined from these experiments and have been used in non-equilibrium flowfield solvers. In addition, these data have been used to calibrate or validate theoretical models of microscopic phenomena in shock heated gases. In shock tube experiments, the test gas is frequently a mixture consisting of the diatomic species of interest highly dilute in a rare gas that acts as a heat bath during the experiment. One can numerically simulate this case in a relatively efficient manner and obtain time dependent solutions of the master equation, e.g., Refs. [2-26]. A linear form of the master equation can be obtained by assuming (i) diatom-diatom collisions are rare and (ii) recombination reactions are negligible [27-30]. We refer to this as the one-way dissociation process. The former assumption is reasonable for highly dilute mixtures while the latter is reasonable only near the shock since recombination becomes important as the diatomic species approach equilibrium further down-
0301-0104/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 030 1-01 0 4 ( 9 5 ) 0 0 0 7 8 - X
D.A. Gonzales, P.L. Varghese/ Chemical Physics 195 (1995) 83-91
84
stream. These approximations lead to an exact solution and have been used to obtain vibrational state populations using efficient eigenvector-eigenvalue routines, e.g., Refs. [14,15,21,23-26]. Simulations of coupled internal state-dissociation-recombination processes of H 2 highly dilute in a rare gas have been performed extensively by Pritchard and co-workers, e.g., Refs. [3-6,9,10,13, 17]. Their master equation studies not only generated a wealth of fundamental insights, but also advanced several numerical techniques for solving the master equation; these include matrix methods and asymptotic expansion techniques, see e.g., Refs. [10,17]. An approximate method to solve the non-linear master equation was suggested by Truhlar and coworkers [18,19] and also by Schwenke [22] to arrive at a linear master equation. In this method the master equation is expanded about equilibrium to terms that are first order in the molecular or atomic concentrations. The resulting linear system of equations can thus be solved by standard matrix routines. Recently we simulated the shock heating of a dilute diatomic gas to arrive at a new model for state-specific dissociation [23-26] by solving the linear vibrational master equation for one-way dissociation. The gas mixtures consisted of N 2, O z, CO or H 2 highly dilute in At. Shock heating was simulated for temperatures in the range 6000-15000 K. In the present work, we consider the shock heating of N 2 highly dilute in Ar for coupled vibration-dissociation-recombination (VDR) processes. The approximate method of Truhlar and co-workers and Schwenke is applied to solve the non-linear master equation. As will be shown below, these approximations provide an excellent reproduction of the relaxation process for times far from (final) equilibrium.
dn a dt
-
('bl t, p
where n, is the number density of N~ in the vibrational state ~, n, is the number density of nitrogen atoms N, and the summations are over all bound vibrational levels (t, = 0 to t' = 59). 2.1. Rate coefficients In Eq. (1) the coefficients k are state-specific pseudo first order rate coefficients for the transition denoted by the appropriate subscripts (c denotes the continuum), i.e., ~:ij kijP, where k,j is a statespecific rate coefficient for the corresponding transition and p is the argon density, assumed to be 1.28 × 1017 molecules cm 3. The set {k,,,} consists of rotationally averaged rate coefficients determined from the surprisal synthesis method of Procaccia and Levine [31]. The rotationally averaged rate coefficient for dissociation from specific vibrational states, k, c, is given by an expression we developed [23,25, 26] using information theoretic ideas [14,31,32], assuming a non-linear (in general) surprisal for dissociation: =
Nm~(v )
k,.c(T)
=
E
yv(L,; r ) k , , v c ( T ) ,
where
Xp ..... + A d .....
Xexp[-A (J YN(C; T) = [Q,ot(t,; T)]
2. Rate equations
dt
-
z
,,~,
],
'g~(N)(2N+
×exp(--E,N/KT),
This study focuses on highly dilute mixtures of N 2 in Ar for which N 2 - N 2 collisions are neglected. The resulting master equation for coupled VDR processes is well known and is given by dn I
(2)
N - Nmi n
-
2
no,
(la)
Orot(l:; T ) =
Nm~x(V) E gs(N)( 2N+ N
-
(3) 1) (4)
1)
Nmi n
Xexp( - - E , . N / K T ) .
(5)
Here, Ad = / 3 ( D 0 - E ' , N ) , Adm,~ = ~ ( D o - E,um,, ), = I / K T , Z = O-Hs(8KT/'utz) 1/2, and O-HS= ~r(d i + dAr)2/4. N is the rotational quantum number in the vibrational manifold v, YN(V; T) is the equilib-
D.A. Gonzales, P.L. Varghese/ Chemical Physics 195 (1995) 83-91 rium fraction of molecules in state (c, N) at the temperature T, E,N is the energy of state (v, N), D~ is the dissociation energy, g~(N) is the nuclear spin degeneracy, Xp and Xpmax a r e variables which arise from thermal averaging approximations [25], Z is the hard sphere rate coefficient, crHS is the hard sphere cross section, and d~ and dAr are the effective diameters of the diatom i and the Ar atom, respectively. The reduced mass of the diatom-Ar pair is /,t and K is the Boltzmann constant. The subscript " m a x " refers to the highest bound, but not quasibound, rotational level in a given vibrational manifold. The surprisal parameter for dissociation AD is assumed to be temperature independent and is set to unity in this study. This value for AD gave the best fit of master equation results to available shock tube data for mixtures of N 2 highly dilute in Ar [26]. The function e ( c ) is required for closure in the derivation of Eq. (3) and represents the dissociation rate for transitions from state (c, Nm~×). We approximate this function by [25,26] ~:(t,) = (1
--
O'0/O'HS)(U//2max) 2 q- O'0/O'Hs ,
(6)
where o-0 is the cross section for dissociative transitions from state (0, Nm, ×) and Cmax is the quantum number of the highest lying vibrational level. We assume tr 0 is small for N e - A r collisions and make the arbitrary assignment [25] ¢r0 = 0 . 0 1 X2. The product e ( v ) Z biases the dissociative rate coefficients to high lying internal states rather than to high lying vibrational states, as is done in other models [8,11,14,33]. Eq. (6) was motivated by the quasiclassical trajectory studies of Haug et al. [21] for H 2 dilute in At. Their study showed t, to be more effective than N in causing dissociation when E N = D 0. Eq. (3) provided good agreement [25] when compared to the state-specific dissociative data published by Haug et al., differing by as little as 10% and by as much as a factor of 4, with an average deviation of 74%. This performance is considered encouraging for engineering purposes. It should be mentioned that closed form expressions were developed by Pritchard and co-workers to compute probabilities for dissociative transitions from specific vibrational [6] and vibrational-rotational [9] states. However these empirically based models were derived for H z - H e mixtures only and therefore are not appropriate for studies of other diatomic systems.
85
Finally, the state-specific rate coefficient for recombination k,,,. is determined from detailed balance assuming a Maxwellian velocity distribution characterized by the heat bath temperature.
2.2. Linear master equation We perform the typical, e.g., Refs. [10,17,34-36], shock heating simulation of a spatially homogeneous dilute thermal gas at the initial temperature T~. The N 2 gas is contained in an imaginary isothermal box containing Ar at the temperature T. The simulation begins at t = 0 +, when the N 2 translational-rotational temperature Ttr is raised to the bath gas temperature T, and Tvi b = Ti, where Tvih is the vibrational temperature. The simulation models the approach to vibrational equilibrium and ends when T~ib = T. Eq. (1) can be linearized by expanding about equilibrium and keeping terms that are first order in n~ and n a [18,19,22], i.e.,
dn i
dn i
~ O(dni/dt) + Z, d t n~ j= o Onj
dt
ne
i, j = 0, 1 . . . . . 59, a,
(7)
where in general n e is the equilibrium number density of state s (including atoms) and n ~ = {n]}. Note that Eq. (7) includes all vibrational levels c and the atomic state a. We normalize Eqs. (1) and (5) in terms of the variable wi, n i --
wi -
-
n~
(8)
-
Similar normalized variables have been used by other authors, e.g. Refs. [2,4,19,22]. Substituting Eq. (8) into Eq. (7), one obtains
dwi.~ ~ dt
j=(I
O(dwi/dt) aWj
w ~Wj,
i=0,
1 . . . . . 59, a, (9)
where w c implies {w~ = 0}. To determine [a(dwi/ dt)/Owj],.~, one substitutes Eq. (8) into Eq. (1), performs the appropriate differentiation, and then evaluates at equilibrium. After rearranging, Eq. (9) can be written in matrix form as d w / d t = A w , or dw i ~- -
~ A i j w j, j - 0
i=0,
1 ,. . . , 5 9 , a,
(10a)
D.A. Gonzales, P.L. Varghese / Chemical Physics 195 (1995) 83-91
86
e
e
3.1. Vibrational population distributions
-
Aij = ( n j / n i )kji(1 - 6ij)(1 - 6jo) i4=a, (lOb)
Aij = 2(n~/n~)[~jc(l - 6j~) - 4n e ~
~:~.,,,6j~,
Ut g=a
i = a,
(10c)
where 3ij is the Kronecker delta. The general solution to Eq. (10) is
wi(t ) = Y'amUim exp(t0mt),
(11)
m
where (.Dmand U,m are the eigenvalues and eigenvectors respectively that satisfy Au
m =
(.DmUm,
01 m
<~O.
(12)
The constants of integration a m are determined from the initial conditions, which we take to be equilibrium conditions at the initial temperature T i. Thus we solve,
wi(O ) = EamUim.
(13)
m
Time dependent vibrational probabilities are then given by
n~,[1 + w,,(t)] P , ( t ) = Y'~n~[1 + w , ( t ) ] '
Fig. 1 shows unit normalized vibrational populations for N 2 - A r collisions with Ti = 4000 K and T = 15000 K. We have chosen this value for T~ in order to compare with our previous results [23]. A high value for T~ was required in that work due to the well known ill-conditioning of the resulting eigenvalue equation for low temperatures (e.g. 300 < Ti < 3000 K). These complications have been discussed in the literature, e.g., Refs. [4,10,13,17] and an efficient method to overcome the ill-conditioning problem was developed by McElwain and Pritchard [6,10]. These problems do not exist in the present formulation (Eq. (10)) and initial temperatures as low as 300 K may be used. In Fig. 1 the dashed line represents the initial Boltzmann distribution, P,!. The distribution labeled 'ss' is the steady state condition and corresponds to t = r~ = 21.9 Ixs. This time corresponds to the final state one obtains for the one-way dissociation case described earlier and is defined as r~ = ] / ~ 1 ] 1 where ] Xj ] 4:0 and is the eigenvalue of the one-way dissociation rate coefficient matrix of smallest absolute value [14]. In contrast, the smallest eigenvalue of the rate coefficient matrix Aij in Eq. (10) must necessarily be zero so that the system returns to
(14) '
~
'
~
'
'
'
I
'
'
'
I
'
'
r
[
,
,
,
u
10-1
while time dependent number densities for atoms and specific molecular states are computed by inverting Eq. (8).
3. Results and discussion
10-3 .g
o e~
104
ss
10 -7 10 -9
In order to solve Eq. (10), a low temperature gas composition must be specified. In this work we assume 5% N 2 at 300 K. This value was chosen to reflect the concentrations studied experimentally by Appleton et al. [37] who examined 2, 3, 5, 10 and 20% N 2 mixtures in Ar. Their results in the temperature range 13000-15000 K are primarily for 5% N 2 mixtures. The energies E,. u w e r e obtained by solving the radial Schr/Sdinger equation using the ground-state potential given by Huxley and Murrell
[381.
pui
10"11 10-13 0.0
-- "~
, L , L , J , , ~ , I J , , I , , 0.2 0.4 0.6 0.8
1.0
E•/D o Fig. 1. T i m e e v o l u t i o n o f n o r m a l i z e d vibrational populations for N 2 collisions in Ar. T i = 4 0 0 0 K, T = 15000 K. D a s h e d line: initial ( B o l t z m a n n ) distribution, labeled P~i. Distribution labeled ' s s ' c o r r e s p o n d s to t = r~. Solid lines: p r e v i o u s w o r k [23] at "r~, 10 2z,~, 10 3r~, and 10-4%~; circles: present results at times c o r r e s p o n d i n g to p r e v i o u s work; squares: present results at lO.3T~; triangles: present results at 34.27~.
D.A. Gonzales, P.L. Varghese / Chemical Physics 195 (1995) 83-91
equilibrium at long times [10,36]. This requirement has been confirmed by our calculations: we observed ] (211 [ = O(10-a), ] O)2 [ = O ( 1 0 4) and [ O)max [ = O(108), where ]toll < [ w21 < ... < [ (.Omax [. These values are similar to those obtained by Pritchard and co-workers in their early studies of H 2 [10]. ( W e computed ] o) 1 ] = O(10 20) when using the double precision option on a Cray Y-MP.) The solid lines in Fig. 1 represent distributions for one-way dissociation reproduced from Ref. [23] and correspond to four different times: 10 4']'SS, 10 3,r~s, 10-2zss, and ~-ss. (The distribution at t = 10- 1~. was essentially identical to that at t = T~ and has been omitted for clarity.) Distributions for the same times were computed using Eq. (10) and are shown as open circles in Fig. 1. The distributions given by the squares and triangles also were computed using Eq. (10) and correspond to t = 10.3%~ and t = 34.2z~s, respectively. We identify the longest time shown in Fig. 1 as that corresponding to equilibrium conditions, i.e., teq = 34.2%~ (750 Ixs). In Fig. 1 we observe excellent agreement with our previous calculations for t = 0.01Z~s and t = r~s, with each distribution having a maximum error of 4% and 3%, respectively. For the earliest times, we observe excellent agreement for all vibrational levels lying below 80% of the dissociation threshold. However, significant departures do occur at high lying vibrational levels. At these early times, populations could not be computed for the highest lying levels due to
'"""1
~ ''"'"1
0.14
' ''"'"l
.......
/
'1
'
~'"'"1
' ''""1
'
''"'"
,
numerical instabilities and thus are not presented. These instabilities are, of course, a consequence of the linearization of Eq. (10). In any case, it is surprising that a master equation derived from an expansion about final equilibrium effectively reproduces distributions very close to the initial (equilibrium) condition. Qualitatively, one observes expected behavior: smooth relaxation to equilibrium, with overpopulation of the high lying states at early times [5,6,39] followed by underpopulation of the high lying states at later times [5-7,21]. 3.2. Temporal evolution o f average vibrational energy
In Fig. 2 we present relaxation profiles for the average vibrational energy per molecule V(t) of the N 2 gas, where
v(t)
= EP,.(t)E,.
(15)
L'
P , ( t ) is given by Eq. (9) and E,. is the energy of vibrational state v. In Fig. 2 we normalize V ( t ) by
the N 2 dissociation energy D o = 79899.7 cm-1. The present results are compared to those obtained by the LRR model and the one-way dissociation computations of Ref. [23]. The LRR model is given by
dr(t)
V° - v ( t ) -
dt
,
0.10
0.06
. . ~
+
ref. 23
o
present
0.02 .......i ........ i ........ i ........ i ........ , ........ i ..... 10-5 10-3 10-1 101 t ] 'Css Fig. 2. Average vibrational energy normalized by dissociation energy D o . Solid line: linear relaxation rate; crosses: previous work; circles: present results.
(16)
Tvib
where V e is the equilibrium vibrational energy and Tvib is the vibrational relaxation time. We compute Tvib by the formula of Millikan and White [40]. One can integrate Eq. (16) to get V ( t ) = V e + [V(0) - V e] exp(--t/Zvib) ,
...
87
(17)
where V(0) corresponds to the vibrational energy at 4000 K. Eq. (17) is shown as a solid line in Fig. 2. The crosses are the one-way dissociation results and the circles are those computed by the present work. Here, as in Fig. 1, excellent agreement is observed between the present results and those of Ref. [23]. However, in Fig. 2 the excellent agreement extends to times much earlier than those presented in Fig. 1. The earliest time we show in Fig. 2 is 0.310 ns (1.42 × 10 5zss); the final time is teq (given above). We observe agreement at these early times because the contribution to the average vibrational energy is
88
D.A. Gonzales, P.L. Varghese / Chemical Physics 195 (1995) 83-91
weighted toward the low and mid lying vibrational levels that are relatively more populated. Thus, the errors in the computation of high ~, populations at very early times (see Fig. 1) do not contribute significantly to the sum in Eq. (15). In Fig. 2 one clearly sees a two-step relaxation process, first to the steady state value and then to equilibrium. Similar behavior was observed by Schwenke [22] in his master equation study of H 2, although to a lesser degree than what we show. The LRR model predicts pure exponential relaxation and thus does not capture this behavior. In the steady state region, the LRR model overpredicts the master equation results by 13%. However, it has been observed [24] that if one replaces V e with the steady state value V ~ in Eqs. (16) and (17), one obtains an excellent reproduction of the master equation results for t ~< r,~. This implies that the approach to steady state in Fig. 2 is essentially exponential. In any case, one might consider the LRR model a reasonable approximation for the VDR process we consider.
3.3. Temporal euolution of t,ibrational relaxation rate For VDR processes in dilute gas mixtures, we show in the Appendix that the instantaneous vibrational relaxation rate is given by
dV
/
dt - (~'> + V~:d- 7~,uk'] na
+
E ~ , - (e,~:,,> ,
(18)
where
<. = Y'~ (E,.,-E,.)k:,.,,,
(19)
I'
led(t) = Y'.P,(t)~:,,.,
(20)
t:
~ ( T ) = E~:c,,,
(21)
l'
n vt(t) = Y'~n,(t).
(22)
L,
The brackets ( . . - ) refer to time dependent averages over P,(t). Eq. (18) is an exact equation for dV/dt for coupled VDR processes of a single di-
atomic species dilute in a rare gas heat bath. At equilibrium, Eq. (18) requires ( S , } = 0 since the terms inside each of the square brackets sum identically to zero. Further, in the absence of dissociation and recombination, Eq. (18) reduces to dV -
~P,(t)S,..
dt
(23)
t
Eq. (23) can be related to the LRR equation by substitution of the following " s u m rule" [17,31,36, 41-43], S,
V ~- E , -
-
-
(24)
%db
Eq. (24) is a well known necessary and sufficient condition for pure exponential relaxation of the average vibrational energy [17,31,36,41-43]. We compute dV/dt from the output of the master equation using Eq. (18) and show these results in Fig. 3. We also present results using the LRR model of Eq. (16) and the recent expression suggested by Park [1]. The Park relation is a modification of the LRR model to account for the assumed diffusive nature of vibrational energy transfer in VDR processes [2] and is given by
dt - \
re
T~--;r~hk
'
(25)
r e = r v w + rc,
(26)
re = 1/(ncO-(U}) ,
(27)
s = 3.5 exp(--5000/T~hk).
(28)
The subscript 'shk' corresponds to conditions immediately downstream of a normal shock. The time rMw is the vibrational relaxation time computed by the Millikan and White formula [40]. Eq. (26) modifies rMw for high temperatures by the addition of the time r e where o-= 10 -16 cm 2, ( u ) is the average relative speed of the colliding partners, and n~ is the number density of the colliding species. Since our simulation is limited to N 2 - A r collisions only, we set n~ equal to the Ar concentration p. Thus % is a constant in the simulation. For convenience we normalize Eqs. (16), (18) and (25) by the factor V~/rMW and instead show d u / d r = - d ( V / V Q / d(t/rMw) in Fig. 3. At T = 15000 K, we compute
D.A. Gonzales, P.L. Varghese / Chemical Physics 195 (1995) 83-91 ' '""T'I
' ' '"'"l
' '"'1"1
' ''''"'
' I'"'"1
10°
0.14
10-3 0.12 <
10-6
0.10 o~ 10-9 ,,'
10-12
--6---- Park [1]
,
• - presem ,~ ....... present fief.) ~x
10q5 10-18 !0-2
10-3
0.08
~
10-1
100
101
0.06 0.04 102
t/'t'ss
10 10"3 ~
..................... ~
'
0.14 0.12
~
0.10
10-
t
10-18 . . . .
5
~ - . . . . . present (ref.)
I .... 5
I,,,,I 10
.... 15
0.06 I, ,,llll
20
25
J ,
0.04
30
t/'~ss
Fig. 3. Normalized time rate of change of average vibrational energy as predicted by several models. The average vibrational energy computed by the present method is included for reference. (a) L o g - l o g scale. (b) Semi-log scale.
V ~ = 11384.0 cm " i, ZMw = 0.457 ItS and "rc = 0.178 ItS.
For reference, Fig. 3 also shows the average vibrational energy V(t) as computed by the present work. Here one can see that only the present results (Eq. (18)) capture the changes in dV/dt due to the two-step behavior of V(t). Eq. (18) predicts a continuous decrease of dV/dt at early times followed by a temporary minimum plateau region through the steady state period. This plateau corresponds to a decreasing growth rate of V(t) through the steady state regime. This behavior is followed by an increase in dV/dt as V(t) begins its ascent to the equilibrium value. Finally, near equilibrium, there is
89
little change in V(t) and dV/dt approaches zero. It is clear that the Park model predicts slower changes in V(t) than the LRR model for the entire simulation. This is consistent with the intended diffusion description of dV/dt. However, the diffusive effect predicted by the Park model is too small in the steady state region and too large near equilibrium. The time interval shown in Fig. 3a is 61.4 ns ~< t ~< teq. The lower limit t = 61.4 ns (2.80 × 10 3"r,s) was chosen because at times less than this value, the present calculations predicted an average recombination flux coefficient k r = k r / p that deviated from its constant value of 4.71 × 10 34 c m 6 molecules 2 s J. This departure is most likely due to the failure of the linear approximation of Eq. (10) at times very far from final equilibrium. This is to be expected. Eq. (18) shows that dV/dt depends on kr(T) in general, but parametrically in this study since T is a constant in the simulation. Thus the present results for dV/dt are unreliable for times when k r departs from a constant. We again find it remarkable that Eq. (10) provides such reasonable results for early times on the order of 10 3~.. It is instructive to examine the results of Fig. 3a on a linear scale, and this is done in Fig. 3b. Here one can observe that the non-equilibrium condition of the gas is primarily described by the steady state condition, as suggested several years ago by Widom [44]. That is, the gas reaches steady state conditions almost instantaneously and then slowly relaxes to equilibrium. If one assumes the velocity downstream of the shock, U, is constant in the relaxation region, then Fig. 3b also corresponds to a spatial distribution through the relation x = Ut, where x is the spatial variable. This would support the notion that one is more likely to observe steady state phenomena in an actual shock tube experiment involving dilute mixtures. Fig. 3b shows that in the steady state region the Park model is more accurate than the LRR model, but is still several orders of magnitude from the value predicted by Eq. (18). As mentioned above, the LRR model provides a better prediction of dV/dt near equilibrium than the Park model. It would be interesting to investigate whether modifications to the Park model could be introduced that would result in a more faithful reproduction of Eq. (18) during the relaxation from steady state to equilibrium conditions.
D.A. Gonzales, P.L. Varghese / Chemical Physics 195 (1995) 83-91
90
4. Conclusions
Then
A numerical study of coupled vibration-dissociation-recombination processes has been performed for a shock heated mixture of N 2 dilute in At. The master equation for this process was linearized about equilibrium and solved using standard eigen-matrix techniques. It was found that the linear approximation was accurate to times very near the initial conditions of the simulation. The objective of this study was to examine the validity of two vibrational relaxation models commonly used in state-of-the-art flow codes by solving the master equation. Our calculations revealed a two-step relaxation behavior for the average vibrational energy, a feature not reproduced by the linear relaxation rate relation. An exact expression for the time rate of change of vibrational energy, dV/dt, was presented for dilute gas mixtures undergoing coupled vibration-dissociation-recombination processes. We observed that the linear relaxation rate relation and the diffusion model of Park do not accurately predict dV/dt during most of the relaxation simulation.
dV(t)
aV
dt
~ On, dt
Acknowledgements This work was supported in part by NASA Minority Graduate Engineering Grant NGT-70032, the NASA Center for Hypersonic Research and Training Grant NAGW-964, and the Texas Advanced Research Program. Computing resources were provided by the University of Texas System Center for High Performance Computing.
Appendix The derivation of dV/dt (Eq. (18)) is presented below. We begin with Eqs. (14) and (15) of the text, which we rewrite here as
V[n,.(t)] = Y'~P,.(t)E~.,
(A.1)
/'
(A.3)
Upon differentiation of Eq. (A.1) and substitution of Eq. (la) for dn,./dt, Eq. (A.3) becomes
dV(t) dt =(riM)
'[~
L,'
Y'. E,~-, n,- ~ ~_. E,~:,,,n,
2]
- EEt. krcflr~-
EEt, icrfl a
z'1~z
l'
l'
'
(A.4) Eq. (A.4) may be simplified in several ways. First, one may switch the indices in the first term inside the first set of square brackets. This term and the remaining double-summed term inside the same set of square brackets may then be combined, i.e.,
E E E,L,,,,,,.,- Z E E,L,,,,,,. t:
t ' ~- l'
l:
l't ~
~'
(A.5) = Eg,,,,,
(A.6)
where we identify the term inside the square brackets of Eq. (A.5) as S,. Second, one may observe that the double-summed terms inside the second set of square brackets in Eq. (A.4) sum identically to zero. Third, one may make use of Eq. (A.2) and use the following notation,
~0_,) =- ~P,(t)O,,.,
t'
,,,(t) P , ( t ) - 2",.('--------5- nM(t "
dn,
(A.7)
t,
(A.2)
where Q,. is a e-dependent quantity. Finally, using Eqs. (15), (20) and (21), we arrive at Eq. (18).
D.A. Gonzales, P.L. Varghese / Chemical Physics 195 (1995) 83-91
References [1] C. Park, J. Thermophys. Heat Transfer 3 (1989) 233. [2] J. Keck and G. Carrier, J. Chem. Phys. 43 (1965) 2284. [3] D.G. Rush and H.O. Pritchard, in: Eleventh Symposium (International) on Combustion (Combustion Institute, Pittsburgh, 1967)p. 13. [4] V.A. LoDato, D.L.S. McElwain and H.O. Pritchard, J. Am. Chem. Soc. 91 (1969) 7688. [5] D.L.S. McElwain and H.O. Pritchard, J. Am. Chcm. Soc. 91 (1969) 7693. [6] D.L.S. McElwain and H.O. Pritchard, in: Thirteenth Symposium (International) on Combustion (Combustion Institute, Pittsburgh, 1971) p. 37. [7] J.E. Dove and D.G. Jones, J. Chem. Phys. 55 (1971) 1531. [8] tt. Johnston and J. Birks, Accounts Chem. Res. 5 (1972) 327. [9] T. Ashton, D.L.S. McElwain and H.O. Pritchard, Can. J. Chem. 51 (1973) 237. [10] H.O. Pritchard, React. Kin. 1 (1975) 243. [11] M. Ramakrishna and S.V. Babu, J. Chem. Phys. 68 (1978) 163. [12] J.E. Dove and S. Raynor, J. Phys. Chem. 83 (1979) 127. [13] A.W. Yau and H.O. Pritchard, J. Phys. Chem. 83 (1979) 134. [14] J.H. Kiefer and J.C. Hajduk, Chem. Phys. 38 (1979) 329. [15] D.G. Truhlar, N.C. Blais, J.C. Hajduk and J.H. Kiefer, Chem. Phys. Letters 63 (1979) 337. [16] C. Lira and D.G. Truhlar, J. Phys. Chem. 87 (1983) 2683. [17] H.O. Pritchard, The quantum theory of unimolecular reactions (Cambridge Univ. Press, Cambridgc, 1984). [18] C. Lim and D.G. Truhlar, Chem. Phys. Letters 114 (1985) 253. [19] K. Haug and D.G. Truhlar, J. Phys. Chem. 89 (19851 3198. [20] H. ltoh, M. Koshi, T. Asaba and H. Matsui, J. Chem. Phys. 82 (1985) 4911. [21] K. Haug, D.G. Truhlar and N.C. Blais, J. Chem. Phys. 86 ([987) 2697; 96 (1992) 5556. [22] D.W. Schwenkc, J. Chem. Phys. 92 (1990) 7267.
91
[23] D.A. Gonzales and P.L. Varghese, AIAA Paper No. 93-0481 (1993). [24] D.A. Gonzales, Ph.D. Dissertation, The University of Texas at Austin, USA (1993). [25] D.A. Gonzales and P.L. Varghese, J. Phys. Chem. 97 (1993) 7612. [26] D.A. Gonzales and P.L. Varghese, J. Thermophys. Heat Transfer 8 (1994) 236. [27] B. Widom, J. Chem. Phys. 34 (196l) 20511. [28] C.A. Brau, J.C. Kcck and G.F. Carrier, Phys. Fluids 9 (19661 1885. [29] N.S. Snider, J. Chem. Phys. 45 (1966) 3299. [30] C.A. Brau, J. Chem. Phys. 47 (1967) 1153. [31] 1. Procaccia and R.D. Levine, J. Chem. Phys. 63 (1975) 4261. [32] A. Kafri and R.D. Levine, J. Chem. Phys. 64 (1976) 5320. [33] J.H. Kiefer, H.P.G. Joosten and W.D. Breshears, Chem. Phys. Letters 30 (1975) 424. [34] R.J. Rubin and K.E. Shuler, J. Chem. Phys. 25 (1956) 59. [35] E.W. Montroll and K.E. Shuler, J. Chem. Phys. 26 (19571 454. [36] 1. Oppenheim, K.E. Shuler and G.H. Weiss, Advan. Mol. Relaxation Processes 1 (1967) 13. [37] J.P. Appleton, M. Steinberg and D. Liquornik, J. Chem. Phys. 48 (1968) 599. [38] P. Huxley and J.N. Murrell, J. Chem. Soc. Faraday Trans. II 79 (1983) 323. [39] H.O. Pritchard, J. Phys. Chem. 65 (1961) 5t14. [40] R.C. Millikan and D.R. White, J. Chem. Phys. 39 (1963) 3209. [41] K.E. Shuler, G.H. Weiss and K. Anderson, J. Math. Phys. 3 ( 19621 550. [42] 1. Procaccia, Y. Shimoni and R.D. gevine, J. Chem. Phys. 65 (1976) 3284. [43] A.W. Yau and H.O. Pritchard, Can. J. Chem. 55 (1977) 737. [44] B. Widom, Science 148 (19651 1555.