Transcript Document

State to state models for
dissociation behind shock waves
Minelli Pierpaolo
IMIP-CNR
Bari, Italy
Boundary conditions for the kinetic
simulation of a shock wave
Rankine-Hugoniot equations
r0
r1
p0
p1
u0
u1
h0
h1
upstream
downstream
r1u1  r0 u0
Perfect gas equation
p

rkT
p1  r1u12  p0  r0 u02
m


u12
u02
h1   h0 
2
2
We use the Rankine-Hugoniot conditions, obtained numerically,
to implement the

Euler code (which use the full state-to-state kinetic scheme). These results allow
us to solve the shock wave problem neglecting transport phenomena (viscosity,
heat conduction, etc.). These last outcomes, obtained in a region where transport
coefficients are negligible and atomic recombination is irrelevant, are assumed
as boundary conditions in the data file of our DSMC code.
Collision Models
•
•
•
Elastic collisions: VSS model
RT collisions: Larsen-Borgnakke method
Vibrational and dissociative collision

A-M collisions


M-M collisions

A2  v  A  A2 v  A


A2  v  A  2 A  A
mono- quantum transitions

A2  v  A2  v  A2  v 1 A2  v 

A2  v  A2  v  A2  v 1 A2  v 1


multi - quantum transitions (dissociation)


A2  v  A2  v  2 A  A2  v 
Detailed balance principle
When a couple of selected particles has been accepted for the collision, the probability
associated to the particular state, described by ’A and ’B, will be given by:
  A ,  B  A , B ;g
pA , B 
 T OT
For definition of probability must be valid the normalization condition:
 p ,   1
A

B
All states
Aand B
In the case of study of the relaxation phenomena, the principle of the detailed balance
leads to the following expression:

 i  f ;gg 2   f  i; gg2
where i and f represents initial and final states of the molecules.
Parker Model
Parker has used an empirical non-impulsive potential model that incorporates
a small degree of asymmetry to derive an expression for the rotational relaxation
time trot. The following approximate expression is obtained:
Z rot 
Z 

1 
 2

3
2
rot 
T    2
T 
      
T   4
 T

1
2
where T * is the characteristic temperature of the intermolecular potential and
(Zrot)∞ is the limiting value. The value employed in our simulation are:

Nitrogen : T   91.5K, Zrot   16
Oxygen :
T *  90K, Zrot   14.4

QCT model for A-M collisions (1)
In the DSMC simulation, rotation and vibration have been uncoupled. QCT
cross sections must be mediated on the rotational spectrum in such a way to
obtain, at every fixed reference average rotational temperature, cross sections
that depend only from the vibrational level of the molecule. If we consider a
generic atom-molecule transition:
A2  v , j  A  A2 v , j A
We can write:
jmax
 j j    v , j  v , j

and
j  0
jmax
  v  
v  
 w j j j
j 0
jmax
w j
j 0
where
 E v , j kTrot
wj e
2 j 1
QCT model for A-M collisions (2)
Let consider the following atom-molecule transition:
A2  v  A  A2 v  A
with ’v = 0,1,…, max, max+1.
Generally, it is valid
 the law of conservation of energy, which, in this case, can
be written as:
ET OT  Ekin  Ev  Ekin  Ev 
From this equation, it is possible to obtain the only unknown variable E’kin,
and so g’ modulus.

Pseudo-level (max+1) models the continuum: if the molecule reaches the
pseudo-level (max+1), then the two bounded atoms dissociate according to:
A2  v  A  A2  max  1 A  3A
QCT and Detailed balance principle
We define the following two types of transitions according to the
corresponding cross sections:
Direct transitions:  i  i  k  with k  1,...,i
Indirect transitions:  i  i  k  with k  1,...,v max  i
But:
 i  f QCT g 2   f  iQCT g2

For this reason, we take into account only QCT results regarding direct
collisions and dissociation processes. For cross sections regarding indirect
processes, weimpose the principle of detailed balance:
 i  i  k; gg2   i  k  i; gQCT g 2
 i  i  k; Ekin 
Ekin
 i  k  i; Ekin QCT
Ekin
Molecule-Molecule collision models (1)
Two difficulties:
•
There are not detailed database for these kind of processes!
•
It is preferred to publish rate coefficients or fit rather than cross section data!
It is well known that the rate constants are related to the reaction cross sections
through the Laplace transforms in situations where thermal equilibrium may be
assumed for the translational degrees of freedom of the reactant molecules.
k  mr  2 kT 
 12
3
2

  EEe
 E kT
dE
0
It is not simple to invert this Laplace transform, because this operation is illconditionedand, moreover, when a mathematical solution is obtained, it
presents, sometimes, strong oscillation or non-physical behaviour.
Molecule-Molecule collision models (2)
In order to overcome these problems, we have used two types of approach, one
for nitrogen case study, and another one for oxygen and hydrogen.
Nitrogen
Despite the potential capabilities of state-to-state direct description of the
chemical kinetics offered by DSMC, phenomenological model for internal state
kinetics have been applied to this method as a rule (e.g. Larsen-Borgnakke
model). We therefore have developed a simple model for collision between
vibrationally excited molecules. Following the approach of Anderson, we use a
simple flexible statistical model for translational and internal energy exchange
that strictly satisfies detailed balance and incorporates some features of the
real cross section.
Molecule-Molecule: Nitrogen (1)
In the simulation we take into account two kinds of mono-quantum energy
exchange collisions:
• Vibration-translation (VT) energy exchange
 v   v 1
 
 cVT  v Ekin
 v   v

 v   v 1
 
 cVT  v 1Ekin
 v   v

Detailed balance principle
A2  v  A2  v   A2  v  1 A2  v 


 22 
 2
1
1
2
2
  
, VTv 
cVT v1,
 v mr gv ,
gv; gg2
 v 1
gg  

v 1,mrvg;g
vc
v  1
2
2




Molecule-Molecule: Nitrogen (2)
•
Vibration-vibration (VV) energy exchange
 v   v 1

 
 cVV  v 1 v Ekin
 v   v 1
For example, in this process

A2  v  A2  v   A2  v  1 A2  v  1
cross sections are proportional, through an adjustable factor cVV, to the product
between
the two higher vibrational numbers in the considered transition and the
kinetic energy of the couple of colliding particles after hit.
We have adjusted the model parameters to reproduce the rate coefficients
calculated by Billing and Fisher in the temperature range of interest.
Molecule-Molecule: Oxygen and Hydrogen
For oxygen and hydrogen, we have used a different approach. Practically, the
diatom-diatom energy transfer processes are modelled by fitting directly the
rate coefficients, of Billing and Kolesnick for oxygen and of Matveyev and
Silakov for hydrogen, by downhill simplex method.
This method is due to Nelder and Mead and it requires only function
evaluations. A simplex is the geometrical figure consisting, in N dimensions, of
N+1 points (or vertices) and all their interconnecting line segments, polygonal
faces, etc.
For multidimensional minimization, the best we can do is give our algorithm a
starting guess, that is, an N-vector of independent variables as the first point to
try. The algorithm is then supposed to make its own way downhill through the
unimaginable complexity of an N-dimensional topography, until it encounters a
(local, at least) minimum. The downhill simplex method must be started not just
with a single point, but with N+1 points, defining an initial simplex.
Downhill simplex method (2)
Mono-quantum transitions (VT) are included along with quasi-resonant VV
energy transfer.
The cross sections for these processes are assumed to be of the form:
  EexpE
E being the collision energy and ,  and  the fitting parameter for each
transition. The crosssections for the backward transitions are determined by
application of the detailed balance principle.
We decide to cut the cross sections at Emax=3eV.
The best fit parameters are able to reproduce the calculated rate coefficients
very well in the temperature range (300-10000K).
Downhill simplex method results (1)
10
2
10
0
10
2
10
0
10
0
H ( v)+H ( w+1) ---> H ( v+1)+H ( w)
-2
10
-4
10
v=1
v=5
v=10
v=14
-6
10
-8
10
0
0,5
1
1,5
E (eV)
(a)
2
2,5
3
2
v=1 ; w=0
v=5 ; w=4
v=10 ; w=9
v=14 ; w=13
-2
2
2
H2( v)+H2( 0) ---> H2( v-1)+H2( 0)
-4
2
10
 ( )
2
 ( )
10
2
10
 ( )
-2
10
2
-1
H2( v)+H2( 5) ---> H2( v-1)+H2( 5)
10
-6
10
-8
-3
10
v=1
v=5
v=10
v=14
-4
10
-5
10
0
0,5
1
1,5
E (eV)
(b)
2
2,5
3
0
0,5
1
1,5
E (eV)
2
2,5
(c)
(a) VT cross sections obtained by Downhill Simplex Method for a process in which molecule, that do not
change vibrational level, is at the ground state.
(b) VT cross sections obtained by Downhill Simplex Method for a process in which molecule, that do not
change vibrational level, is not at the ground state.
(c) VV cross sections obtained by Downhill Simplex Method for various processes.
3
Downhill simplex method results (2)
Downhill Sim plex Method
Matvey ev and Silakov
-14
10
Downhill Sim plex Method
Matvey ev and Silakov
-13
6 10
10
Downhill Sim plex Method
Matvey ev and Silakov
-12
-15
2
2
-15
2
2
2
2
-15
7 10
H ( 14)+H ( 14) ---> 2H+H( 13)
-13
3
3
8 10
H ( 10)+H ( 10) ---> H ( 11)+H ( 9)
-13
5 10
9 10
2
2
2
-1
2
-1
2
K(cm *s )
H ( 1)+H ( 1) ---> H ( 2)+H ( 0)
K(cm *s )
3
-1
K(cm *s )
9 10
-13
8 10
-13
4 10
-13
7 10
-15
6 10
-13
6 10
-13
-15
3 10
5 10
0
2000
4000
6000
T(K)
(a)
8000
10000
0
2000
4000
6000
T(K)
(b)
8000
10000
0
2000
4000
6000
T(K)
8000
10000
(c)
(a), (b), (c) Comparison between rate coefficients obtained by Downhill Simplex Method and that
obtained by fit from Matveyev and Silakov.
2
10
1
10
-8
10
-9
3
2
 ( )
-1
10
K(cm *s )
Downhill simplex method results (3)
Downhill Simplex Method
QCT
10
Downhill Simplex Method
QCT
0
H (14)+H ---> 3H
H (14)+H ---> 3H
2
10
2
-1
10
0
0,5
1
1,5
2
2,5
3
-10
0
2000
4000
E (eV)
(a)
6000
8000
10000
T(K)
(b)
(a) Comparison between cross sections obtained by Downhill Simplex Method and by QCT
(b) Comparison between rate coefficients, relative to the same process of (a), obtained by Downhill
Simplex Method and by QCT.
Downhill simplex method results (4)
-3
10
10
-14
10
-15
-4
3
2
 ( )
-1
K(cm *s )
10
Downhill Simplex Met hod
QCT
-5
10
Downhill Simplex Method
QCT
10
H (0)+H ---> 3H
-16
H (0)+H ---> 3H
2
2
-6
10
10
0
0,5
1
1,5
2
2,5
3
-17
0
2000
4000
8000
10000
T(K)
E (eV)
(a)
6000
(b)
(a) Comparison between cross sections obtained by Downhill Simplex Method and by QCT
(b) Comparison between rate coefficients, relative to the same process of (a), obtained by Downhill
Simplex Method and by QCT.
Downhill simplex method results (5)
10
2
Downhill Simplex Method
QCT
Downhill Simplex Method
QCT
-10
5 10
H (14)+H ---> H (13)+H
2
H (14)+H ---> H (13)+H
2
2
2
10
10
1
3
2
 ( )
-1
K(cm *s )
-10
4 10
-10
3 10
-10
0
2 10
0
0,5
1
1,5
E (eV)
(a)
2
2,5
3
0
2000
4000
6000
8000
10000
T(K)
(b)
(a) Comparison between cross sections obtained by Downhill Simplex Method and by QCT
(b) Comparison between rate coefficients, relative to the same process of (a), obtained by Downhill
Simplex Method and by QCT.
Molecule-Molecule: dissociation
•
Multi-quantum VT transition
Dissociation cross sections from lower levels are equal to mono-quantum cross
section from the last bounded level vmax to the pseudo-level vdiss, calculated at
the same total energy, dampened from a negative exponential factor.
  v ,  v   diss ,  v ; Etot    max ,  v   diss ,  v ; Etot  e

E max E v
kU
Such damping factor is called vibrational favouring.
Finally,
A2  max 
A2  v  multi-quantum
2 A  A2  v 
we can
write forVTdissociative
 through
If we consider
the following
reaction: reactions,
transition, the following expression of cross section:
Cross section can be written as:


  max v




,  v   diss ,  v ; Etot  E  E  cVTE diss eEkin kU
 v ,max
v   diss ,  v ; Etot  Ekin  E v  cmax
VT diss
kin


kin

   
E
E
Shock wave produced in pure Nitrogen (I)
16
1 10
-1
15
8 10
15
6 10
1 10
-2
-2
at
6 10
(a)
15
4 10

density , cm
-3
8 10
-2
4 10
(b)
-2
2 10
15
2 10
0
0 10
0
0 10
-200
0
200
400 600
x/
800
1000 1200
-200
0
200
400 600
x/
800
1000 1200
In Fig. (a) and (b) we report, the total number density and the atomic molar
fraction. The spatial scale is always expressed using the upstream mean free
path  as a unit. We can observe the effect of the thermal relaxation on the gas
number density. The compression ratio behind the shock waves reaches rapidly
the value of 6 valid for strong shocks in a gas of rigid rotors. From here it
grows slowly to the value of 8, following the slow vibrational relaxation.
Further increase of the compression ratio is due to the increase of the atom
fraction.
Shock wave produced in pure Nitrogen (II)
4
4
1,2 10
1,2 10
T
T
3
8,0 10
3
8,0 10
T
T,K
01
(c)
3
4,0 10
T,K
rot
3
T
T
4,0 10
(d)
rot
T
01
0
0,0 10
-20
0
0
20
40
x/
60
80
0,0 10
-200
0
200
400 600
x/
800
1000 1200
In Fig. (c) and (d) we report the three ‘temperature’ which are relevant for this
study: i.e. the static translational temperature T, the rotational temperature Trot
and the T01 vibrational temperature
1
T01  01k ln n1 n0 
It can be seen that the sudden jump in T is rapidly followed by slower jump for
Trot. The vibrational relaxation follows with a longer relaxation time.
Shock wave produced in pure Nitrogen (III)
0
0
10
10
x/=-138
x/=17
x/=26
x/=47
x/=74
x/=145
x/=850
-1
10
-2
10
-3
10
10
-2
10
(a)
(a)
-3
10
-4
10
-5
10
(b)
-4
10
10
x/=-138
x/=5
x/=11
x/=17
x/=23
x/=32
-1
-5
0
2
4
Energy , eV
6
8
0
1
2
3
4
Energy , eV
5
6
7
This slower relaxation is also affected by its mesoscopic structure in terms of
the vibrational distribution functions (vdf). In Fig. (a) we report some of these
vdf’s at different positions along the flow. It can be noticed that the vdf is
sensibly non-Boltzmann due to the highly effective VV/VT processes.
In Fig. (b) is plotted the distribution of (classical) rotational energy. The
different curves refer to different positions along the shock front and during the
rotational relaxation. It can be seen the strong deviation from the equilibrium
distribution during relaxation.
Shock wave produced in pure Nitrogen (IV)
0
5
10
5 10
5
v , cm/s
4 10
x/=-138
x/=5
x/=11
x/=14
-1
10
x/=17
x/=20
x/=32
x/=1146
5
3 10
-2
(c)
(a)
5
2 10
10
(d)
-3
10
5
1 10
-4
0
0 10
-200
0
200
400 600
x/
800
1000 1200
10
0
0 10
5
2 10
5
4 10
v , cm/s
5
6 10
5
8 10
In Fig. (c) we report the velocity along the flow, while in Fig. (d) we report the
translational distribution function at different positions along the shock front.
At the two extreme points the distribution are equilibrium distributions at the
upstream and downstream temperatures, respectively. At the points in between
a sensible deviation from the Maxwell one is observed, in the form of a ‘tail’.
Even if the effect of such a deviation on macroscopic quantities is probably
small under typical conditions, these plots show the degree of detail which our
model can achieve in describing the shock thickness.
Conclusions
In our work, we have studied by direct simulation the vibrational and dissociation
kinetics in a strong normal shock wave in molecular gases. The main achievement
of our work in comparison with previous studies is the self-consistent modelling
of translational, rotational and vibrational kinetics of a diatomic gas which is
undergoing dissociation in the shock wave. Furthermore, we used only data in
form of cross sections from the literature, or calculated by molecular dynamics
methods. Results show that the vibrational, rotational and translational
distribution functions can significantly deviate from the equilibrium one in the
shock region. Such deviations affect in turn the rate of elementary processes.
Improvements
Critical points:
•
Computational cost
•
There are few data for molecule-molecule interaction
In order to overcome the first problem we are led to consider a parallelization
of the code. This will also permit to extend the calculations to 2D.
This is necessary in order to couple the DSMC with fluid-dynamics codes for
the simulation of near-continuum flows.
A hybrid particle-continuum computational method would be more accurate
in order to describe re-entry phenomena, paying attention to the information
exchange at the interface between the particle and the continuum domain.
Finally we are studying how to couple DSMC with a PIC method in order to
reproduce plasma phenomena or shock waves in ionized gases.