Non-equilibrium Effects in Viscous Reacting Gas Flows Elena Kustova - - PowerPoint PPT Presentation

Issues in Solving the Boltzmann Equation for Aerospace ICERM, Brown University, Providence June 3–7, 2013

Non-equilibrium Effects in Viscous Reacting Gas Flows

Elena Kustova Saint Petersburg State University

The Boltzmann equation (1872)

More than 140 years

• f studying the Boltzmann equation

∂f ∂t +

u · ∂f

∂ r +

F · ∂f

∂ u = △collf

Still unsolved Still plenty of surprises and contradictions Still inspires new studies

Outline

Introduction Reduced-order non-equilibrium fluid dynamic models derived from the Boltzmann equation General idea State-to-state model Multi-temperature models One-temperature models Limitations of models commonly used in CFD Reaction rates and normal mean stress in one-temperature viscous flows Vibrationally non-equilibrium flows. Rate of vibrational relaxation On different contributions to the heat transfer Some features of transport in gases with electronic excitation Conclusions

• Introduction. Methods for solving the Boltzmann equation
• D. Hilbert, S. Chapman, D. Enskog, L. Waldmann, H. Grad, G. Bird,
• M. Kogan, C. Cercignani, S. Vallander, R. Brun, V. Zhdanov,
• E. Nagnibeda, and many others

Linearized Boltzmann equation Model equations (BGK, ES and other modifications) Integral form of the Boltzmann equation Using of variational principle Moment methods (Grad’s method and its generalization) Discrete velocities method Asymptotic methods (Hilbert, Chapman–Enskog and its generalizations) Numerical solution of the Boltzmann equation Direct simulations Monte Carlo

• Introduction. Gases with internal degrees of freedom

Different ways of description: Classical: both translational and internal degrees of freedom are described classically (Taxman, Kagan) Quantum mechanical: both translational and internal degrees

• f freedom are quantized (Waldmann, Snider).

Quasi-classical: while translational degrees of freedom are treated classically, the internal modes are quantized (Wang Chang, Uhlenbeck)

• Introduction. Gases with internal degrees of freedom

Let fcij( r, u, t) be a distribution function of c particles over velocity

• uc, vibrational and rotational energies εc

i , εci j .

The Wang Chang–Uhlenbeck equation (1951): ∂fcij ∂t + uc · ∂fcij ∂ r + F · ∂fcij ∂ uc = Jcij =

• γ

J(γ)

cij ,

J(γ)

cij is specified by the cross section of a microscopic process γ

Dimensionless form (in the absence of mass forces): ∂fcij ∂t + uc · ∂fcij ∂ r =

• γ

1 εγ Jγ

cij,

εγ ∼ τγ θ ≪ 1.

Reduced-order models for fluid dynamics

Weak and strong non-equilibrium flows Weakly non-equilibrium flows τγ ≪ θ ∂fcij ∂t + uc · ∂fcij ∂ r = 1 εJtotal

cij

, ε ∼ τfp θ ≪ 1. Strongly non-equilibrium flows ∃γ : τγ ∼ θ ∂fcij ∂t + uc · ∂fcij ∂ r = 1 εJrap

cij + Jsl cij,

ε ∼ τrap τsl ≪ 1.

Modified Chapman–Enskog method

Distribution functions depend on r and t only through macroscopic parameters and their gradients: fcij ( r, u, t) = fcij ( uc, ρλ ( r, t) , ∇ρλ ( r, t) , ...) , Characteristic times of physical–chemical processes differ essentially, some of them proceed on the gas-dynamic time scale θ. The basis

• f the method is to establish the hierarchy of characteristic times

τrap ≪ τsl ∼ θ Collision operators are divided into two groups: operators of rapid and slow processes: Jrap

cij ,

Jsl

cij,

Modified Chapman–Enskog method. Collision invariants

Collision invariants ψcij + ψdkl = ψc′i′j′ + ψd′k′l′ Collision invariants for all processes ψ(λ)

cij ,

λ = 1, ..., 5 : mc, mc uc, mcu2

c

2 + εc

ij

Additional collision invariants for the most frequent collisions. ˜ ψ(µ)

cij ,

µ = 1, ..., M Number of additional invariants depend on the deviation from equilibrium.

Modified Chapman–Enskog method

Fluid dynamic variables corresponding to the collision invariants of all processes ρλ =

• cij
• ψ(λ)fcijd

uc, λ = 1, ..., 5 λ = 1: density ρ λ = 2, 3, 4: velocity v λ = 5: specific energy U Macroscopic variables corresponding to additional invariants of rapid processes ˜ ρµ =

• cij
• ˜

ψ(µ)

cij fcijd

uc, µ = 1, ..., M Can be different depending on flow conditions

Modified Chapman–Enskog method

Governing equations

Conservation equations correspond to the invariants of all processes ∂ρ(λ) ∂t +

• cij
• ψ(λ)

cij

uc · ∂fcij ∂ r d uc = 0, λ = 1, ..., 5 M relaxation equations ∂˜ ρ(µ) ∂t +

• cij
• ˜

ψ(µ)

cij

uc · ∂fcij ∂ r d uc =

• cij
• ˜

ψ(µ)

cij Jsl cijd

uc, µ = 1, ..., M Production term in the right-hand side is specified by slow processes.

Modified Chapman–Enskog method

Zero- and first-order solutions

The solution is sought in the form fcij =

∞

• n=0

ǫnf (n)

cij

Zero-order solution Jrap

cij

• f (0), f (0)

= 0 Zero-order distribution function is not local equilibrium First-order solution f (1)

cij

= f (0)

cij (1 + ϕcij)

The first order correction is found from the integral equations I rap

cij (ϕ) = Jsl(0) cij

− Df (0)

cij

I rap

cij

is the linearized operator of rapid processes, Df (0)

cij

is the streaming operator

State-to-state model

Time hierarchy τtr < τrot ≪ τvibr < τreact ∼ θ Macroscopic variables vibrational state populations velocity temperature Macroscopic equations ρdαci dt = −∇ · Jmci +

• r

˙ ξrνr,ciMc, c = 1, .., L, i = 0, ..., Lc ρdv dt = ∇· P ρdu dt = −∇· q + P : ∇v

Multi-temperature models

Time relation τtr < τrot < τVV ≪ τVT < τreact ∼ θ Macroscopic variables chemical species mass fractions velocity temperature vibrational temperatures Macroscopic equations ρdαc dt = −∇ · Jmc +

• r

˙ ξrνr,cMc, c = 1, .., L, ρdαcEv,c dt = −∇· qv + ˙ Ev,c ρdv dt = ∇· P ρdu dt = −∇· q + P : ∇v

One-temperature model

Time relation τtr < τint ≪ τreact ∼ θ Macroscopic variables chemical species mass fractions velocity temperature Macroscopic equations ρdαc dt = −∇ · Jmc +

• r

˙ ξrνr,cMc, c = 1, .., L, ρdv dt = ∇· P ρdu dt = −∇· q + P : ∇v

Comparison of models. Compressive flows

N2 flow behind a shock wave, M0 = 15, T0 = 293 K, p0 = 100 Pa

5 10 15 20 10-7 10-6 10-5 10-4 10-3 10-2 10-1 ni / n i

3' 2' 1' 3 2 1

0.5 1.0 1.5 2.0

• 500
• 400
• 300
• 200
• 100

q, kW/m

2

x, cm

3 2 1

Figure : Vibrational populations (a) and heat flux (b) behind the shock front.

Comparison of models. Expanding flows

N2 flow in a conic nozzle, T∗ = 7000 K, p∗ = 1 atm

10 20 30 40 10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 1 ni / n

2

i

4 3 1

10 20 30 40 50 103 104 105 106 q, W/m

2

x/R

4 3 2 1

Figure : Vibrational populations (a) and heat flux (b) along the nozzle axis.

Limitations of commonly used models

CFD, common practice: Using the Law of Mass Action (LMA) in viscous flow solvers Using the Landau-Teller expression for the rate of vibrational relaxation Neglecting the bulk viscosity and non-equilibrium reaction contributions in the normal mean stress Neglecting thermal diffusion in heat and mass transfer Neglecting electronic excitation

Normal stress and reaction rates in one-temperature viscous flows

Integral operators of rapid and slow processes Jrap

ci

= Jtr

ci + Jint ci ,

Jsl

ci = Jreact ci

, ǫ = τtr τreact Governing equations ρdαc dt + ∇ · J mc =

• r

˙ ξrνrcMc, c = 1, .., L, ρdv dt + ∇· P = 0, ρdu dt + ∇· q + P : ∇v = 0.

Transport and production terms

Mass diffusive flux Jmc(r, t) = mc

• i
• ccfci(r, u, t) duc

Pressure tensor P(r, t) =

• ci
• mcccccfci(r, u, t) duc

Energy flux q(r, t) =

• ci

mcc2

c

2 + εc

i

• ccfci(r, u, t) duc

Production term

• r

˙ ξrνrcMc = mc

• i
• Jsl

ci duc,

Jsl

ci = J2⇋2 ci

+ J2⇋3

ci

Integral operators for slow processes

Exchange reactions J2⇋2

ci

=

• dd′c′
• kk′i′

fc′i′ m3

c′sc′ i′

fd′k′ m3

d′sd′ k′

− fci m3

csc i

fdk m3

dsd k

• ×

m3

c′m3 d′sc′ i′ sd′ k′ gσc′i′d′k′ cidk

dΩd ud Dissociation reactions J2⇋3

ci

=

• d
• kk′

f ′

dk′fc′ff ′h3sc i

• mc

mc′mf ′ 3 − fcifdk

• ×

×gσdiss

ci, ddudduc′duf ′du′ d,

σc′i′d′k′

cidk

, σdiss

ci, d are the reaction cross sections

Production terms

Reaction rates ˙ ξr = kf , r

L

• c=1

ρc Mc ν(r)

rc

− kb, r

L

• c=1

ρc Mc ν(p)

rc

Reaction rate coefficients

kf , r = NA

• iki′k′

fcifdk ncnd g σf , r d2Ω dud duc, r = ex, di, kb, ex = NA

• iki′k′

fc′i′fd′k′ nc′nd′ g ′ σb, ex d2Ω dud′ duc′, kb, di = N 2

A

• ik

fc′ff ′f ′dk nc′nf ′nd σb, di duc dud duc′ duf ′ du′

d.

σf , r is the cross section of rth reaction

Zero-order approximation

Maxwell-Boltzmann distribution function f (0)

ci

= mc 2πkT 3/2 nc Z int

c (T)sc i exp

• −mcc2

c

2kT − εc

i

kT

• Transport terms:

J mc = 0, q = 0, P = pU Governing equations: ρdαc dt =

• r

˙ ξ(0)

r

νrcMc, c = 1, .., L, ρdv dt + ∇p = 0, ρdu dt + p ∇· v = 0.

Zero-order production terms

Zero-order reaction rates: ˙ ξ(0)

r

= k(0)

f , r L

• c=1

ρc Mc ν(r)

rc

− k(0)

b, r L

• c=1

ρc Mc ν(p)

rc

Zero-order reaction rate coefficient: k(0)

f , r = NA

• iki′k′

f (0)

ci f (0) dk

ncnd g σf , r d2Ω dud duc, If the reaction cross sections σf , r are known ⇒ the zero order rate coefficients can be easily calculated. Alternatively, the Arrhenius law can be applied for the k(0)

f , r

calculation, and the equilibrium constant for the backward reaction rate coefficient.

Using the detailed balance principle, we obtain: ˙ ξ(0)

r

= ωr k(0)

f ,r L

• c=1

nc NA ν(r)

rc

= ωr k(0)

f ,r L

• c=1

ρc Mc ν(r)

rc

ωr is the chemical reaction characteristics: ωr = 1 − exp Ar RGT

• Ar is the affinity of a chemical reaction r.

Summary on the zero-order (Euler) approximation

In the zero-order approximation, the stress tensor takes the diagonal form pU and does not depend on chemical reactions The rate of chemical reaction r depends only on the affinity of the appropriate reaction Ar and does not depend on the affinities of other reactions Therefore, no cross effects between different chemical reactions and between normal mean stress and chemical reactions appear in inviscid gas flows

Summary on the zero-order (Euler) approximation

The law of mass action is valid in the zero-order (inviscid) flow approximation of the Chapman-Enskog method The zero-order chemical-reaction rate coefficient can be calculated by averaging the corresponding cross section over the Maxwell-Boltzmann distribution The main problem in the modeling of inviscid reacting flows is the correct determination of the zero-order rate coefficients of chemical reactions k(0)

f ,r , k(0) b,r . No other uncertainties occur in

the Euler equations Existing CFD models used for one-temperature inviscid flows are more or less self-consistent

First-order solution

In the first-order approximation, the distribution function takes the form fci = f (0)

ci

• 1 + φ(1)

ci

• The first-order correction to the distribution function:

φ(1)

ci =−1

nAci·∇ ln T−1 n

• d

Dd

ci·dd−1

nBci :∇v−1 nFci∇ · v − 1 n

• r

G r

ciωr

functions Aci, Dd

ci, Bci, Fci and G r ci are found from the linear

integral equations. The first-order normal mean stress and chemical reaction rates are determined by the scalar functions Fci, G r

ci, velocity

divergence ∇ · v and chemical reaction characteristics ωr

First-order stress tensor

P = πU + 2µ (∇v)0

S

The normal mean stress

−(π + p) = RGT

• r

lvr ωr − lvv ∇ · v,

coefficients lvr, lvv are determined by the bracket integrals

lvr = −[F, G r] NA , lvv = kT [F, F] .

There is a connection between coefficients lvr, lvv, relaxation pressure prel and bulk viscosity coefficient ζ:

−RGT

• r

lvrωr = prel, lvv = ζ

First-order reaction rates

The chemical-reaction rate in the first-order approximation takes the form ˙ ξr = ˙ ξ(0)

r

+ ˙ ξ(1)

r

The first-order correction ˙ ξ(1)

r

= −lrv ∇ · v + RGT

• s

lrs ωs The kinetic coefficients lrv = −[G r, F] NA , lrs = 1 RGT [G r, G s] NA

First order reaction rate coefficients

First order reaction rate coefficients: kf , r = k(0)

f , r (T) − ¯

k(1)

f , r (α1, ..., αL, ρ, T) − ˜

k(1)

f , r (α1, ..., αL, ρ, T, ∇ · v) .

First order corrections: ¯ k(1)

f , r = NA

n

• iki′k′

f (0)

ci f (0) dk

ncnd (Gci + Gdk) g σf , r d2Ω duc dud ˜ k(1)

f , r = ∇ · v NA

n

• iki′k′

f (0)

ci f (0) dk

ncnd (Fci + Fdk) g σf , rd2Ω ducdud

¯ k(1)

f , r are due to deviations from Maxwell-Boltzmann

distributions ˜ k(1)

f , r are due to spatial non-homogeneity

If there is no internal degrees of freedom, coefficients ˜ k(1)

f , r = 0

Summary on the first-order (Navier–Stockes) approximation

In the first-order approximation, we found the cross-coupling effects between chemical reactions and normal mean stress The reaction rates depend on the velocity divergence The rate of the rth reaction is affected by other reactions, therefore the existence of the cross-coupling effects between various chemical reactions in viscous flows becomes evident The law of mass action does not hold in the first-order approximation The ratio kf , r/kb, r is not equal to the equilibrium constant Keq Symmetry of coefficients lvr, lrv, and lrs can be proved based on the symmetry properties of bracket integrals. Therefore, the Onsager–Casimir reciprocity relations are valid

Summary on the first-order (Navier–Stockes) approximation

Existing CFD models for viscous flows do not account for these effects and are not completely self-consistence. This can lead to loss

• f accuracy in strongly non-equilibrium flow simulations

The effects can be small but apriori we do not know this. We encourage to check above contributions before neglecting them To implement the non-equilibrium effects we should modify the commonly used Navier–Stockes governing equations in order to take into account the bulk viscosity, the contribution of chemical reactions to the normal mean stress, and the first-order corrections to the reaction rates Developing the efficient numerical algorithms for the calculation of the first-order contributions: ζ, prel, ¯ k(1)

f , r, ˜

k(1)

f , r, ¯

k(1)

b, r, and ˜

k(1)

b, r is of

particular importance

Compressive N2/N flow behind a shock wave

The initial conditions are M0 = 15, T0 = 293 K, p0 = 100 Pa.

Approximate study: First, the flow parameters and their derivatives are found in the inviscid flow approximation; then these quantities are used as input parameters for the calculation of the first-order normal mean stress and reaction rate coefficients.

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 0.2 0.3 0.4 0.5 0.6 0.7

(a)

shear viscosity coefficient bulk viscosity coefficient

η, ζ · 10

3, Pa · s

x, cm

0.0 0.5 1.0 1.5 2.0

• 100
• 80
• 60
• 40
• 20

(b) prel

ζ div · v

prel , ζ div · v, Pa x, cm

Figure : Shear and bulk viscosity coefficients (a); first-order corrections to the

normal mean stress (b) behind the shock front as functions of x.

Compressive N2/N flow behind a shock wave

0.0 0.5 1.0 1.5 2.0 0.262 0.264 0.266 0.268 0.270 0.272

(a)

π

p p, π, atm x,cm

0.0 0.5 1.0 1.5 2.0 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45

(b)

δ, %

x, cm

Figure : Pressure p and normal mean stress π (a); per cent first-order

correction to the normal mean stress (b) behind the shock front. the bulk viscosity coefficient exceeds the shear viscosity coefficient prel and the term ζ∇ · v associated to bulk viscosity are of the same order; moreover, the absolute value of prel is higher in a compressive flow. In order to stay self-consistent, one should take into account simultaneously both effects. The total contribution of the correction terms to the normal mean stress is weak: π is slightly higher compared to p; the correction δ < 0.5% .

Compressive N2/N flow behind a shock wave

0.0 0.5 1.0 1.5 2.0 10000 20000 30000 40000 50000 60000 70000 80000

(a) kf , m

3/mole/s

x, cm kf

(1) (div · v)

kf

(1) (prel)

kf kf

(0)

0.0 0.5 1.0 1.5 2.0

• 3000
• 2000
• 1000

1000 2000 3000 4000 5000 6000

(b) kb

(1) (prel)

kb

(1) (div · v)

kb kb

(0)

x, cm kb , m

6/mole 2/s

Figure : Dissociation (a) and recombination (b) rate coefficients behind the

shock front as functions of x. the contribution of the terms ˜ k(1)

f

, ˜ k(1)

b

associated to the bulk viscosity effect, is weak the correction terms ¯ k(1)

f

, ¯ k(1)

b

connected to the relaxation pressure play an important role: in the beginning of the relaxation zone, their values are of the same order (and even higher for ¯ k(1)

b ) as the zero-order rate

coefficients k(0)

f

, k(0)

b

Compressive N2/N flow behind a shock wave

0.0 0.5 1.0 1.5 2.0 20 40 60 80 100 120 140 160 180

(a)

δ kf=δ ξ· δ kb δ, %

x, cm

0.0 0.5 1.0 1.5 2.0 20 40 60 80 100 120 140 160 180

(b)

ξ·(0) ξ· ξ·, mole/m

3/s

x, cm

Figure : First-order corrections to the reaction rate coefficients (a) and

reaction rate (b) behind the shock front. Taking into account the first-order effects leads to a considerable decrease of the total reaction rate coefficient close to the shock front the per cent correction reaches 80% for the dissociation and 170% for the recombination rate coefficients The difference between ˙ ξ and ˙ ξ(0) is significant; the value of the first-order correction to ˙ ξ is approximately equal to that for the dissociation rate coefficient (up to 80%).

Expanding N2/N flow in a nozzle

A conic nozzle with an angle 21◦. The low pressure throat conditions are T∗ = 7000 K, p∗ = 1 atm.

10 20 30 40 50 0.00 0.05 0.10 0.15 0.20 0.25 0.30

(a)

shear viscosity coefficient bulk viscosity coefficient

η, ζ · 10

3, Pa · s

x/R 10 20 30 40 50 1E-3 0.01 0.1 1 10 100

(b)

prel

ζ div · v

prel , ζ div · v, Pa x/R

Figure : Shear and bulk viscosity coefficients (a); first-order corrections to the

normal mean stress (b) along the nozzle axis as functions of x/R. T∗ = 7000 K, p∗ = 1 atm.

Expanding N2/N flow in a nozzle

1 2 3 0.05 0.10 0.15 0.20

(a)

π

p x/R p, π, atm

10 20 30 40 50 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

(b)

δ, %

x/R

Figure : Pressure p and normal mean stress π (a); per cent first-order

correction to the normal mean stress (b) along the nozzle axis. shear and bulk viscosity coefficients are of the same order of magnitude the relaxation pressure is one to two orders lower compared to the term ζ∇ · v The first-order effects lead to a slight decrease of the normal mean stress; the per cent correction is about 2% close to the throat and decreases rapidly with x/R.

Expanding N2/N flow in a nozzle

1.0 1.5 2.0 2.5 3.0 0.00 0.01 0.02 0.03 0.04 0.05 0.06

(a) x/R kf kf

(1) (div · v)

kf

(1) (prel)

kf

(0)

kf , m

3/mole/s

10 20 30 40 50

• 1x10

6

1x10

6

2x10

6

3x10

6

4x10

6

5x10

6

(b) kb kb

(1) (div · v)

kb

(1) (prel)

kb

(0)

kb , m

6/mole 2/s

x/R

Figure : Dissociation (a) and recombination (b) along the nozzle axis.

the contribution of the corrections ¯ k(1)

f

, ¯ k(1)

b

determined by the relaxation pressure is considerably lower compared to the role of the terms ˜ k(1)

f

, ˜ k(1)

b

governed by the velocity divergence: the term ¯ k(1)

f

is negligible, and ¯ k(1)

b

≪ ˜ k(1)

b

Expanding N2/N flow in a nozzle

10 20 30 40 50 20 40 60 80 100 120

(a)

δ kb=ξ· δ kf

x/R

δ, %

1 2 3 4 5 6 7 8 9 10

• 30
• 25
• 20
• 15
• 10
• 5

(b)

ξ· ξ·(0)

x/R

ξ·, mole/m

3/s

Figure : First-order corrections to the reaction rate coefficients (a) and

reaction rate (b) along the nozzle axis. For the dissociation rate coefficient, the first-order corrections are important close to the throat; their contribution is about 30% near the critical section and decreases to 0.3-0.5% in the expanding part For the recombination rate coefficient, the first-order effects are significant in the whole flow range For the reaction rate, the first-order effects are found to be important; the difference is determined mainly by recombination reaction; the first-order correction is approximately equal to δkb.

Expanding N2/N flow in a nozzle

The high pressure throat conditions: T∗ = 7000 K, p∗ = 100 atm. Qualitative results are similar:

bulk and viscosity coefficients take approximately the same values the contribution of relaxation pressure is much less compared to the term ζ∇ · v

Quantitatively, the first-order corrections are very low:

the correction to the normal mean stress does not exceed 0.05% the correction to the dissociation rate coefficient is less than 0.25% the correction to the recombination rate coefficient as well as to the reaction rate achieves the maximum value 2.5%.

Thus one can conclude, that in a high pressure expanding flow, the first-order effects appear to be negligible.

Summary on the first order corrections in

• ne-temperature flows

Although the bulk and shear viscosity coefficients take close values, the effect of bulk viscosity on the normal mean stress is found to be weak; the same holds for the contribution of prel to π. The maximum correction to the normal mean stress is about 2% in a low-pressure expanding flow and is much less for other cases. Alternatively, the contribution of the first-order corrections to the reaction rates and rate coefficients is much more important: up to 170% for the recombination rate coefficient and 120% for the reaction rates. Exception is for the high-pressure expanding flow, where all first-order corrections are found to be small.

Vibrationally non-equilibrium flows

State-to-state model

Relaxation equations (equations of detailed vibrational-chemical kinetics for the level populations) ρdαci dt = −∇ · Jmci +

• r

˙ ξrνr,ciMc, c = 1, .., L, i = 0, ..., Lc Rate of transitions/reactions ˙ ξr = 1 NA

• j
• Jr

cij duc,

Jr

cij = Jr, vibr cij

+ Jr, react

cij

Vibrationally non-equilibrium flows

Multi-temperature models

Relaxation equations ρdαc dt = −∇ · Jmc +

• r

˙ ξrνr,cMc, c = 1, .., L, ρdαcEv,c dt = −∇· qv + ˙ Ev,c, c = 1, ..., Lmol Rate of vibrational energy relaxation ˙ Ev,c = NA

• r
• i

εc

i ˙

ξrνr,ci For anharmonic oscillators, vibrational energy εc

i is not a

collisional invariant, and relaxation equation should be written for the specific numbers of vibrational quanta Wc, ρcWc =

i inci

State-to-state model. Zero-order approximation

The distribution function f (0)

cij

is given by Maxwell-Boltzmann distribution over velocity and rotational energy The rate of non-equilibrium process r NA ˙ ξ(0)

r

=

• 1 − exp

Ar kT

jj′ll′

• f (0)

cij f (0) dkl ˜

σf , r(g)ducdud r stands for any vibrational energy transition or state-specific chemical reaction

State-to-state model. Zero-order approximation

Ar are generalized affinities of state-specific chemical reactions and vibrational transitions: Ar, ex = 3 2kT ln mcmd mc′md′ + kT ln Z rot

ci Z rot dk

Z rot

c′i′Z rot d′k′

− kT ln ncindk nc′i′nd′k′ + + (εc′ + εd′ − εc − εd) +

• εc′

i′ + εd′ k′ − εc i − εd k

• ,

Ar, dr = 3 2kT ln mc mc′mf ′ − 3 2 ln(2πkT)+3kT ln h +kT ln Z rot

ci −

−kT ln nci nc′nf ′ + (εc′ + εf ′ − εc − εc

i ) .

State-to-state model. Zero-order approximation

For the particular case of VV transitions within the same species c, Ar takes the simplified form Ai + Ak = Ai′ + Ak′ Ar, VV = −kT ln nink ni′nk′ + (εi′ + εk′ − εi − εk) For VT transitions Ai + M = Ai′ + M Ar, VT = −kT ln ni ni′ + (εi′ − εi)

State-to-state model. Zero-order approximation

If we introduce ωr = 1 − exp Ar kT

• then the zero-order rate of vibrational transitions and

state-specific chemical reactions can be written as a linear function of ωr: ˙ ξ(0)

r

= ωrk(0)

f ,r (T)ΠL c=1ΠLc i=0

ρci Mc ν(r)

r,ci

k(0)

f ,r is the zero-order rate coefficient for the rth vibrational

transition or chemical reaction. The last expression is the generalized mass action law for coupled chemical reactions and vibrational transitions

Limit transitions for generalized affinities

Multi-temperature model, harmonic oscillators

Boltzmann vibrational distribution nci = nc Z vibr

c

(Tv,c) exp

• −

εc

i

kTv,c

• Generalized affinity

Ar, ex = 3 2kT ln mcmd mc′md′ + kT ln Z rot

c Z rot d

Z rot

c′ Z rot d′

Z vibr

c

(Tv,c)Z vibr

d

(Tv,d) Z vibr

c′ (Tv,c′)Z vibr d′ (Tv,d′)−

−kT ln ncnd nc′nd′ + (εc′ + εd′ − εc − εd) +

• εc′

i′ + εd′ k′ − εc i − εd k

• −

−kT

• εc′

i′

kTv,c′ + εd′

k′

kTv,d′ − εc

i

kTv,c − εd

k

kTv,d

Limit transitions for generalized affinities

Multi-temperature model, harmonic oscillators

VT relaxation in a single-component gas Ar, VT = (εi′ − εi)

• 1 − T

Tv

• = ∆εii′
• 1 − T

Tv

• For harmonic oscillators, only single-quantum jumps are allowed.

Therefore, for VT relaxation in a single-component system there are

• nly two types of reactions

ω1 = 1−exp hν kT

• 1 − T

Tv

• ;

ω2 = 1−exp

• − hν

kT

• 1 − T

Tv

• Zero-order rate of vibrational relaxation ˙

E (0)

v

˙ E (0)

v

= nhν 2NA

• i

ni

• r=1,2

ωrk(0)

f ,r

Limit transitions for generalized affinities

Multi-temperature model, anharmonic oscillators

Treanor vibrational distribution nci = nc Z vibr

c

(T, T1,c) exp

• −εc

i − iεc 1

kT − iεc

1

kT1,c

• Generalized affinity

Ar, ex = 3 2kT ln mcmd mc′md′ +kT ln Z rot

c Z rot d

Z rot

c′ Z rot d′

Z vibr

c

(T, T1,c)Z vibr

d

(T, T1,d) Z vibr

c′ (T, T1,c′)Z vibr d′ (T, T1,d′)−

−kT ln ncnd nc′nd′ + (εc′ + εd′ − εc − εd) + +iεc

1

T T1,c − 1

• +kεd

1

T T1,d − 1

• −i′εc′

1

T T1,c′ − 1

• −k′εd′

1

T T1,d′ − 1

Limit transitions for generalized affinities

Multi-temperature model, anharmonic oscillators

VT relaxation in a single-component gas Ar, VT =

• i′ − i
• ε1
• 1 − T

T1

• = ∆iε1
• 1 − T

T1

• Multi-quantum (M) transitions are allowed:

1-quantum : A11 = ε1

• 1 − T

T1

• ;

A12 = −ε1

• 1 − T

T1

• ,

2-quantum : A21 = 2ε1

• 1 − T

T1

• ;

A22 = −2ε1

• 1 − T

T1

• ,

... M-quantum : AM1 = Mε1

• 1 − T

T1

• ;

AM2 = −Mε1

• 1 − T

T1

• Then the number of possible VT reactions is 2M.

Zero-order rate of vibrational relaxation ˙ W (0) is obtained in the similar form as a linear function of ωr but it depends on ∆i rather than ∆εi and summation is taken over 2M reactions

Limit transitions for generalized affinities

One-temperature model

Vibrational distribution nci = nc Z vibr

c

(T) exp

• − εc

i

kT

• Generalized affinity

Ar, ex = 3 2kT ln mcmd mc′md′ +kT ln Z int

c Z int d

Z int

c′ Z int d′

−kT ln ncnd nc′nd′ +(εc′ + εd′ − εc − εd) Ar coincides with the classical definition of the affinity of exchange reaction If only vibrational transitions take place in the mixture (and no chemical reactions) then the Ar are identically zero, and we have the case of complete thermodynamic equilibrium

State-to-state model. First-order solution

First-order distribution function f (1)

cij = f (0) cij

• − 1

nAcij ·∇ ln T − 1 n

• dk

Ddk

cij ·ddk − 1

nBcij : ∇v − − 1 nFcij∇ · v − 1 n

• r

G r

cijωr

• Scalar fluxes are specified by the terms Fcij∇ · v and

r G r cijωr

This representation of the last term becomes possible because the zero order rates of transitions are expressed as linear functions of the scalar force ωr.

Cross-coupling effects for the state-to-state model

Stress tensor P = πU + 2µ(∇v) s

• µ is the shear viscosity, π is the normal mean stress

π = −p−kT

• r

lvr ωr +lvv ∇ · v, lvr = − [F, G r] , lvv = kT [F, F] First-order rate of non-equilibrium vibrational transitions and chemical reactions ˙ ξr = ˙ ξ(0)

r

+ ˙ ξ(1)

r

NA ˙ ξ(1)

r

= kT

• s

lrs ωs−lrv ∇ · v, lrv = − [G r, F] , lrs = [G r, G s] kT Normal mean stress and rates of transitions-reactions are the linear functions of the same scalar forces ∇ · v and ωr and are strongly coupled. The kinetic coefficients are symmetric lvr = lrv, lrs = lsr due to symmetry properties of bracket integrals. Therefore the Onsager-Casimir reciprocity relations are verified for the case of strong vibrational-chemical coupling.

Cross-coupling effects for the state-to-state model

In the linearized case ωr ≈ − Ar kT Scalar fluxes become linear functions of generalized affinities, which is consistent with the results of linear irreversible thermodynamics The generalized mass action law does not work in a viscous flow since ˙ ξr depends on the velocity divergence and affinities

• f all transitions/reactions

For self-consistent CFD simulations of viscous compressible flows, it is necessary to include all the first-order correction

• terms. Including only the bulk viscosity in the fluid-dynamics

equations is not self-consistent if the corresponding term lrv ∇ · v in the reaction-rate expressions is neglected.

Cross-coupling effects for multi-temperature models

Single-component gas

Distribution function f (1)

ij

=f (0)

ij

• − 1

n Aij·∇ ln T− 1 n A(1)

ij ·∇ ln Tv− 1

n Bij : ∇v− 1 n Fij∇·v− 1 n

• r=1,2

G r

ijωr

• For anharmonic oscillators, we have ∇T1 instead of ∇T, the summation

is taken over r = 1, 2, ..., 2M, and ωr are calculated differently Normal mean stress is obtained in a similar form The rate of vibrational energy relaxation for harmonic oscillators ˙ Ev = ˙ E (0)

v

+ ˙ E (1)

v ,

˙ E (1)

v

=

• r=1,2

∆εr kT ˙ E (1)

v,r ,

˙ E (1)

v,r = kT s=1,2

lrs ωs − lrv ∇ · v

• ,

lrv = − [G r, F] , lrs = [G r, G s] ∆ε1 = hν, ∆ε2 = −hν. For anharmonic oscillators, the summation is taken from 1 to 2M, and ∆εr should be replaced by mhν, where m is the number of quanta transferred in the collision.

Cross-coupling effects for multi-temperature models

Cross-coupling terms are written in the form similar to that in the state-to-state model The kinetic coefficients are symmetric lvr = lrv, lrs = lsr Therefore the reciprocity Onsager–Casimir relations are valid for the multi-temperature case The use of Landau-Teller expression for the rate of vibrational relaxation ˙ Ev = ρEv(T) − Ev(Tv) τvibr is not justified for viscous flows because of cross-coupling terms

Numerical example

4000 4500 5000 5500 6000

• 80
• 60
• 40
• 20

20 40 60 80 100

T=5000 K N2 (a)

δ, %

Tv, K p=100 Pa, div v=1000s

• 1

p=1000 Pa, div v=1000s

• 1

p=1000 Pa, div v=2000s

• 1

4000 4500 5000 5500 6000 20 40 60 80 100

T=5000 K O2 (b)

δ, %

Tv, K p=100 Pa, div v=1000s

• 1

p=100 Pa, div v=2000s

• 1

p=1000 Pa, div v=1000s

• 1

p=1000 Pa, div v=2000s

• 1

Contribution of the first-order correction to the total rate δ as a function of Tv for N2 (a) and O2 (b).

For low pressure and large velocity divergence ˙ E (1)

v

may be of the same

• rder as ˙

E (0)

v . For O2 the effect is weaker, the mean contribution of ˙

E (1)

v

is within 1-2%.

Numerical example

4000 4500 5000 5500 6000

• 600
• 400
• 200

200 400

(a) T=5000 K, p=100 Pa, div v=1000 s

• 1

uv

(1), J/m 3/c

Tv, K N2 O2

5 10 15 20

• 12
• 10
• 8
• 6
• 4
• 2

2

(b)

δ, %

x/R N2, p0=1atm, T0=7000 K O2, p0=1atm, T0=4000 K

First-order correction ˙ E (1)

v

as a function of Tv (a) and contribution of the first-order correction to the total rate δ as a function of x/R in a nozzle.

Close to the throat (particularly for nitrogen), the first-order effects can influence noticeably the rate of vibrational relaxation, whereas with rising x/R (R is the throat radius), the contribution of the first-order correction decreases.

Summary on vibrational non-equilibrium flows

A self-consistent kinetic model relating the rates of non-equilibrium processes and the normal mean stress to the velocity divergence and chemical reaction/transition affinities is proposed. In the inviscid approximation, cross effects between reaction rates and diagonal elements of the viscous stress tensor do not appear. Cross effects between reaction/transition rates and diagonal elements of the viscous stress tensor exist in viscous gas flows; the rate of each reaction is affected by other reactions and flow compressibility; the law of mass action is violated for viscous flows; the Landau-Teller expression for the rate of vibrational relaxation does not hold.

State-to-state model. Heat and mass transfer features

State-specific diffusion velocity Vci = −

• dk

Dcidkddk − DTc∇ ln T = VDVE

ci

+ VMD

c

+ VTD

c

VDVE

ci

is the contribution of vibrational energy diffusion, characteristic feature of the state-to-state approach VMD

c

is the mass diffusion VTD

c

is the thermal diffusion

In CFD, basically the Fick’s law is used Vc = −Dc∇xc either with the effective species diffusion coefficients or constant Schmidt number. Thermal diffusion is systematically neglected.

State-to-state model. Heat and mass transfer features

Heat flux in the state-to-state model q = −λ∇T−p

• ci

DTci+

• ci

nciVci 5 2kT+ <εci>rot +εc

v + εc

• Taking into account state-specific diffusion velocity:

q = qHC + qMD + qTD + qDVE

qHC is the contribution due to heat conduction (Fourier flux)

Thermal diffusion and diffusion of vibrational energy is usually neglected

Contributions to the heat flux

Different contributions to the heat flux in compressive and expanding flows

0.5 1.0 1.5 2.0

• 600
• 400
• 200

200 400 q, kW/m

2

x, cm

5 4 3 2 1

10 20 30 40 50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 qα /q x/R

4 3 2 1

Shock wave. 1: total heat flux, 2: contri- bution of heat conduction, 3: thermal dif- fusion, 4: mass diffusion, 5: diffusion of vibrational energy.

• Nozzle. 1: contribution of heat con-

duction, 2: mass diffusion, 3: diffu- sion of vibrational energy, 4: thermal diffusion

Contributions to the heat flux in a boundary layer

(a) Mixture 25%CO2, 25%CO, 25%(O2+O), 25%C, Tw = 500 K, Te = 3000 K, pe = 1000 Pa, β = 500 s−1. Non-catalytic surface (b) Mixture 10.4% CO2, 57.6% CO, 31.7% (O2+O), 0.3% C, Tw = 1500 K, Te = 6000 K, pe = 888 Pa, β = 2708 s−1. Non-catalytic surface

2 104 4 104 6 104 8 104 1 105 1 2 3 4 5 6 7 8 1 - q

F

2 - q

TD

3 - q

MD

4 - q

DVE

5 - q q (W m-2) 1 2 3 4 5 2

e

2 104 4 104 6 104 8 104 1 105 1 2 3 4 5 6 7 8 1 - q

F

2 - q

TD

3 - q

MD

4 - q

DVE

5 - q q (W m-2) 1 2 3 4 5

f

Contributions to the diffusion velocity in a boundary layer

Mixture 10.4% CO2, 57.6% CO, 31.7% (O2+O), 0.3% C, Tw = 1500 K, Te = 6000 K, pe = 888 Pa, β = 2708 s−1. Non-catalytic surface

1 101 2 101 3 101 4 101 5 101 1 2 3 4 5 6 7 8 1- - V

TD-CO2(0,0,0)

2- - V

MD,DVE-CO2(0,0,0)

• VTD-CO2(0,0,0) and - VMD,DVE-CO2(0,0,0) (m s-1)

 1 2

• 5 10-4

5 10-4 1 10-3 1.5 10-3 2 10-3 2.5 10-3 3 10-3 1 2 3 4 5 6 7 2

a

• 1.2 100
• 1 100
• 8 10-1
• 6 10-1
• 4 10-1
• 2 10-1

2 10-1 4 10-1 1 2 3 4 5 6 7 8 1- - V

TD-O

2- - V

MD,DVE-O

• VTD-O and - VMD,DVE-O (m s-1)

 1 2

b

Contributions to the heat flux in a boundary layer

78.58% (N2+N), 21.38% (O2+O), 0.04% NO, Tw = 1000 K, Te = 7000 K, pe = 1000 Pa, β = 5000 s−1. (a) Non-catalytic surface (b) Partially catalytic surface

• 5 104

5 104 1 105 1.5 105 2 105 2.5 105 3 105 1 2 3 4 5 6 7 8 1 - q

F

2 - q

TD

3 - q

MD

4 - q

DVE

5 - q q (W m-2) 1 2 3 4 5

a

• 2 105
• 1 105

1 105 2 105 3 105 4 105 5 105 6 105 1 2 3 4 5 6 7 8 1 - q

F

2 - q

TD

3 - q

MD

4 - q

DVE

5 - q q (W m-2) 1 2 3 4 5

Summary on heat flux contributions

For a non-catalytic surface, the main contribution to the heat transfer is given by thermal conductivity and thermal diffusion; the contribution of the vibrational energy diffusion varies depending on the deviation of the flow from thermal equilibrium; the mass diffusion process is negligible. The mass flux near the surface is specified mainly by thermal diffusion. For a catalytic surface, mass diffusion is the main process responsible for the heat transfer; the contribution of thermal diffusion is found to be small; diffusion of vibrational energy can be important close to the wall. For shock waves and nozzle flows, the contribution of thermal diffusion is small. Diffusion of vibrational energy is of importance close to the shock front

Some features of transport in gases with electronic excitation

For temperatures lower than 8000-10000 K electronic excitation practically does not contribute to the heat capacities and transport coefficients. Therefore it is usually neglected At higher temperatures the role of electronic excitation occurs important for both molecules and atoms

Capitelli, first works in 1970-s, then interrupted Galkin, Zhdanov, some private discussions Capitelli (Bari), 2000-s, equilibrium plasma of argon and hydrogen atoms with electronic excitation SPbSU, 2000-s, non-equilibrium mixtures of atoms and molecules with electronic excitation

Contribution of electronic excitation

Contribution of translational and internal degrees of freedom to the heat conductivity coefficient % T, K N N2 λtr,N λint,N λtr,N2 λrv,N2 λint,N2 500 100 72.8 27.2 27.2 1000 100 66.1 33.9 33.9 5000 89.8 10.2 56.3 43.6 43.7 10000 70.1 29.9 47.9 40.5 52.1 15000 46.7 53.3 32.4 28.6 67.6 20000 21.4 78.6 33.0 27.8 67.0 25000 23.5 76.5 42.4 31.7 57.6 30000 38.3 61.7 53.9 34.2 46.1

Contribution of electronic excitation

New effects: Bulk viscosity in mixtures of atoms/ions/electrons Internal heat conductivity coefficients in atomic plasmas Prandtl number

500 10000 20000 30000 40000 50000 0,66 0,67 0,68 0,69 0,70 0,71 0,72

• Prandtl number

Temperature [K] N N

2

O O

2

Prandtl number for atomic and molecular species with electronic excitation.

Conclusions

Many non-equilibrium effects are systematically neglected in the CFD, among them

bulk viscosity and relaxation pressure viscous contributions to the reaction/transition rates electronic excitation thermal diffusion

We encourage to estimate these non-equilibrium effects before taking decision that they are negligible