April 2014 

NASA/TP–2014-218254 

Extension of a Kinetic Approach to Chemical 

Reactions to Electronic Energy Levels and 

Reactions Involving Charged Species with 

Application to DSMC Simulations 

Derek Shane Liechty 

Langley Research Center, Hampton, Virginia 
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page 

at http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443-757-5803 

• Phone the NASA STI Information Desk at 

443-757-5802 

• Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

April 2014 

NASA/TP–2014-218254 

Extension of a Kinetic Approach to Chemical 

Reactions to Electronic Energy Levels and 

Reactions Involving Charged Species with 

Application to DSMC Simulations 

Derek Shane Liechty 

Langley Research Center, Hampton, Virginia 

.
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

Acknowledgments 

This document was submitted as a dissertation to the Faculty of the University of 

Maryland in partial fulfillment of the requirements for the degree of Doctor of 

Philosophy in Aerospace Engineering. 

The author would like to thank Dr. Christopher Johnston for his assistance in postprocessing the DSMC simulations presented herein to obtain the radiative heating to the 

surface of the Fire II vehicle. 

The thesis advisor was Professor Mark J. Lewis, Willis Young Professor and Department 

Chair, Department of Aerospace Engineering, University of Maryland. Advisory 

Committee members were Professor Mark J. Lewis, Chair, Associate Professor James 

Baeder, Assistant Professor Ray Sedwick, Associate Professor Kenneth Yu, and 

Professor Gregory Jackson. 
iii 

 

Abstract 

The ability to compute rarefied, ionized hypersonic flows is becoming more important as 

missions such as Earth reentry, landing high mass payloads on Mars, and the exploration 

of the outer planets and their satellites are being considered. Recently introduced 

molecular-level chemistry models that predict equilibrium and nonequilibrium reaction 

rates using only kinetic theory and fundamental molecular properties are extended in the 

current work to include electronic energy level transitions and reactions involving 

charged particles. These extensions are shown to agree favorably with reported transition 

and reaction rates from the literature for near-equilibrium conditions. Also, the 

extensions are applied to the second flight of the Project FIRE flight experiment at 1634 

seconds with a Knudsen number of 0.001 at an altitude of 76.4 km. In order to 

accomplish this, NASA’s direct simulation Monte Carlo code DAC was rewritten to 

include the ability to simulate charge-neutral ionized flows, take advantage of the 

recently introduced chemistry model, and to include the extensions presented in this 

work. The 1634 second data point was chosen for comparisons to be made in order to 

include a CFD solution. The Knudsen number at this point in time is such that the 

DSMC simulations are still tractable and the CFD computations are at the edge of what is 

considered valid because, although near-transitional, the flow is still considered to be 

continuum. It is shown that the inclusion of electronic energy levels in the DSMC 

simulation is necessary for flows of this nature and is required for comparison to the CFD 

solution. The flow field solutions are also post-processed by the nonequilibrium 

radiation code HARA to compute the radiative portion of the heating and is then 

compared to the total heating measured in flight. 
iv 

Table of Contents 

Section 1 Introduction ............................................................................................ 2

1.1 Objectives of this Work ................................................................................ 3

1.2 Outline of Thesis .......................................................................................... 3

Section 2 Simulation of Rarefied and Weakly Ionized Reentry Flows ................ 5

2.1 Theoretical Basis of Particle Simulation ....................................................... 5

2.1.1 The Boltzmann Equation ....................................................................... 5

2.2 Traditional Direct Simulation Monte Carlo Methods ................................... 11

2.2.1 Collision Algorithm .............................................................................. 11

2.2.2 Internal Degrees-of-Freedom .............................................................. 16

2.2.3 Chemical Reactions ............................................................................ 23

2.2.4 Electric Field ....................................................................................... 27

2.3 Recent Direct Simulation Monte Carlo Advances ...................................... 29

2.3.1 Vibrational Relaxation ......................................................................... 29

2.3.2 Chemical Reactions – The Quantum-Kinetic Chemistry Model .......... 30

Section 3 Extending the Quantum Kinetic Chemistry Model ............................ 38

3.1 Detailed Electronic Energy Level Model .................................................... 38

3.1.1 Equilibrium Sampling .......................................................................... 38

3.1.2 Post-Collision Sampling ...................................................................... 39

3.1.3 Electronic Energy Level Transitions .................................................... 40

3.2 Chemical Reactions Involving Charged Species ........................................ 41

3.2.1 Ionization Reactions ........................................................................... 41

3.2.2 Associative Ionization Reactions ........................................................ 42

3.2.3 Charge Exchange Reactions .............................................................. 44

3.2.4 Other Exchange Reactions ................................................................. 45

3.3 Selecting a Reaction .................................................................................. 45

Section 4 Application of New Models to Upper-Atmosphere Hypersonic 

 Flows ..................................................................................................... 47

4.1 Measurement of Reaction/Transition Rates ............................................... 47

4.2 Definition of Translational, Rotational, Vibrational, and Overall 

Temperatures ....................................................................................................... 48

4.3 Electronic Energy Level Transitions and Comparison to Equilibrium ......... 50

4.3.1 Electron-Impact Excitation .................................................................. 50

4.3.2 Heavy Particle Impact Excitation ........................................................ 53

4.3.3 Equilibrium Conditions ........................................................................ 57

4.4 Associative Ionization Reactions ................................................................ 57

4.5 Ionization Reactions ................................................................................... 59

4.6 Charge Exchange Reactions ..................................................................... 59

4.7 Exchange Reactions Involving Charged Species ....................................... 60

Section 5 Modeling of the FIRE II Reentry Flow with Ionization and Electronic 

Energy Levels with DSMC ................................................................... 62
v 

5.1 The Project FIRE Flight Experiment ........................................................... 62

5.2 Details of the Continuum Flow Model and Non-Equilibrium Radiation Model

 63

5.3 DSMC Simulation Parameters ................................................................... 64

5.3.1 Collision Modeling ............................................................................... 64

5.3.2 Internal Energy Modeling .................................................................... 64

5.3.3 Chemistry Modeling ............................................................................ 65

5.3.4 Electric Field Modeling ........................................................................ 65

5.3.5 Gas-Surface Interactions .................................................................... 65

5.4 FIRE II Simulations .................................................................................... 66

5.4.1 Baseline Solution ................................................................................ 66

5.4.2 Addition of Electronic Energy Levels .................................................. 67

5.4.3 Quantum-Kinetic Solution ................................................................... 69

5.4.4 Comparison of Surface Heating .......................................................... 71

Section 6 Conclusions .......................................................................................... 73

6.1 Summary of Results ................................................................................... 73

6.2 Contributions .............................................................................................. 74

Section 7 References ............................................................................................ 75

Appendices ............................................................................................................. 83

Appendix A: Species and Chemistry Data ........................................................... 83

Appendix B: Quantum-Kinetic Electronic Energy Level Transition Rates. ....... 97

Appendix C: Quantum-Kinetic Reaction Rates.................................................. 113

 
vi 

List of Tables 

Table 5.1: List of DSMC FIRE II cases. ............................................................... 62 

Table 5.2: Free stream and surface parameters for FIRE II entry condition. ....... 62 

Table 5.3: FIRE II surface heating at r = 0.02 m. ................................................. 72 

Table A.1: Parameters used in the VHS molecular model and characteristic 

temperatures. ...................................................................................... 83 

Table A.2: Internal energy relaxation parameters. ............................................... 83 

Table A.3: Parameters used in new vibrational relaxation model. ........................ 83 

Table A.4: Electronic energy levels of O2. ............................................................ 83 

Table A.5: Electronic energy levels of N2. ............................................................ 84 

Table A.6: Electronic energy levels of O. ............................................................. 84 

Table A.7: Electronic energy levels of N. ............................................................. 85 

Table A.8: Electronic energy levels of NO. ........................................................... 87 

Table A.9: Electronic energy levels of O2

+

. .......................................................... 87 

Table A.10: Electronic energy levels of N2

+

. ........................................................... 88 

Table A.11: Electronic energy levels of O+

. ............................................................ 88 

Table A.12: Electronic energy levels of N+

. ............................................................ 89 

Table A.13: Electronic energy levels of NO+. ......................................................... 89 

Table A.14: Endothermic chemical reaction rate coefficients (m3/molecule/s) used 

in the TCE model (the first rate listed for each reaction is the rate used 

in the simulations). kf = aTbexp(-Ea/kT) ............................................. 89 

Table A.15: Exothermic chemical reaction rate coefficients (m3

/molecule/s) used in 

the TCE model. kr = aTbexp(-Ea/kT) .................................................. 93 

Table A.16: Coefficients required for the collision volume in Q-K recombination 

reactions (Vcoll = aTbVref). .................................................................... 94 

Table A.17: Coefficients required for the collision volume in Q-K exchange 

reactions (Ea = a(T/Tref)

b

|Eexch|; Ea = Ea – Eexch if Eexch < 0.0). ............. 95 

Table B.1: Sampled temperatures and their departure from the set equilibrium 

temperature. ..................................................................................... 103 

Table B.2: Sampled and Boltzmann distributions for O2 at 10,000 K. ................ 104 

Table B.3: Sampled and Boltzmann distributions for O2 at 20,000 K. ................ 104 

Table B.4: Sampled and Boltzmann distributions for O2 at 30,000 K. ................ 104 

Table B.5: Sampled and Boltzmann distributions for N2 at 10,000 K. ................ 105 

Table B.6: Sampled and Boltzmann distributions for N2 at 20,000 K. ................ 105 

Table B.7: Sampled and Boltzmann distribution for N2 at 30,000 K. .................. 105 

Table B.8: Sampled and Boltzmann distributions for O at 10,000 K. .................. 105 

Table B.9: Sampled and Boltzmann distributions for O at 20,000 K. .................. 106 

Table B.10: Sampled and Boltzmann distributions for O at 30,000 K. .................. 106 

Table B.11: Sampled and Boltzmann distributions for N at 10,000 K. .................. 106 

Table B.12: Sampled and Boltzmann distributions for N at 20,000 K. .................. 107 

Table B.13: Sampled and Boltzmann distributions for N at 30,000 K. .................. 107 

Table B.14: Sampled and Boltzmann distributions for NO at 10,000 K. ............... 107 

Table B.15: Sampled and Boltzmann distributions for NO at 20,000 K. ............... 107 

Table B.16: Sampled and Boltzmann distributions for NO at 30,000 K. ............... 108 

Table B.17: Sampled and Boltzmann distributions for O2

+

 at 10,000 K. ............... 108 

Table B.18: Sampled and Boltzmann distributions for O2

+

 at 20,000 K. ............... 108 
vii 

Table B.19: Sampled and Boltzmann distributions for O2

+

 at 30,000 K. ............... 108 

Table B.20: Sampled and Boltzmann distributions for N2

+

 at 10,000 K. ............... 109 

Table B.21: Sampled and Boltzmann distributions for N2

+

 at 20,000 K. ............... 109 

Table B.22: Sampled and Boltzmann distributions for N2

+

 at 30,000 K. ............... 109 

Table B.23: Sampled and Boltzmann distributions for O+

 at 10,000 K. ................ 110 

Table B.24: Sampled and Boltzmann distributions for O+

 at 20,000 K. ................ 110 

Table B.25: Sampled and Boltzmann distributions for O+

 at 30,000 K. ................ 110 

Table B.26: Sampled and Boltzmann distributions for N+ at 10,000 K. ................ 110 

Table B.27: Sampled and Boltzmann distributions for N+ at 20,000 K. ................ 110 

Table B.28: Sampled and Boltzmann distributions for N+ at 30,000 K. ................ 111 

Table B.29: Sampled and Boltzmann distributions for NO+ at 10,000 K. .............. 111 

Table B.30: Sampled and Boltzmann distributions for NO+

 at 20,000 K. .............. 111 

Table B.31: Sampled and Boltzmann distributions for NO+ at 30,000 K. .............. 112 
viii 

List of Figures 

Figure 2.1: Selection methodology for inelastic and elastic collisions. .................. 22

Figure 2.2: Selection methodology for inelastic collisions if the sum of the 

probabilities is greater than unity. ....................................................... 23

Figure 2.3: Comparison of Q-K dissociation rate of the reaction N2 + N2 → N + N + 

N2 with those from literature. ............................................................... 32

Figure 2.4: Comparison of Q-K recombination rate of atomic nitrogen with those 

from literature. ..................................................................................... 34

Figure 2.5: Comparison of Q-K exchange rate for the endothermic reaction ........ 36

Figure 2.6: Comparison of Q-K exchange rate for the exothermic reaction .......... 37

Figure 3.1: Sampled and analytic reaction rates when dissociation competes with 

ionization. ............................................................................................ 46

Figure 4.1: Comparison of sampled and reported electronic energy level transition 

rates from the first to fifth energy levels of atomic nitrogen. ................ 50

Figure 4.2: Comparison of sampled and reported electronic energy level transition 

rates from the second to the fifth energy levels of atomic nitrogen. .... 51

Figure 4.3: Comparison of sampled and reported electronic energy level transition 

rates from the first to fifth energy levels of atomic oxygen. ................. 51

Figure 4.4: Comparison of sampled and reported electronic energy level transition 

rates from the second to the fifth energy levels of atomic oxygen. ..... 52

Figure 4.5: Comparison of sampled and reported electronic energy level transition 

rates from the second to the third energy levels of molecular nitrogen.52

Figure 4.6: Comparison of sampled and reported electronic energy level transition 

rates from the first to the fourth energy levels of nitric oxide. .............. 53

Figure 4.7: Comparison of sampled and reported electronic energy level transition 

rates from the first to the third energy levels of the molecular nitrogen 

ion. ...................................................................................................... 53

Figure 4.8: Comparison of sampled and reported electronic energy level transition 

rates for the molecular nitrogen ion when colliding with molecular 

nitrogen. .............................................................................................. 54

Figure 4.9: Comparison of sampled and reported electronic energy level transition 

rates for molecular nitrogen when colliding with atomic nitrogen. ....... 55

Figure 4.10: Comparison of sampled and reported electronic energy level transition 

rates for nitric oxide when colliding with molecular nitrogen. .............. 55

Figure 4.11: Comparison of sampled and reported electronic energy level transition 

rates for the nitric oxide when colliding with molecular nitrogen. ........ 56

Figure 4.12: Comparison of sampled and reported electronic energy level transition 

rates for molecular oxygen when colliding with atomic oxygen. .......... 56

Figure 4.13: Endothermic associative reaction rates for 

compared to rates from the literature. ................................................. 58

Figure 4.14: Exothermic associative reaction rates for compared 

to rates from the literature. .................................................................. 58

Figure 4.15: Endothermic ionization reaction rates for 

compared to rates from the literature. ................................................. 59

N + O → NO+ + e−

NO+ + e− → N + O

O + e

− →O+ + e− + e−
ix 

Figure 4.16: Endothermic charge exchange reaction rates for to 

rates from the literature. ...................................................................... 60

Figure 4.17: Exothermic charge exchange reaction rates for 

compared to rates from the literature. ................................................. 60

Figure 4.18: Endothermic reaction rate for the exchange reaction involving charged 

species with comparisons to rates from the 

literature. ............................................................................................. 61

Figure 4.19: Exothermic reaction rate for the exchange reaction involving charged 

species with comparisons to rates from the 

literature. ............................................................................................. 61

Figure 5.1: Schematic of FIRE II test vehicle. ....................................................... 63

Figure 5.2: Translational temperature comparison between CFD and Case 1 (no 

electronic energy levels). .................................................................... 67

Figure 5.3: Comparison of temperature profiles between the CFD solution and 

Case 2 (electronic energy levels added). ............................................ 68

Figure 5.4: Comparison of neutral species mole fractions between the CFD 

solution and Case 2. ........................................................................... 69

Figure 5.5: Comparison of charged species mole fractions between the CFD 

solution and Case 2. ........................................................................... 69

Figure 5.6: Comparison of temperature profiles between the CFD solution and 

Case 3 (with electronic energy levels and Q-K chemistry). ................. 70

Figure 5.7: Comparison of neutral species mole fractions between the CFD 

solution and Case 3. ........................................................................... 71

Figure 5.8: Comparison of charged species mole fractions between the CFD 

solution and Case 3. ........................................................................... 71 

Figure B.1: Transitions for O2 + O2. ....................................................................... 97 

Figure B.2: Transitions for O2 + N2. ....................................................................... 97 

Figure B.3: Transitions for O2 + O. ........................................................................ 97 

Figure B.4: Transitions for O2 + N. ........................................................................ 97 

Figure B.5: Transitions for O2 + NO....................................................................... 98 

Figure B.6: Transitions for N2 + O2. ....................................................................... 98 

Figure B.7: Transitions for N2 + N2. ....................................................................... 98 

Figure B.8: Transitions for N2 + O. ........................................................................ 98 

Figure B.9: Transitions for N2 + N. ........................................................................ 98 

Figure B.10: Transitions for N2 + NO. ...................................................................... 98 

Figure B.11: Transition from level 1-2 for N2 + e. .................................................... 99 

Figure B.12: Transition from level 2-3 for N2 + e. .................................................... 99 

Figure B.13: Transition from level 1-2 for O + e. ..................................................... 99 

Figure B.14: Transition from level 1-3 for O + e. ..................................................... 99 

Figure B.15: Transition from level 1-4 for O + e. ..................................................... 99 

Figure B.16: Transition from level 1-5 for O + e. ..................................................... 99 

Figure B.17: Transition from level 2-3 for O + e. ................................................... 100 

Figure B.18: Transition from level 2-5 for O + e. ................................................... 100 

Figure B.19: Transition from level 1-2 for N + e. .................................................... 100 

Figure B.20: Transition from level 1-3 for N + e. .................................................... 100 

N+ + N2 → N2

+ + N

N2

+ + N → N+ + N2

NO+ + N ↔ N2

+ + O

N2

+ + O ↔ NO+ + N
x 

Figure B.21: Transition from level 1-4 for N + e. .................................................... 100 

Figure B.22: Transition from level 1-5 for N + e. .................................................... 100 

Figure B.23: Transition from level 2-3 for N + e. .................................................... 101 

Figure B.24: Transition from level 2-4 for N + e. .................................................... 101 

Figure B.25: Transition from level 2-5 for N + e. .................................................... 101 

Figure B.26: Transition from level 3-4 for N + e. .................................................... 101 

Figure B.27: Transition from level 3-5 for N + e. .................................................... 101 

Figure B.28: Transitions for NO + O2. .................................................................... 101 

Figure B.29: Transitions for NO + N2. .................................................................... 102 

Figure B.30: Transitions for NO + O. ..................................................................... 102 

Figure B.31: Transitions for NO + N. ..................................................................... 102 

Figure B.32: Transitions for NO + NO. .................................................................. 102 

Figure B.33: Transition from level 1-3 for NO + e. ................................................. 102 

Figure B.34: Transition from level 1-4 for NO + e. ................................................. 102 

Figure B.35: Transition from level 3-4 for NO + e. ................................................. 103 

Figure B.36: Transition from level 1-2 for N2

+

 + e. ................................................. 103 

Figure B.37: Transition from level 1-3 for N2

+

 + e. ................................................. 103 

Figure B.38: Transition from level 2-3 for N2

+

 + e. ................................................. 103 

Figure C.1: Q-K reaction rates for reaction 1. ...................................................... 113 

Figure C.2: Q-K reaction rates for reaction 2. ...................................................... 113 

Figure C.3: Q-K reaction rates for reaction 3. ...................................................... 113 

Figure C.4: Q-K reaction rates for reaction 4. ...................................................... 113 

Figure C.5: Q-K reaction rates for reaction 5. ...................................................... 114 

Figure C.6: Q-K reaction rates for reaction 6. ...................................................... 114 

Figure C.7: Q-K reaction rates for reaction 7. ...................................................... 114 

Figure C.8: Q-K reaction rates for reaction 8. ...................................................... 114 

Figure C.9: Q-K reaction rates for reaction 9. ...................................................... 114 

Figure C.10: Q-K reaction rates for reaction 10. .................................................... 114 

Figure C.11: Q-K reaction rates for reaction 11. .................................................... 115 

Figure C.12: Q-K reaction rates for reaction 12. .................................................... 115 

Figure C.13: Q-K reaction rates for reaction 13. .................................................... 115 

Figure C.14: Q-K reaction rates for reaction 14. .................................................... 115 

Figure C.15: Q-K reaction rates for reaction 15. .................................................... 115 

Figure C.16: Q-K reaction rates for reaction 16. .................................................... 115 

Figure C.17: Q-K reaction rates for reaction 17. .................................................... 116 

Figure C.18: Q-K reaction rates for reaction 18. .................................................... 116 

Figure C.19: Q-K reaction rates for reaction 19. .................................................... 116 

Figure C.20: Q-K reaction rates for reaction 20. .................................................... 116 

Figure C.21: Q-K reaction rates for reaction 21. .................................................... 116 

Figure C.22: Q-K reaction rates for reaction 22. .................................................... 116 

Figure C.23: Q-K reaction rates for reaction 23. .................................................... 117 

Figure C.24: Q-K reaction rates for reaction 24. .................................................... 117 

Figure C.25: Q-K reaction rates for reaction 25. .................................................... 117 

Figure C.26: Q-K reaction rates for reaction 26. .................................................... 117 

Figure C.27: Q-K reaction rates for reaction 27. .................................................... 117 

Figure C.28: Q-K reaction rates for reaction 28. .................................................... 117 
xi 

Figure C.29: Q-K reaction rates for reaction 29. .................................................... 118 

Figure C.30: Q-K reaction rates for reaction 30. .................................................... 118 

Figure C.31: Q-K reaction rates for reaction 31. .................................................... 118 

Figure C.32: Q-K reaction rates for reaction 32. .................................................... 118 

Figure C.33: Q-K reaction rates for reaction 33. .................................................... 118 

Figure C.34: Q-K reaction rates for reaction 34. .................................................... 118 

Figure C.35: Q-K reaction rates for reaction 35. .................................................... 119 

Figure C.36: Q-K reaction rates for reaction 36. .................................................... 119 

Figure C.37: Q-K reaction rates for reaction 37. .................................................... 119 

Figure C.38: Q-K reaction rates for reaction 38. .................................................... 119 

Figure C.39: Q-K reaction rates for reaction 39. .................................................... 119 

Figure C.40: Q-K reaction rates for reaction 40. .................................................... 119 

Figure C.41: Q-K reaction rates for reaction 41. .................................................... 120 

Figure C.42: Q-K reaction rates for reaction 42. .................................................... 120 

Figure C.43: Q-K reaction rates for reaction 43. .................................................... 120 

Figure C.44: Q-K reaction rates for reaction 44. .................................................... 120 

Figure C.45: Q-K reaction rates for reaction 45. .................................................... 120 

Figure C.46: Q-K reaction rates for reaction 46. .................................................... 120 

Figure C.47: Q-K reaction rates for reaction 47. .................................................... 121 

Figure C.48: Q-K reaction rates for reaction 48. .................................................... 121 

Figure C.49: Q-K reaction rates for reaction 49. .................................................... 121 

Figure C.50: Q-K reaction rates for reaction 50. .................................................... 121 

Figure C.51: Q-K reaction rates for reaction 51. .................................................... 121 

Figure C.52: Q-K reaction rates for reaction 52. .................................................... 121 

Figure C.53: Q-K reaction rates for reaction 53. .................................................... 122 

Figure C.54: Q-K reaction rates for reaction 54. .................................................... 122 

Figure C.55: Q-K reaction rates for reaction 55. .................................................... 122 

Figure C.56: Q-K reaction rates for reaction 56. .................................................... 122 

Figure C.57: Q-K reaction rates for reaction 57. .................................................... 122 

Figure C.58: Q-K reaction rates for reaction 58. .................................................... 122 

Figure C.59: Q-K reaction rates for reaction 59. .................................................... 123 

Figure C.60: Q-K reaction rates for reaction 60. .................................................... 123 

Figure C.61: Q-K reaction rates for reaction 61. .................................................... 123 

Figure C.62: Q-K reaction rates for reaction 62. .................................................... 123 

Figure C.63: Q-K reaction rates for reaction 63. .................................................... 123 

Figure C.64: Q-K reaction rates for reaction 64. .................................................... 123 

Figure C.65: Q-K reaction rates for reaction 65. .................................................... 124 

Figure C.66: Q-K reaction rates for reaction 66. .................................................... 124 

Figure C.67: Q-K reaction rates for reaction 67. .................................................... 124 

Figure C.68: Q-K reaction rates for reaction 68. .................................................... 124 

 
xii 

Nomenclature 

a constant in inverse power potential and Arrhenius equation

B magnetic field (T) 

b constant in Arrhenius equation

C1,C2 constants in the Millikan-White vibrational relaxation model and TCE 

chemistry model

D non-collisional operator defined for the Boltzmann equation

E electric field (N/C), energy (J) 

e electron charge (C)

F force acting on a particle (N)

f distribution function 

g relative velocity (m/s), degeneracy 

h Planck’s constant 

i energy level

J collisional operator defined for the Boltzmann equation 

k Boltzmann constant 

kf, kr forward/reverse reaction or transition rate

M normalized collision mass

m mass (kg) 

mr reduced mass 

N number of particles in a computational cell 

n number density (m-3) 

P probability 

p pressure (N/m2

)

QA partition function of species A

Q(v) probability density function for the post-collision velocity v 

Q(a,b) incomplete Gamma function 

q charge on a particle (C) 

R random number, Collision rate parameter 

T temperature (K)

t time (s) 

W number of real particles represented by a simulated particle 

V volume (m3)

u,v,w velocity in the x-,y-,z-directions (m/s) 

Z collision number

x,y,z Cartesian coordinates (m) 

α exponent in inverse power potential 

χ,ψ collision angles (rad) 

δ Dirac delta function 

ε symmetry factor (=1 if dissimilar, =2 if same), specific energy (J/kg) 

ε 0 permittivity of free space 

φ electric potential 

Υ fraction of collisions that meet conditions for chemical reaction 

Γ( ) a gamma function of a 
xiii 

μ viscosity 

Θ characteristic temperature (K) 

σ collision cross-section (m2

) 

τ mean free time (s) 

Ω solid angle in a collision (strad) 

ω exponent in VHS interaction model 

ζ degrees-of-freedom 

Subscripts 

AI associative ionization

a activation 

ambipolar ambipolar electric field

chex charge exchange

coll collision 

diss dissociation

e electron 

el electronic energy mode 

eq equilibrium 

exch exchange 

f forward 

int internal 

ion ionization

max maximum value 

r reverse 

react reaction 

rec recombination 

ref reference

rot rotational energy mode 

T total 

trans translational energy 

vib vibrational energy mode 

Superscripts 

AI associative ionization

chex charge exchange

cont continuum value 

DSMC DSMC value 

ion ionization
1 

Section 1 

Introduction

The aerothermodynamic environments experienced by hypersonic vehicles are characterized 

by strong shock waves, which produce high temperatures resulting in high heating of the 

vehicle. The severity of the environment is dependent upon the mission profile and the 

vehicle configuration[1]. A high-speed planetary-entry blunt-body vehicle, such as the 

Galileo probe entering Jupiter or a high speed Earth return capsule, encounters an 

environment dominated by high enthalpy effects, which include elevated levels of ionization 

and may result in significant radiative heating. The aerodynamic and aerothermal loads on a 

hypersonic vehicle set the performance requirements for various sub-systems, such as the 

thermal protection system (TPS) and the control system. The design and safe operation of 

these vehicles, therefore, require adequate definition of the flight environment.

During the design phase, the aerothermodynamic performance of a hypersonic vehicle is 

defined by use of either experimental or computational tools. However, as in most flight 

regimes, these definitions cannot be solely obtained from ground test facilities, as no facility 

can reproduce all aspects of the flight environment. This limitation is particularly true for the 

hypersonic flight regime because of the high energy requirement needed to create a true 

hypersonic environment at a reasonable scale on the ground, and especially applicable to the 

transitional/rarefied portion of the trajectory because of the very low density requirement. 

Therefore, designers rely on computational tools for the prediction of the aerothermal loads 

experienced by vehicles entering or re-entering the atmosphere from space.

In rarefied flows, because of the breakdown of the continuum assumption, kinetic methods 

are employed. Since its introduction in the early 1960s, the direct simulation Monte Carlo

(DSMC) method of Bird[2] has become the de facto standard for simulating low-density gas 

dynamics. For these types of flows, traditional computational fluid dynamics (CFD) methods 

are invalid because assumptions made in developing the differential equations on which they 

are based break down under rarefied conditions. In contrast to continuum CFD methods, the 

DSMC method performs a direct physical simulation of the gas at the molecular level that 

closely mimics the Boltzmann equations. In the simulation, particles are tracked in space and 

time, accounting for both gas surface interactions and intermolecular collisions.

In the upcoming years, the ability to compute rarefied, ionized hypersonic flows will become 

increasingly important. Missions such as Earth (re)entry, landing high mass payloads on 

Mars, and the exploration of the outer planets and their satellites will continue to push the 

state of the art in DSMC simulations. These missions may result in shock-layer temperatures 

reaching into the tens-of-thousands of degrees, where the treatment of electronic energy 

levels and ionization processes becomes increasingly important. The proper treatment of 

these mechanisms is needed to accurately predict overall shock wave thickness and structure, 
2 

as well as heat transfer to the vehicle. Even when the reactants are minor species, the overall 

flow field is not affected by details of the energy modes and reaction model, but a particular 

radiation signature may be strongly affected[3]. Either way, the accurate treatment and 

prediction of electronic energy distribution and ionization is an important component of the 

design of vehicle sub-systems.

1.1 Objectives of this Work

Although there has been much recent interest in the modeling of ionization processes[4-9], 

the inclusion of electronic energy levels in DSMC simulations has largely been overlooked 

and they are not included in many DSMC implementations. However, recent advances in 

DSMC algorithms[10-14], specifically the Quantum-Kinetic (Q-K) Chemistry Model, have 

allowed the prediction of equilibrium and non-equilibrium reaction rates using only kinetic 

theory and fundamental molecular properties (i.e., no macroscopic reaction-rate information). 

Models such as these are important especially when reaction rate data are limited or the gas is 

far from equilibrium and an engineering approximation is required. Thus far, these new 

models have only covered neutral-species reactions. It is the purpose of this thesis to expand 

the Q-K model to include electronic energy level transitions and also to include reactions 

involving charged species.

In addition to ionization, radiation can become an important factor depending on the mission 

and thermal protection system. For example, inflatable ballutes can be used to decelerate 

spacecraft while still high in the atmosphere. The material used for the skin of the ballute 

would have a very low tolerance for heating, necessitating the knowledge of what the 

radiative heating load would be. Several methods of calculating radiative heating for 

rarefield flows[15-18] have been developed. For the current study, however, an uncoupled 

methodology used in many computational fluid dynamics radiation evaluations[19, 20] was 

adopted. While this study does not improve on the methodologies used for calculating the 

radiative heat transfer, having separate species non-equilibrium electronic temperatures 

would make radiative predictions much more accurate.

The models developed herein will be validated against reaction rates reported in the literature 

and applied to the second flight of the Flight Investigation Reentry Environment (FIRE II) 

Project. This was an Apollo-era flight experiment to measure the radiative and convective 

heating on spacecraft material during atmospheric reentry at lunar return speeds. The 

trajectory point at 76.42 km was chosen for the computations because, while the flow is 

nearly in the transitional regime and DSMC is more appropriate, computational fluid 

dynamics (CFD) has frequently been applied at this altitude with success. Therefore, the 

DSMC simulations presented in this thesis will be compared to CFD solutions and flight data 

from this trajectory point. However, while the methods developed herein are completely 

applicable to higher altitudes, the use of CFD is not recommended because of the break down 

of the assumption of continuum flow.

1.2 Outline of Thesis

The remainder of this thesis is organized as follows. The second Section provides an 

overview of the governing equations that are the basis for the simulation of fluid flows from 
3 

the continuum through rarefied regimes. Also, a description of the methodology and 

structure of the DSMC algorithm, both “traditional” DSMC and the more recent algorithmic 

advances, is presented as well as previously existing electronic energy level treatments and 

ionized gas models. Section Three describes the extensions to the Q-K model and their 

implementation. In Section Four, these extensions are applied to an adiabatic heat bath and 

sampled transition/reaction rates are compared to those rates found in the literature. In 

Section Five, these extensions are then compared against previous chemical reaction models 

as applied to the second Flight Investigation of Reentry Environment (FIRE II) 

experiment[21]. Finally, Section Six summarizes the results of the thesis, the contributions 

of the thesis to the state-of-the-art, and recommended future research directions.
4 

Section 2 

Simulation of Rarefied and Weakly Ionized Reentry Flows

In this Section, the theoretical basis of particle simulation is first summarized. Next, a 

description of the direct simulation Monte Carlo numerical method is given, including the 

traditional methods and more recent algorithmic advances, as well as previous attempts to 

include electronic energy levels and models of ionized flow. Throughout this description, 

mention will be given of any modifications that were required to be performed on the DSMC 

algorithm of choice.

2.1 Theoretical Basis of Particle Simulation

The particle simulation method has been developed independently by plasma physicists and 

aerothermodynamisists. Plasma physicists, who simulate charged particles, have developed 

the particle-in-cell (PIC) method based on the theory of charge assignment and force 

interpolation. This theory has made it possible to obtain smooth distributions of charge 

density and current density in Maxwell’s equation of electromagnetic fields. The details of 

the PIC method can be found in the classical works of Birdsall and Langdon[22] and 

Hockney and Eastwood[23]. On the other hand, aerothermodynamicists have been 

developing a particle simulation method of neutral species beginning with the work of 

Bird[2, 24]. This method is called the direct simulation Monte Carlo (DSMC) method.

Particle simulation of both charged and neutral species is necessary to examine the structures 

of high-speed, rarefied, reacting flows. However, for flows where the plasma density is 

considered low (~10%), the assumption of charge neutrality is generally assumed and the 

charged particles are treated as neutral particles[25-30]. More will be discussed on these 

techniques in Section 2.2.4. The current section describes how the probability laws used in 

particle simulation may be derived[31] from a relevant Boltzmann equation for a weakly 

ionized plasma and simplifications will be pointed out for our assumption of charge 

neutrality.

2.1.1 The Boltzmann Equation

The theoretical basis of particle simulation can be explained using a weakly ionized plasma 

consisting of electrons (e), ions (A), and molecules (B). If Coulomb collisions e-e, e-A, and 

A-A can be disregarded, only the remaining e-B, A-B, and B-B collisions need to be 

considered. The particle motions of species e, A, and B are governed by the respective 

Boltzmann equation for each of these species. Except for Coulomb collisions, the Boltzmann 

equation for electrons or ions is linear, but the Boltzmann equation for molecules is nonlinear 

because of the B-B collisions. Although the PIC method was developed for charged species

and the DSMC method for neutral species, the idea behind the two methods is similar. The

Boltzmann equation of ions, A, is written as:

∂ fA

∂t

+ vj

∂ fA

∂xj

+

1

mA

∂

∂vj



Fj f ( ) A = fA



( ) v′ f ( ) w′ − fA



( ) v f  ⎡⎣ ( ) w ⎤⎦ gσ AB dΩd 

w ∫∫ (2.1.1)
5 

Here, Coulomb collisions are neglected and only A-B collisions are considered, is the force 

acting on particle A with a mass of mA, fA( ) and fB( ) are the velocity distribution functions 

of species A and B, and are the post-collision velocities, is the relative 

speed, σAB is the differential cross section, and dΩ (=sinχ dχ dψ) is the solid angle. In the 

remainder of the derivation, and are expressed as and in 

the right hand side of the equation for simplicity. The differential cross section σAB is a 

function of g and the deflection angle χ. The force is given by:



F = qA



E + v ×  ( ) B (2.1.2)

where qA is the charge of the particle A and and are the electric and magnetic fields, 

respectively. However, for our purposes, the force experienced by the charged particles is 

effectively zero since there is an assumption of no charge separation (i.e., charge neutrality). 

For the sake of completeness, though, this term will be kept in the remainder of the 

derivation.

If the method to obtain is given, the solution of fA can be found at any time by 

step-by-step calculations started from an initial velocity distribution function. For a small Δt, 

we have:

fA



v,

 ( ) x,t + Δt = fA

v,

 ( ) x,t + Δt

∂ fA

∂t

⎛

⎝

⎜ ⎞

⎠

⎟ (2.1.3)

Substituting of Equation (2.1.1) into Equation (2.1.3) gives us:

fA



v,

 ( ) x,t + Δt = ( ) 1+ ΔtJ ( ) 1+ ΔtD fA

v,

 ( ) x,t (2.1.4)

where D and J are operators defined by

D f ( ) A = −vj

∂ fA

∂xj

− 1

mA

∂

∂vj

Fj f ( ) A (2.1.5)

J f ( ) A = fA



( ) v′ fB



( ) w′ − fA



( ) v fB

 ⎡⎣ ( ) w ⎤⎦ gσ AB dΩd 

w ∫∫ (2.1.6)

D and J denote non-collisional and collisional operators, respectively, and therefore, 

Equation (2.1.4) decouples the non-collisional and collisional parts of the motion. Equation 

(2.1.4) is equal to Equation (2.1.3) if the term of order (Δt)

2 is disregarded such that the 

solution obtained from Equation (2.1.4) is first-order accurate. In general, Δt should be much 

less than the mean free time τ (the mean free path divided by the mean speed of the 

particles). This condition should be strictly satisfied in the case when particle A is charged 

[32]. Note, however, that if particle A is a molecule, the condition can be replaced by a 

weaker one of Δt < τ. The reason why the weaker condition can be applied is that molecular 



Fj



v w 

v′ w′ g = 

v − w

fA



v,

 ( ) x,t fB

v,

 ( ) x,t fA

( ) v fB



( ) v



F



E 

B

fA



v,

 ( ) x,t

∂fA / ∂t
6 

trajectories are straight and, therefore, accurate for any Δt and that the state of the gas 

changes little in the mean free time.

Equation (2.1.4) suggests that a method of solving the Boltzmann equation is to first obtain

(non-collisional) and then obtain (collisional). The first 

procedure is equivalent to solving:

∂ fA

∂t

+ vj

∂ fA

∂xj

+

1

mA

∂

∂vj

Fj f ( ) A = 0 (2.1.7)

by use of the initial distribution function . If the solution at t+Δt is , then 

the solution is , and therefore

fA

* 

v + Δv,



x + Δ ( x,t) = fA



v,

 ( ) x,t (2.1.8)

where and . Equation (2.1.8) implies that the particles of species 

A move according to the equation of motion

mA

d



v

dt = qA



E + v ×  ( ) B

d



x

dt = 

v

(2.1.9)

These equations of motion are used in PIC simulations which means that the Boltzmann 

equation without the collision term is that solved by the PIC method. When Equations 

(2.1.9) are solved for all particles, which is trivial for DSMC since there is no net force, we 

have . The second procedure is equivalent to solving

∂ fA

∂t = fA



( ) v′ fB



( ) w′ − fA



( ) v fB

 ⎡⎣ ( ) w ⎤ ∫∫ ⎦gσ ABdΩd 

w (2.1.10)

using the initial distribution function . The solution of Equation (2.1.10) at t+Δt is

.

Equation (2.1.10), which does not include the term , represents the collisional 

relaxation of the velocity distribution in a spatially uniform state of species A. Equation

(2.1.10) can be applied only to a small cell where the gradient of density, temperature, and 

flow velocity of species A can be regarded as uniform. Because of this requirement, the cell

dimension should be of the order of the mean free path of species A, and is the basis for the 

cell size of the DSMC method. 

A solution of equation (2.1.10) can be obtained by an equivalent stochastic process[33]. 

First, in particle simulation, it is sufficient to consider the motions of a set of particles 

fA

* ≡ (1+ ΔtD)fA ( ) 1+ ΔtJ fA

*

fA



v,

 ( ) x,t fA

* 

v,

 ( ) x,t

fA

* 

v,

 ( ) x,t = fA

v − 

( ) F / mA Δt,



x −  ( vΔt,t)

Δv = 

( ) F / mA Δt Δ

x = vΔt

fA

* 

v,

 ( ) x,t

fA

* 

v,

 ( ) x,t

fA



v,

 ( ) x,t + Δt

∂fA / ∂xj
7 

randomly sampled from all particles. The sampled particles are referred to as simulated 

particles. The number of simulated particles in a cell is much smaller than is the number of 

real particles. The weight, W, is defined as the number of real particles represented by each 

simulated particle. After moving all simulated particles by Δt, we identify the simulated 

particles in each cell by examining the positions of the particles. Let NA be the number of 

simulated particles in a cell and be the particle velocities. The velocity 

distribution function for a cell located at is expressed as

fA

* 

v,

 ( ) x,t = nA

NA

δ 3 

v − vAi ( )

i=1

NA

∑ (2.1.11)

where nA is the number density of the cell, δ3( ) is the Dirac delta function, and is the 

solution of Eq. (2.1.9) at time t+Δt. The function is the probability density 

function for the velocity of particle Ai. Similarly, the velocity distribution function of 

species B at the end of the particle’s motion is

fB

* 

w,

 ( ) x,t = nB

NB

δ 3 

w − v ( ) Bj

j=1

NB

∑ (2.1.12)

where is the velocity of simulated particle Bj. The solution of can be 

obtained from

fA



v,

 ( ) x,t + Δt = ( ) 1+ ΔtJ fA

* 

v,

 ( ) x,t

= fA

* 

v,

 ( ) x,t + Δt fA

( ) v′ fB



( ) w′ − fA



( ) v fB

 ⎡⎣ ( ) w ⎤⎦ gσ AB dΩd 

w ∫∫ (2.1.13)

Substitution of Equations (2.1.11) and (2.1.12) into Equation (2.1.13) yields

fA



v,

 ( ) x,t + Δt = nA

NA

fAi



( ) v

i=1

NA

∑ (2.1.14)

where is the probability density function for the velocity of particle Ai at t+Δt. It is 

given by

fAi



( ) v = 1− PAi ( )δ 3 

v − vAi ( ) + PAi

QAi



( ) v (2.1.15)

The structure of Eq. (2.1.15) states that PAi is the collisional probability of particle Ai. Note 

that if PAi = 0, , and there is no change in velocity because of a collision. 

The probability PAi is given by

PAi = PAi ,Bj

j=1

NB

∑ (2.1.16)



vA1,



vA2 ,...

fA

* 

v,

 ( ) x,t 

x



vAi

δ3 

v − v ( ) Ai



vAi



vBi fA



v,

 ( ) x,t + Δt

fAi



( ) v

fAi



( ) v = δ3 

v − v ( ) Ai
8 

where

PAi ,Bj = NB

−1

nBgAi ,Bjσ T

AB

gAi ( ) ,Bj Δt (2.1.17)

Here, is the relative speed and is the total collision cross section. 

Equation (2.1.16) reveals that PAi,Bj is the probability of particle Ai colliding with particle Bj

in Δt. The function is the probability density function for the post-collision velocity

. It is written as:

QAi



( ) v = PAi ,Bj

PAi

QAi ,Bj



( ) v

j=1

NB

∑ (2.1.18)

Here, the ratio PAi,Bj/PAi represents the conditional probability that the collision partner of 

particle Ai is particle Bj under the condition that particle Ai collides in Δt. Here, is 

the probability density function for the post-collision velocity of particle Ai when its collision 

partner is particle Bj, which is given by:

QAi ,Bj



( ) v =

σ ABgAi ,Bj ( ) ,χ δ 

g′ − gAi ( ) ,Bj

gAi ( ) ,Bj

2

MB

3

σ T

AB

gAi ( ) ,Bj

(2.1.19)

where



g′ = 1

MB



v − MA



vAi − MB



v ( Bj) (2.1.20)

MA = mA/(mA+mB), MB = mB/(mA+mB), χ is the angle between and , and 

σAB(gAi,Bj,χ) is the differential cross section. From Eq. (2.1.19), we have

QAi ,Bj



( ) v d



v = δ 

g′ − gAi ( ) ,Bj dg′

σ ABgAi ,Bj ( ) ,χ sin χdχdψ

σ T

AB

gAi ( ) ,Bj

(2.1.21)

where . Note that is the Jacobian of the 

transformation from to , and then Eq. (2.1.20) yields the post-collision velocity

of particle Ai



v′

Ai = 1

mA + mB

mA



vAi + mB



vBj + mB

 ( g′) (2.1.22)

gAi,Bj = 

vAi − vBj σT

AB

QAi



( ) v 

v

QAi,Bj



( ) v



g′ 

vAi − 

vBj

d



v = MB

3

d



g′ = MB

3

g′

2

dg′sin χdχdψ MB

3



v 

g′ 

v′

Ai ≡  ( ) v
9 

where is the post-collision relative velocity and . The probability 

that is included in the solid angle is . 

The direction of can be sampled using this probability.

Equations (2.1.16) and (2.1.17), which are derived from the Boltzmann equation, can also be 

obtained by use of the elementary free-path theory. The collision frequency of particle Ai is

vAi = nB gAi ,Bjσ T

AB

gAi ( ) ,Bj

where gAi,B is the relative speed between Ai and some particle B, and the bar denotes the 

average of particles B. Therefore, the average of all simulated particles B in a cell and the 

collision probability of particle Ai take the form

PAi = vAiΔt = nBΔt

NB

gAi ,Bjσ T

AB

gAi ( ) ,Bj

j=1

NB

∑

This equation agrees with Equation (2.1.16). We cannot obtain the probability laws to 

determine the post-collision velocities from the free-path theory, however. In many cases, 

we can write such probability laws intuitively because the laws derived from the Boltzmann 

equation agree with the physical image of the collision.

The probability laws used to determine the velocity of particle Ai at t+Δt are then: 

1) Particle Ai collides in Δt with a probability of PAi.

2) Under the condition that particle Ai collides, its collision partner is particle Bj with a 

probability of PAi,Bj/PAi.

3) If the collision partner is Bj, the post-collision velocity is given by Equation

(2.1.22).

Therefore, the particle simulation method can be derived from the Boltzmann equations of 

ions and is true for any species. Although the Boltzmann equation for a simple gas is 

nonlinear because of collisions between like molecules, the simulation method was derived 

in 1980[34]. Before this work, those in the field of rarefied gas dynamics had regarded the 

particle simulation method as a kind of numerical experiment. By use of the measure theory, 

mathematicians[35] verified that Nanbu’s particle simulation method is an appropriate 

method for solving the Boltzmann equation.

The most important ideas of the particle simulation method can be summarized as: 1) 

decoupling of collisionless motion and collisions; and 2) calculating collisions independently 

in each cell. These basic concepts have been shown to naturally result from the operator 

splitting in the form of (1+ΔtJ)(1+ΔtD). Now consider the conditions for choosing Δt. The 

operator splitting is allowed if and only if and . The first 



g′ = 

v′

Ai − 

v′ ( ) Bj



g′ = gAi,Bj



g′ sin χdχdψ σAB

gAi,Bj ( ) , χ sin χdχdψ / σT

AB

g( ) Ai,Bj

( ) χ,ψ 

g′



v′

Ai

Δt J f ( ) A  fA Δt D f ( ) A  fA
10 

condition states that the effect of collisions in Δt must be small; therefore the collision 

probability is small, where τ is the mean free time. The time step must satisfy 

the condition . The second condition requires that and in Eq. 

(2.1.8). That is, and , where is the mean particle speed and 

is the mean acceleration. The condition that is usually replaced by the 

Courant condition , where Δc denotes the cell size. Therefore, a particle’s 

displacement in Δt must be less than the cell size. The condition requires 

that the change in the particle’s speed in Δt must be much less than the original speed for PIC 

simulations. Therefore, a smaller Δt should be used for an external field.

2.2 Traditional Direct Simulation Monte Carlo Methods

The DSMC technique was developed primarily by Graeme Bird, now Emeritus Professor at

the University of Sydney, Australia, nearly a half century ago. Since its inception, this 

technique has been used in many different applications ranging from the simulation of microelectromechanical devices (MEMs), to various industrial and pharmacological deposition 

simulations, to the simulation of (re)entry and other rarefied flow fields.

In the work outlined herein, the DSMC method is used to simulate the Boltzmann equation 

for rarefied/transitional reentry flows in which the Navier-Stokes equations are not valid 

throughout the entire flow field. Specifically, the algorithm employed is the DSMC Analysis 

Code (DAC)[36, 37], which was developed by a team led by LeBeau beginning at the NASA 

Langley Research Center and continuing at the NASA Johnson Space Center. The DSMC 

algorithms in the release version of DAC are based on those developed by Bird[2].

Regardless of the DSMC code used, the basic structure of the algorithm is the same. The 

following is included to briefly describe the flow of a generic DSMC algorithm.

2.2.1 Collision Algorithm

In each computational cell, particles are chosen to collide in such a way as to reproduce the 

required macroscopic collision rate. The simplest implementation of collision pair selection 

is to disregard particle position within the cell when selecting the pair. This methodology 

restricts the size of a computational cell in that it must be smaller than the local average mean 

free path of the flow field. This was the first selection methodology and is still utilized 

today. The next selection methodology was to select the first particle randomly in the cell 

and then select the particle closest to the first excluding its previous collision partner, or the 

so-called virtual sub-cell (VSC) algorithm of LeBeau[38]. This methodology allows for a 

larger collision cell size as long as the resulting mean collision distance divided by mean free 

path remains at least as small as the resulting value from the original selection methodology.

Other selection methodologies have been proposed[39, 40], but the two discussed above are 

the most prevalent. The VSC algorithm is that which was used in the current study.

PAi  Δt / τ

Δt  τ Δx  x Δv  v

vΔt  x ( ) F / mA Δt  v v

F / mA vΔt  x

vΔt < Δc

( ) F / mA Δt  v
11 

In many modern DSMC codes, Bird’s No Time Counter (NTC) method is used to select 

colliding pairs[2]. At each time step and in each computational cell, the number of particle 

pairs to be tested for collisions is calculated as follows:

Ncoll = 1

2 nN σ T ( ) g max Δt (2.2.1)

where n is the number density of gas particles in the cell, is the average number of 

simulator particles in the cell, Δt is the simulation time step, g is the relative velocity of a 

colliding pair of particles, and σT is the total collision cross section of the colliding pair. 

Each of these pairs are then selected to collide if the probability:

Pcoll = ( ) σ T g

( ) σ T g max

(2.2.2)

is less than a random number uniformly distributed between zero and one. This method 

assumes that the collision rates between species will be implicitly reproduced, even though 

all species are lumped into the same group for the pair selection process. This method can 

become inefficient if the product of the collision cross section and relative velocity becomes 

very much different between species pairs, as is the case when electrons are introduced into 

the flow, because the maximum value of is very large and results in many null 

collisions. This is the basic algorithm used by the DAC code.

In order to address this issue, the collision algorithm in DAC was modified as follows[31].

The collision algorithm can be divided into two groups. One includes collisions between 

similar particles, and the other collisions between dissimilar particles. The collision term of 

the Boltzmann equation is nonlinear for like collisions and linear for unlike collisions. 

Consequently, the collision algorithm for like collisions is different from that for the unlike 

collisions. 

There are five classes of collisions involving charged particles: electrons with neutrals, ions 

with neutrals, electrons with ions, ions with ions, and electrons with electrons. The 

interactions of the last three classes obey the Coulomb force law:

F = 1

4πε 0

q1q2

r

2 (2.2.3)

where ε 0 is the permittivity of free space and r is the distance between the particles. 

Bird[41] and Taylor et. al. [25], however, used the variable hard sphere (VHS) model of Bird 

to describe the first four classes of collisions by setting the diameter of the heavy ions equal 

to the diameter of their neutral counterparts. The reference diameter of the electrons is set 

equal to 1e-10 m, which is orders of magnitude larger than the classical value of 6e-15 m. 

This choice was somewhat arbitrary, but Bird did mention that the electron diameter was 

varied in his simulations by a factor of three with no significant effect on the flow field 

parameters. Neither modeled electron-electron collisions. This is the approach that will be 

followed herein.

N

σT g
12 

In general, then, one can divide the collision algorithm into two for weakly ionized plasmas, 

neglecting Coulomb collisions: 1) collisions between unlike particles; and 2) collisions 

between like particles. The algorithm for type (1) is derived from the linear Boltzmann 

equation as in Section 2.1.1, and that for type (2) is from the nonlinear Boltzmann equation. 

Let us consider the elastic collisions:

A v( ) A + B v( ) B → A v′ ( ) A + B v′ ( ) B

If gas B is at rest and in equilibrium at temperature TB, the probability that particle A

elastically collides with some particle B is:

Pcoll = nB gσ T ( ) g Δt (2.2.4)

where nB is the number density of gas B, is the relative velocity, and the average 

is taken over the Maxwellian distribution of . The average in Equation(2.2.4) is for all 

simulated particles in the cell, which gives:

Pcoll = nB

NB

gjσ T g( )j Δt

j=1

NB

∑ (2.2.5)

where NB is the number of simulated B particles. If the collision is treated by the hard-sphere 

model, the scattering is isotropic. The post-collision velocities are then:



v′

A = 1

mA + mB

mA



vA + mB



vB + mB



vA − 

vB

 ( R)



v′

B = 1

mA + mB

mA



vA + mB



vB − mA



vA − 

vB

 ( R)

(2.2.6)

where the two R’s are the same unit vector with a random direction. However, the hard 

sphere model is generally insufficient for most applications because it does not match the 

gas’s viscosity or diffusion rate as a function of temperature. On the other hand, the most 

accurate Lennard-Jones potential is not suited for the DSMC method because the collision 

calculation is very time-consuming. But, the simple inverse-power potential is 

sufficiently accurate and very efficient for collision calculations. The differential cross 

section σ for the inverse-power potential takes the form:

σ = Ag−4 α (2.2.7)

where g is the relative speed and A is a function of the deflection angle χ. If A is independent 

of χ, the probability that the post-collision relative velocity is directed into the 

solid angle dΩ becomes dΩ/4π, or the scattering is isotropic and the post-collision velocities 

are the same as that of the hard-sphere model. The variable hard-sphere (VHS) model was 

g = vA − 

vB



vB

V r( ) = a rα

σ dΩ σT
13 

proposed by Bird[41] to avoid the divergence of the total collision cross section and using 

isotropic scattering for the inverse-power potential. In the VHS model, the coefficient A is 

constant and, once a gas is selected, A and α can be chosen to best fit the Chapman-Enskog 

viscosity to experimental data. This leads to a power law temperature dependence of the 

coefficient of viscosity such that:

μ ∝T ω (2.2.8)

where

ω = α + 4

2α

(2.2.9)

and ω is bounded by 0.5 and 1.0 at the hard-sphere and Maxwell molecule limits, 

respectively. The VHS parameters used in the current study are discussed in Section 5.3.

2.2.1.1 Simple Gas

For a single species gas in a computational cell, let be the velocities of the 

simulated particles. The probability Pi that molecule i collides in Δt with some other particle 

is given by Equation (2.2.5). For a simple gas, this is:

Pi = n

N gijσ T g( )ij Δt = Pij

j=1

N

∑

j=1

N

∑ (2.2.10)

where and

Pij = N−1ngijσ T ( ) gij Δt (2.2.11)

Note that Pii = 0 because gii = 0. The term for j = i should therefore be omitted in Equation

(2.2.10).

If gij is then replaced by gmax, the resulting probability is:

( ) Pij max = N−1

ngmaxσ T g( ) max Δt (2.2.12)

This probability is greater than Pij. For Pij = (Pij)max, Equation (2.2.10) takes the form:

( ) Pi max = nN−1

( ) N −1 gmaxσ T g( ) max Δt (2.2.13)

The number of collisions in Δt for (Pi)max is

Nmax = 1

2 ( ) Pi max = 12 n N( ) −1 gmaxσ T g( ) max Δt

i=1

N

∑ (2.2.14)



v1,



v2 ,...,



vN

gij = 

vi − vj
14 

The factor of 1⁄2 on the right hand side of Equation (2.2.14) is to avoid double counting 

collisions between like particles. Nmax is called the maximum collision number. The 

computing time necessary to obtain the real collision number

Ncoll = 1

2 Pij

j=1

N

∑

i=1

N

∑ (2.2.15)

is proportional to N2. This is not insignificant because collisions in all cells in a 

computational domain must be calculated at each time step. In the null-collision method, the 

task is avoided by regarding Nmax as the collision number in Δt. To compensate for excess 

collisions, Nmax-Nc, some of the Nmax collisions are regarded as null collisions. We now 

define:

Pij = ( ) Pij max qij (2.2.16)

where . This means that the pair (i,j) first makes a tentative 

collision with a probability of (Pij)max. The tentative collision results in a real collision with a 

probability of qij and a null collision with a probability of 1-qij. The number of tentative 

collisions is (Pij)max times the number of pairs N(N-1)/2; the product is equal to Nmax. A 

tentative collision pair can be chosen randomly because (Pij)max is common to all (i,j). So, to 

summarize, we first obtain Nmax. In general, Nmax is not an integer. Let , 

where is the integral part of Nmax. We then set with a probability 

of p and with a probability of 1-p. The following two steps are then repeated 

Nmax times:

1. Choose a pair (i,j) of i ≠ j randomly and obtain qij.

2. Call a random number U. If U < qij, the pair collides. If they do collide, the precollision velocities are replaced by the post-collision velocities given in Equation (2.2.6)

where the particles A and B are of the same kind.

2.2.1.2 Gas Mixture

For the sake of simplicity, consider a binary mixture of species A and B. The types of 

collisions are A-A, B-B, and A-B. The A-B collision is a collision between unlike particles. 

The other two types are treated as in the gas of a simple gas above. Let NA and NB be the 

numbers of simulated particles A and B in the cell. The probability that molecule Ai collides 

with some molecule B in Δt is

PAi = nB gσ T

ABΔt (2.2.17)

where g is the relative speed and is the total cross section for the A-B collision. The 

average is for all simulated molecules B in the cell, so the probability takes the form:

qij = gijσT g( )ij gmaxσT g( ) max

p = Nmax − Nmax ⎢⎣ ⎥⎦

Nmax ⎢⎣ ⎥⎦ Nmax = Nmax ⎢⎣ ⎥⎦ + 1

Nmax = Nmax ⎢⎣ ⎥⎦

σT

AB
15 

PAi = PAi ,Bj

j=1

NB

∑ (2.2.18)

where

PAi ,Bj = NB

−1

nBgAi ,Bjσ T

ABΔt (2.2.19)

If we again introduce the idea of the maximum collision number where gAi,Bi is replaced by 

, we have the upper bound of PAi,Bi:

Pmax

AB = NB−1

nBgmax

AB σ TAB gmaxAB ( )Δt (2.2.20)

The maximum number of A-B collisions in Δt for is

Nmax

AB = PmaxAB = PmaxAB

j=1

NB

∑

i=1

NA

∑

i=1

NA

∑ (2.2.21)

Note that the factor of 1⁄2 is not necessary for A-B collisions. Substitution of Eq. (2.2.20) into 

Eq. (2.2.21) gives us:

Nmax

AB = nBNAgmaxAB σ TAB gmaxAB ( )Δt (2.2.22)

The number density nB is expressed as:

nB = WNB

Vc

(2.2.23)

where Vc is the cell volume and W is the weight. As discussed previously, one simulated 

particle represents W real particles. Presently, assume that all the particles in a cell have a 

common weight. The collision procedure in the previous section is then repeated for all

species combinations A-A, A-B, and B-B using the updated definitions of number of 

collisions and probabilities for like and unlike collision pairs as required.

2.2.2 Internal Degrees-of-Freedom

A typical application of the DSMC method involves the calculation of many billions of 

intermolecular collisions, and this leads to the requirement for computationally simple 

molecular models. This has led to the development of “phenomenological” models that 

focus on the generation of accurate results for the observed gas properties but do not attempt 

to provide the most realistic model of individual collisions at the microscopic level[13]. The 

developers of classical kinetic theory faced a similar problem in that they required 

mathematically tractable molecular models that permitted the derivation of analytical 

expressions for the macroscopic transport properties. The classical models were the obvious 

choice for early DSMC work, but were found to have serious shortcomings that have been 

gmax

AB

Pmax

AB
16 

overcome by phenomenological models that have been introduced in the context of the 

DSMC method.

We have already seen a phenomenological model be applied to elastic collisions by way of 

the variable hard sphere model, which combines the unrealistic but simple isotropic 

scattering of a hard sphere with a realistic variation of diameter with the relative speed of the 

collision pair. This model leads to accurate results for viscosity and heat conduction 

coefficients.

In the same way, classical kinetic theory did not lead to any useful models for inelastic 

collisions. For example, the rough sphere model retains the deficiencies of the hard sphere 

model with regard to the transport properties, and it is unable to deal with the quantum 

effects that cause most gases to have fewer than three rotational degrees of freedom and has a 

fixed and unrealistically fast rotational relaxation time. This problem was solved in the 

context of the DSMC method by the introduction of the Larsen-Borgnakke, or L-B, 

model[42] for rotation which can be regarded as the archetypal phenomenological procedure. 

It is an add-on to the VHS model and, for a fraction of the collisions that are chosen to match 

measured rotational relaxation rates, the post-collision rotational energies are selected from 

the equilibrium distribution that correspond to the collision energy. The L-B model was 

extended to the vibrational modes, but it assumed that the internal energies were 

continuously distributed and was found lacking. The situation was corrected by the 

introduction[43] of the quantum version of the L-B model. 

DSMC methods test colliding particles for possible internal energy exchanges or chemical 

reactions. The collision selection rules use probabilities that must be defined such that they 

result in the desired macroscopic relaxation rate behavior. This typically involves identifying 

a macroscopic rate equation that the simulation attempts to match through appropriate 

definitions of selection probabilities P. These rate equations often contain parameters such 

as collision numbers Z, characteristic times τ, or rate coefficients that are estimated from 

experiments or analysis. Ambiguity regarding specific definitions of the targeted 

macroscopic rates and coefficients has often led to confusion in defining the appropriate 

selection probabilities in DMSC methods.

In studies of vibrational relaxation, Haas[44] developed a relationship between the 

appropriate DSMC selection probability P and the resulting Landau-Teller macroscopic 

relaxation behavior. This relationship was later generalized, allowing for rotational 

relaxation, by Lumpkin et al.[45] in the form of a “correction factor” for collision numbers Z. 

It was observed that the DSMC probabilities P, in general, cannot be specified by the 

assumption P=1/Z, which has been employed throughout the DSMC community. Dissimilar 

applications of the selection probability and incompatible definitions of collision number led 

to disagreement among DSMC researchers regarding these correction factors.

Haas et al.[46] went on to clarify the definitions of macroscopic rate parameters employed in 

rotational relaxation, including differing interpretations of Jeans’ equation[47]. It was shown 

that relaxing energies, when plotted against the cumulative number of molecular collisions, 

lead to a single solution, regardless of the assumed intermolecular potential in an adiabatic 

bath. Differing definitions of collision time τc between experiments and DSMC were 
17 

clarified. The two most widely used selection methodologies were discussed, and DSMC 

rotational relaxation probabilities P, specific to each methodology, were derived to match the 

macroscopic rate equations exactly. Prevalent selection methods cannot reproduce Jeans’ 

equation exactly for the case of heteromolecular collisions, nor can they reproduce Jeans’ and 

Landau-Teller[48] rate behavior for simultaneous rotational and vibrational relaxation. Thus, 

a new selection methodology was introduced that was better suited to the general case of gas 

mixtures experiencing multimode internal energy relaxation. Although constant Z values 

were used, the developments were valid for variable Z as well.

Gimelshein et al.[49, 50] proposed a simplified inelastic collision algorithm using particle 

selection prohibiting multiple relaxations. They felt that the selection methodology proposed 

by Haas[46] lead to an unnecessarily complicated functional dependence between the 

probability and the collision number. To avoid the complexity, a straightforward selection 

methodology was developed which also resulted in relaxation behavior matching the LandauTeller equation when the correct correction factor was used.

Gimelshein et al.[50] took a closer look at vibrational relaxation and derived an exact 

relationship that connects the vibrational relaxation number, Zvib

DSMC, used in the DSMC 

method and that employed in continuum simulations. An approximate expression for 

Zvib

DSMC was also derived that is cost-effective and applicable when translational temperature 

is larger than the vibrational temperature.

2.2.2.1 Rotational and Vibrational Energy Transfer

Models for rotational and vibrational energy transfer with translational energy in a collision 

have been explored previously. The Larsen-Borgnakke (L-B) model[42] is usually used to 

describe the rotational and vibrational energy transfer in the DSMC method. A fraction of 

the total number of collisions, 1/Zvib or 1/Zrot, lead to the energy exchange, where Zvib and Zrot

are the vibrational and rotational collision numbers, respectively. Originally, the internal 

energy modes were treated as continuous. These energy modes are, of course, discrete. To 

address this issue, Bergemann and Boyd introduced the algorithms of the L-B model with 

discrete vibrational and rotational energy levels for diatomic molecules[43, 51]. Gimelshein 

et al. then extended the L-B model for the case of polyatomic gas mixtures with discrete 

vibrational and rotational energies[49].

The temperature dependent nature of rotational and vibrational relaxation rates has been 

developed by Parker[52] and Millikan and White[53], which provides the characteristic 

relaxation times as functions of the temperature for a wide range of pure and mixed diatomic 

gases. In the L-B model, the vibrational and rotational post-collision energies are modeled 

according to the local equilibrium distribution functions. The derivation of these models for 

application in the DSMC method can be viewed in Reference [2], but the vibrational 

relaxation model resulting from Millikan and White, which is used in this study, is:

Zvib = C1

T ω

⎛

⎝

⎜ ⎞

⎠

⎟ exp C2T −1

3

( ) (2.2.24)

where C1 and C2 are constants used to fit the experimental data. 
18 

However, these lead to the incorrect macroscopic relaxation rate. Lumpkin and Haas et 

al.[45, 46] introduced a correction factor for the rotational relaxation rate in the DSMC 

method. Most continuum analyses and computational codes use Jeans’ equation[47] to 

model rotational relaxation. In a single species gas, Jeans’ equation uses the instantaneous 

equilibrium rotational energy Erot*(t) and a time constant, which can be written as the product 

of the instantaneous collision time τc and a characteristic collision number Z as:

dErot

dt = Erot

* ( )t − Erot ( )t

τ cZr

cont (2.2.25)

where Erot is the mean rotational energy per molecule in the ensemble, and t represents time. 

Denoting the rotational degrees of freedom by ζrot, Erot* is defined by:

Erot

* = ζ rot

2

kT (2.2.26)

and will change with time as the translational temperature T varies during relaxation under 

adiabatic conditions. Similarly, τc may vary with T, depending upon the intermolecular 

potential assumed for the gas. Most DSMC methods employ the VHS model, which assumes 

rigid-sphere collision mechanics, but establishes the collision rate in a manner consistent 

with an inverse-power potential. The relative translational energy of colliding particles near 

equilibrium, as a result of biasing due to collision selection, is distributed with ζtrans degrees 

of freedom given by:

ζ trans = 4 − 4

α = 5 − 2ω (2.2.27)

Note that this can differ from ζT = 3 associated with the translational energy of an 

equilibrium ensemble of particles in a reservoir. Biasing from collision selection plays a 

significant role in the relaxation of the energy in a reservoir and becomes apparent in the 

derivation of the relaxation probability. Haas[46] derived the relationship for the particle 

selection method for a single species gas (allowing for multiple relaxation) as:

Ppart

DSMC = ζ T + 2ζ rot

2ζ rot

⎛

⎝

⎜

⎞

⎠

⎟ 1− 1− ζ rot

ζ T

⎛

⎝

⎜

⎞

⎠

⎟

ζ T +ζ rot

ζ T + 2ζ rot

⎛

⎝

⎜

⎞

⎠

⎟

4

Zrot

cont

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ (2.2.28)

that correctly reproduced the relaxation rate given by Jeans’ equation. He goes on to say that 

Ppart = 1/Z could be a fair assumption for adiabatic relaxation, but only if one redefined Z by 

using adiabatic relaxation data that has been cast into the isothermal form of Jeans’ equation 

as suggested by Bird[2]. As such for adiabatic flows, Ppart = (1/ZBird) = (1/Z)(1+ζrot/ζT) can 

be a fair assumption to the exact expression. Better yet, the probabilities were shown to 

approach asymptotic limits described by:
19 

Ppart

DSMC = 1+ ζ rot

ζ T

⎛

⎝

⎜

⎞

⎠

⎟

1

Zrot

cont (2.2.29)

which represents an excellent approximation when Z > 10 for the particle-selection 

methodology.

Gimelshein[50] then analyzed the vibrational relaxation rates in the DSMC method. The 

main difference between rotation-translation and vibration-translation relaxation is that while 

the rotational mode can be describes as having a constant number of degrees-of-freedom at 

all temperatures except cryogenic, the number of vibrational degrees-of-freedom varies with 

the temperature. For a simple harmonic oscillator, the number of vibrational degrees-offreedom ζvib is:

ζ vib = 2Θvib T

exp( ) Θvib T −1 (2.2.30)

where Θvib is the characteristic vibrational temperature of the mode. This means that the 

relationship for rotational relaxation probability given above are not applicable to the 

vibration-translation energy transfer. Gimelshein derives the relationships that allow one to 

match the vibrational relaxation rate in DSMC with that predicted by the Landau-Teller 

equation[48]:

dEvib

dt = Evib

* ( )t − Evib ( )t

τ collZvib

cont (2.2.31)

where t is again time, Evib is the mean vibrational energy per molecule in the ensemble, and 

Evib* is the instantaneous equilibrium vibrational energy, defined as:

Evib

* = ζ vib ( ) Ttrans

2

kTtrans (2.2.32)

He goes on to derive the expression for the DSMC collision number for vibration as:

Pvib

DSMC = ζ vib ( ) Ttrans Ttrans −ζ vib ( ) Tvib Tvib

ζ T ( ) Ttrans − T ′

1

Zvib

cont (2.2.33)

The probability of a vibrationally inelastic collision in the DSMC method is then 

. Although exact, Equation (2.2.33) is costly for direct use in DSMC computations because it 

involves the calculation of the vibrational and translational temperatures, and then the 

subsequent numerical determination of T’ (post collisional temperature). However, the exact 

expression can be significantly simplified when the vibrational and translational temperatures 

are close:

P = 1 / Zvib

cont
20 

Pvib

DSMC ≈ ζ T + A

ζ T

1

Zvib

cont

A = 1

2 ζ vib

2 ( ) T exp( ) Θvib T

(2.2.34)

The above equation is applicable only when the difference between the vibrational and 

translational temperatures is small. However, it was found that if the translational 

temperature is used as T, Equation (2.2.34) reproduced well the correct expression in 

Equation (2.2.33) for all cases where the vibrational temperature was lower than the 

translational one.

2.2.2.2 Electronic Energy Levels

In the past, Bird[15] and Carlson et al.[17, 25] modeled electronically excited levels and 

electronic transitions in DSMC. In their work, each excited particle was assigned a 

distribution over all available states, assuming that the electronic energy mode was in 

equilibrium. Gallis and Harvey[16, 54] modeled non-equilibrium thermal radiation in 

DSMC as well as electronic excitation processes. However, for their method, accurate cross 

sections for electronic excitation and absorption were required and the uncertainty in the 

cross sections is not well known for all the species. Because of these deficiencies, a much 

more accessible non-equilibrium electronic energy model is developed in this thesis in

Section 3.1.

2.2.2.3 Selection Methodology

Haas[46] derives the relationship for a versatile L-B selection methodology that is applicable 

to gas mixtures and is adapted from the particle-selection method; probability, P, is applied 

to each molecule individually within the colliding pair, but the event in which both particles 

relax is prohibited (called particle-selection prohibiting double relaxation). This ensures that 

relaxation of A molecules is independent of the energy of B molecules, and the L-B 

mechanics can therefore reproduce the relaxation rate behavior for each species 

independently. The appropriate functional form for the relaxation probability was defined as:

1

2 PA|B

2 − ζ TA|B +ζ rA

2ζ TA|B

1

ZA|B

−

ζ TA|B +ζ rB

2ζ TA|B

1

ZB|A

+1 ⎛

⎝

⎜

⎞

⎠

⎟ PA|B +

ζ TA|B +ζ rA

ζ TA|B

1

ZA|B

= 0 (2.2.35)

where the Z values are from continuum methods. This expression involves only rotational 

energy exchange. The expression for multimode energy relaxation is presented by Haas[46], 

and the reader is referred to the reference for the complete derivation.

Gimelshein[49, 50] believed that, although correct, the above selection methodology leads to 

an unnecessarily complicated functional dependence between the probability and the 

collision number. This is especially true when multi-mode relaxation is considered. To 

avoid the complexity, a straightforward selection methodology was used that is illustrated in 

Figure 2.1. The first step is to select a random number Rn uniformly between 0 and 1. 

Then, using the value of Rn, the possible relaxation events are tested, so that either rotational 
21 

or vibrational relaxation may occur for particle A or particle B. Multiple relaxation events 

are prohibited.

Figure 2.1: Selection methodology for inelastic and elastic collisions.

Similar to that suggested by Haas, the simplified selection methodology allows only one 

relaxation event to take place in a collision. Vibrational and rotational relaxation rates of 

different species in the gas mixture are independent, and the necessary relaxation rates may 

also be accurately captured. It can be easily implemented and is characterized by a simpler 

dependence between the collision numbers and the relaxation probabilities. This method 

should not be applied when the total probability of internal-translational energy exchange is 

larger than unity. If the sum of the probabilities is greater than unity, an alternative 

methodology is implemented as shown in Figure 2.2. These selection methodologies are 

easily updated to include electronic energy exchange.
22 

Figure 2.2: Selection methodology for inelastic collisions if the sum of the 

probabilities is greater than unity.

2.2.3 Chemical Reactions

Once the collision partners have been selected as described in Section 2.2.1, they are then 

tested to see if a chemical reaction will occur. A typical bimolecular reaction may be written 

as:

A + B ↔ C + D

As long as the reaction takes place in a single step with no species other than the reactants 

present, it has been found that the rate equation for the change of concentration of species A

may be written:

− dnA

dt = kf ( ) T nAnB − kr ( ) T nCnD (2.2.36)

The rate coefficients are functions of temperature and are independent of the number 

densities and time. A further result that originally had an empirical basis is that the rate 

coefficients are of the form:
23 

kf = aT b exp − Ea

kT

⎛

⎝

⎜ ⎞

⎠

⎟ (2.2.37)

where a and b are constants and Ea is called the activation energy of the reaction.

The Total Collision Energy (TCE) model[55] of Bird was developed in order to obtain a rate 

coefficient based on a chemical reaction cross section which takes the form:

σ react =

0 if Ecoll < Ea

σ TC1 ( ) Ecoll − Ea

C2

1− Ea

Ecoll

⎛

⎝

⎜

⎞

⎠

⎟

ζ +3

2−ω AB

if Ecoll > Ea

⎧

⎨

⎪

⎪

⎩

⎪

⎪

(2.2.38)

where C1 and C2 are constants, Ecoll is the collision energy available to the reaction, and is 

the average number of internal degrees of freedom that contribute to the collision energy. 

The resulting rate coefficient becomes:

k T( ) = 2C1

σ ref

π1

2

ε

2kTref

mr

⎛

⎝

⎜

⎞

⎠

⎟

1

2 Γ ζ + C2 + 5 ( 2 −ω AB )

Γ ζ + b + 3 ( )2

kC2T C2+1−ω AB

Tref

1−ω AB

exp − Ea

kT

⎛

⎝

⎜ ⎞

⎠

⎟ (2.2.39)

where the constants C1 and C2 are:

C1 = π1

2

εaTref

b

2σ ref

Γ ζ + 5 ( ) 2 −ω AB

Γ ζ + b + 3 ( )2

mr

2kTref

⎛

⎝

⎜

⎞

⎠

⎟

1

2

Tref

1−ω AB

kb−1+ω AB

C2 = b −1+ω AB

(2.2.40)

The probability of the reaction then occurring at each collision with Ec > Ea between powerlaw VHS particles is then:

Preact = σ R

σ T

= π1

2

εaTref

b

2σ ref ( ) kTref

b−1+ω AB

Γ ζ + 5 ( ) 2 −ω AB

Γ ζ + b + 3 ( )2

mr

2kTref

⎛

⎝

⎜

⎞

⎠

⎟

1

2

( ) Ecoll − Ea

b+ζ +1

2

Ecoll

ζ +3

2−ω AB

(2.2.41)

The same can be done for termolecular, dissociation/recombination reactions represented as:

AB + T ↔ A + B + T

where AB is the dissociating molecule and T is the third body particle. The rate of formation 

of either species A or B in the forward reaction can be written:

dnA

dt = kf nABnT (2.2.42)

ζ
24 

while the loss of species A in the reverse reaction is

dnA

dt = krnAnBnT (2.2.43)

The corresponding reaction probability for the forward dissociation reaction is the same as 

that in Equation (2.2.41) and the probability of a recombination can be written:

Preact = σ react

σ T

= π1

2

nT εaTref

b

2σ ref

Γ 5 ( ) 2 −ω AB

Γ b + 3 ( )2

mr

2kTref

⎛

⎝

⎜

⎞

⎠

⎟

1

2

Ecoll

kTref

⎛

⎝

⎜

⎞

⎠

⎟

b−1+ω AB

(2.2.44)

where the activation energy and, therefore, contributing internal degrees of freedom are taken 

as zero.

These are the typically implemented set of probabilities for both forward and reverse 

reactions in DSMC codes. They require available data for each forward and reverse reaction. 

However, these data are not always readily available for the reverse reactions. Also, many 

CFD codes make use of statistical mechanics in order to specify the reverse reaction rate. If 

comparisons are to be made to CFD results, it follows that the DSMC code should use 

statistical mechanics to obtain the reverse reaction rate where data are not available. This 

methodology is not used in general where data are available because the reaction 

probabilities are functions of the macroscopic temperature, which is less physically realistic 

than the dependence on the collision energy for the forward reactions. Another use of 

statistical mechanics is to calculate the reverse reaction rate from the forward reaction data, 

and then choose the relevant coefficients from the available reverse reaction data to most 

closely match the rate computed from statistical mechanics. The following modifications 

were implemented in the DAC software for those reverse reactions where data were

unavailable.

Chemical equilibrium requires a balance between the forward and reverse reactions such that 

the concentrations of all species remain constant with time. Therefore, for a bimolecular 

reaction,

Keq ≡ kf

kr

= nCnD

nAnB

(2.2.45)

where Keq is the equilibrium constant. Similarly, the equilibrium constant for the 

termolecular dissociation-recombination reaction is written as:

Keq ≡ kf

kr

= nAnB

nAB

(2.2.46)

Equilibrium statistical mechanics relates the numbers of the various species in a system to the 

partition functions. For the bimolecular reaction,
25 

NC ND

NANB

= QC

QD

QAQB exp − EakT

⎛

⎝

⎜ ⎞

⎠

⎟ (2.2.47)

where Q is the partition function for the species that is indicated by the superscript, and Ea

should be interpreted as the difference between the forward and reverse activation energies. 

The number N is related to the volume V of the system by N = nV, which results in the 

bimolecular equilibrium constant written as:

Keq ≡ kf

kr

= nCnD

nAnB

= QC

QD

QAQB exp − EakT

⎛

⎝

⎜ ⎞

⎠

⎟ (2.2.48)

For the termolecular reaction, the system volume does not cancel and the corresponding 

equation is:

Keq ≡ kf

kr

= nAnB

nAB

= QA

QB

VQAB exp − EakT

⎛

⎝

⎜ ⎞

⎠

⎟ (2.2.49)

The partition function for each species can be written as the product of the separate 

translational, rotational, vibrational, and electronic partition functions:

Q = QtransQrotQvibQel (2.2.50)

First, the translational partition function is:

Qtrans = V 2πmkT h2 ( )3

2

(2.2.51)

The rotational partition function is:

Qrot = T

ε′Θrot

(2.2.52)

where ε’ is a symmetry factor equal to 2 for homonuclear molecules and 1 for heteronuclear 

molecules. The harmonic oscillator model leads to a vibrational partition function of:

Qvib = 1 1( − exp( ) −Θvib T ) (2.2.53)

Finally, the electronic partition function is summed over the j levels as:

Qel = gi exp − Eel

i

kT

⎛

⎝

⎜

⎞

⎠

⎟

i=0

j

∑ (2.2.54)

where gi are the degeneracies. Finally, it can be shown[2] that, for a bimolecular reverse 

reaction given the data for the forward reaction (a’s and b’s), the reverse reaction probability 

is:
26 

Preact = σ react

σ T

= afT bf

( ) 2 ε π−1

2

σ ref ( ) T Tref

1−ω

( ) 2kTref mr

1

2

QAQB

QCQD (2.2.55)

Also, for the termolecular reverse reaction:

Preact = σ react

σ T

= afT bf

( ) 2 ε π−1

2

σ ref ( ) T Tref

1−ω

( ) 2kTref mr

1

2

nTVQAB

QAQB (2.2.56)

For the case where the reverse reaction has finite activation energy, the probability is then:

Preact = σ react

σ T

= afT bf exp −Ea ( ) kT

( ) 2 ε π−1

2

σ ref ( ) T Tref

1−ω

( ) 2kTref mr

1

2

QAQB

QCQD

Γ ζ + 5 ( ) 2 −ω

Γ ζ + 5

2 −ω,Ea ( kT ) (2.2.57)

2.2.4 Electric Field

Various researchers have modeled the electric field structure in a hypersonic shock layer in 

the past. The first such calculation was by Bird[26]. In his work, the electrons were moved 

through the computational grid with the ions that they were created with. No explicit electric 

field calculation was made, so the ions and electrons were not accelerated, and the simulation 

time step was that of the heavy particles. Later, Bird modified this method[27] to include a 

calculation of the ambipolar electric field using a form of the equation derived by Langmuir 

and Tonks[56]

Eambipolar = − kTe

e

d ln n( )e ⎡

⎣ ⎤⎦

dx (2.2.58)

In this equation, Te and ne are the average electron translational temperature and number 

density, macroscopic quantities derived from the DSMC solution. This equation can be 

derived from the macroscopic conservation of momentum for the electron species using the 

assumption of negligible inertial effects, negligible friction due to collisions, zero magnetic 

field and constant translational temperature. In his work, Bird calculated Eambipolar using the 

results from a previously converged calculation. He then re-converged the DSMC 

calculation, applying the prescribed acceleration to the charged particles at each time step 

from Eambipolar. The field Eambipolar was calculated from the re-converged DMSC solution, and 

the entire process repeated until there was no change in the resulting flow field parameters. 

The movement of the electrons was still tied to that of the ions in this method.

In the work by Gallis and Harvey[28], the electrons and ions were again moved together, and 

the electric field was calculated in a similar manner to equation (2.2.58). Here, however, the 

assumption of isothermal electrons was not made:

Eambipolar = − kTe

e

d ln p( )e ⎡

⎣ ⎤⎦

dx (2.2.59)
27 

and the gradient of the electron pressure was instead used to calculate the ambipolar electric 

field. Gallis and Harvey used a simulation time step corresponding to the cell crossing time 

of the heavy particles, and accelerated the heavy particles according to the field Eambipolar. 

They then computed the average ion velocity in a given cell, and adjusted the average 

electron velocity in the same cell to match. In this way, the electron energies were affected 

by the electric field in an average sense.

Carlson, Hassan and Taylor[25, 29] devised a method in which assumptions of zero net 

current and charge neutrality were used to calculate the electric field without the use of any 

macroscopic quantities. They wrote down Newton’s law for a charged particle of species s

moving under the influence of an electric field:

ms

d



v

dt = qs



E (2.2.60)

and thus the expressions for average electron and ion velocity by summing over individual 

particles:



ve =



v ∑ e,0

Ne

− e



EΔt

2me



vi =



v ∑ i,0

Ni

+

e



EΔt

2mi

(2.2.61)

where is the average electric field over a simulation time step, N is the number of 

simulator particles, is a velocity vector at the start of a simulation iteration, and the 

subscripts e and i refer to electron and ion species, respectively. By enforcing charge 

neutrality, Ne = Ni, they then solved for the ambipolar electric field required to make the 

average ion and electron velocities equal, while enforcing the assumption of zero net current:

e



EΔt

2 =



ve,0 − 

v ∑∑ i,0

Ne

1

mi

− 1

me

∑∑⎛

⎝

⎜

⎞

⎠

⎟

(2.2.62)

The electric field given above was used to accelerate the charged particles during the move 

portion of the DSMC algorithm. Additionally, a simple model of the particle acceleration in 

the plasma sheath at the vehicle surface was implemented. The electrons were reflected 

specularly from the vehicle surface to simulate their reflection in the strong negative 

potential gradient that would exist in the plasma sheath. An additional energy increment, eφ, 

was added to the surface heat transfer measurement for each ion impacting the surface. The 

expected potential drop across the sheath, φ, was calculated from one-dimensional 

collisionless sheath theory.

Boyd[30] developed a model for the electric field of a weakly ionized plasma that is similar 

to the one originally proposed by Bird. The model was first implemented in a simulation of 



E 

v0
28 

the plasma flow through an arcjet type thruster. Each electron is moved with the velocity of 

an ion throughout the domain, and the charged particles do not receive any velocity 

increment because of their response to an electric field. However, in this model, the average 

ion velocity of the cell in which the electron is located is used to move the electron particles, 

and it is not associated with a specific ion particle. This model maintains charge neutrality in 

an approximate sense, and results in lower computational overhead than the first model 

proposed by Bird.

In order to calculate the self-induced electric field without any assumption of charge 

neutrality, zero net current, or ambipolar diffusion, the electrostatic Poisson equation must be 

solved. This was carried out for very rarefied flows surrounding spacecraft in low Earth orbit 

by two groups of researchers[57, 58]. While interesting, the preliminary results presented by 

these researchers in the cited references are for flow fields of significantly different structure 

than the shock layers being considered in this work. Gallis et al.[59] presented preliminary 

results in which they solved the electrostatic Poisson equation in conjunction with the DSMC 

method for a hypersonic helium flow field. Unfortunately, the results presented in this work 

were sparse, and the physical processes associated with reacting air were not included in the 

analysis, and a systematic study of the effect of rigorously including the electric field in the 

computation on the flow field parameters was not presented.

Presently, a hybrid of Bird’s original model and Boyd’s model[30] was implemented. The 

electrons were kept separate but retained their own velocities as was done in Boyd’s model; 

however, they are moved along with their associated ion as was done in Bird’s model. By 

doing this, all particles are moved at the heavy particle time scale, but collisions occur at the 

electron time scale.

2.3 Recent Direct Simulation Monte Carlo Advances

Molecular-level chemistry models that predict equilibrium and non-equilibrium reaction rates 

using only kinetic theory and fundamental molecular properties (i.e., no macroscopic 

reaction-rate information) have recently been proposed[10, 11, 13, 14]. In addition, a new 

vibrational relaxation model has also been introduced that only relies on fundamental 

molecular properties and a single measure of vibrational collision number at a reference 

temperature[10]. Preliminary evaluations[14] indicate that the resultant reaction rates are in 

very good agreement with measured Arrhenius rates for near-equilibrium conditions and with 

both measured rates and other theoretical models for far-from-equilibrium conditions.

2.3.1 Vibrational Relaxation

A new vibrational relaxation model has recently been introduced by Bird[10]. The data for 

the vibrational collision number, Zvib, is generally presented as a function of temperature in a 

Millikan-White plot. A convenient way of using this information in a DSMC procedure is 

through an expression that employs the values of Zvib at two temperatures. The temperatures 

may be set to the characteristic temperature of vibration, Θvib, and the characteristic 

temperature of dissociation, Θdiss. Then, if the corresponding collision numbers are Zc and 

Zdiss, the variation of the collision number with temperature T is given by:
29 

Zvib = Zdiss

Θdiss

T

⎛

⎝

⎜ ⎞

⎠

⎟

ω Zc

Zdiss

Θvib

Θdiss

⎛

⎝

⎜

⎞

⎠

⎟

ω ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

Θdiss

T

⎛

⎝

⎜ ⎞

⎠

⎟

1

3

−1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

Θdiss

Θvib

⎛

⎝

⎜

⎞

⎠

⎟

1

3

−1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

(2.3.1)

The problem with the Millikan-White data for most gases is that it leads to collision numbers 

less than unity at the dissociation temperature. This is physically impossible and the most 

effective way to deal with the problem is to set the collision number to unity at the 

dissociation temperature. Also, it is preferable to set the second temperature, Tref, where the 

collision number is Zref. The collision number then becomes:

Zvib = Θdiss

T

⎛

⎝

⎜ ⎞

⎠

⎟

ω

Zref

Θvib

Tref

⎛

⎝

⎜

⎞

⎠

⎟

ω ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

Θdiss

T

⎛

⎝

⎜ ⎞

⎠

⎟

1

3

−1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

Θdiss

Tref

⎛

⎝

⎜

⎞

⎠

⎟

1

3

−1

⎡

⎣

⎢

⎢

⎢

⎤

⎦

⎥

⎥

⎥

(2.3.2)

For a procedure based on this equation, the required data for each vibrational mode is the 

vibrational collision number at the reference temperature. The one-third-power law is 

preserved, but there is a decrease in the rate at higher temperatures. This is essentially a 

more precise alternative to the “Park correction factor”[60] that is commonly employed in 

CFD.

Consider a collision in which the sum of the relative translational energy and the precollision vibrational energy of the molecule under consideration is Ecoll. The “quantized 

collision temperature” is defined by:

Tcoll = imaxΘvib

7

2 −ω

imax = Int

Ecoll

kΘvib

⎡

⎣

⎢ ⎤

⎦

⎥

(2.3.3)

where the collision energy has been quantized through the truncation of imax to an integer. 

The relevant parameters for the new vibrational model used are listed in Table A.3.

2.3.2 Chemical Reactions – The Quantum-Kinetic Chemistry Model

Bird has proposed[10, 11, 13] a set of molecular-level chemistry models based solely on 

fundamental properties of the two colliding particles: total collision energy, quantized 

vibrational levels, and molecular dissociation energies. These models link chemical-reaction 

cross-sections to the energy-exchange process and the probability of transition between 

vibrational energy states. The Larsen-Borgnakke[42] procedures and the principle of 

microscopic reversibility are then used to derive simple models for recombination reactions 

and for the reverse (exothermic) reactions. These models do not require any macroscopic 

data, and they balance the fluxes into and out of each state, thus satisfying microscopic 

reversibility. As mentioned earlier, preliminary evaluations[14] indicate that the resultant 

reaction rates are in very good agreement with measured Arrhenius rates for near-equilibrium 
30 

conditions and with both measured rates and other theoretical models for far-from 

equilibrium conditions.

2.3.2.1 Dissociation Reactions

Bird suggests[10] that a complete and consistent description of vibrational energy exchange 

between colliding molecules must include dissociation. This concept is embodied in the 

energy-threshold model of Bird[2]. According to this model, during the energy-exchange 

process of a diatomic molecule, if the energy content of the vibrational mode exceeds the 

energy threshold of the dissociation or exchange reaction, the reaction

occurs.

In his critique of Bird’s model, Lord[61] pointed out that a “continuum” of states (i.e., an 

infinite number) exists above the dissociation limit. Therefore, Lord argues, if a dissociation 

reaction is energetically possible, it occurs. Thus, assuming a particular collision between 

two simulators, where at least one of them is a diatomic molecule, the serial application of 

the Larsen-Borgnakke energy exchange model would make the collision energy:

Ecoll = Etrans,pair + Evib,mol (2.3.4)

available to the vibrational mode of the molecule in question. In this expression of collision 

energy, Etrans,pair is the translational energy of the pair, and Evib,mol is the vibrational energy of 

the molecule considered for dissociation. Assuming a harmonic oscillator model, the 

maximum vibrational level that can be reached is that given in Equation (2.3.3), where the 

vibrational collision number has been defined as unity. If this level is higher than the 

dissociation level (id = Θd / Θv), or

imax >

Θdiss

Θvib

(2.3.5)

a dissociation occurs.

When the serial application of the quantum Larsen-Borgnakke model is considered, potential 

post-collision states i

* are chosen uniformly from the states equal to or below imax and are 

selected through an acceptance-rejection routine[2] with probability

Pvib

i

*

= 1− i

*

kΘvib

Ecoll

⎛

⎝

⎜

⎞

⎠

⎟

3

2−ω

(2.3.6)

where ω is the temperature exponent of the coefficient of viscosity. Then, the corresponding 

rate coefficient in an equilibrium VHS gas at temperature T is:

kf

diss ( ) T = gRcollAB,TΥ i( ) max

AB,T

(2.3.7)

The degeneracy g of the reaction is included to cover the case of a polyatomic molecule in 

which the reaction can occur through more than one mode with the same characteristic

AB + T → A + B + T
31 

vibrational temperature. The collision rate parameter is the collision rate for collisions 

between gas species AB and T divided by the number density product. It can be written for 

an equilibrium VHS gas as

Rcoll

AB,T = 2π1

2

ε

⎛

⎝

⎜

⎞

⎠

⎟ rref

AB + rrefT

( )2 T

Tref

⎛

⎝

⎜

⎞

⎠

⎟

1−ωAB,T

2kTref

mr

AB,T

⎛

⎝

⎜

⎞

⎠

⎟

1

2

(2.3.8)

where rref is the molecular radius at temperature Tref and mr is the reduced mass. For 

dissociations, it is necessary to add the symmetry factor ε that is 1 for homonuclear and 2 for 

heteronuclear molecules. The parameter is the fraction of collisions between AB 

and T that have sufficient energy to meet the condition of equation (2.3.5). With i as the precollision vibrational state of AB, the result for an equilibrium VHS gas is

Υ i( ) max

AB,T

= Q 5

2 −ωAB,T ( ), i( ) max − i Θvib

AB

T

⎡

⎣

⎢ ⎤

⎦

⎥exp − iΘvib

AB

T

⎛

⎝

⎜

⎞

⎠

⎟

⎧

⎨

⎪

⎩⎪

⎫

⎬

⎪

i=0 ⎭⎪

imax−1

∑ zvib ( ) T AB (2.3.9)

Here, is a form of the incomplete Gamma function and 

is the vibrational partition function. As an example, the 

dissociation rate of molecular nitrogen is presented in Figure 2.3 for the Q-K model as 

compared to several rates from the literature, the first of which is the value used in the 

DSMC simulations using the TCE chemistry model.

Figure 2.3: Comparison of Q-K dissociation rate of the reaction N2 + N2 → N + 

N + N2 with those from literature.

Rcoll

AB,T

Υ i( ) max

AB,T

Q a( ) , x = Γ( ) a, x Γ( ) a

zv ( ) T = 1 1⎡ − exp( ) −Θv T ⎣ ⎤

⎦
32 

2.3.2.2 Recombination Reactions

The condition for recombination in a collision between species A and B is that another atom 

or molecule is within the “ternary collision volume” Vcoll. The probability of recombination 

is therefore

Prec = nVcoll (2.3.10)

where n is the number density. The corresponding rate equation is

kr ( ) T = Rcoll

A,B

nVcoll (2.3.11)

The ratio from Equations (2.3.7) and (2.3.11) may be equated to the 

equilibrium constant from statistical mechanics in order to determine the ternary collision 

volume. It has been found that a volume in the form

Vcoll = aT bVref (2.3.12)

allows nearly an exact fit of this ratio to the equilibrium constant. The coefficients used for 

these adjustments can be viewed in Table A.16. The reference volume Vref is most 

conveniently set to the volume of the sphere with radius equal to the sum of the reference 

radii of the three particles.

The values of the constants a and b are determined in the context of the macroscopic 

temperature, but the evaluation of the ternary collision volume in the DSMC implementation 

is based on the quantities associated with the collision. The “collision temperature” based on 

the relative translational energy in a collision with relative speed g is:

Tcoll = mr

A,B

g

2

2k 5

2 −ωA,B ( ) (2.3.13)

As an example, the cumulative recombination rate of two atomic nitrogen atoms is presented 

in Figure 2.4 as compared to values from the literature, the first of which is the rate 

computed by statistical mechanics from the dissociation of molecular nitrogen.

kf ( ) T / kr ( ) T
33 

Figure 2.4: Comparison of Q-K recombination rate of atomic nitrogen with 

those from literature.

2.3.2.3 Exchange Reactions

Exchange reactions are binary reactions written as in which A and C are the 

molecules that split, while B and D are the atoms. The phenomenological DSMC procedure 

for this reaction is analogous to that for the dissociation model. For both the forward and 

reverse reactions, the reaction probability is equal to that of Equation (2.3.6) with i

* set to the 

vibrational level corresponding to the activation energy of the reaction, ia

ia = Int

Ea

kΘvib

⎡

⎣

⎢ ⎤

⎦

⎥ (2.3.14)

This is a relative probability that must be normalized through division by the sum of 

probabilities from the ground state to imax, resulting in

Pexch = 1− ia ( kΘvib Ecoll)

3

2−ω

(1− ikΘvib Ecoll)

3

2−ω

i=0

imax

∑

(2.3.15)

This is sufficient if the level ia is much larger than unity, but some endothermic reactions and 

all exothermic reactions have either zero or small activation energies. The use of discrete 

levels then leads to unacceptably abrupt changes in the rate coefficients and it is impossible 

to obtain a near exact agreement between the ratio of forward to reverse rate coefficients and 

the equilibrium constant. The level ia is the highest level with energy below the activation 

energy Ea and a smooth curve is obtained for the rate coefficients if the energy iakΘv in the 

numerator is replaced by Ea. The DSMC probability of an exchange reaction in a collision 

with Ec > Ea is then

10-48

10-47

10-46

10-45

10-44

0 5000 10000 15000 20000 25000 30000 35000

Q-K Analytic

Q-K Sampled

Park[65]*

Farbar[6]

Gupta[78]

Rate Constant (m3/molecule/s)

Temperature (K)

N + N + T > N2 + T

A + B ↔ C + D
34 

Pexch = ( ) 1− Ea Ecoll

3

2−ω

(1− ikΘvib Ecoll)

3

2−ω

i=0

imax

∑

(2.3.16)

An exact expression can be written down for the fraction of collisions that satisfy the 

condition of Equation (2.3.16) in an equilibrium VHS gas. With j as the pre-collision 

vibrational level, this is

dN

N = 1− exp( ) −Θvib T

Γ 5 ( ) 2 −ω

× exp − jΘvib

T

⎛

⎝

⎜ ⎞

⎠

⎟ mrg

2

2kT

1− Ea

1

2 mrg

2 + jkΘvib

⎛

⎝

⎜

⎞

⎠

⎟

⎡

⎣

⎢ ⎤

⎦

⎥

3

2−ω

Ea ( ) k− jΘvib T if >0,else0

∞

∫

j=0

∞

∑

× 1− ikΘvib

1

2 mrg

2 + jkΘvib

⎛

⎝

⎜

⎞

⎠

⎟

3

2−ω

i=0

1

2mrg2 ( + jkΘvib ) ( ) kΘvib ⎡

⎣ ⎤⎦

∑

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

−1

exp − mrg

2

2kT

⎛

⎝

⎜

⎞

⎠

⎟ d − mrg

2

2kT

⎛

⎝

⎜

⎞

⎠

⎟

(2.3.17)

and the corresponding rate coefficient is

k T( ) = gRcoll( ) dN N (2.3.18)

This is a complex expression not easily programmed, but it is possible to develop an 

alternative expression that is based on physical reasoning rather than formal kinetic theory. 

The selection of possible post-collision vibrational levels in Equation (2.3.6) is from an 

equilibrium distribution so that the rate coefficient is simply the product of the degeneracy g, 

the collision rate parameter, and the fraction of the vibrational level ia in the Boltzmann 

distribution. The term iakΘv in the exponential is replaced by Ea in order to maintain 

consistency with Equation (2.3.16), and the normalizing factor is simply the vibrational 

partition function, thus

kf ( ) T = g

A

Rcoll

A,B exp − Ea

A,B

kT

⎛

⎝

⎜

⎞

⎠

⎟ zvib ( ) T A (2.3.19)

and

kr ( ) T = gC Rcoll

C,D exp − Ea

C,D

kT

⎛

⎝

⎜

⎞

⎠

⎟ zvib ( ) T A (2.3.20)

The default activation energy for the forward reaction is the heat of reaction and that for the 

reverse reaction is zero. The ratio of the forward to the reverse rate is compared with the 

equilibrium constant and the two are brought into the best possible agreement by the upward 

adjustment of either the forward or reverse activation energy. These adjustments are 
35 

temperature dependent through a relation similar to that in Equation (2.3.12) for the ternary 

collision factor:

Ea = a

Tcoll

Tref

⎛

⎝

⎜

⎞

⎠

⎟

b

Eexch

Ea = Ea − Eexch if Eexch < 0

(2.3.21)

For the endothermic reaction, Eexch is the negative of the activation energy, which would then 

result in the second part of the equation adding the prescribed energy to the activation 

energy. For the exothermic reaction, Eexch is equal to the activation energy of the 

endothermic reaction, so the second part of the equation above does not come into effect. 

The coefficients used in these adjustments are presented in Table A.17. In addition to these 

phenomenological adjustments, there may be experimental or theoretical evidence for 

activation energies above the default values. Should a phenomenological adjustment fall 

short of the desired value, both the forward and reverse activation energies may be adjusted 

to keep the forward to reverse rate ratio in agreement with the equilibrium constant. As an 

example, the endothermic and exothermic reaction rates for the exchange reaction N2 + O →

NO + N are compared to rates from the literature in Figure 2.5 and Figure 2.6, 

respectively. The Q-K results presented in these figures are in good agreement with those 

presented from the literature.

Figure 2.5: Comparison of Q-K exchange rate for the endothermic reaction

N2 + O → NO + N with those from literature.

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Bird[15]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

N2 + O > NO + N
36 

Figure 2.6: Comparison of Q-K exchange rate for the exothermic reaction 

NO + N → N2 + O with those from literature.

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]*

Bird[15]

Gimelshein[67]

Boyd[69]

Rate Constant (m

Carlson[80]

3/molecule/s)

Temperature (K)

NO + N > N2 + O
37 

Section 3 

Extending the Quantum Kinetic Chemistry Model

For missions where the shock-layer temperature reaches tens-of-thousands of degrees Kelvin, 

the importance of including electronic energy levels and ionization reactions becomes more 

pronounced. Although the treatment of electronic energy levels and ionization using the 

DSMC method is not a new development, most of these methods are based on measured, 

equilibrium rates, which are always questionable at high temperature, especially when 

applied to non-equilibrium problems.

In Section 2, the Quantum-Kinetic Chemistry Model of Bird was introduced. One of the 

attractive properties of this model is the fact that it does not rely on any macroscopic 

information, leveraging instead off of kinetic theory and physical data corresponding to the 

molecules undergoing each individual collision. The purpose of this Section is to introduce 

new methods of treating electronic energy level transitions and chemical reactions involving 

charged species following the Q-K methodologies utilizing not only the vibrational energy 

levels, but also the electronic energy levels of the atoms and molecules.

3.1 Detailed Electronic Energy Level Model

To model the distribution of electronic energy in the DSMC technique, three procedures must 

be defined. First, when a particle is introduced into a simulation, it is necessary to obtain a 

new electronic energy from a given distribution through statistical sampling. A new 

electronic energy is also required when a particle changes its electronic energy after a 

collision with a surface. Under equilibrium conditions, the distribution has the well-known 

Boltzmann form. Second, it is necessary to statistically sample a new electronic energy 

following a collision, which involves electronic energy transfer. Third, a method to 

reproduce the electronic energy transition rate must be defined. Detailed examples of 

applying this model at hypersonic (re)entry conditions in air are presented subsequently in 

Section 4.3.

3.1.1 Equilibrium Sampling

Each electronic energy level j has a distinct energy, Ej, and degeneracy, gj. The Boltzmann 

distribution for the electronic energy levels at a given temperature T gives the following 

result for the fraction of particles in level j:

fel( )j = Nj

N = gj exp −Ej ( ) kT

gi exp −Ei ( ) kT

i=0

imax

∑

(2.3.22)

This distribution is used when either creating a new particle or after a surface collision in the 

DSMC simulation at a boundary specified at the temperature T. However, it is not possible 

to sample an electronic energy level j directly from the distribution. Therefore, an 
38 

acceptance-rejection procedure is employed by selecting values for j from the following 

distribution:

f ′

el( )j = gj exp −Ej ( ) kT

gJ exp −EJ ( ) kT (2.3.23)

where J is the value of j for which Equation (2.3.22) is a maximum. Unfortunately, since the 

degeneracy of each level is a variable specific to each species, there is no way to determine 

the maximum level a priori, so the maximum level must be searched for at each 

implementation or saved and recalled for a constant boundary temperature. The sampling of 

a new electronic energy level then proceeds as follows:

1) select at random an electronic energy level evenly distributed between 0 and Jmax

where Jmax is the maximum possible energy level;

2) determine the value J for which Equation (2.3.22) is a maximum;

3) accept the value of j if fel’(j) > RAND, a random number;

4) if the value of j is not accepted, return to step 1.

3.1.2 Post-Collision Sampling

A phenomenological approach is usually adopted in the DSMC method when a collision 

occurs that involves energy transfer. The Larsen-Borgnakke method[42] samples a postcollision state from a combined distribution of the translational and electronic collision 

energies of the colliding particles. Based on the approach of Bergemann and Boyd[43] and 

Boyd[51], the Dirac delta function is employed to define the distribution of energies in 

Equation (2.3.22) in the following continuous form:

fel Eel ( ) , j = gj exp −Ej ( ) kT

gi exp −Ei ( ) kT

i=0

imax

∑

δ Eel

kT − Eel

j

kT

⎛

⎝

⎜

⎞

⎠

⎟ (2.3.24)

Consideration must now be given to the distribution of translational energy of the colliding 

particles. This distribution is naturally affected by the intermolecular model used since this 

determines the collision probability. For the present study, the variable hard sphere (VHS) 

collision model of Bird[41] is used. However, it is straightforward to develop the 

formulation for an alternative collision model such as the variable soft sphere of Koura[62]. 

For the VHS model, the distribution of the translational energies is:

ftrans ( ) Etrans = 1

Γ 5 ( ) 2 −ω

Etrans

kT

⎛

⎝

⎜ ⎞

⎠

⎟

3

2−ω

exp − Etrans

kT

⎛

⎝

⎜ ⎞

⎠

⎟ (2.3.25)
39 

Using Equations (2.3.24) and (2.3.25), the combined distribution for sampling a postcollision electronic energy level j’ from the total collision energy Ecoll = Etrans+Eel = 

Etrans’+Eel’ is:

f ( ) j′;Ecoll = 1

Γ 5 ( ) 2 −ω

gj′

gi exp −Eel

i ( ) kT

i=0

imax

∑

Ecoll − Eel

j′

kT

⎛

⎝

⎜

⎞

⎠

⎟

3

2−ω

exp − Ecoll

kT

⎛

⎝

⎜ ⎞

⎠

⎟ (2.3.26)

In applying the general Larsen-Borgnakke scheme, it is assumed that local thermodynamic

equilibrium prevails. Therefore, the temperature T in Equation (2.3.26) is constant. Also, the 

total collision energy is constant, so it is only necessary to perform sampling of the postcollision state from the following distribution form:

g( ) j′;Ecoll ∝ gj′ Ecoll − Eel

j′ ( )3

2−ω

(2.3.27)

Again, an acceptance-rejection procedure is used. The normalized distribution that is 

required is obtained by finding the value of j’ for which Equation (2.3.27) is maximum, J’, 

which is different for each value of Ecoll. Therefore, the following distribution is obtained:

g′( ) j′;Ecoll = gj′ Ecoll − Eel

j′ ( )3

2−ω

gJ′ Ecoll − Eel

J′ ( )3

2−ω (2.3.28)

The procedures are implemented as follows:

1. given a pair of particles with total collision energy Ecoll which undergoes electronic 

energy exchange, determine J’;

2. determine the maximum allowable electronic energy level obtainable from Ecoll, J’’;

3. take J

* to be the smaller of J’ and J’’;

4. as described in the previous section on equilibrium sampling, perform an acceptancerejection procedure to sample j’ from:

g′( ) j′;Ecoll = gj′ Ecoll − Eel

j′ ( )3

2−ω

gJ* Ecoll − Eel

J*

( )3

2−ω (2.3.29)

3.1.3 Electronic Energy Level Transitions

A transition from electronic energy level i to level j for species A can be written in the form 

of a chemical reaction as . Following the Q-K model, the simplest model is 

to assume that the transition occurs if the electronic energy level of A after a trial LarsenBorgnakke redistribution of the collision energy, Ecoll = Etrans+Eel

A

, is j

* = j. Since the 

Ai + B → Aj + B
40 

transition is treated as a chemical reaction, there are Jmax transitions to consider where Jmax is 

the maximum energy level energetically possible.

In the DSMC code, the process for determining if an energy level transition occurs, and if so, 

to which transition level, is as follows:

1. determine the value of Jmax and compute the denominator of Equation (2.3.28);

2. perform acceptance-rejection according to Equation (2.3.28) until a test level is 

accepted.

This selection process is analogous to the selection process used in selecting a chemical 

reaction from a list of possible reactions.

3.2 Chemical Reactions Involving Charged Species

As outlined in Section 2.3.2, Bird[10, 11, 13], supported by Gallis[14], proposes a set of 

molecular-level chemistry models based solely on fundamental properties of the two 

colliding molecules: total collision energy, quantized vibrational energy levels, and 

molecular dissociation energies. These models link chemical cross-sections to the energy 

exchange process and the probability of transition between vibrational energy states. The 

Larsen-Borgnakke procedures and the principle of microscopic reversibility are then used to 

derive simple models. The models do not require any macroscopic data and they function by 

seeking to balance fluxes into and out of each state, thus satisfying microscopic reversibility. 

Since reactions involving only neutral species have been covered thus far, processes 

involving charged species are presented in the current section. Again, more detailed 

comparisons for these reactions pertaining to hypersonic (re)entry in are be presented in 

Sections 4.4 through 4.7.

3.2.1 Ionization Reactions

The procedure for ionization is analogous to that for dissociation. The difference is that 

ionization is linked to the electronic energy level instead of the vibrational energy level as

with dissociation. For ionization, the collision energy is defined as Ecoll=Etrans,pair+Eele,particle

where the first term on the right is the relative translational energy of the colliding pair and 

the second term is the electronic energy of the particle being considered for ionization. If the 

collision energy is greater than the ionization energy of the particle, ionization occurs. The 

corresponding rate coefficient in an equilibrium VHS gas at temperature T is:

kf

ion ( ) T = RcollA,B

Υ J( ) max

A,B

(2.4.1)

where the collision rate parameter is identical to that of Equation (2.3.8). The fraction 

of collisions between A and B that have sufficient energy to meet the requirement of 

ionization is then:

Rcoll

A,B
41 

Υ J( ) max

A,B

= Q 5

2 −ωA,B

,

Eion − Eel

i,A

kT

⎛

⎝

⎜

⎞

⎠

⎟

i=0

Jmax

∑ gi

A exp −Eeli,A ( ) kT

zel( ) T A (2.4.2)

where zel(T)A is the electronic partition function for particle A.

3.2.2 Associative Ionization Reactions

The forward (endothermic) associative ionization reaction can be thought of occurring as 

where AB* is an intermediate, excited molecule that may lead to 

the ionized species AB+ plus an electron. The first step is similar to the recombination 

procedure described in Ref. [10]. If the AB* molecule is in the ground vibrational state after 

a trial Larsen-Borgnakke redistribution of the collision energy, the intermediate molecule is 

“formed.” In the second step, the dissociation energy is first added to the original collision 

energy. The new collision energy is then distributed among the internal energy modes of the 

AB* molecule and then this molecule is tested for ionization as described above. If the 

remaining collision energy after redistribution plus the electronic energy of the molecule is 

greater than the ionization energy of the molecule, the associative ionization reaction occurs. 

The corresponding rate coefficient can then be written as:

kf

AI = RcollA,B Υ f ,AIAB*

zvib ( ) T AB*

(2.4.3)

where the collision rate parameter is again identical to that of Equation (2.3.8). The 

vibrational partition function is included to account for the probability that the AB* molecule 

will be in the ground state. The fraction of collisions with sufficient energy for the molecule 

AB* to ionize are then calculated by summing over all possible initial internal energy states of 

atoms A and B (electronic energy states) and all possible vibrational, rotational, and 

electronic energy states of the intermediate AB* molecule as follows:

Υ f ,AI

AB*

= fel iel

A ( ) fel ielB ( )

iel

AB*

=0

iel,max

AB*

∑

irot

AB*

=0

irot ,max AB*

∑

ivib

AB*

=0

ivib,max AB*

∑

iel

B =0

iel,max

B

∑

iel

A =0

iel,max

A

∑

×Pvib,trans

ivib

AB*

Ecoll

0 ( )Prot,transirot

AB*

Ecoll

1 ( )Pel,transiel

AB*

Ecoll

2 ( )Q 5

2 −ωAB*

,

Eion − Eel

AB*

kT

⎛

⎝

⎜

⎞

⎠

⎟

(2.4.4)

The sums over the atomic electronic energy levels is over all possible levels. For the 

molecular internal energies, the maximum level is that which is the highest energetically 

possible. To start, an initial translational energy is assumed to be the average relative 

translational energy, or . This is an approximation and may result in a 

discrepancy between the above analytical expression and the actual sampled rate. It is 

possible to integrate over all possible relative translational energies, but it makes an already 

long calculation (summing over all possible internal states) much longer. As will be shown 

subsequently, this formulation results in an acceptable comparison with the sampled values 

in Section 4.4. The first two expressions after the summations are simply the probabilities of 

atoms A and B being in electronic levels iele

A and ieleB

, respectively. The next three terms are 

A + B → AB* → AB+ + e−

Rcoll

A,B

5 ( ) 2 −ω kT
42 

analogous to Equation (2.3.29) in that they are the probabilities of the respective internal 

mode of energy being at the specified quantum level at the specified collision energy. The 

collision energies are defined as:

Ecoll

0 = 5 ( ) 2 −ω kT + EdissAB*

+ Eel

A + EelB

Ecoll

1 = Ecoll0 − EvibAB*

Ecoll

2 = Ecoll1 − ErotAB*

(2.4.5)

where each successive collision energy subtracts the previous internal energy mode. Then 

the last term, as in the dissociation and ionization rate equations, is the fraction of collisions 

for which the translational energy is greater than the difference between the ionization energy 

of AB* and the electronic energy of AB*.

The reverse (exothermic) associative ionization reaction is similar to the forward reaction. If 

the AB* molecule is in the ground electronic state after a trial Larsen-Borgnakke 

redistribution of the collision energy, an intermediate molecule is “formed.” In the second 

step, the ionization energy is first added to the original collision energy. The new collision 

energy is again distributed among the internal energy modes of the AB* molecule, which is

then tested for dissociation as described above. If the remaining collision energy after 

redistribution plus the vibrational energy of the molecule is greater than the dissociation 

energy of the molecule, the reverse associative ionization reaction occurs. The 

corresponding rate coefficient can then be written as:

kr

AI = RcollAB+ ,e−

Υr,AI

AB* g0

AB*

exp −Eel

0 ( ) kT

gi

AB*

exp −Eel

i ( ) kT

i=0

imax

∑

(2.4.6)

where the collision rate parameter is again identical to that of Equation (2.3.8). The 

last term is included to account for the probability that the AB* molecule will be in the 

ground electronic state. The fraction of collisions with sufficient energy for the molecule 

AB* to ionize are then calculated by summing over all possible initial internal energy states of 

AB+ (the electron has none) and all possible vibrational, rotational, and electronic energy 

states of the intermediate AB* molecule as follows:

Υr,AI

AB*

= fel iel

AB+

( ) frot irot

AB+

( ) fvib ivib

AB+

( )

iel

AB*

=0

iel,max

AB*

∑

irot

AB*

=0

irot ,max AB*

∑

ivib

AB*

=0

ivib,max AB*

∑

ivib

AB+

=0

ivib,max AB+

∑

irot

AB+

=0

irot ,max AB+

∑

iel

AB+

=0

iel,max

AB+

∑

×Pvib,trans

ivib

AB*

Ecoll

0 ( )Prot,transirot

AB*

Ecoll

1 ( )Pel,transiel

AB*

Ecoll

2 ( )Q 5

2 −ωAB*

,

Ediss − Evib

AB*

kT

⎛

⎝

⎜

⎞

⎠

⎟

(2.4.7)

The sums over the internal energy levels of AB+ are over all possible levels. For the internal 

energies of AB*, the maximum level is that which is the highest energetically possible. In 

order to do this, an initial translational energy is again assumed to be the average relative 

Rcoll

A,B
43 

translational energy, or , along with the aforementioned cautions. The first three 

expressions after the summations are simply the probabilities of the AB+ molecule being in 

the respective states. The next three terms are again analogous to Equation (2.3.29) in that 

they are the probabilities of the respective internal mode of energy being at the specified 

quantum level at the specified collision energy. The collision energies are defined as:

Ecoll

0 = 5 ( ) 2 −ω kT + EionAB*

+ Eel

AB*

+ Erot

AB*

+ Evib

AB*

Ecoll

1 = Ecoll0 − EvibAB+

Ecoll

2 = Ecoll1 − ErotAB+

(2.4.8)

where each successive collision energy subtracts the previous internal energy mode. Then 

the last term, as in the dissociation and ionization rate equations, is the fraction of collisions 

for which the translational energy is greater than the difference between the dissociation 

energy of AB* and the vibrational energy of AB*

3.2.3 Charge Exchange Reactions

The charge exchange reaction is analogous to the exchange reaction outlined in Section 

2.3.2.3. The binary exchange reaction can be written as . This reaction 

takes place when the electronic energy level of the colliding neutral particle, A, after a trial 

Larsen-Borgnakke redistribution of the sum of the relative translational energy and electronic 

energy of the neutral particle is equal to the electronic energy level with energy just below 

the activation energy of the reaction. The probability of the reaction occurring is then:

Pchex = ga Ecoll − Eel

a ( )3

2−ω

gi Ecoll − Eel

i ( )3

2−ω

i=0

imax

∑

(2.4.9)

The corresponding rate equation for the forward reaction is then:

kf

chex

( ) T = Rcoll

A,B+ ga exp −Ea

A,B+

( ) kT

gi

A exp −Eeli,A ( ) kT

i=0

imax

A

∑

(2.4.10)

and for the reverse reaction is:

kr

chex

( ) T = Rcoll

A+ ,B ga exp −Ea

A+ ,B

( ) kT

gi

B exp −Eeli,B ( ) kT

i=0

imax

B

∑

(2.4.11)

5 ( ) 2 −ω kT

A + B+ ↔ A+ + B
44 

As with the exchange reactions in Section 2.3.2.3, the ratio of the forward and reverse 

reactions can be compared with the equilibrium constant and adjustments can be made to the 

phenomenological model to bring the model in closer agreement with statistical mechanics.

3.2.4 Other Exchange Reactions

In other exchange reactions involving charged species, the predominant mechanism is the 

breaking of a bond of one species (either neutral or charged) in order to exchange one atom 

to bond with the original atom, just as with an exchange reaction involving only neutral 

particles. Therefore, the logic for performing this reaction involving charged species is 

similar to that involving only neutral species.

3.3 Selecting a Reaction

Next, the Quantum-Kinetic chemistry model requires the selection of which reaction occurs. 

For many collision pairs, there may be more than one possible reaction. In the TCE 

chemistry model, the reaction probability is computed for each of the possible reactions and 

summed. This sum is then compared to a random number between zero and one. If the 

random number is less than the sum of the probabilities, a reaction is to take place. The 

selection of the reaction then proceeds by normalizing the reaction probabilities, and another 

random number is chosen to select which reaction is allowed. In the Q-K method, there is no 

clear analogy, and previous implementations have utilized a hierarchical approach where the 

reactions are considered by type first. For the neutral reactions, the dissociation reaction was 

always considered first and allowed to proceed even if an exchange reaction was a 

possibility. However, it was found that this approach tended to be biased in favor of 

dissociation reactions, especially at high temperatures where most of the collisions have 

larger relative translational energies and the vibrational level distribution is shifted to the 

higher levels. This method becomes even more undesirable when charged species reactions 

are included in the reaction set. Specifically, what happens when both dissociation and 

ionization reactions occur simultaneously? This event doesn’t actually occur for the current 

reaction set under consideration (only electron impact ionization of atomic species is 

considered herein), but this is obviously an issue that needs to be resolved.

An alternate approach to the hierarchical method was therefore implemented. This approach 

is, like the TCE model, based on assigning reaction probabilities to all of the possible 

reactions, summing them, then comparing the sum to a random number. If the random 

number is less than the sum of the probabilities, then the probabilities are normalized and 

then another random number is generated to select the reaction. The reaction probabilities 

for exchange and charge exchange reactions have already been given in Equations (2.3.16)

and (2.4.10), respectively. For the remaining reactions that occur when the collision energy 

is such that the reaction occurs (dissociation, ionization, and associative ionization), the 

probability is set equal to one if it is determined that the reaction is allowed to occur. These 

probabilities are then summed and compared to a random number between zero and one 

(obviously, if a dissociation, ionization or associative ionization reaction are possible, a 

reaction will occur). If a reaction is to occur, Pr,i is defined as the reaction probability for the 

i

th reaction, and the cumulative distribution for all reactions up to and including the jth

reaction is:
45 

Pr,i =

Pr,n

n=1

j

∑

Pr,n

n=1

NR

∑

(2.5.1)

where NR is the total number of reactions possible with the current reactants. A random 

number is then chosen to determine which reaction proceeds. 

As shown subsequently in Section 4, this selection methodology allows the sampled reaction 

rate to accurately reproduce the analytical rates shown in this Section. There is one item that 

requires mentioning at this point that the chemistry set utilized in this thesis does not include. 

That is the possibility of a dissociation and ionization reaction occurring for the same set of 

particles. For example, when molecular nitrogen collides with an electron, the molecule can 

either dissociate or ionize. Under the rules for these processes listed above, when the 

respective criterion is met for either reaction, the probability of that reaction occurring is 

unity. What happens when both have a probability of unity? Using the above reaction 

selection methodology, it is shown in Figure 3.1 that the sampled rates still accurately 

reproduce the analytical rates. This occurs because the energy required for dissociation is 

less than that for ionization, so will always occur more frequently as shown in the figure. 

The highest reaction rate physically possible would be equal to the collision rate of the 

particles. As the limit where all of the collisions have sufficient energy for both reactions is 

approached, a departure from the analytical rates occurs since both are equal, resulting in 

sampling rates one half the analytic rates. 

Figure 3.1: Sampled and analytic reaction rates when dissociation competes 

with ionization.

10-29

10-27

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical: Dissociation

Q-K Analytical: Ionization

Q-K Sampled: Dissociation

Rate Constant (m

Q-K Sampled: Ionization

3/molecule/s)

Temperature (K)
46 

Section 4 

Application of New Models to Upper-Atmosphere Hypersonic Flows

In this Section, the new electronic energy level transition and chemistry models are validated

by examining their ability to reproduce rates of chemical reactions in the literature[5-7, 9, 15, 

29, 60, 63-78]. The DSMC models are applied in a two-dimensional code where the results 

are accumulated over all cells to approximate a zero-dimensional simulation. The test gas is 

comprised of eleven species: O, N, O2, N2, NO, O+, N+, O2

+

, N2

+

, NO+, and e-. The 

simulations were run at various temperatures between 10,000 K and 60,000 K and at a 

continuum number density of 1.0e23 /m3 in an adiabatic box 0.002 m on a side with 2 million 

particles. During the simulations, relaxation of the particles (including rotational, vibrational 

and electronic) was allowed to proceed as usual, but when a reaction was determined to take 

place, the number of reactions was advanced by one but the simulators were left unchanged, 

to maintain the thermodynamic energy balance with energy neither added nor removed from 

the gas. 

4.1 Measurement of Reaction/Transition Rates

For a bi-molecular reaction , the forward reaction rate is calculated as:

kf =

dnA

dt

nAnB

(3.1.1)

where the numerator is the change in number density of species A due to the reaction in the 

sampled time and the denominator is the product of the number densities of species A and B. 

The corresponding reverse reaction rate is:

kr =

dnC

dt

nCnD

(3.1.2)

For a termolecular reaction , the forward reaction rate is calculated as:

kf =

dnAB

dt

nABnT

(3.1.3)

as the bi-molecular reaction rate is, but the reverse reaction rate is now:

kr =

dnA

dt

2nAnBnT

(3.1.4)

The capability of these techniques to reproduce measured equilibrium reaction rates will now 

be investigated. However, the uncertainty of the currently available measured rates for 

A + B ↔ C + D

AB + T ↔ A + B + T
47 

atmospheric hypersonic flows can be large. Most of the experimental data on which reaction 

rate sets are based were obtained for low temperatures (2,000 – 7,000 K). Even though 

uncertainty may be large, Park[60] suggests that the experimentally obtained reaction rates 

can be extrapolated up to 30,000 K.

The uncertainty of these measured rates is not always reported and in some cases is 

unknown. The problem becomes even more severe in the extrapolated portion of the 

temperature range (T ≥ 10,000 K). Park[60] reports that some rates measured for the same 

reaction differ by more than one order of magnitude, which will be apparent in the figures 

that follow. Since this work is focused on the extension of the Q-K model to electronic 

energy levels and reactions involving charged particles, only these transition and reaction 

rates will be discussed in this Section. However, comparisons of the analytical Q-K rate 

expressions, the sampled Q-K rates, and rates found throughout the literature for every 

reaction included in this study are presented in Appendix C.

4.2 Definition of Translational, Rotational, Vibrational, and Overall

Temperatures

There are various ways to define temperatures with regards to the DSMC methodology. In 

this work, we will refer to statistical mechanics for the definition of temperatures by defining 

the energy as:

ε = 1

2ζ RT (3.2.1)

where ε is the specific energy of the quantity of interest, ζ is the number of degrees-offreedom, R is the universal gas constant, and T is the temperature. By finding the energy, 

recognizing that , and rearranging Equation (3.2.1), one has the temperature of 

interest. For the translational temperature, we define εtrans as:

ε trans = 1

2 c′

2 (3.2.2)

where is the average square of the thermal velocity defined as:

c′

2 = c2 − 2

c



c0 + c0

2 (3.2.3)

where is the mean molecular velocity and is the macroscopic stream or mean velocity. 

The translational temperature can then be defined as:

Ttrans = 2ε trans

3k (3.2.4)

since there are 3 translational degrees of freedom.

For the internal degrees of freedom of a specific species n, we define the average internal 

energy of a single particle as:

k = mR

c′

2



c c0
48 

Eint

n = mnε intn (3.2.5)

and the temperature is then:

Tint

n = 2Eint

n

kζ int

n (3.2.6)

For the overall internal temperature averaged over all of the N species, we then define:

Tint =

ζ int

n Tintn

n=1

N

∑

ζ int

n

n=1

N

∑

(3.2.7)

For rotational energy, the average energy of each species is found by summing the rotational 

energies of the particles in the computational cell and dividing by the number of particles. 

For the rotational degrees of freedom, it is assumed that the mode is fully excited and is set to 

2 for diatomic molecules.

For vibrational energy, the average energy of each species is again found by summing the 

vibrational energies of the particles and dividing by the number of particles. To then find the 

number of vibrational degrees of freedom for the species, using the definition of the average 

vibrational energy level in = Evib

n kΘvib and some manipulation, one gets:

ζ vib

n = 2in ln 1+1 i ( ) n (3.2.8)

The electronic temperature is treated in a similar manner, but some modifications must be 

made to account for the discretized spacing of the electronic energy levels. To begin, the 

species electronic temperature is defined as:

Tele

n = Eele

1,n − Eele0,n

k ln n0

n

g1

n

n1

n

g0

n

⎛

⎝

⎜

⎞

⎠

⎟

(3.2.9)

where the 0 and 1 denote the ground and first excited electronic energy levels, g is the 

degeneracy of the level, and n can be taken as the number density or total number of samples 

(result is the same since the ratio is being taken). This definition comes directly from the 

ratio of the Boltzmann distribution of the ground and first excited electronic energy levels. 

Then, the effective number of degrees of freedom is defined as:

ζ ele

n = 2Eele

n

kTele

n (3.2.10)
49 

as we have defined the effective degrees of freedom of the other internal modes.

4.3 Electronic Energy Level Transitions and Comparison to 

Equilibrium

The following section applies the methodologies developed in Section 3.1 to the 

aforementioned equilibrium heat bath. There are excitation and quenching processes 

between each of the atomic/molecular energy levels, resulting in a large number of rates 

being required for previous electronic energy level transition models, both for CFD and 

DSMC. The new Q-K based transition model negates the necessity for each of these 

transition rates (most of which use approximate analytic formulas). In the following subsections, only a portion of the transitions presented in Appendix C are discussed for 

exemplary purposes only.

4.3.1 Electron-Impact Excitation

The process of electron-impact excitation is significant for both atoms and molecules in 

hypersonic shock-layer applications. To begin, the electron-impact transitions of atomic 

nitrogen are examined. The sampled transition rates are compared to values reported by 

Tayal[79], Park[60], and Frost[80]. The transitions from the first to the fifth level and from 

the second to the fifth level are presented in Figure 4.1 and Figure 4.2, respectively. It is 

important to note that the 1-5 transition is considered a “forbidden” transition, while the 2-5 

transition is an “allowed” transition. Even so, the phenomenological Q-K model compares 

very well to the reported rates from the literature.

Figure 4.1: Comparison of sampled and reported electronic energy level 

transition rates from the first to fifth energy levels of atomic 

nitrogen.

10-25

10-24

10-23

10-22

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N1 + e- > N5 + e-
50 

Figure 4.2: Comparison of sampled and reported electronic energy level 

transition rates from the second to the fifth energy levels of 

atomic nitrogen.

For the electron-impact transitions of atomic oxygen, the sampled transition rates are 

compared to values reported by Zatsarinny[81], Park[60], and Bhatia[82]. Once again, even 

though one transition is allowed (level 1-5, Figure 4.3) and one is forbidden (level 2-5, 

Figure 4.4), the comparisons are quite favorable.

Figure 4.3: Comparison of sampled and reported electronic energy level 

transition rates from the first to fifth energy levels of atomic 

oxygen.

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Park[62]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N2 + e- > N5 + e10-24

10-22

10-20

10-18

10-16

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Zatsarinny[83]

Park[62]

Rate Constant (m

Bhatia[84]

3/molecule/s)

Temperature (K)

O1 + e- > O5 + e-
51 

Figure 4.4: Comparison of sampled and reported electronic energy level 

transition rates from the second to the fifth energy levels of 

atomic oxygen.

The electron-impact transition rates of the molecular species are compared to rates presented 

by Cherni[83], Telulet[84], and Gorelov[85]. Example transition rates for molecular 

nitrogen, nitric oxide, and the molecular nitrogen ion are shown in Figure 4.5, Figure 4.6, 

and Figure 4.7, respectively. Again, good agreement is shown between the sampled rates 

and the reported rates.

Figure 4.5: Comparison of sampled and reported electronic energy level 

transition rates from the second to the third energy levels of 

molecular nitrogen.

10-24

10-22

10-20

10-18

10-16

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Zatsarinny[83]

Rate Constant (m

Bhatia[84]

3/molecule/s)

Temperature (K)

O2 + e- > O5 + e10-18

10-17

10-16

10-15

10-14

10-13

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Cherni[85]

Telulet[86]

Rate Constant (m

Gorelov[87]

3/molecule/s)

Temperature (K)

N2

2

 + e- > N2

3

 + e-
52 

Figure 4.6: Comparison of sampled and reported electronic energy level 

transition rates from the first to the fourth energy levels of nitric 

oxide.

Figure 4.7: Comparison of sampled and reported electronic energy level

transition rates from the first to the third energy levels of the 

molecular nitrogen ion.

4.3.2 Heavy Particle Impact Excitation

Transitions due to heavy particle collisions are similar to electron-impact transitions, but they 

are significantly less important for most applications, as discussed by Johnston[86], where he 

shows that the process is negligible for atomic and molecular species in lunar-return shock 

layers. The comparisons shown in this section are, for the most part, as favorable as those in 

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Telulet[86]

Rate Constant (m

Gorelov[87]

3/molecule/s)

Temperature (K)

NO1 + e- > NO4 + e10-20

10-19

10-18

10-17

10-16

10-15

10-14

10-13

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Gorelov[87]

Rate Constant (m

Telulet[86]

3/molecule/s)

Temperature (K)

N2

+ 1 + e-

 > N2

+ 3 + e-
53 

the previous section for electron-impact transition rates. Comparisons are made with values 

reported in References [85, 87-93]. For example, the transitions of the molecular nitrogen

ion when colliding with molecular nitrogen, as shown in Figure 4.8, is generally within an 

order of magnitude of the reported values (note the 3-4 order of magnitude in the spread of 

data for the level 1-2 transition). There is also good agreement for collisions of molecular 

nitrogen and nitric oxide for most transitions, as exemplified in Figure 4.9 and Figure 

4.10, respectively.

Figure 4.8: Comparison of sampled and reported electronic energy level

transition rates for the molecular nitrogen ion when colliding with 

molecular nitrogen.

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

1 - 2 Flagan[91]

Rate Constant (m

1 - 2 Nagy[92]

3/molecule/s)

Temperature (K)
54 

Figure 4.9: Comparison of sampled and reported electronic energy level 

transition rates for molecular nitrogen when colliding with atomic 

nitrogen.

Figure 4.10: Comparison of sampled and reported electronic energy level 

transition rates for nitric oxide when colliding with molecular 

nitrogen.

However, there seem to be certain combinations of heavy-impact partners where the current

phenomenological model can’t accurately reproduce the measured/calculated transition rates. 

These combinations generally include transitions where the collision partner is molecular 

nitrogen, as shown in Figure 4.11, and the transition is from the lower levels (i.e., 1-2 or 

even some times 2-3).

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

Rate Constant (m

1 - 2 Gorelov[89]

3/molecule/s)

Temperature (K)

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

Rate Constant (m

3 - 4 Park[90]

3/molecule/s)

Temperature (K)
55 

Figure 4.11: Comparison of sampled and reported electronic energy level 

transition rates for the nitric oxide when colliding with molecular 

nitrogen.

The current model also appears to have some difficulty in reproducing electronic energy 

level transitions in collisions of molecular oxygen, as shown in Figure 4.12. However, 

even though these comparisons have been less than ideal, they are not completely unreliable 

considering the four to five order of magnitude spread in the rates previously reported from 

the literature.

Figure 4.12: Comparison of sampled and reported electronic energy level 

transition rates for molecular oxygen when colliding with atomic 

oxygen.

10-27

10-25

10-23

10-21

10-19

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

NO-O2 3-4

1 - 2 Park[90]

2 - 3 Park[90]

Rate Constant (m

3 - 4 Park[90]

3/molecule/s)

Temperature (K)

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

Rate Constant (m3/molecule/s)

Temperature (K)
56 

4.3.3 Equilibrium Conditions

In order for the electronic energy level procedures to be acceptable, they need to 

reproduce equilibrium conditions after the simulation has been allowed to run for 

a sufficient amount of time to ensure that detailed balance has been enforced. The 

first requirement is for the sampled translational and electronic temperatures to 

remain very near the input equilibrium temperature. The sampled translational, 

electron, and electronic temperatures for all of the species are presented in

Table B.1 at 10,000, 20,000, and 30,000 K. All of the temperatures depart from the set 

equilibrium value by less than 0.5%.

The resultant population distribution for the first seven energy levels of all atomic and 

molecular species (where applicable) for a range of temperatures is compared to the 

Boltzmann distribution in Table B.2 through Table B.31. Again, the simulation represents 

the equilibrium state very well, with sampled departures from the Boltzmann distribution less 

than 1% for most levels where the number fraction in the state is above 1.0x10-3. The larger 

departure from the equilibrium distribution is expected for levels with a very small 

population due to statistical noise.

4.4 Associative Ionization Reactions

Because of their low threshold energies, associative ionization reactions play an important 

role in determining the level of ionization and structure of the flow field at hypersonic 

reentry conditions. There exists very little experimental or computational data in the 

literature regarding the associative ionization reactions in air. The associative ionization 

reaction has the lowest energy threshold, however the associative 

ionization reactions of and are also very important 

mechanisms in the current application. The sampled DSMC reaction rates of endothermic

and exothermic reactions of the N + O reaction are compared to measured rates from the 

literature in Figure 4.13 and Figure 4.14, respectively.

N + O ↔ NO+ + e−

N + N ↔ N2

+ + e− O + O ↔ O2+ + e−
57 

Figure 4.13: Endothermic associative reaction rates for 

compared to rates from the literature.

Figure 4.14: Exothermic associative reaction rates for 

compared to rates from the literature.

In these figures, the black lines are the analytical solutions from Equations (2.4.3) and (2.4.6)

, the black symbols are the sampled rates from the DSMC code, and the first of the values 

from literature is the curve used by LAURA. As discussed earlier, there is a wide spread in 

the reaction rates reported in the literature. Given the assumed measurement uncertainty, the 

Q-K reaction rates are in good agreement with the rates from the literature.

10-23

10-22

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Park[62]

Kang[66]

Park[71]

Rate Constant (m

Gupta[74]

3/molecule/s)

Temperature (K)

N + O > NO+ + eN + O → NO+ + e−

10-18

10-17

10-16

10-15

10-14

0 5000 10000 15000 20000 25000 30000 35000

Q-K Analytical

Q-K Sampled

Park[65]*

Park[62]

Bird[15]

Ozawa[7]

Boyd[79]

Boyd[69]

Kang[66]

Gupta[74]

Rate Constant (m3/molecule/s)

Temperature (K)

NO+ + e- > N + O

NO+ + e− → N + O

59 

Figure 4.16: Endothermic charge exchange reaction rates for 

 to rates from the literature.

Figure 4.17: Exothermic charge exchange reaction rates for 

compared to rates from the literature.

4.7 Exchange Reactions Involving Charged Species

The remainder of relevant chemical reactions involves various exchange reactions involving 

charged species. As an example, the endothermic and exothermic reaction rates for the 

reaction are shown in Figure 4.18 and Figure 4.19, respectively, 

compared to values from the literature as shown.

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Rate Constant (m

Phelps[75]

3/molecule/s)

Temperature (K)

N+ + N2 > N2

+

 + N

N+ + N2 → N2

+ + N

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]*

Bird[15]

Kang[66]

Rate Constant (m

Boyd[79]

3/molecule/s)

Temperature (K)

N2

+

 + N > N+ + N2

N2

+ + N → N+ + N2

NO+ + N ↔ N2

+ + O
60 

Figure 4.18: Endothermic reaction rate for the exchange reaction involving 

charged species with comparisons to rates 

from the literature.

Figure 4.19: Exothermic reaction rate for the exchange reaction involving 

charged species with comparisons to rates 

from the literature.

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]

Rate Constant (m

Bird[15]

3/molecule/s)

Temperature (K)

NO+ + N > N2

+

 + O

NO+ + N ↔ N2

+ + O

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]*

Bird[15]

Boyd[79]

Rate Constant (m3/molecule/s)

Temperature (K)

N2

+

 + O > NO+ + N

N2

+ + O ↔ NO+ + N
61 

Section 5 

Modeling of the FIRE II Reentry Flow with Ionization and 

Electronic Energy Levels with DSMC

In this Section, the newly developed DSMC models presented in previous Sections are 

exercised on a demanding near-continuum reenty flow field with ionization. These results 

are compared with results from traditional DSMC models, a solution of the FIRE II flow 

field obtained using current best practices in computational fluid dynamics (CFD), and

values inferred from flight. Both flow field and surface parameters are examined. To this 

end, three DSMC cases are summarized in Table 5.1. 

Table 5.1: List of DSMC FIRE II cases.

Case Electronic Energy Model Chemistry Model 

Case 1 None TCE

Case 2 Present TCE

Case 3 Present Q-K 

5.1 The Project FIRE Flight Experiment

The Project FIRE (Flight Investigation Reentry Environment) was an Apollo-era flight 

experiment to measure the radiative and convective heating on spacecraft material during 

atmospheric reentry at lunar return speeds (11.35 km/s)[21]. Project FIRE consisted of two 

sub-orbital flight experiments launched April 14, 1964 and May 22, 1965 from Cape 

Kennedy, Florida. For the current study, computations are performed for the reentry 

conditions of the second flight at 1634 seconds, which corresponds to an altitude of 76.42 

km, and conditions for the simulation are provided in Table 5.2. The geometry and 

dimensions of the second vehicle in the experimental series, FIRE II, are presented in Figure 

5.1. The forebody of the FIRE II reentry vehicle consisted of three phenolic-asbestos heat 

shields sandwiched between beryllium calorimeters. The first two heat shield and 

calorimeter packages were designed to be ejected after the onset of melting, yielding heating 

data free of the effects of ablation and three separate data gathering periods. Calorimeter 

plugs and radiometers were located at various positions on the heat shields. The flight 

condition examined in this study occurred during the first data collection period.

Table 5.2: Free stream and surface parameters for FIRE II entry condition.

Time (s) Altitude (km) ρ∞ (kg/m3

) n∞ (m-3) V∞ (m/s) T∞ (K) KnD Twall (K) XN2 XO2

1634 76.42 3.72E-05 7.73e20 11360 195 2.59E-03 615 0.76 0.24
62 

Figure 5.1: Schematic of FIRE II test vehicle.

5.2 Details of the Continuum Flow Model and Non-Equilibrium 

Radiation Model

Two different numerical approaches are applied: a CFD solution of the Navier-Stokes 

equations using the Langley Aerothermodynamic Upwind Relaxation Algorithm code and 

particle-based DSMC computations using the DAC code, which is the focus of this study. 

All computations presented herein were steady state and axisymmetric with the x-direction 

along the axis. In each case, an 11-species, 21-reaction, thermochemical non-equilibrium 

approach is adopted. The surface catalytic boundary condition used in the CFD solution was 

“fully-catalytic,” which means that all atomic species recombine to their molecular 

counterpart (i.e., ) and also charged heavy species recombine with electrons to 

form their charge neutral counterpart (i.e., ).

The Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) is a high fidelity, 

structured grid computer code, specialized for hypersonic re-entry physics, utilizing state-ofthe-art algorithms for computational fluid dynamic simulations[94, 95]. Key elements of 

LAURA include Roe’s averaging[96] and Yee’s Symmetric Total Variation Diminishing 

(STVD)[97] formulation of second-order, inviscid flux. Yee’s STVD formulation has been 

found to be exceptionally robust and Courant number independent using first point-implicit 

and then line-implicit relaxation for hypersonic flow simulations.

A two-temperature thermochemical non-equilibrium model[95] is applied. The chemical 

reaction rates are compiled from previous studies (listed as the references for the first 

reaction rate parameters for each reaction in Table A.14). The thermophysical properties 

are taken from the work of Mcbride et al.[98] and Gupta et al.[76]. Multicomponent 

diffusion is approximated using Sutton and Gnoffo’s[99] approximate-corrected approach.

0.6716&m

R&0.9347&m

R&0.0102&m

0.5258&m

0.4122&m



N + N → N2

N+ + e− → N
63 

The shock-layer radiation is modeled with the HARA (High-temperature 

Aerothermodynamic RAdiation) code. The details of this code for treating air species are 

presented by Johnson et al.[19, 20]. Briefly, it is based on a set of atomic levels and lines 

obtained from the National Institute of Standards and Technology (NIST) online 

database[100] and the Opacity Project[101], as well as atomic bound-free cross sections from 

the TOPbase[102]. The negative nitrogen and oxygen ions are treated using cross sections 

suggested by Soon and Kunc[103] and Chauveau et al.[104], respectively. The molecular 

band systems are treated using a smeared-rotational band (SRB) model[105], which was 

shown by Johnston et al.[19] to be sufficient for treating VUV absorbing and optically-thin 

emitting band systems in air. The molecular data for modeling these band systems are 

obtained from Laux[106], except for the VUV N2 systems, which are obtained from various 

other sources[107-109]. The non-Boltzmann modeling of the atomic and molecular 

electronic states is based on a set of electron-impact excitation rates compiled from the 

literature and presented in detail by Johnston et al.[20]. Following the work of Park[110], 

the quasi-steady state assumption is made when solving the Master Equation. The tangentslab approximation is applied to calculate radiative flux and the divergence of the radiative 

flux.

5.3 DSMC Simulation Parameters

In this section, the parameters used in the simulation of the FIRE II flight experiment are 

discussed.

5.3.1 Collision Modeling

The baseline collision model implemented in DAC, for both neutral and charged species 

interactions, is that originally published by Bird[41], in which all particle interactions are 

treated using the VHS collision model. The reference diameters of the heavy ions O2

+

, N2

+

, 

O+, N+, and NO+ are set equal to the diameters of their associated neutral particles, and the 

reference diameter of the electrons is set to 1e-10 m. This value of the electron reference 

diameter is much larger than the classical value of 6e-15 m, and is used to produce electron 

collision cross sections of the required magnitude within the mathematical framework of the 

VHS collision model. A complete list of the VHS reference diameters for all species and the 

parameters Tref and ω is given in Table A.1. As shown previously, the present choice of 

electron reference diameter does not introduce significant error in the DSMC results.

5.3.2 Internal Energy Modeling

For the current thesis, the rotational collision number is set to unity. The reason for this is 

because the solutions are being compared to a computational fluid dynamics (CFD) solution 

where the translational and rotational temperatures are assumed equal. Therefore, for DSMC 

to best match CFD, the rotational energy is sampled from the equilibrium distribution at each 

collision (Zrot = 1.0). The vibrational energy relaxation model used is that based on the work

of Millikan and White and the parameters used for each of the molecular species are 

presented in Table A.2.
64 

One goal of this work is to examine the effects of adding an electronic energy transition 

algorithm to DSMC solutions of hypersonic (re)entry flows. To that end, the electronic 

energy levels of all relevant simulated species (O2, N2, O, N, NO, O2

+

, N2

+

, O+, N+, and NO+) 

are included in Cases 2 and 3. The energy levels included can be viewed in Table A.4

through Table A.13. The energy levels for O and N were taken from Johnston[86], while 

the remaining species were taken from Hash[111]. A more comprehensive set of energy 

levels could be compiled for O+ and N+, but these are the energy levels used by the CFD code 

to be compared with, as well as the nonequilibrium radiation code.

5.3.3 Chemistry Modeling

The baseline chemistry model is the TCE model using the reaction rates used in LAURA. 

These rates are listed in Table A.14 as the first rate listed for each reaction. As discussed in 

Section 2.2.3, the exothermic (reverse) reaction rates are then determined by considering the 

statistical mechanics equilibrium constant, and then the coefficients are chosen from those in 

the literature to most closely match the rate as computed by statistical mechanics.

The other chemistry model, which is one of the main contributions, is the Q-K chemistry 

model and the extensions to charged species. The characteristic temperatures and reference 

parameters for vibrational relaxation can be found in Table A.1 and Table A.3, 

respectively. Additionally, the parameters required to account for the statistical mechanics 

equilibrium constant and adjust the activation energies for recombination and 

exchange/charge exchange can be found in Table A.16 and Table A.17, respectively. 

Comparisons of the Q-K reaction set to various rates found in the literature can be viewed in 

Appendix C.

5.3.4 Electric Field Modeling

The electric field model implemented in the current study, which was also added to DAC, is 

a hybrid of Bird’s model[26] and that of Boyd[30]. In this model, the electron particles are 

moved along with their associated ion during the movement portion of the DSMC algorithm. 

However, the electrons retain their individual velocity components, which are used in the 

collision portion of the DSMC algorithm. Neither the electrons nor ions actually experience 

an electric field, in that their velocity components are not adjusted to account for an electric 

field as one is not computed.

5.3.5 Gas-Surface Interactions

The surface catalycity must be treated differently with the present DSMC simulations 

because the CFD solutions use the fully-catalytic surface boundary condition to account for

the higher surface temperatures in the continuum regime. However, there is evidence[25]

that, for surface temperatures such as those at the higher altitudes for the FIRE II experiment, 

the catalycity of neutral atomic species recombining is much closer to zero. Therefore, the 

interaction of particles with the vehicle surface is managed by assuming both a fully catalytic 

surface as was used with the CFD solution for comparison and also a partially catalytic 

surface that only allows charged species to become neutralized. The fully catalytic 

simulations will be referred to simply by their case number (i.e., Case 1) and the partially

catalytic simulations will have an extra notation indicating that it is partially catalytic (i.e., 
65 

Case 1pc). The remaining collisions (neutral molecules), along with the recombined particles, 

are then allowed to reflect diffusely from the surface, after thermally accommodating to the 

specified wall temperature.

5.4 FIRE II Simulations

In this section, a progression of DSMC solutions is discussed and compared with CFD 

beginning with a traditional DSMC solution without electronic energy levels (Case 1), 

followed by the addition of electronic energy levels (Case 2), concluding with the Q-K

chemistry model (Case 3). At the end of the section, a discussion of the surface heating, both 

convective and radiative, is included where the DSMC rates are compared to both the CFD 

solution and values inferred from flight.

5.4.1 Baseline Solution

The baseline solution, referred to as Case 1 in Table 5.1, is the FIRE II simulation at 1634 

seconds using the TCE chemistry model without the addition of electronic energy levels. 

Various temperatures for this case are compared to that from the CFD solution in Figure 

5.2. Clearly, the exclusion of electronic energy levels from the DSMC solution changes the 

shock-layer profile. The shock standoff distance of the DSMC solution is nearly 40% greater 

than that of the CFD solution, while the shock-layer temperature is also greater than that 

predicted by CFD. It can be seen that, inside the shock layer of the DSMC solution, the 

temperatures are not in equilibrium. This has been reported in Refs. [25] and [4] as an effect 

of the addition of ionization physics at this altitude. There is some scatter present in the data 

for both the rotational and vibrational temperatures, which may be caused by the low number 

density of molecular species in the flow field. The spike in vibrational temperature outside 

of the shock layer is the result of trace species of vibrationally excited molecules that have 

escaped from the shock layer. This is a phenomenon typical of flows such as this in DSMC 

simulations. Since the physics are clearly not the same, comparisons of species 

concentrations are not provided for this case.
66 

Figure 5.2: Translational temperature comparison between CFD and Case 1

(no electronic energy levels).

5.4.2 Addition of Electronic Energy Levels

Electronic energy levels were included in DSMC Case 2 and are more appropriate for 

comparison to the CFD solution. The various temperature profiles from CFD and Case 2 are 

compared in Figure 5.3. The DSMC shock standoff distance has been reduced as compared 

to Case 1 and now the location and magnitude of the peak translational temperatures of the 

CFD and DSMC Case 2 solutions compare quite well. The differences between the solutions 

are the larger equilibrium temperature within the shock layer and the increased electronic 

temperature that persists upstream of the location of maximum translational temperature, 

which is discussed below.

0

1 104

2 104

3 104

4 104

5 104

-0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0

CFD: Ttrans, Trot

CFD: Tvib, Tele

Case 1: Ttrans

Case 1: Trot

Case 1: Tvib

Case 1: Telectron

Temperature (K)

X (m)
67 

Figure 5.3: Comparison of temperature profiles between the CFD solution 

and Case 2 (electronic energy levels added).

Comparisons of species mole fractions between the CFD solution and Case 2 are presented in 

Figure 5.4 and Figure 5.5 for neutral and charged species, respectively. In general, the 

species fractions compare well. As with the vibrational and rotational temperatures, there is 

noise in the mole fractions extracted from the DSMC simulation below around 0.05% 

because of the low sample size of these species. The DSMC solution contains trace species 

further upstream of the shock front as compared to the CFD solution. This is most likely 

caused by rarefaction effects where there are insufficient collisions to keep all particles 

within the shock layer. This also accounts for the increased electronic temperatures upstream 

of the maximum translational temperature as seen in Figure 5.3, and is analogous to what 

occurs with the vibrational temperature upstream of the shock, except with the electronic 

energy levels, not only do the molecules carry energy upstream of the shock, but also the 

atomic species. Even for a small population of excited particles, the electronic temperature 

can be large.

0

5000

1 104

1.5 104

2 104

2.5 104

3 104

3.5 104

4 104

-0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0

CFD: Ttrans, Trot

CFD: Tvib, Tele

Case 2: Ttrans

Case 2: Trot

Case 2: Tvib

Case 2: Tele

Case 2: Telectron

Temperature (K)

X (m)
68 

Figure 5.4: Comparison of neutral species mole fractions between the CFD 

solution and Case 2.

Figure 5.5: Comparison of charged species mole fractions between the CFD 

solution and Case 2.

5.4.3 Quantum-Kinetic Solution

As was discussed earlier, DSMC Case 3 switches chemistry modules from the traditional 

TCE model to the new Q-K model. The various temperature profiles from CFD and Case 3 

are compared in Figure 5.6. For this case, the peak translational temperature is slightly 

lower and moved upstream as compared to the CFD solution. This discrepancy may

attributable both to the difference in the resultant equilibrium reaction rates discussed in 

Section 4 and also because the Q-K model does not rely on equilibrium-based information. 

10-8

10-7

10-6

10-5

0.0001

0.001

0.01

0.1

1

-0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0

CFD: XO2

CFD: XN2

CFD: XO

CFD: XN

CFD: XNO

Case 2: XO2

Case 2: XN2

Case 2: XO

Case 2: XN

Case 2: XNO

Mole Fraction

X (m)

10-8

10-7

10-6

10-5

0.0001

0.001

0.01

0.1

1

-0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0

CFD: XO2+

CFD: XN2+

CFD: XO+

CFD: XN+

CFD: XNO+

CFD: XeCase 2: XO2+

Case 2: XN2+

Case 2: XO+

Case 2: XN+

Case 2: XNO+

Case 2: XeMole Fraction

X (m)
69 

If the Q-K reaction rates are compared to the reaction rates used in the CFD solution (the first 

referenced rate in each figure) in Appendix C, there can be a large discrepancy, even though 

the Q-K rate is generally still well within the spread of reported reaction rates.

Figure 5.6: Comparison of temperature profiles between the CFD solution 

and Case 3 (with electronic energy levels and Q-K chemistry).

Comparisons of species mole fractions are presented between the CFD solution and DSMC 

Case 3 in Figure 5.7 and Figure 5.8 for neutral and charged species, respectively. Again, 

the mole fractions generally compare well between the CFD and Case 2 solutions with some 

notable differences. One difference is again the presence of trace species well upstream of 

the shock layer.

0

5000

1 104

1.5 104

2 104

2.5 104

3 104

3.5 104

4 104

-0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0

CFD: Ttrans, Trot

CFD: Tvib, Tele

Case 3: Ttrans

Case 3: Trot

Case 3: Tvib

Case 3: Tele

Case 3: Telectron

Temperature (K)

X (m)
70 

Figure 5.7: Comparison of neutral species mole fractions between the CFD 

solution and Case 3.

Figure 5.8: Comparison of charged species mole fractions between the CFD 

solution and Case 3.

5.4.4 Comparison of Surface Heating

The FIRE II calorimeter data presented by Cornette[21] provide heating values that contain 

the convective heating plus the contribution of the radiative flux absorbed by the beryllium 

calorimeter. In order to compare to these data, the code HARA was used to post-process the 

flow field solutions to obtain the non-equilibrium radiation heating to the surface. It was 

determined[21] that the fraction of the radiative flux absorbed by the calorimeter, α, was 

0.72.

10-8

10-7

10-6

10-5

0.0001

0.001

0.01

0.1

1

-0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0

CFD: XO2

CFD: XN2

CFD: XO

CFD: XN

CFD: XNO

Case 3: XO2

Case 3: XN2

Case 3: XO

Case 3: XN

Case 3: XNO

Mole Fraction

X (m)

10-8

10-7

10-6

10-5

0.0001

0.001

0.01

0.1

1

-0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0

CFD: XO2+

CFD: XN2+

CFD: XO+

CFD: XN+

CFD: XNO+

CFD: XeCase 3: XO2+

Case 3: XN2+

Case 3: XO+

Case 3: XN+

Case 3: XNO+

Case 3: XeMole Fraction

X (m)
71 

The computed heat flux to the surface of the vehicle is compared to the values inferred from 

the flight data at a radial location of approximately 0.02 m in Table 5.3. When comparing 

the convective heating component between the CFD and fully catalytic DSMC solutions

including electronic energy levels, very good agreement is found. The solutions agree to 

within 1% of each other. However, DSMC Case 1, without electronic energy levels, has a 

convective heating rate approximately 15 W/cm2 higher, again highlighting the fact that 

electronic energy levels are required for the simulation. When the cases with electronic 

energy levels are compared to the inferred flight value, the convective heating is overpredicted by 34%. If the partially catalytic solutions are then considered, they under-predict 

the convective heating by around 15%. Therefore, the surface of the flight vehicle was likely 

at least partially catalytic to atomic recombination.

The radiative heating component is also listed in Table 5.3 in the third column. Only the 

DSMC simulations with electronic energy levels modeled have been post-processed for 

radiative heating. Because of the differences in molar fractions between the DSMC 

simulations and the CFD solution, there is about a 10 W/cm2 increase in radiative heating for 

the DSMC simulations. There are generally more species present in the shock layer that 

radiate in the DSMC simulations, such as N, O, and N2

+

. The increase in radiative heating 

for the partially catalytic DSMC simulations is accounted for by the increased levels of N 

and O near the surface since these species do not recombine upon colliding with the wall.

Table 5.3: FIRE II surface heating at r = 0.02 m.

Case Convective Heating (W/cm2) Radiative Heating (W/cm2) qconv + αqrad

(W/cm2

)

CFD 217 25 235

Case 1 230 --- ---

Case 1pc 148 --- ---

Case 2 219 36 245

Case 2pc 133 95 201

Case 3 217 35 242

Case 3pc 137 72 189

Flight 162 18 175
72 

Section 6 

Conclusions

This section contains a summary of the results of this thesis, the original contributions made 

to the field, and identifies future research directions.

6.1 Summary of Results

Extensions to the recently introduced kinetic-theory approach for computing chemical 

reaction rates (the Q-K method) have been developed and validated for DSMC for electronic 

energy level transitions and reactions involving charged species. These new models make 

use of the principles of microscopic reversibility and molecular-level energy exchange to 

predict the probability that an energy level transition or a chemical reaction occurs during a 

collision between two particles. Procedures have been defined that are required to implement 

these models in the DSMC technique.

For each of these new Q-K reaction-rate models, comparisons have been made between the 

predicted analytic curves, the sampled transition or reaction rates from DSMC simulations, 

and rates obtained from the literature. The analytic and sampled values compare favorably to 

the rates from the literature; the historical rates often differ by several orders of magnitude 

for the same reactions. It should be noted again, however, that the Q-K model is 

phenomenological in nature and there may be cases for which it leads to rates that are outside 

of the range of uncertainty of the measured rates. Even so, the usefulness of these models 

has been shown as engineering models, especially where there are no data or the uncertainty 

in reported values is quite large. If there are reactions for which there is a high degree of 

certainty in the reported rate, they can be specified to use the TCE model.

The proposed extensions to the Q-K model were implemented in the DSMC Analysis Code 

(DAC) of NASA and tested on the Project Fire II flight experiment. The addition of 

electronic energy levels to the standard chemistry model resulted in better agreement 

between DSMC and CFD. Both the location and magnitude of the maximum translational 

temperature agree. The DSMC solution without electronic energy levels had predicted a 

higher peak translational temperature and shock thickness approximately 40% greater than 

the CFD solution. However, because of rarefaction effects, trace species with elevated 

electronic energies were found upstream of the shock layer, resulting in an elevated 

electronic temperature. Similar results were found when the chemistry model was changed 

from the TCE model to the Q-K model. The only major difference was a slightly thicker 

shock layer and lower peak translational temperature because of the differences in resultant 

reaction rates. Most likely, the Q-K model better represents the reaction rates because it is a 

phenomenological approach, whereas for the other models, the reaction rates are extrapolated 

from low temperature experiments.

Once surface convective heating had been obtained, the flow field was post-processed using 

the HARA non-equilibrium radiation solver to compare DSMC, CFD and values obtained 

from the flight experiment. The convective heating of the fully catalytic DSMC simulations 

with electronic energy levels included compared to within 1% of values obtained using CFD. 

These values are, however, around 34% higher than the value inferred from flight. When a 

partially catalytic surface boundary condition is employed, allowing only electron/ion 
73 

recombinations, the convective heating dropped to levels on the order of 15% lower than the 

flight value. Therefore, the fully catalytic and partially catalytic simulations have been 

shown to bound the flight value of convective heating. The radiative heating computed from 

the DSMC simulations is around 10 W/cm2 higher than that of the CFD simulation, which is 

caused by slightly larger amounts of N, O, and N2

+

. However, the nonequilibrium radiation 

solver is required to make certain assumptions in its calculations, so a direct input of the 

sampled non-Boltzmann distribution of electronic energy levels from the DSMC simulation 

would negate the necessity of the assumptions.

6.2 Contributions

The ability to simulate rarefied, ionized, hypersonic flows has been added to the suite of 

simulation capabilities for the National Aeronautics and Space Administration. Specifically, 

physical models for the treatment of charged particles using the DSMC method have been 

implemented in DAC. It has been shown that the inclusion of electronic energy levels, which 

few other DSMC solvers have the capability to do, is necessary for lunar-reentry velocity 

missions for the accurate prediction of shock shape and surface heating. The ability to 

predict electronic energy level transitions and resulting distributions using only energy level 

data has been developed as a result of this thesis. Also, the QK chemistry model was 

extended to include reactions involving charged species. This model allows the inclusion of 

chemical reactions for which specific reaction rate data is not known and may be more 

appropriate for non-equilibrium conditions. Both electronic energy level transition rates and 

chemical reaction rates have been shown to agree well with rates from the literature.
74 

Section 7 

References

1. Bose, D., et al., Uncertainty Assessment of Hypersonic Aerothermodynamics 

Prediction Capability, in 42nd AIAA Thermophysics Conference2011, AIAA Paper 

No. 2011-3141: Honolulu, HI.

2. Bird, G.A., Molecular Dynamics and the Direct Simulation of Gas Flows. 1994, 

Oxford, UK: Oxford University Press.

3. Ozawa, T., et al., Quasi-Classical Trajectory Modeling of OH Production in Direct 

Simulation Monte Carlo. Journal of Thermophysics and Heat Transfer, 2005. 19(2): 

p. 235–244.

4. Farbar, E.D., Kinetic Simulation of Rarefied and Weakly Ionized Hypersonic Fow 

Fields, in Aerospace Engineering2010, The University of Michigan: Ann Arbor, 

MI. p. 170.

5. Farbar, E.D. and I.D. Boyd, Simulation of Reactions Involving Charged Particles in 

Hypersonic Rarefied Flows, 2009, AIAA Paper No. 2009-267.

6. Farbar, E.D. and I.D. Boyd, Simulation of FIREII II Reentry Flow Using the Direct 

Simulation Monte Carlo Method, 2008, AIAA Paper No. 2008-4103.

7. Ozawa, T., et al., DSMC-CFD Comparison of a High Altitude, Extreme-Mach 

Number Reentry Flow, 2008, AIAA Paper No. 2008-1216.

8. Ozawa, T., J. Zhong, and D.A. Levin, Development of kinetic-based energy 

exchange models for noncontinuum, ionized hypersonic flows. Physics of Fluids, 

2008. 20.

9. Ozawa, T., et al., Modeling of the Stardust Reentry Flows with Ionization in DSMC, 

in 45th AIAA Aerospace Sciences Meeting and Exhibit2007, AIAA Paper No. 2007-

611: Reno, NV.

10. Bird, G.A. A Comparison of Collision Energy-based and Temperature-based 

Procedures in DSMC. in 26th Symposium on Rarefied Gas Dynamics. 2009. Kyoto, 

Japan: American Institue of Physics.

11. Bird, G.A. Chemical Reactions in DSMC. in 27th Symposium on Rarefied Gas 

Dynamics. 2010. Pacific Grove, CA: American Institute of Physics.

12. Bird, G.A. The Quantum-Kinetic Chemistry Model. in 27th International 

Symposium on Rarefied Gas Dynamics. 2010. Pacific Grove, CA.

13. Bird, G.A., The Q-K Model for Gas-Phase Chemical Reaction Rates. Physics of 

Fluids, 2011. 23(10): p. 106101.

14. Gallis, M.A., R.B. Bond, and J.R. Torczynski, A kinetic-theory approach for 

computing chemical-reaction rates in upper-atmosphere hypersonic flows. The 

Journal of Chemical Physics, 2009. 131(12): p. 124311.

15. Bird, G.A., Nonequilibrium Radiation During Re-Entry at 10 km/s, 1987, AIAA 

Paper No. 1987-1543.
75 

16. Gallis, M.A. and J.K. Harvey, Nonequilibrium Thermal Radiation from Air Shock 

Layers Modeled with Direct Simulation Monte Carlo. Journal of Thermophysics 

and Heat Transfer, 1994. 8(4): p. 765–772.

17. Carlson, A.B. and H.A. Hassan, Radiation Modeling with Direct Simulation Monte 

Carlo. Journal of Thermophysics and Heat Transfer, 1992. 6: p. 631–636.

18. Ozawa, T., et al., Development of Coupled Particle Hypersonic Flowfield-Photon 

Monte Carlo Radiation Methods. Journal of Thermophysics and Heat Transfer, 

2010. 24(3): p. 612–622.

19. Johnston, C.O., B.R. Hollis, and K. Sutton, Spectrum Modeling for Air Shock-Layer 

Radiation at Lunar-Return Conditions. Journal of Spacecraft and Rockets, 2008.

45(5): p. 865–878.

20. Johnston, C.O., B.R. Hollis, and K. Sutton, Non-Boltzmann Modeling for Air 

Shock-Layer Radiation at Lunar-Return Conditions. Journal of Spacecraft and 

Rockets, 2008. 45(5): p. 879–890.

21. Cornette, E.S., Forebody Temperature and Calorimeter Heating Rates Measured 

During Project Fire II Reentry at 11.35 km/s, 1966, NASA TM X-1305.

22. Birdsall, C.K. and A.B. Langdon, Plasma Physics via Computer Simulation. 1985, 

New York: McGraw-Hill.

23. Hockney, R.W. and J.W. Eastwood, Computer Simulation Using Particles. 1988, 

Bristol, U.K.: Inst. of Phys. Publishing.

24. Bird, G.A., Molecular Gas Dynamics. 1976, Oxford, U.K.: Clarendon Press.

25. Taylor, J.C., A.B. Carlson, and H.A. Hassan, Monte Carlo Simulation of Radiating 

Re-Entry Flow. Journal of Thermophysics and Heat Transfer, 1994. 8(3): p. 478– 485.

26. Bird, G.A., Low Density Aerothermodynamics, 1985, AIAA Paper No. 1985-0994.

27. Bird, G.A., Direct Simulation of Typical AOTV Entry Flows, 1986, AIAA Paper 

No. 1986-1310.

28. Gallis, M.A. and J.K. Harvey, Ionization Reactions and Electric Fields in Plane 

Hypersonic Shock Waves. Progess in Astronautics and Aeronautics, 1992. 160: p. 

234–244.

29. Carlson, A.B. and H.A. Hassan, Direct Simulation of Reentry Flows with 

Ionization. Journal of Thermophysics and Heat Transfer, 1992. 6(3): p. 400–404.

30. Boyd, I.D., Monte Carlo Simulation of Nonequilibrium Flow in a Low-power 

Hydrogen Arcjet. Physcis of Fluids, 1997. 9(10): p. 4575–4584.

31. Nanbu, K., Probability Theory of Electron-Molecule, Ion-Molecule, MoleculeMolecule, and Coulomb Collisions for Particle Modeling of Materials processing 

Plasmas and Gases. IEEE Transactions of Plasma Science, 2000. 28(3): p. 971–

990.
76 

32. Nanbu, K., Stochasitc Theory of Motion and Collision of Charged Particle in a 

Uniform Electric Field. Journal of the Physical Society of Japan, 1994. 63: p. 979–

983.

33. Nanbu, K., Direct Simulation Scheme Derived from the Boltzmann Equation.-II. 

Multicomponent Gas Mixtures. Journal of the Physical Society of Japan, 1980. 49:

p. 2042–2049.

34. Nanbu, K., Direct Simulation Scheme derived from the Boltzmann Equation-I: 

Monocomponent Gases. Journal of the Physical Society of Japan, 1980. 49: p. 

2042–2049.

35. Babovsky, H. and R. Illner, A Convergence Proof for Nanbu's Simulation Method 

for the Full Boltzmann Equation. SIAM J. Num. Anal., 1989. 26: p. 45–65.

36. Wilmoth, R.G., G.J. LeBeau, and A.B. Carlson, DSMC Grid Methodologies for 

Computing Low-Density, Hypersonic Flows About Reusable Launch Vehicles, 1996, AIAA 1996-1812.

37. LeBeau, G.J. and F.E. Lumpkin III, Application Highlights of the DSMC Analysis 

Code (DAC) Software for Simulating Rarefied Flows. Computer Methods in 

Applied Mechanics and Engineering, 2001. 191(6-7): p. 595–609.

38. LeBeau, G.J., K.A. Boyles, and F.E. Lumpkin III, Virtual Sub-Cells for the Direct 

Simulation Monte Carlo Method, 2003, AIAA 2003-1031.

39. Gallis, M. and J.R. Torczynski, Effect of Collision-Partner Selection Schemes on 

the Accuracy and Efficiency of the Direct Simulation Monte Carlo Method.

International Journal for Numerical Methods in Fluids, 2010.

40. Burt, J.M., E. Josyula, and I.D. Boyd, Techniques for Reducing Collision 

Separation in Direct Simulation Monte Carlo Calculations, in 42nd AIAA 

Thermophysics Conference2011, AIAA Paper No. 2011-3312: Honolulu, HI.

41. Bird, G.A., Monte Carlo simulation in an engineering context. Progess in 

Astronautics and Aeronautics, 1981. 74(1).

42. Borgnakke, C. and P.S. Larsen, Statistical collision model for Monte Carlo 

simulation of polyatomic gas mixtures. Journal of Computational Physics, 1975. 18.

43. Bergemann, F. and I.D. Boyd, New Discrete Vibrational Energy Model for Direct 

Simulation Monte Carlo Method. Progess in Astronautics and Aeronautics, 1994. 

160: p. 174–183.

44. Haas, B.L., Thermochemistry Models Applicable to a Vectorized Particle 

Simulation, in Department fo Aeronautics and Astronautics1990, Stanford 

University.

45. Lumpkin III, F.E., B.L. Haas, and I.D. Boyd, Resolution of Differences Between 

Collision Number Definitions in Particle and Continuum Simulations. Physics of 

Fluids A, 1991. 3: p. 2282.
77 

46. Haas, B.L., et al., Rates of Thermal Relaxation in Direct Simulation Monte Carlo 

Methods. Physics of Fluids, 1994. 6(6): p. 2191–2201.

47. Jeans, J.H., The Dynamical Theory of Gases. 2 ed. 1916, London: Cambridge 

University Press.

48. Landau, L. and E. Teller, Zur Theorie der Schalldispersion. Phys. Z. Sowjetunion, 

1936. 10: p. 14.

49. Gimelshein, S.F., I.D. Boyd, and M.S. Ivanov, Modeling of Internal Energy 

Transfer in Plume Flows of Polyatomic Molecules by the DSMC Method, 1999, 

AIAA Paper No. 1999-738.

50. Gimelshein, N.E., S.F. Gimelshein, and D.A. Levin, Vibrational Relaxation Rates 

in the Direct Simulation Monte Carlo Method. Physics of Fluids, 2002. 14(12): p. 

4452–4455.

51. Boyd, I.D., Relaxation of Discrete Rotational Energy Distributions Using a Monte 

Carlo Method. Physics of Fluids A, 1993. 5(9): p. 2278–2286.

52. Parker, J.G., Rotational and Vibrational Relaxation in Diatomic Gases. Physics of 

Fluids, 1959. 2(4): p. 449–462.

53. Millikan, R.C. and D.R. White, Systematics of Vibrational Relaxation. Journal of 

Chemical Physics, 1963. 39(12): p. 3209–3213.

54. Gallis, M.A. and J.K. Harvey, Atomic Species Radiation from Air Modeled with 

Direct Simulation Monte Carlo Method. Journal of Thermophysics and Heat 

Transfer, 1995. 9(3): p. 456–463.

55. Bird, G.A. Simulation of Multi-Dimensional and Chemically Reacting Flows. in 

11th International Symposium on Rarefied Gas Dynamics. 1978. Cannes, France: 

Commissariat a Lenergie Atomique.

56. Tonks, L. and I. Langmuir, General Theory of a Plasma Arc. Physical Review, 

1929. 34(1): p. 876–922.

57. Bartel, T.J. and C.R. Justiz, DSMC Simulation of Ionized Rarefied Flows, in 24th 

AIAA Fluid Dynamics Conference1993, AIAA Paper 1993-3095: Orlando, FL.

58. Justiz, C.R. and C. Dalton, A Hybrid Flow Model for Charged and Neutral 

Particles Around Spacecraft in Low Earth Orbit, in 27th AIAA Thermophysics 

Conference1992, AIAA Paper 1992-2935: Nashville, TN.

59. Gallis, M., R. Prasad, and J.K. Harvey, The Effect of Plasmas on the Aerodynamic

Performance of Vehicles, 1998, AIAA 1998-2666.

60. Park, C., Nonequilibrium Hypersonic Aerothermodynamics. 1990, New York: 

Wiley.

61. Lord, R.G., Modeling Vibrational Energy Exchange of Diatomic Molecules Using 

the Morse Interatomic Potential. Physics of Fluids, 1998. 10(3): p. 4.
78 

62. Koura, K. and H. Matsumoto, Variable soft sphere molecular model for inversepower-law or Lennard-Jones potential. Physics of Fluids A, 1991. 3.

63. Park, C., R.L. Jaffe, and H. Partridge, Chemical-Kinetic Parameters of Hyperbolic 

Earth Entry. Journal of Thermophysics and Heat Transfer, 2001. 15(1): p. 76–90.

64. Kang, S.W. and M.G. Dunn, Theoretical and Measured Electron-Density 

Distributions for the RAM Vehicle at High Altitudes, 1972, AIAA Paper No. 1972- 689.

65. Gimelshein, N.E. and D.A. Levin, Modeling of Glow Radiation in the Rarefied 

Flow About an Orbiting Spacecraft. Journal of Thermophysics and Heat Transfer, 

2000. 14(4): p. 471–479.

66. Johnston, H.S., Gas Phase Reaction Kinetics of Neutral Oxygen Species, U.S.D.o. 

Commerce, Editor 1968, National Bureau of Standards: Washington, D.C.

67. Boyd, I.D. and T. Gokcen, Computation of Axisymmetric and Ionized Hypersonic 

Flows Using Particle and Continuum Methods. AIAA Journal, 1994. 32(9): p. 

1828–1835.

68. Bourdon, A. and P. Vervisch, Study of a Low-Pressure Nitrogen Plasma Boundary 

Layer Over a Metallic Plate. Physics of Plasmas, 1997. 4: p. 4144–4157.

69. Park, C., Chemical-Kinetic Problems of Future NASA Missions, 1991, AIAA Paper 

No. 1991-464.

70. Losev, S.A., et al., Thermochemical Nonequilibrium Kinetic Models in Strong 

Shock Waves on Air, 1994, AIAA Paper No. 1994-1990.

71. Fujita, K., T. Yamada, and N. Ishii, Impacts of Ablation Gas Kinetics on Hyperbolic 

Entry Radiative Heating, in 44th AIAA Aerospace Sciences Meeting2006, AIAA 

Paper No. 2006-1185: Reno, NV.

72. Gupta, R.N., Navier-Stokes and Viscous Shock-Layer Solutions for Radiating 

Hypersonic Flows, 1987, AIAA Paper No. 1987-1576.

73. Phelps, A.V., Cross Sections and Swarm Coefficients for Nitrogen Ions and 

Neutrals in N2 and Argon Ions and Neutrals in Ar for Energies from 0.1 eV to 10 

keV. J. Phys. Chem. Ref. Data, 1991. 20: p. 557.

74. Park, C., Review of Chemical-Kinetic Problems of Future NASA Missions, I: Earth 

Entries. Journal of Thermophysics and Heat Transfer, 1993. 7: p. 385–398.

75. Park, C., J.T. Howe, and R.L. Jaffe, Review of Chemical-Kinetic Problems of 

Future NASA Missions, II: Mars Entries. Journal of Thermophysics and Heat 

Transfer, 1994. 8(1): p. 9–23.

76. Gupta, R.N., et al., A Review of Reaction Rates and Thermodynamic and Transport 

Properties for an 11-Species Air Model for Chemical and Thermal Nonequilibrium 

Calculations to 30,000 K, 1990, NASA RP-1232.

77. Boyd, I.D., Modeling of Plasma Formation in Rarefied Hypersonic Entry Flows, 

2007, AIAA Paper No. 2007-206.
79 

78. Carlson, A.B., Direct Simulation Monte Carlo with Ionization and Radiation, in 

Department of Mechanical and Aerospace Engineering1990, North Carolina Sate 

University: Raleigh. p. 82.

79. Tayal, S.S., Effective Collision Strengths for Electron Impact of N I. Atomic Data 

and Nuclear Data Tables, 2000. 76: p. 191–212.

80. Frost, R.M., et al., Calculated cross sections and measured rate coefficients for 

electron-impact excitation of neutral and singly ionized nitrogen. Journal of 

Applied Physics, 1998. 84(6): p. 2989–3003.

81. Zatsarinny, O. and S.S. Tayal, Electron Collisional Excitation Rates for O I using 

the B-Spline R-Matrix Approach. The Astrophysical Journal Supplement Series, 

2003. 148: p. 575–582.

82. Bhatia, A.K. and S.O. Kastner, The Neutral Oxygen Spectrum. I. Collisionally 

Excited Level Populations and Line Intensities under Optically Thin Conditions.

The Astrophysical Journal Supplement Series, 1995. 96: p. 325–341.

83. Chernyi, G.G. and S.A. Losev, Development of Thermal Protection Systems for 

Interplanetary Flight, R.I.o. Mechanics, Editor 1999, International Science and 

Technology Center: Moscow.

84. Telulet, P., J.P. Sarrette, and A.M. Gomes, Calculation of Electron-Impact Inelastic 

Cross-Sections and Rate Coefficients for Diatomic Molecules. Application to Air 

Molecules. Journal of Qualitative Spectroscopy and Radiative Transfer, 1999. 62: p. 

549–569.

85. Gorelov, V.A., et al., Experimental and Numerical Study of Nonequilibrium 

Ultraviolet NO and N2+ Emission in a Shock Layer. Journal of Thermophysics and 

Heat Transfer, 1998. 12(2): p. 172–179.

86. Johnston, C.O., Nonequilibruim Shock-Layer Radiative Heating for Earth and Titan 

Entry, in Aerospace Engineering2006, Virginia Polytechnic Institute and State 

University: Blacksburg, VA. p. 245.

87. Gorelov, V.A., et al. Nonequilibrium Ionization and Radiation Behind Shock Wave 

in Martian Atmosphere. in Third European Symposium on Aerothermodynamics for 

Space Vehicles. 1998. Noordwijk, The Netherlands: European Space Agency.

88. Park, C., Rate Parameters for Electronic Excitation of Diatomic Molecules II. 

Heavy Particle-Impact Processes, 2008, AIAA Paper No. 2008-1446.

89. Flagan, R.C. and J.P. Appleton, Excitation Mechanisms of the Nitrogen First- Positive and First-Negative Radiation at High-Temperature. Journal of Chemical 

Physics, 1972. 56(3): p. 1163–1173.

90. Nagy, O., Excitation Cross-Sections of N2+ Molecular Ion by Electron Impact and 

the Vibrational Energy Levels of the Three Target States. Chemical Physics, 2003. 

286(1): p. 109–114.
80 

91. Guerra, V. and J. Loureiro, Electron and heavy particle kinetics in a low-pressure 

nitrogen glow discharge. Plasma Sources Science & Technology, 1997. 6: p. 361–

372.

92. Raynaud, E., et al., Huygens Aerothermal Environment: Radiative Heating, in 3rd 

International Planetary Probe Workshop2005, ESA: Anavyssos, Greece.

93. Magin, T.E., et al., Nonequilibrium Radiation Modeling for Huygens Entry, in 3rd 

International Planetary Probe Workshop2005: Anavyssos, Greece.

94. Cheatwood, F.M. and P. Gnoffo, User's Manual for the Langley 

Aerothermodynamic Upwind Relaxation Algorithm (LAURA), April 1996, NASA 

TM 4674.

95. Gnoffo, P., R.N. Gupta, and J.L. Shinn. Conservation Equations and Physical 

Models for Hypersonic Air Flows in Thermal and Chemical Nonequilibrium. in 

NASA TP 2867. Feb. 1989.

96. Roe, P.L., Approximate Riemann Solvers, Parameter Vectors, and Difference 

Schemes. Journal of Computational Physics, 1981. 43(2): p. 357–372.

97. Yee, H.C., On Symmetric and Upwind TVD Schemes, in NASA TM 883251986.

98. Mcbride, B.J., M.J. Zehe, and S. Gordon, NASA Glenn Coefficients for Calculating 

Thermodynamic Properties of Individual Species, Sept. 2002, NASA TP 2002- 211556.

99. Sutton, K. and P. Gnoffo, Multi-Component Diffusion with Application to 

Computational Aerothermodynamics, 1998, AIAA Paper 98-2575.

100. Ralchenko, Y. and e. al. NIST Atomic Spectra Database, Version 3.1.0. July 2006; 

Available from: http://physics.nist.gov/PhysRefData/ASD/index.html.

101. Team, T.O.P., The Opacity Project, Vol. 1. 1995, Bristol and Philedelphia: Institute 

of Physics Publishing.

102. Cunto, W. and e. al, TOPbase at the CDS. Astronomy and Astrophysics, 1993. 275: 

p. L5–L8.

103. Soon, W.H. and J.A. Kunc, Nitrogen Plasma Continuum Emission Associated with 

N-(3P) and N-(1D) Ions. Physical Review A, 1990. 41: p. 4531–4533.

104. Chauveau, S., et al., Radiative Transfer in LTE Air Plasmas for Temperatures up to 

15,000 K. Journal of Qualitative Spectroscopy and Radiative Transfer, 2003. 77: p. 

113–130.

105. Chambers, L.H., Predicting Radiative Heat Transfer in Thermochemical 

Nonequilibrium Flow Fields, 1994, NASA TM-4564.

106. Laux, C.O., Optical Diagnostics and Radiative Emission of Air Plasmas, 1993, 

High Temperature Gas Dynamics Lab, Mechanical Engineering Dept., Rept. T-288, 

Stanford University.
81 

107. Whang, T.J., et al., Franck-Condon Factors of the ..... Transitions of N2. Journal of 

Quantitative Spectroscopy and Radiative Transfer, 1996. 55: p. 335–344.

108. Stahel, D., M. Leoni, and K. Dresslar, Nonadiabatic Representations of the 1Σu and 

1Πu States of the N2 Molecule. Journal of Chemical Physics, Sept. 1983. 79: p. 

2541–2558.

109. Chauveau, S., et al., Contributions of Diatomic Molecular Electronic Systems to 

Heated Air Radiation. Journal of Quantitative Spectroscopy and Radiative Transfer, 

2002. 72: p. 503v530.

110. Park, C., Radiation Enhancement by Nonequilibrium in Earth's Atmosphere.

Journal of Spacecraft and Rockets, 1985. 22: p. 27–36.

111. Hash, D., et al., FIRE II Calculations for Hypersonic Nonequilibrium 

Aerothermodynamics Code Verification: DPLR, LAURA and US3D, 2007, AIAA 

2007-605.
82 

Appendices

Appendix A: 

Species and Chemistry Data

Table A.1: Parameters used in the VHS molecular model and characteristic 

temperatures.

Species Tref (K) dref (m) ω Θrot (K) Θvib (K) Θdiss (K) Θion (K)

O2 288.0 3.96E-10 0.77 2.080 2256.0 59500.0 140066.5

N2 288.0 4.07E-10 0.74 2.874 3371.0 11355.0 180798.3

O 288.0 3.00E-10 0.75 - - - 158053.5

N 288.0 3.00E-10 0.75 - - - 168613.6

NO 288.0 4.00E-10 0.79 2.080 2719.0 75500.0 107457.8

O2

+ 288.0 3.96E-10 0.77 2.080 2256.0 59500.0 - 

N2

+ 288.0 4.07E-10 0.74 2.874 3371.0 11355.0 - 

O+ 288.0 3.00E-10 0.75 - - - - 

N+ 288.0 3.00E-10 0.75 - - - - 

NO+ 288.0 4.00E-10 0.79 2.080 2719.0 75500.0 - 

e- 288.0 1.00E-10 0.70 - - - - 

Table A.2: Internal energy relaxation parameters.

Species ZrZot

C1 C2

O2 1 56.5 153.5

N2 1 9.1 220

NO 1 9.1 220

O2

+ 1 56.5 153.5

N2

+ 1 9.1 220

NO+ 1 9.1 220

Table A.3: Parameters used in new vibrational relaxation model.

O2 N2 NO

Zref 17900.0 52500.0 123000.0

Tref (K) 2256.0 3371.0 2791.0

Θd (K) 59500.0 113500.0 75500.0

Table A.4: Electronic energy levels of O2.

Level g Eel (J) State

1 3 0.0000E+01

2 2 1.5728E-19

Zvi

b = C1 Tω ( )expC 2

T −1

3 ( )
83 

Level g Eel (J) State

3 1 2.6211E-19

4 1 6.5663E-19

5 6 6.8912E-19

6 3 7.0306E-19

7 3 9.9267E-19

Table A.5: Electronic energy levels of N2.

Level g Eel (J) State

1 1 0.0000E+01

2 3 9.9727E-19

3 6 1.1843E-18

4 6 1.1881E-18

5 3 1.3165E-18

6 1 1.3538E-18

7 2 1.3763E-18

8 2 1.4483E-18

9 5 1.5415E-18

10 1 1.6925E-18

11 6 1.7242E-18

12 6 1.7707E-18

13 10 1.8474E-18

14 6 1.9388E-18

15 6 2.0778E-18

Table A.6: Electronic energy levels of O.

Level g E (J) Core Config. Term

1 9 1.55E-21 2s2 2p4 3P 

2 5 3.15E-19 2s2 2p4 1D 

3 1 6.71E-19 2s2 2p4 1S 

4 5 1.47E-18 2s32p2 3s 5So 

5 3 1.53E-18 2s32p2 3s 3So 

6 15 1.72E-18 2s 3p 5P 

b 1

Σ

g

+
84 

Level g E (J) Core Config. Term

7 9 1.76E-18 2s32p2 3p 3P 

8 5 1.90E-18 2s32p2 4s 5So 

9 3 1.91E-18 2s32p2 4s 3So 

10 25 1.94E-18 2s32p2 3d 5Do 

11 15 1.94E-18 2s32p2 3d 3Do 

12 15 1.97E-18 2s32p2 4p 5P 

13 9 1.98E-18 2s32p2 4p 3P 

14 15 2.01E-18 2s32p2 3s 3Do 

15 5 2.03E-18 2s32p2 5s 5So 

16 3 2.03E-18 2s32p2 5s 3So 

17 5 2.04E-18 2s32p2 3s 1Do 

18 25 2.04E-18 2s32p2 4d 5Do 

19 15 2.04E-18 2s32p2 4d 3Do 

20 35 2.05E-18 2s32p2 4f 5F 

21 21 2.05E-18 2s32p2 4f 3F 

22 15 2.06E-18 2s32p2 5p 5P 

23 9 2.06E-18 2s32p2 5p 3P 

24 25 2.09E-18 2s32p2 5d 5Do 

25 15 2.09E-18 2s32p2 5d 3Do 

26 35 2.09E-18 2s32p2 5f 5F 

27 21 2.09E-18 2s32p2 5f 3F 

28 288 2.12E-18 2s32p2 n=6 - 

29 392 2.14E-18 2s32p2 n=7 - 

30 512 2.15E-18 2s32p2 n=8 - 

31 648 2.15E-18 2s32p2 n=9 - 

32 800 2.16E-18 2s32p2 n=10 - 

Table A.7: Electronic energy levels of N.

Level g Eel (J) Core Config. Term

1 4 0.0000E+01 2s2 2p3 4So 

2 10 3.8195E-19 2s2 2p3 2Do 

3 6 5.7287E-19 2s2 2p3 2Po 

4 12 1.6554E-18 2s22p2 3s 4P 
85 

Level g Eel (J) Core Config. Term

5 6 1.7122E-18 2s22p2 3s 2P 

6 12 1.7507E-18 2s 2p4 4P 

7 2 1.8589E-18 2s22p2 3p 2So 

8 20 1.8839E-18 2s22p2 3p 4Do 

9 12 1.8973E-18 2s22p2 3p 4Po 

10 4 1.9219E-18 2s22p2 3p 4So 

11 10 1.9235E-18 2s22p2 3p 2Do 

12 6 1.9426E-18 2s22p2 3p 2Po 

13 10 1.9798E-18 2s22p2 3s 2D 

14 12 2.0598E-18 2s22p2 4s 4P 

15 6 2.0698E-18 2s22p2 4s 2P 

16 6 2.0784E-18 2s22p2 3d 2P 

17 28 2.0802E-18 2s22p2 3d 4F 

18 14 2.0828E-18 2s22p2 3d 2F 

19 12 2.0828E-18 2s22p2 3d 4P 

20 20 2.0859E-18 2s22p2 3d 4D 

21 10 2.0884E-18 2s22p2 3d 2D 

22 2 2.1151E-18 2s22p2 4p 2So 

23 20 2.1220E-18 2s22p2 4p 4Do 

24 12 2.1258E-18 2s22p2 4p 4Po 

25 10 2.1300E-18 2s22p2 4p 2Do 

26 4 2.1343E-18 2s22p2 4p 4So 

27 6 2.1377E-18 2s22p2 4p 2Po 

28 90 2.1912E-18 2s22p2 4d2P,4d4F,4d4D,4d4F,4d4P,4d2D 

29 126 2.1946E-18 2s22p2 4f4D,4f4F,4f4G,4f2D,4f2F,4f2G 

30 450 2.2368E-18 2s22p2 n=5

31 648 2.2703E-18 2s22p2 n=6

32 822 2.2864E-18 2s22p2 n=7

33 1152 2.2968E-18 2s22p2 n=8

34 1458 2.3040E-18 2s22p2 n=9
86 

Level g Eel (J) Core Config. Term

35 1800 2.3091E-18 2s22p2 n=10

Table A.8: Electronic energy levels of NO.

Level g Eel (J) State

1 4 0.0000E+01

2 8 7.5485E-19

3 2 8.7218E-19

4 4 9.1115E-19

5 4 9.5349E-19

6 4 9.7336E-19

7 4 1.0343E-18

8 2 1.0533E-18

9 4 1.1979E-18

10 2 1.2032E-18

11 4 1.2269E-18

12 4 1.2401E-18

13 2 1.2410E-18

14 2 1.2485E-18

15 2 1.2515E-18

16 4 1.2580E-18

Table A.9: Electronic energy levels of O2

+

.

Level g Eel (J) State

1 4 0.0000E+01

2 8 6.5380E-19

3 4 8.0594E-19

4 6 8.0650E-19

5 4 8.6013E-19

6 2 9.2966E-19

7 4 9.8329E-19

8 4 1.0568E-18

X 2Πr
87 

Level g Eel (J) State

9 4 1.2177E-18

10 4 1.2276E-18

11 8 1.3011E-18

12 4 1.3111E-18

13 2 1.3243E-18

14 2 1.3786E-18

15 4 1.4302E-18

Table A.10: Electronic energy levels of N2

+

.

Level g Eel (J) State

1 2 0.0000E+01

2 4 1.8211E-19

3 2 5.0578E-19

4 4 5.0654E-19

5 8 8.2636E-19

6 8 9.1377E-19

7 4 1.0492E-18

8 4 1.0528E-18

9 4 1.1323E-18

10 4 1.1621E-18

11 8 1.1919E-18

12 8 1.2316E-18

13 4 1.2713E-18

14 4 1.2733E-18

15 2 1.2831E-18

16 2 1.3309E-18

17 4 1.4302E-18

Table A.11: Electronic energy levels of O+.

Level g Eel (J) Core Config. Term

1 4 0.0000E+01 2s2 2p3 4So 

2 10 5.3270E-19 2s2 2p3 2Do 

3 6 8.0386E-19 2s2 2p3 2Po 

D 2

Πu
88 

Table A.12: Electronic energy levels of N+.

Level g Eel (J) Core Config. Term

1 9 0.0000E+01 2s2 2p2 3P 

2 5 3.0424E-19 2s2 2p2 1D 

3 1 6.4929E-19 2s2 2p2 1S 

4 15 1.8324E-18 2s 2p3 3Do 

5 9 2.1696E-18 2s 2p3 3Po 

6 5 2.8642E-18 2s 2p3 1Do 

7 12 2.9610E-18 2s22p 3s 3Po 

Table A.13: Electronic energy levels of NO+.

Level g Eel (J) State

1 1 0.0000E+01

2 3 1.0367E-18

3 6 1.1771E-18

4 6 1.2293E-18

5 3 1.3457E-18

6 1 1.3814E-18

7 2 1.4194E-18

8 2 1.4595E-18

Table A.14: Endothermic chemical reaction rate coefficients (m3/molecule/s) 

used in the TCE model (the first rate listed for each reaction is the 

rate used in the simulations). kf = aTbexp(-Ea/kT)

Number Reaction a b Ea (J) Erx (J) Reference

1 O2 + O2 → O + O + O2 3.3211E-09 -1.50 1.5 8.1995E-19 -8.1995E-19 Park[63] 

4.5800E-11 -1.00 8.1970E-19 -8.1970E-19 Bird[15] 

5.9780E-12 -1.00 8.2150E-19 -8.2150E-19 Kang[64] 

5.3930E-11 -1.00 8.1970E-19 -8.1970E-19 Gimelshein[65] 

2 O2 + N2 → O + O + N2 3.3211E-09 -1.50 1.5 8.1995E-19 -8.1995E-19 Park[63] 

4.5800E-11 -1.00 8.1970E-19 -8.1970E-19 Bird[15] 

5.9780E-12 -1.00 8.2150E-19 -8.2150E-19 Kang[64] 

1.1980E-11 -1.00 8.1970E-19 -8.1970E-19 Gimelshein[65] 

3 O2 + O → O + O + O 1.6605E-08 -1.50 1.0 8.1995E-19 -8.1995E-19 Park[63] 

1.3750E-10 -1.00 8.1970E-19 -8.1970E-19 Bird[15] 

5.9780E-12 -1.00 8.2150E-19 -8.2150E-19 Kang[64] 

3.3200E-09 -1.50 8.2140E-19 -8.2140E-19 Ozawa[9] 

1.4980E-10 -1.00 8.1970E-19 -8.1970E-19 Gimelshein[65] 

4.5700E-11 -1.00 8.1970E-19 -8.1970E-19 Johnston[66] 

4 O2 + N → O + O + N 1.6605E-08 -1.50 1.0 8.1995E-19 -8.1995E-19 Park[63] 

1.3750E-10 -1.00 8.1970E-19 -8.1970E-19 Bird[15] 

5.9780E-12 -1.00 8.2150E-19 -8.2150E-19 Kang[64] 

3.3200E-09 -1.50 8.2140E-19 -8.2140E-19 Ozawa[9] 

5.9930E-12 -1.00 8.1970E-19 -8.1970E-19 Gimelshein[65] 

5 O2 + NO → O + O + NO 3.3211E-09 -1.50 1.0 8.1995E-19 -8.1995E-19 Park[63] 

ζ

X1

Σ+
89 

Number Reaction a b 

Ea (J)Erx (J) Reference

4.5800E-11-1.00 8.1970E-19-8.1970E-19 Bird[15] 

5.9780E-12-1.00 8.2150E-19-8.2150E-19 Kang[64] 

5.9930E-12-1.00 8.1970E-19-8.1970E-19 Gimelshein[65] 

6 O

2 + O2

+

→ O + O + O

2

+ 3.3211E

-09-1.50 1.5 8.1995E-19-8.1995E-19 Park[63] 

4.5800E-11-1.00 8.1970E-19-8.1970E-19 Bird[15] 

5.9780E-12-1.00 8.2150E-19-8.2150E-19 Kang[64] 

5.3930E-11-1.00 8.1970E-19-8.1970E-19 Gimelshein[65] 

7 O

2 + N2

+

→ O + O + N

2

+ 3.3211E

-09-1.50 1.5 8.1995E-19-8.1995E-19 Park[63] 

4.5800E-11-1.00 8.1970E-19-8.1970E-19 Bird[15] 

5.9780E-12-1.00 8.2150E-19-8.2150E-19 Kang[64] 

1.1980E-11-1.00 8.1970E-19-8.1970E-19 Gimelshein[65] 

8 O2 + O

+

→ O + O + O

+ 1.6605E

-08-1.50 1.0 8.1995E-19-8.1995E-19 Park[63] 

1.3750E-10-1.00 8.1970E-19-8.1970E-19 Bird[15] 

5.9780E-12-1.00 8.2150E-19-8.2150E-19 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

1.4980E-10-1.00 8.1970E-19-8.1970E-19 Gimelshein[65] 

4.5700E-11-1.00 8.1970E-19-8.1970E-19 Johnston[66] 

9 O

2 + N

+

→ O + O + N

+ 1.6605E

-08-1.50 1.0 8.1995E-19-8.1995E-19 Park[63] 

1.3750E-10-1.00 8.1970E-19-8.1970E-19 Bird[15] 

5.9780E-12-1.00 8.2150E-19-8.2150E-19 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

5.9930E-12-1.00 8.1970E-19-8.1970E-19 Gimelshein[65] 

10O

2 + NO

+

→ O + O + NO

+ 3.3211E

-09-1.50 1.0 8.1995E-19-8.1995E-19 Park[63] 

4.5800E-11-1.00 8.1970E-19-8.1970E-19 Bird[15] 

5.9780E-12-1.00 8.2150E-19-8.2150E-19 Kang[64] 

5.9930E-12-1.00 8.1970E-19-8.1970E-19 Gimelshein[65] 

11N

2 + O2

→ N + N + O

2 1.1620E

-08-1.60 0.5 1.5630E-18-1.5630E-18 Park[63] 

6.1700E-09-1.60 1.5610E-18-1.5610E-18 Bird[15] 

3.1550E-13-0.50 1.5600E-18-1.5600E-18 Kang[64] 

4.9800E-08-1.60 1.5630E-18-1.5630E-18 Ozawa[9] 

7.9700E-13-0.50 1.5610E-18-1.5610E-18 Boyd[67] 

3.1870E-13-0.50 1.5610E-18-1.5610E-18 Gimelshein[65] 

12N

2 + N2

→ N + N + N

2 1.1620E

-08-1.60 1.0 1.5630E-18-1.5630E-18 Park[63] 

6.1700E-09-1.60 1.5610E-18-1.5610E-18 Bird[15] 

3.1550E-13-0.50 1.5600E-18-1.5600E-18 Kang[64] 

1.1600E-08-1.60 1.5630E-18-1.5630E-18 Ozawa[9] 

7.9700E-13-0.50 1.5610E-18-1.5610E-18 Boyd[67] 

13N

2 + O 

→ N + N + O 4.9980E-08-1.60 0.5 1.5630E-18-1.5630E-18 Park[63] 

1.8500E-08-1.60 1.5610E-18-1.5610E-18 Bird[15] 

3.1550E-13-0.50 1.5600E-18-1.5600E-18 Kang[64] 

7.1400E-08-1.50 1.5630E-18-1.5630E-18 Boyd[67] 

3.1870E-13-0.50 1.5610E-18-1.5610E-18 Gimelshein[65] 

14N

2 + N 

→ N + N + N 4.9980E-08-1.60 1.0 1.5630E-18-1.5630E-18 Park[63] 

1.8500E-08-1.60 1.5610E-18-1.5610E-18 Bird[15] 

6.7830E-08-1.50 1.5600E-18-1.5600E-18 Kang[64] 

7.1400E-08-1.50 1.5630E-18-1.5630E-18 Boyd[67] 

6.9000E-08-1.50 1.5610E-18-1.5610E-18 Gimelshein[65] 

15N

2 + NO 

→ N + N + NO 1.1620E-08-1.60 0.5 1.5630E-18-1.5630E-18 Park[63] 

6.1700E-09-1.60 1.5610E-18-1.5610E-18 Bird[15] 

3.1550E-13-0.50 1.5600E-18-1.5600E-18 Kang[64] 

4.9800E-08-1.60 1.5630E-18-1.5630E-18 Ozawa[9] 

7.9700E-13-0.50 1.5610E-18-1.5610E-18 Boyd[67] 

3.1870E-13-0.50 1.5610E-18-1.5610E-18 Gimelshein[65] 

16N

2 + O2

+

→ N + N + O

2

+ 1.1620E

-08-1.60 0.5 1.5630E-18-1.5630E-18 Park[63] 

6.1700E-09-1.60 1.5610E-18-1.5610E-18 Bird[15] 

3.1550E-13-0.50 1.5600E-18-1.5600E-18 Kang[64] 

4.9800E-08-1.60 1.5630E-18-1.5630E-18 Ozawa[9] 

ζ
90 

Number Reaction a b 

Ea (J)Erx (J) Reference

7.9700E-13-0.50 1.5610E-18-1.5610E-18 Boyd[67] 

3.1870E-13-0.50 1.5610E-18-1.5610E-18 Gimelshein[65] 

17N

2 + N2

+

→ N + N + N

2

+ 1.1620E

-08-1.60 1.0 1.5630E-18-1.5630E-18 Park[63] 

6.1700E-09-1.60 1.5610E-18-1.5610E-18 Bird[15] 

3.1550E-13-0.50 1.5600E-18-1.5600E-18 Kang[64] 

1.1600E-08-1.60 1.5630E-18-1.5630E-18 Ozawa[9] 

7.9700E-13-0.50 1.5610E-18-1.5610E-18 Boyd[67] 

18N

2 + O

+

→ N + N + O

+ 4.9980E

-08-1.60 0.5 1.5630E-18-1.5630E-18 Park[63] 

1.8500E-08-1.60 1.5610E-18-1.5610E-18 Bird[15] 

3.1550E-13-0.50 1.5600E-18-1.5600E-18 Kang[64] 

7.1400E-08-1.50 1.5630E-18-1.5630E-18 Boyd[67] 

3.1870E-13-0.50 1.5610E-18-1.5610E-18 Gimelshein[65] 

19N

2 + N

+

→ N + N + N

+ 4.9980E

-08-1.60 1.0 1.5630E-18-1.5630E-18 Park[63] 

1.8500E-08-1.60 1.5610E-18-1.5610E-18 Bird[15] 

6.7830E-08-1.50 1.5600E-18-1.5600E-18 Kang[64] 

7.1400E-08-1.50 1.5630E-18-1.5630E-18 Boyd[67] 

6.9000E-08-1.50 1.5610E-18-1.5610E-18 Gimelshein[65] 

20N

2 + NO

+

→ N + N + NO

+ 1.1620E

-08-1.60 0.5 1.5630E-18-1.5630E-18 Park[63] 

6.1700E-09-1.60 1.5610E-18-1.5610E-18 Bird[15] 

3.1550E-13-0.50 1.5600E-18-1.5600E-18 Kang[64] 

4.9800E-08-1.60 1.5630E-18-1.5630E-18 Ozawa[9] 

7.9700E-13-0.50 1.5610E-18-1.5610E-18 Boyd[67] 

3.1870E-13-0.50 1.5610E-18-1.5610E-18 Gimelshein[65] 

21N

2 + e

-

→ N + N + e

- 9.9670E

-27 2.60 0.5 1.5630E-18-1.5630E-18 Bourdon[68] 

4.9800E-06-1.60 1.5610E-18-1.5630E-18 Ozawa[9] 

1.4900E-05-1.60 1.5610E-18-1.5610E-18 Park[69] 

1.3450E-10-1.28 1.5630E-18-1.5630E-18 Losev[70] 

22 NO + O

2

→ N + O + O

2 3.8300E

-13-0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

8.3000E-15 0.00 1.0420E-18-1.0420E-18 Park[69] 

6.5900E-10-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

23 NO + N

2

→ N + O + N

2 3.8300E

-13-0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

8.3000E-15 0.00 1.0420E-18-1.0420E-18 Park[69] 

6.5900E-10-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

24 NO + O → N + O + O 7.6600E-13-0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

8.3000E-15 0.00 1.0420E-18-1.0420E-18 Ozawa[9] 

1.8300E-13 0.00 1.0420E-18-1.0420E-18 Park[69] 

1.3180E-08-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

25 NO + N → N + O + N 7.6600E-13 -0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

1.8300E-13 0.00 1.0420E-18-1.0420E-18 Park[69] 

1.3180E-08-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

26 NO + NO → N + O + NO 3.8300E-13-0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

8.3000E-15 0.00 1.0420E-18-1.0420E-18 Park[69] 

1.3180E-08-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

27 NO + O

2

+

→ N + O + O

2

+ 3.8300E

-13-0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

3.8300E-13-0.50 1.0430E-18-1.0430E-18 Bird[15] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

8.3000E-15 0.00 1.0420E-18-1.0420E-18 Park[69] 

ζ
91 

Number Reaction a b 

Ea (J)Erx (J) Reference

6.5900E-10-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

28 NO + N

2

+

→ N + O + N

2

+ 3.8300E

-13-0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

3.8300E-13-0.50 1.0430E-18-1.0430E-18 Bird[15] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

8.3000E-15 0.00 1.0420E-18-1.0420E-18 Park[69] 

6.5900E-10-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

29 NO + O

+

→ N + O + O

+ 7.6600E

-13-0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

7.6600E-13-0.50 1.0430E-18-1.0430E-18 Bird[15] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

8.3000E-15 0.00 1.0420E-18-1.0420E-18 Ozawa[9] 

1.8300E-13 0.00 1.0420E-18-1.0420E-18 Park[69] 

1.3180E-08-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

30 NO + N

+

→ N + O + N

+ 7.6600E

-13-0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

7.6600E-13-0.50 1.0430E-18-1.0430E-18 Bird[15] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

1.8300E-13 0.00 1.0420E-18-1.0420E-18 Park[69] 

1.3180E-08-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

31 NO + NO

+

→ N + O + NO

+ 3.8300E

-13-0.50 1.0 1.0420E-18-1.0420E-18 Fujita[71] 

3.8300E-13-0.50 1.0430E-18-1.0430E-18 Bird[15] 

6.4760E-10-1.50 1.0420E-18-1.0420E-18 Kang[64] 

3.3200E-09-1.50 8.2140E-19-8.2140E-19 Ozawa[9] 

8.3000E-15 0.00 1.0420E-18-1.0420E-18 Park[69] 

1.3180E-08-1.50 1.0430E-18-1.0430E-18 Gimelshein[65] 

32 N + e

-

→ N

+ + e-

+ e

- 5.8100E

-08-1.00 0.0 2.3220E-18-2.3220E-18 Park[63] 

4.1500E+03-3.82 2.3300E-18-2.3300E-18 Park[60] 

9.0000E-11-1.00 1.5000E-18-1.5000E-18 Ozawa[7] 

2.9890E-17 0.60 2.3300E-18-2.3300E-18 Losev[70] 

5.8100E-08-1.00 2.3300E-18-2.3300E-18 Boyd[67] 

1.0000E-14 0.00 2.3300E-18-2.3300E-18 Bird[15] 

1.8270E+02-3.14 2.3330E-18-2.3300E-18 Kang[64] 

8.4340E-14 0.00 2.3300E-18-2.3300E-18 Farbar[5] 

33 O + e

-

→ O

+ + e-

+ e

- 1.5900E

-08-1.00 0.0 2.1880E-18-2.1880E-18 Park[63] 

6.4800E+03-3.78 2.1900E-18-2.1900E-18 Park[60] 

3.0000E-10-1.00 1.4000E-18-1.4000E-18 Ozawa[7] 

8.6350E-18 0.68 2.1810E-18-2.1810E-18 Losev[70] 

1.5900E-08-1.00 2.1810E-18-2.1810E-18 Boyd[67] 

3.0000E-12 0.00 2.1810E-18-2.1810E-18 Bird[15] 

5.9780E+01-2.91 2.1810E-18-2.1810E-18 Kang[64] 

1.0540E-14 0.00 2.1810E-18-2.1810E-18 Farbar[5] 

34N

2 + O 

→ NO + N 9.9630E-17 0.10 0.0 5.2470E-19-5.2470E-19 Fujita[71] 

5.3000E-17 0.10 5.1750E-19-5.1750E-19 Bird[15] 

1.1620E-16 0.00 5.2460E-19-5.2460E-19 Kang[64] 

9.4500E-18 0.42 5.9280E-19-5.9280E-19 Ozawa[9] 

1.0600E-12-1.00 5.1750E-19-5.1750E-19 Park[69] 

1.1200E-16 0.00 5.1750E-19-5.1750E-19 Gimelshein[65] 

35 NO + O → O

2 + N 1.3950E

-17 0.00 0.0 2.6790E-19-2.6790E-19 Park[63] 

3.6000E-22 1.29 2.2380E-19-2.2380E-19 Bird[15] 

5.3140E-21 1.00 2.7200E-19-2.7200E-19 Kang[64] 

5.2790E-21 1.00 2.7190E-19-2.7190E-19 Gimelshein[65] 

36 N + N → N

2

+ + e- 7.3060E

-23 1.50 0.0 9.3190E-19-9.3190E-19 Park[63] 

3.3870E-17 0.00 9.3200E-19-9.3200E-19 Park[60] 

2.9800E-20 0.77 9.3400E-19-9.3400E-19 Bird[15] 

2.3250E-17 0.00 9.3610E-19-9.3610E-19 Kang[64] 

37 N + O → NO

+ + e- 8.8010E

-18 0.00 0.0 4.4040E-19-4.4040E-19 Park[63] 

8.7660E-18 0.00 4.4000E-19-4.4000E-19 Park[60] 

ζ
92 

Number Reaction a b Ea (J) Erx (J) Reference

1.4600E-21 1.00 4.4000E-19 -4.4000E-19 Park[69] 

2.5500E-20 0.37 4.4220E-19 -4.4220E-19 Bird[15] 

2.3250E-24 1.50 4.4040E-19 -4.4040E-19 Kang[64] 

1.4990E-20 0.50 4.4730E-19 -4.4730E-19 Gupta[72] 

38 O + O → O2

+ + e- 1.1790E-27 2.70 0.0 1.1130E-18 -1.1130E-18 Fujita[71] 

1.8300E-17 0.00 1.1100E-18 -1.1100E-18 Park[60] 

6.4200E-22 0.49 1.1200E-18 -1.1200E-18 Bird[15] 

2.6570E-13 -0.98 1.1160E-18 -1.1160E-18 Kang[64] 

39 N+ + N2 → N2

+ + N 1.6610E-18 0.50 1.0 1.6840E-19 -1.6840E-18 Park[63] 

1.6700E-17 -0.18 1.6700E-19 -1.6700E-19 Bird[15] 

2.3250E-26 2.10 1.6700E-19 -1.6700E-19 Phelps[73] 

3.3540E-19 0.81 1.7950E-19 -1.7950E-19 Kang[64] 

40 O+ + N2 → N2

+ + O 1.5110E-18 0.36 1.0 3.1480E-19 -3.1480E-19 Park[74] 

1.4940E-18 0.36 3.1480E-19 -3.1480E-19 Losev[70] 

1.0550E-16 -0.21 3.0600E-19 -3.0600E-19 Bird[15] 

5.6460E-11 -2.00 3.1750E-19 -3.1750E-19 Kang[64] 

41 O+ + NO → N+ + O2 2.3240E-25 1.90 1.0 3.6730E-19 -3.6730E-19 Park[74] 

42 NO+ + N → N2

+ + O 1.1960E-16 0.00 1.0 4.9010E-19 -4.1090E-19 Park[74] 

2.8300E-17 0.40 4.9010E-19 -4.9010E-19 Bird[15] 

43 NO+ + N → O+ + N2 5.6460E-17 -1.08 1.0 1.7670E-19 -1.7670E-19 Park[74] 

44 NO+ + O → N+ + O2 1.6610E-18 0.50 1.0 1.0660E-18 -1.0660E-18 Park[74] 

2.2250E-17 0.31 1.0670E-18 -1.0670E-18 Kang[64] 

45 NO+ + O → O2

+ + N 1.1960E-17 0.29 1.0 6.7100E-19 -6.7100E-19 Fujita[71] 

46 NO+ + O2 → O2

+ + NO 3.9850E-17 0.41 1.5 5.4010E-19 -5.4010E-19 Park[74] 

3.9850E-17 0.00 4.5010E-19 -4.5010E-19 Losev[70] 

1.7200E-14 -0.17 4.4730E-19 -4.4730E-19 Bird[15] 

2.9890E-15 0.17 4.5560E-19 -4.5560E-19 Kang[64] 

47 O2

+ + N → N+ + O2 1.4450E-16 0.14 1.0 3.9490E-19 -3.9490E-19 Park[74] 

48 O2

+ + N2 → N2+ + O2 1.6440E-17 0.00 1.5 5.6190E-19 -5.6190E-19 Park[74] 

49 O2

+ + O → O+ + O2 6.6420E-18 -0.09 1.0 2.4850E-19 -2.4850E-19 Park[75] 

1.8900E-16 -0.52 2.5900E-19 -2.5900E-19 Bird[15] 

4.8490E-12 -1.11 3.8650E-19 -3.8650E-19 Kang[64] 

Table A.15: Exothermic chemical reaction rate coefficients (m3/molecule/s) 

used in the TCE model. kr = aTbexp(-Ea/kT)

Number Reaction a b Ea (J) Erx (J) Reference

50 O + O + T → O2 + T 6.3050E-44 -0.50 0.0 0.0000E+01 8.1995E-19 Farbar[6] 

8.3000E-45 -0.50 0.0000E+01 8.1995E-19 Gupta[76] 

51 N + N + T → N2 + T 5.6910E-40 -1.50 0.0 0.0000E+01 1.5630E-18 Farbar[6] 

3.0060E-44 -0.50 0.0000E+01 1.5630E-18 Gupta[76] 

52 N + O + T → NO + T 1.5830E-43 -0.50 0.0 0.0000E+01 1.0420E-18 Farbar[6] 

2.7850E-40 -1.50 0.0000E+01 1.0420E-18 Gupta[76] 

53 NO + N → N2 + O 2.4900E-17 0.10 0.0 0.0000E+01 5.1750E-19 Carlson[78] 

2.0200E-17 0.10 0.0000E+01 5.9280E-19 Bird[15] 

4.0600E-12 -1.36 0.0000E+01 5.1750E-19 Boyd[67] 

2.4900E-17 0.00 0.0000E+01 5.1750E-19 Gimelshein[65] 

54 O2 +N → NO + O 5.2000E-22 1.29 0.0 4.9680E-20 2.2380E-19 Bird[15] 

4.1300E-21 1.18 0.0000E+01 5.5300E-20 Ozawa[9] 

4.6000E-15 -0.55 0.0000E+01 2.2380E-20 Boyd[67] 

1.5980E-18 0.50 4.9680E-20 2.7190E-19 Gimelshein[65] 

55 N2

+ + e-

→ N + N 4.4830E-12 -0.50 0.0 0.0000E+01 9.3400E-19 Park[60] 

8.8800E-10 -1.23 0.0000E+01 9.3400E-19 Bird[15] 

1.0000E-08 -1.43 0.0000E+01 9.3400E-19 Ozawa[7] 

7.2740E-12 -0.65 0.0000E+01 9.3400E-19 Boyd[77] 

ζ

ζ
93 

Number Reaction a b Ea (J) Erx (J) Reference

1.5700E-17 0.85 0.0000E+01 9.3400E-19 Boyd[67] 

2.4910E-08 -1.50 0.0000E+01 9.3400E-19 Kang[64] 

56 NO+ + e-→ N + O 1.4940E-11 -0.65 0.0 0.0000E+01 4.4220E-19 Park[60] 

4.0300E-09 -1.63 0.0000E+01 4.4220E-19 Bird[15] 

2.0000E-06 -2.05 0.0000E+01 4.4220E-19 Ozawa[7] 

1.3210E-09 -1.19 0.0000E+01 4.4220E-19 Boyd[77] 

2.2000E-13 -0.19 0.0000E+01 4.4220E-19 Boyd[67] 

1.3280E-08 -1.50 0.0000E+01 4.4220E-19 Kang[64] 

2.9890E-11 -1.00 0.0000E+01 4.4730E-19 Gupta[72] 

57 O2

+ + e-

→ O + O 2.4900E-12 -0.50 0.0 0.0000E+01 1.1200E-18 Park[60] 

3.8300E-09 -1.51 0.0000E+01 1.1200E-18 Bird[15] 

6.0000E-08 -1.60 0.0000E+01 1.1200E-18 Ozawa[7] 

1.4530E-04 -2.41 0.0000E+01 1.1200E-18 Boyd[77] 

9.2200E-15 0.29 0.0000E+01 1.1200E-18 Boyd[67] 

1.3280E-08 -1.50 0.0000E+01 1.1160E-18 Gupta[72] 

58 N2

+ + N → N+ + N2 2.3700E-18 -0.52 1.0 0.0000E+01 1.6700E-19 Bird[15] 

2.3430E-14 -0.61 0.0000E+01 1.6700E-19 Boyd[77] 

4.6500E-19 0.50 0.0000E+01 1.6700E-19 Kang[64] 

59 N2

+ + O → O+ + N2 1.7700E-17 -0.21 1.0 0.0000E+01 3.1480E-19 Bird[15] 

1.9780E-18 0.11 0.0000E+01 3.1480E-19 Boyd[77] 

4.1180E-11 -2.20 0.0000E+01 3.1480E-19 Kang[64] 

60 N+ + O2 → O+ + NO 2.4430E-26 2.10 1.0 0.0000E+01 2.1120E-19 Boyd[77] 

61 N2

+ + O → NO+ + N 4.1000E-18 0.40 1.0 0.0000E+01 4.9010E-19 Bird[15] 

1.7440E-18 0.30 0.0000E+01 4.9010E-19 Boyd[77] 

62 O+ + N2 → NO+ + N 3.9700E-18 -0.71 1.0 0.0000E+01 1.7670E-19 Boyd[77] 

63 N+ + O2 → NO+ + O 3.0400E-18 0.29 1.0 0.0000E+01 1.0660E-18 Boyd[67] 

1.6610E-16 0.00 0.0000E+01 1.0670E-18 Kang[64] 

64 O2

+ + N → NO+ + O 8.9180E-13 -0.97 1.0 0.0000E+01 6.7090E-19 Boyd[77] 

65 O2

+ + NO → NO+ + O2 3.9900E-17 0.41 1.5 0.0000E+01 4.5010E-19 Ozawa[9] 

6.1950E-16 -0.05 0.0000E+01 4.5010E-19 Boyd[77] 

4.4700E-15 -0.17 0.0000E+01 4.5010E-19 Carlson[78] 

2.9890E-17 0.50 0.0000E+01 4.5010E-19 Kang[64] 

66 N+ + O2 → O2

+ + N 1.0

67 N2

+ + O2 → O2+ + N2 4.5890E-18 -0.04 1.5 0.0000E+01 5.6190E-19 Boyd[77] 

68 O+ + O2 → O2

+ + O 1.8900E-16 -0.52 1.0 0.0000E+01 2.5900E-19 Bird[15] 

4.9930E-18 -0.00 0.0000E+01 2.5900E-19 Boyd[77] 

4.6500E-19 0.50 0.0000E+01 2.5900E-19 Kang[64] 

4.9900E-14 0.00 0.0000E+01 2.5900E-19 Boyd[67] 

Table A.16: Coefficients required for the collision volume in Q-K 

recombination reactions (Vcoll = aTbVref).

Number Reaction a b 

50 O + O + O2 → O2 + O2 0.040 –1.30

O + O + N2 → O2 + N2 0.070 –1.20

O + O + O → O2 + O 0.080 –1.20

O + O + N → O2 + N 0.090 –1.20

O + O + NO → O2 + NO 0.065 –1.20

O + O + O2

+ → O2 + O2+ 0.040 –1.30

O + O + N2

+ → O2 + N2+ 0.070 –1.20

O + O + O+ → O2 + O+ 0.080 –1.20

O + O + N+ → O2 + N+ 0.090 –1.20

O + O + NO+ → O2 + NO+ 0.065 –1.20

O + O + e- → O2 + e- 0.000 0.00

ζ
94 

Number Reaction a b 

51 N + N + O2 → N2 + O2 0.150 –2.05

N + N + N2 → N2 + N2 0.090 –2.10

N + N + O → N2 + O 0.160 –2.00

N + N + N → N2 + N 0.170 –2.00

N + N + NO → N2 + NO 0.170 –2.10

N + N + O2

+ → N2 + O2+ 0.150 –2.05

N + N + N2

+ → N2 + N2+ 0.090 –2.10

N + N + O+ → N2 + O+ 0.160 –2.00

N + N + N+ → N2 + N+ 0.170 –2.00

N + N + NO+ → N2 + NO+ 0.170 –2.10

N + N + e- → N2 + e- 0.000 0.00

52 N + O + O2 → NO + O2 0.300 –1.90

N + O + N2 → NO + N2 0.400 –2.00

N + O + O → NO + O 0.300 –1.75

N + O + N → NO + N 0.300 –1.75

N + O + NO → NO + NO 0.150 –1.90

N + O + O2

+ → NO + O2+ 0.300 –1.90

N + O + N2

+ → NO + N2+ 0.400 –2.00

N + O + O+ → NO + O+ 0.300 –1.75

N + O + N+ → NO + N+ 0.300 –1.75

N + O + NO+ → NO + NO+ 0.150 –1.90

N + O + e- → NO + e- 0.000 0.00

Table A.17: Coefficients required for the collision volume in Q-K exchange

reactions (Ea = a(T/Tref)

b

|Eexch|; Ea = Ea – Eexch if Eexch < 0.0).

Number Reaction a b Eexch (J)

34 N2 + O → NO + N 0.150 0.00 –5.2465E-19

35 NO + O → O2 + N 0.100 0.68 –2.6785E-19

39 N+ + N2 → N2

+ + N 0.400 0.10 –1.6844E-19

40 O+ + N2 → N2

+ + O 0.500 0.20 –3.1479E-19

41 O+ + NO → N+ + O2 0.100 0.50 –3.6725E-19

42 NO+ + N → N2

+ + O 0.050 0.00 –4.9013E-19

43 NO+ + N → O+ + N2 0.050 0.70 –1.7672E-19

44 NO+ + O → N+ + O2 0.000 0.00 –1.0659E-19

45 NO+ + O → O2

+ + N 0.000 0.00 –6.7100E-19

46 NO+ + O2 → O2

+ + NO 0.000 0.00 –4.5009E-19

47 O2

+ + N → N+ + O2 0.020 0.10 –3.9487E-19

48 O2

+ + N2 → N2+ + O2 0.050 0.55 –5.6192E-19

49 O2

+ + O → O+ + O2 0.700 0.15 –2.4852E-19

53 NO + N → N2 + O 0.025 0.90 5.2465E-19

54 O2 + N → NO + O 0.100 0.48 2.6785E-19

58 N2

+ + N → N+ + N2 0.400 0.25 1.6844E-19

59 N2

+ + O → O+ + N2 0.600 0.35 3.1479E-19

60 N+ + O2 → O+ + NO 0.150 0.40 3.6725E-19
95 

Number Reaction a b Eexch (J)

61 N2

+ + O → NO+ + N 0.020 0.80 4.9013E-19

62 O+ + N2 → NO+ + N 0.000 0.00 1.7672E-19

63 N+ + O2 → NO+ + O 0.000 0.00 1.0659E-18

64 O2

+ + N → NO+ + O 0.000 0.00 6.7100E-19

65 O2

+ + NO → NO+ + O2 0.070 0.50 5.4009E-19

66 N+ + O2 → O2

+ + N 0.000 0.00 3.9487E-19

67 N2

+ + O2 → O2+ + N2 0.200 0.30 5.6192E-19

68 O+ + O2 → O2

+ + O 0.500 0.10 2.4852E-19
96 

Appendix B: 

Quantum-Kinetic Electronic Energy Level Transition Rates

Figure B.1: Transitions for O2 + O2.

Figure B.2: Transitions for O2 + O.

Figure B.3: Transitions for O2 + N2.

Figure B.4: Transitions for O2 + N.

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

Rate Constant (m3/molecule/s)

Temperature (K)

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

Rate Constant (m3/molecule/s)

Temperature (K)

10-27

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

Rate Constant (m3/molecule/s)

Temperature (K)

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

Rate Constant (m3/molecule/s)

Temperature (K)
97 

Figure B.5: Transitions for O2 + 

NO.

Figure B.6: Transitions for N2 + N2.

Figure B.7: Transitions for N2 + N.

Figure B.8: Transitions for N2 + O2.

Figure B.9: Transitions for N2 + O.

Figure B.10: Transitions for N2 + 

NO.

10-27

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

Rate Constant (m

3 - 4 Park[90]

3/molecule/s)

Temperature (K)

10-28

10-26

10-24

10-22

10-20

10-18

10-16

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

1 - 2 Gorelov[89]

Rate Constant (m

2 - 3 Cherni[85]

3/molecule/s)

Temperature (K)

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

Rate Constant (m

1 - 2 Gorelov[89]

3/molecule/s)

Temperature (K)

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

1 - 2 Gorelov[89]

Rate Constant (m3/molecule/s)

Temperature (K)

10-24

10-22

10-20

10-18

10-16

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

3 - 4 Park[90]

1 - 2 Gorelov[89]

Rate Constant (m3/molecule/s)

Temperature (K)

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

Rate Constant (m

3 - 4 Park[90]

3/molecule/s)

Temperature (K)
98 

Figure B.11: Transition from level 1–

2 for N2 + e.

Figure B.12: Transition from level 1–

2 for O + e.

Figure B.13: Transition from level 1–

4 for O + e.

Figure B.14: Transition from level 2–

3 for N2 + e.

Figure B.15: Transition from level 1–

3 for O + e.

Figure B.16: Transition from level 

1v5 for O + e.

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Cherni[85]

Telulet[86]

Rate Constant (m

Gorelov[87]

3/molecule/s)

Temperature (K)

N2

1

 + e- > N2

2

 + e10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Zatsarinny[83]

Park[62]

Rate Constant (m

Bhatia[84]

3/molecule/s)

Temperature (K)

O1 + e- > O2 + e10-23

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Zatsarinny[83]

Park[62]

Rate Constant (m

Bhatia[84]

3/molecule/s)

Temperature (K)

O1 + e- > O4 + e10-18

10-17

10-16

10-15

10-14

10-13

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Cherni[85]

Telulet[86]

Rate Constant (m

Gorelov[87]

3/molecule/s)

Temperature (K)

N2

2

 + e- > N2

3

 + e10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Zatsarinny[83]

Park[62]

Rate Constant (m

Bhatia[84]

3/molecule/s)

Temperature (K)

O1 + e- > O3 + e10-24

10-22

10-20

10-18

10-16

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Zatsarinny[83]

Park[62]

Rate Constant (m

Bhatia[84]

3/molecule/s)

Temperature (K)

O1 + e- > O5 + e-
99 

Figure B.17: Transition from level 2–

3 for O + e.

Figure B.18: Transition from level 1–

2 for N + e.

Figure B.19: Transition from level 1–

4 for N + e.

Figure B.20: Transition from level 2–

5 for O + e.

Figure B.21: Transition from level 1–

3 for N + e.

Figure B.22: Transition from level 1–

5 for N + e.

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Zatsarinny[83]

Park[62]

Rate Constant (m

Bhatia[84]

3/molecule/s)

Temperature (K)

O2 + e- > O3 + e10-18

10-17

10-16

10-15

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Park[62]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N1 + e- > N2 + e10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Park[62]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N1 + e- > N4 + e10-24

10-22

10-20

10-18

10-16

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Zatsarinny[83]

Rate Constant (m

Bhatia[84]

3/molecule/s)

Temperature (K)

O2 + e- > O5 + e10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Park[62]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N1 + e- > N3 + e10-25

10-24

10-23

10-22

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N1 + e- > N5 + e-
100 

Figure B.23: Transition from level 2–

3 for N + e.

Figure B.24: Transition from level 2–

5 for N + e.

Figure B.25: Transition from level 3–

5 for N + e.

Figure B.26: Transition from level 2–

4 for N + e.

Figure B.27: Transition from level 3–

4 for N + e.

Figure B.28: Transitions for NO + 

O2.

10-16

10-15

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Park[62]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N2 + e- > N3 + e10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Park[62]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N2 + e- > N5 + e10-21

10-20

10-19

10-18

10-17

10-16

10-15

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Park[62]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N3 + e- > N5 + e10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N2 + e- > N4 + e10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Tayal[81]

Rate Constant (m

Frost[82]

3/molecule/s)

Temperature (K)

N3 + e- > N4 + e10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

Rate Constant (m

3 - 4 Park[90]

3/molecule/s)

Temperature (K)
101 

Figure B.29: Transitions for NO + 

N2.

Figure B.30: Transitions for NO + N.

Figure B.31: Transition from level 1–

3 for NO + e.

Figure B.32: Transitions for NO + O.

Figure B.33: Transitions for NO + 

NO.

Figure B.34: Transition from level 1–

4 for NO + e.

10-27

10-25

10-23

10-21

10-19

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

NO-O2 3-4

1 - 2 Park[90]

2 - 3 Park[90]

Rate Constant (m

3 - 4 Park[90]

3/molecule/s)

Temperature (K)

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

Rate Constant (m

3 - 4 Park[90]

3/molecule/s)

Temperature (K)

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Telulet[86]

Rate Constant (m

Gorelov[87]

3/molecule/s)

Temperature (K)

NO1 + e- > NO3 + e10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

Rate Constant (m

3 - 4 Park[90]

3/molecule/s)

Temperature (K)

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

1 - 2 DSMC

2 - 3 DSMC

3 - 4 DSMC

1 - 2 Park[90]

2 - 3 Park[90]

Rate Constant (m

3 - 4 Park[90]

3/molecule/s)

Temperature (K)

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Telulet[86]

Rate Constant (m

Gorelov[87]

3/molecule/s)

Temperature (K)

NO1 + e- > NO4 + e-
102 

Figure B.35: Transition from level 3–

4 for NO + e.

Figure B.36: Transition from level 1–

3 for N2

+ + e.

Figure B.37: Transition from level 1–

2 for N2

+ + e.

Figure B.38: Transition from level 2–

3 for N2

+ + e.

Table B.1: Sampled temperatures and their departure from the set 

equilibrium temperature.

Tequilibrium (K) 10000.00 % Difference 20000.00 % Difference 30000.00 % Difference

Ttr (K) 9949.81 –0.502% 19918.00 –0.410% 29901.50 –0.328%

Telectron (K) 9953.79 –0.462% 19927.70 –0.361% 29942.00 –0.193%

Tel,O2 (K) 10000.40 0.004% 20018.80 0.094% 30301.90 1.006%

Tel, N2 (K) 10011.80 0.118% 20017.00 0.085% 30060.00 0.200%

Tel, O (K) 10019.80 0.198% 19992.40 –0.038% 29967.50 –0.108%

Tel, N (K) 10009.50 0.095% 19996.20 –0.019% 30040.60 0.135%

Tel, NO (K) 10011.20 0.112% 20037.70 0.189% 30119.40 0.398%

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Telulet[86]

Rate Constant (m

Gorelov[87]

3/molecule/s)

Temperature (K)

NO3 + e- > NO4 + e10-20

10-19

10-18

10-17

10-16

10-15

10-14

10-13

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Gorelov[87]

Rate Constant (m

Telulet[86]

3/molecule/s)

Temperature (K)

N2

+ 1 + e-

 > N2

+ 3 + e10-18

10-17

10-16

10-15

10-14

10-13

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Telulet[86]

Gorelov[87]

Rate Constant (m3/molecule/s)

Temperature (K)

N2

+ 1 + e-

 > N2

+ 2 + e10-18

10-17

10-16

10-15

10-14

10-13

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

DSMC

Gorelov[87]

Rate Constant (m

Telulet[86]

3/molecule/s)

Temperature (K)

N2

+ 2 + e-

 > N2

+ 3 + e-
103 

Tel,O2+ (K) 10012.50 0.125% 20046.50 0.233% 30140.30 0.468%

Tel, N2+ (K) 10015.70 0.157% 20031.50 0.158% 30106.20 0.354%

Tel, O+ (K) 10011.30 0.113% 19993.80 –0.031% 29923.40 –0.255%

Tel, N+ (K) 10022.70 0.227% 19989.50 –0.053% 29896.20 –0.346%

Tel, NO+ (K) 10011.20 0.112% 20014.50 0.073% 30062.80 0.209%

Table B.2: Sampled and Boltzmann distributions for O2 at 10,000 K.

Level Sampled Boltzmann % Difference

1 7.7708E-01 7.7719E-01 –0.014%

2 1.6583E-01 1.6585E-01 –0.010%

3 3.8865E-02 3.8807E-02 0.149%

4 2.2448E-03 2.2280E-03 0.754%

5 1.0584E-02 1.0566E-02 0.175%

6 4.8032E-03 4.7751E-03 0.588%

7 5.9326E-04 5.8615E-04 1.212%

Table B.3: Sampled and Boltzmann distributions for O2 at 20,000 K.

Level Sampled Boltzmann % Difference

1 5.5362E-01 5.5315E-01 0.086%

2 2.0892E-01 2.0863E-01 0.138%

3 7.1264E-02 7.1363E-02 –0.138%

4 1.7149E-02 1.7099E-02 0.292%

5 9.0780E-02 9.1209E-02 –0.470%

6 4.3126E-02 4.3358E-02 –0.535%

7 1.5139E-02 1.5191E-02 –0.342%

Table B.4: Sampled and Boltzmann distributions for O2 at 30,000 K.

Level Sampled Boltzmann % Difference

1 4.2546E-01 4.2474E-01 0.170%

2 1.9476E-01 1.9370E-01 0.550%

3 7.5215E-02 7.5191E-02 0.032%

4 2.9078E-02 2.9007E-02 0.245%

5 1.6010E-01 1.6091E-01 –0.506%

6 7.7062E-02 7.7793E-02 –0.939%

7 3.8326E-02 3.8662E-02 –0.869%
104 

Table B.5: Sampled and Boltzmann distributions for N2 at 10,000 K.

Level Sampled Boltzmann % Difference

1 9.9500E-01 9.9505E-01 –0.005%

2 2.1963E-03 2.1776E-03 0.859%

3 1.1353E-03 1.1238E-03 1.026%

4 1.1045E-03 1.0937E-03 0.992%

5 2.1789E-04 2.1573E-04 1.000%

6 5.5589E-05 5.4867E-05 1.316%

7 9.4670E-05 9.3263E-05 1.509%

Table B.6: Sampled and Boltzmann distributions for N2 at 20,000 K.

Level Sampled Boltzmann % Difference

1 7.3201E-01 7.3248E-01 –0.064%

2 5.9496E-02 5.9350E-02 0.246%

3 6.0508E-02 6.0296E-02 0.352%

4 5.9591E-02 5.9482E-02 0.183%

5 1.8720E-02 1.8681E-02 0.211%

6 5.4476E-03 5.4391E-03 0.156%

7 1.0032E-02 1.0029E-02 0.034%

Table B.7: Sampled and Boltzmann distribution for N2 at 30,000 K.

Level Sampled Boltzmann % Difference

1 3.6019E-01 3.6023E-01 –0.011%

2 9.7743E-02 9.7283E-02 0.473%

3 1.2460E-01 1.2387E-01 0.592%

4 1.2323E-01 1.2275E-01 0.391%

5 4.5125E-02 4.5014E-02 0.246%

6 1.3711E-02 1.3711E-02 0.001%

7 2.5906E-02 2.5975E-02 –0.265%

Table B.8: Sampled and Boltzmann distributions for O at 10,000 K.

Level Sampled Boltzmann % Difference

1 9.4529E-01 9.4553E-01 –0.025%

2 5.3797E-02 5.3567E-02 0.430%

3 8.1975E-04 8.1264E-04 0.875%

4 1.3074E-05 1.2912E-05 1.258%
105 

5 5.1522E-06 5.0119E-06 2.799%

6 6.0933E-06 6.0882E-06 0.084%

7 2.7848E-06 2.7387E-06 1.682%

Table B.9: Sampled and Boltzmann distributions for O at 20,000 K.

Level Sampled Boltzmann % Difference

1 7.4519E-01 7.4541E-01 –0.029%

2 1.3215E-01 1.3224E-01 –0.070%

3 7.2810E-03 7.2842E-03 –0.044%

4 2.0615E-03 2.0531E-03 0.409%

5 9.9362E-04 9.9083E-04 0.282%

6 2.4553E-03 2.4419E-03 0.550%

7 1.2762E-03 1.2686E-03 0.597%

Table B.10: Sampled and Boltzmann distributions for O at 30,000 K.

Level Sampled Boltzmann % Difference

1 3.0872E-01 3.0795E-01 0.250%

2 8.0065E-02 7.9930E-02 0.170%

3 6.7504E-03 6.7670E-03 –0.245%

4 4.9942E-03 4.9743E-03 0.399%

5 2.5918E-03 2.5813E-03 0.406%

6 8.1065E-03 8.0535E-03 0.658%

7 4.4226E-03 4.3898E-03 0.748%

Table B.11: Sampled and Boltzmann distributions for N at 10,000 K.

Level Sampled Boltzmann % Difference

1 8.4636E-01 8.4671E-01 –0.042%

2 1.3341E-01 1.3311E-01 0.225%

3 2.0090E-02 2.0036E-02 0.271%

4 1.5808E-05 1.5762E-05 0.292%

5 5.2667E-06 5.2246E-06 0.807%

6 7.9995E-06 7.9046E-06 1.201%

7 6.2563E-07 6.0151E-07 4.011%
106 

Table B.12: Sampled and Boltzmann distributions for N at 20,000 K.

Level Sampled Boltzmann % Difference

1 4.3952E-01 4.3890E-01 0.141%

2 2.7547E-01 2.7515E-01 0.116%

3 8.2188E-02 8.2688E-02 –0.605%

4 3.2850E-03 3.2799E-03 0.155%

5 1.3379E-03 1.3353E-03 0.198%

6 2.3279E-03 2.3227E-03 0.223%

7 2.6346E-04 2.6158E-04 0.719%

Table B.13: Sampled and Boltzmann distributions for N at 30,000 K.

Level Sampled Boltzmann % Difference

1 1.0682E-01 1.0620E-01 0.580%

2 1.0632E-01 1.0558E-01 0.697%

3 3.9934E-02 3.9954E-02 –0.051%

4 5.8929E-03 5.8551E-03 0.646%

5 2.5708E-03 2.5524E-03 0.722%

6 4.6881E-03 4.6517E-03 0.783%

7 6.0134E-04 5.9704E-04 0.720%

Table B.14: Sampled and Boltzmann distributions for NO at 10,000 K.

Level Sampled Boltzmann % Difference

1 9.8591E-01 9.8601E-01 –0.010%

2 8.3771E-03 8.3267E-03 0.605%

3 8.9642E-04 8.8993E-04 0.730%

4 1.3534E-03 1.3421E-03 0.843%

5 9.9780E-04 9.8763E-04 1.030%

6 8.6417E-04 8.5528E-04 1.040%

7 5.5615E-04 5.5025E-04 1.073%

Table B.15: Sampled and Boltzmann distributions for NO at 20,000 K.

Level Sampled Boltzmann % Difference

1 7.3783E-01 7.3899E-01 –0.157%

2 9.6384E-02 9.6040E-02 0.359%

3 1.5806E-02 1.5699E-02 0.684%

4 2.7450E-02 2.7264E-02 0.683%

5 2.3546E-02 2.3388E-02 0.675%
107 

6 2.1886E-02 2.1765E-02 0.557%

7 1.7534E-02 1.7457E-02 0.439%

Table B.16: Sampled and Boltzmann distributions for NO at 30,000 K.

Level Sampled Boltzmann % Difference

1 4.7041E-01 4.7202E-01 –0.342%

2 1.5317E-01 1.5259E-01 0.383%

3 2.9009E-02 2.8737E-02 0.948%

4 5.2855E-02 5.2311E-02 1.040%

5 4.7618E-02 4.7228E-02 0.826%

6 4.5222E-02 4.5016E-02 0.457%

7 3.8913E-02 3.8862E-02 0.132%

Table B.17: Sampled and Boltzmann distributions for O2

+ at 10,000 K.

Level Sampled Boltzmann % Difference

1 9.7135E-01 9.7155E-01 –0.020%

2 1.7156E-02 1.7058E-02 0.577%

3 2.8567E-03 2.8334E-03 0.821%

4 4.2667E-03 4.2331E-03 0.795%

5 1.9302E-03 1.9136E-03 0.867%

6 5.8441E-04 5.7826E-04 1.063%

7 7.9150E-04 7.8423E-04 0.927%

Table B.18: Sampled and Boltzmann distributions for O2

+ at 20,000 K.

Level Sampled Boltzmann % Difference

1 6.6647E-01 6.6768E-01 –0.181%

2 1.2558E-01 1.2512E-01 0.372%

3 3.6258E-02 3.6057E-02 0.557%

4 5.4176E-02 5.3977E-02 0.369%

5 2.9739E-02 2.9632E-02 0.361%

6 1.1570E-02 1.1518E-02 0.451%

7 1.9025E-02 1.8970E-02 0.293%

Table B.19: Sampled and Boltzmann distributions for O2

+ at 30,000 K.

Level Sampled Boltzmann % Difference

1 4.1158E-01 4.1287E-01 –0.313%

2 1.7107E-01 1.7034E-01 0.427%

3 5.9420E-02 5.8989E-02 0.731%
108 

4 8.8882E-02 8.8364E-02 0.586%

5 5.1954E-02 5.1754E-02 0.386%

6 2.1926E-02 2.1879E-02 0.217%

7 3.8413E-02 3.8443E-02 –0.077%

Table B.20: Sampled and Boltzmann distributions for N2

+ at 10,000 K.

Level Sampled Boltzmann % Difference

1 6.1227E-01 6.1286E-01 –0.096%

2 3.2813E-01 3.2776E-01 0.113%

3 1.5790E-02 1.5718E-02 0.457%

4 3.1316E-02 3.1262E-02 0.173%

5 6.2032E-03 6.1665E-03 0.596%

6 3.2954E-03 3.2741E-03 0.649%

7 6.1995E-04 6.1403E-04 0.964%

Table B.21: Sampled and Boltzmann distributions for N2

+ at 20,000 K.

Level Sampled Boltzmann % Difference

1 3.1457E-01 3.1438E-01 0.059%

2 3.2567E-01 3.2514E-01 0.163%

3 5.0346E-02 5.0348E-02 –0.003%

4 1.0017E-01 1.0042E-01 –0.245%

5 6.3131E-02 6.3071E-02 0.096%

6 4.5969E-02 4.5958E-02 0.025%

7 1.4075E-02 1.4073E-02 0.014%

Table B.22: Sampled and Boltzmann distributions for N2

+ at 30,000 K.

Level Sampled Boltzmann % Difference

1 1.8259E-01 1.8224E-01 0.195%

2 2.3563E-01 2.3481E-01 0.349%

3 5.3865E-02 5.3742E-02 0.229%

4 1.0748E-01 1.0728E-01 0.183%

5 9.9579E-02 9.9136E-02 0.447%

6 8.0402E-02 8.0276E-02 0.157%

7 2.8918E-02 2.8946E-02 –0.098%
109 

Table B.23: Sampled and Boltzmann distributions for O+ at 10,000 K.

Level Sampled Boltzmann % Difference

1 9.4566E-01 9.4590E-01 –0.025%

2 5.0110E-02 4.9904E-02 0.414%

3 4.2276E-03 4.2007E-03 0.641%

Table B.24: Sampled and Boltzmann distributions for O+ at 20,000 K.

Level Sampled Boltzmann % Difference

1 6.9240E-01 6.9214E-01 0.037%

2 2.5131E-01 2.5137E-01 –0.023%

3 5.6287E-02 5.6491E-02 –0.361%

Table B.25: Sampled and Boltzmann distributions for O+ at 30,000 K.

Level Sampled Boltzmann % Difference

1 5.2579E-01 5.2459E-01 0.229%

2 3.6206E-01 3.6242E-01 –0.099%

3 1.1216E-01 1.1299E-01 –0.735%

Table B.26: Sampled and Boltzmann distributions for N+ at 10,000 K.

Level Sampled Boltzmann % Difference

1 9.4103E-01 9.4131E-01 –0.030%

2 5.8007E-02 5.7735E-02 0.471%

3 9.5709E-04 9.4864E-04 0.891%

4 2.7276E-06 2.7013E-06 0.974%

5 1.3391E-07 1.4098E-07 -5.016%

6 0.0000E+01 5.1142E-10 –100.000%

7 9.9063E-10 6.0919E-10 62.615%

Table B.27: Sampled and Boltzmann distributions for N+ at 20,000 K.

Level Sampled Boltzmann % Difference

1 8.3493E-01 8.3487E-01 0.008%

2 1.5404E-01 1.5411E-01 –0.046%

3 8.8414E-03 8.8344E-03 0.079%
110 

4 1.8248E-03 1.8258E-03 –0.056%

5 3.2492E-04 3.2310E-04 0.565%

6 1.4622E-05 1.4505E-05 0.810%

7 2.4824E-05 2.4524E-05 1.222%

Table B.28: Sampled and Boltzmann distributions for N+ at 30,000 K.

Level Sampled Boltzmann % Difference

1 7.6006E-01 7.5955E-01 0.067%

2 2.0205E-01 2.0243E-01 –0.188%

3 1.7556E-02 1.7600E-02 –0.251%

4 1.5094E-02 1.5173E-02 –0.519%

5 4.0278E-03 4.0337E-03 –0.145%

6 4.1872E-04 4.1885E-04 –0.031%

7 7.9535E-04 7.9589E-04 –0.068%

Table B.29: Sampled and Boltzmann distributions for NO+ at 10,000 K.

Level Sampled Boltzmann % Difference

1 9.9599E-01 9.9603E-01 –0.004%

2 1.6518E-03 1.6379E-03 0.849%

3 1.1971E-03 1.1854E-03 0.988%

4 8.2190E-04 8.1217E-04 1.198%

5 1.7649E-04 1.7472E-04 1.015%

6 4.5549E-05 4.4971E-05 1.285%

7 6.9164E-05 6.8332E-05 1.218%

Table B.30: Sampled and Boltzmann distributions for NO+ at 20,000 K.

Level Sampled Boltzmann % Difference

1 7.8349E-01 7.8358E-01 –0.012%

2 5.5182E-02 5.5037E-02 0.263%

3 6.6323E-02 6.6215E-02 0.163%

4 5.4691E-02 5.4809E-02 –0.215%

5 1.7924E-02 1.7975E-02 –0.286%

6 5.2648E-03 5.2652E-03 –0.008%

7 9.1870E-03 9.1786E-03 0.091%
111 

Table B.31: Sampled and Boltzmann distributions for NO+ at 30,000 K.

Level Sampled Boltzmann % Difference

1 4.5955E-01 4.5873E-01 0.178%

2 1.1342E-01 1.1263E-01 0.703%

3 1.6075E-01 1.6052E-01 0.145%

4 1.4048E-01 1.4151E-01 –0.728%

5 5.2954E-02 5.3415E-02 –0.862%

6 1.6242E-02 1.6335E-02 –0.568%

7 2.9677E-02 2.9810E-02 –0.446%
112 

Appendix C: 

Quantum-Kinetic Reaction Rates

Figure C.1: Q-K reaction rates for 

reaction 1.

Figure C.2: Q-K reaction rates for 

reaction 3.

Figure C.3: Q-K reaction rates for 

reaction 2.

Figure C.4: Q-K reaction rates for 

reaction 4.

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytic

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

O2 + O2 > O + O + O2

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytic

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Ozawa[9]

Gimelshein[67]

Rate Constant (m

Johnston[68]

3/molecule/s)

Temperature (K)

O2 + O > O + O + O

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytic

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

O2 + N2 > O + O + N2

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Ozawa[9]

Bird[15]

Kang[66]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

O2 + N > O + O + N
113 

Figure C.5: Q-K reaction rates for 

reaction 5.

Figure C.6: Q-K reaction rates for 

reaction 7.

Figure C.7: Q-K reaction rates for 

reaction 9.

Figure C.8: Q-K reaction rates for 

reaction 6.

Figure C.9: Q-K reaction rates for 

reaction 8.

Figure C.10: Q-K reaction rates for 

reaction 10.

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytic

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

O2 + NO > O + O + NO

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytic

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

O2 + N2

+

 > O + O + N2

+

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytic

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Ozawa[9]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

O2 + N+ > O + O + N+

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytic

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

O2 + O2

+

 > O + O + O2

+

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytic

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Ozawa[9]

Gimelshein[67]

Rate Constant (m

Johnston[68]

3/molecule/s)

Temperature (K)

O2 + O+ > O + O + O+

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytic

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

O2 + NO+ > O + O + NO+
114 

Figure C.11: Q-K reaction rates for 

reaction 11.

Figure C.12: Q-K reaction rates for 

reaction 13.

Figure C.13: Q-K reaction rates for 

reaction 15.

Figure C.14: Q-K reaction rates for 

reaction 12.

Figure C.15: Q-K reaction rates for 

reaction 14.

Figure C.16: Q-K reaction rates for 

reaction 16.

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Ozawa[9]

Bird[15]

Kang[66]

Boyd[69]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

N2 + O2 > N + N + O2

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Boyd[69]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

N2 + O > N + N + O

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Ozawa[9]

Boyd[69]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

N2 + NO > N + N + NO

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Ozawa[9]

Bird[15]

Kang[66]

Rate Constant (m

Boyd[69]

3/molecule/s)

Temperature (K)

N2 + N2 > N + N + N2

10-24

10-22

10-20

10-18

10-16

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Boyd[69]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

N2 + N > N + N + N

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Ozawa[9]

Boyd[69]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

N2 + O2

+

 > N + N + O2

+
115 

Figure C.17: Q-K reaction rates for 

reaction 17.

Figure C.18: Q-K reaction rates for 

reaction 19.

Figure C.19: Q-K reaction rates for 

reaction 21.

Figure C.20: Q-K reaction rates for 

reaction 18.

Figure C.21: Q-K reaction rates for 

reaction 20.

Figure C.22: Q-K reaction rates for 

reaction 22.

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Ozawa[9]

Rate Constant (m

Boyd[69]

3/molecule/s)

Temperature (K)

N2 + N2

+

 > N + N + N2

+

10-24

10-22

10-20

10-18

10-16

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Boyd[69]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

N2 + N+ > N + N + N+

10-27

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Bourdon[70]

Ozawa[9]

Park[71]

Rate Constant (m

Losev[72]

3/molecule/s)

Temperature (K)

N2 + e- > N + N + e10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Boyd[69]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

N2 + O+ > N + N + O+

10-25

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Ozawa[9]

Boyd[69]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

N2 + NO+ > N + N + NO+

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + O2 > N + O + O2
116 

Figure C.23: Q-K reaction rates for 

reaction 23.

Figure C.24: Q-K reaction rates for 

reaction 25.

Figure C.25: Q-K reaction rates for 

reaction 27.

Figure C.26: Q-K reaction rates for 

reaction 24.

Figure C.27: Q-K reaction rates for 

reaction 26.

Figure C.28: Q-K reaction rates for 

reaction 28.

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + N2 > N + O + N2

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + N > N + O + N

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + O2

+

 > N + O + O2

+

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + O > N + O + O

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + NO > N + O + NO

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + N2

+

 > N + O + N2

+
117 

Figure C.29: Q-K reaction rates for 

reaction 29.

Figure C.30: Q-K reaction rates for 

reaction 31.

Figure C.31: Q-K reaction rates for 

reaction 33.

Figure C.32: Q-K reaction rates for 

reaction 30.

Figure C.33: Q-K reaction rates for 

reaction 32.

Figure C.34: Q-K reaction rates for 

reaction 34.

10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + O+ > N + O + O+

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + NO+ > N + O + NO+

10-30

10-28

10-26

10-24

10-22

10-20

10-18

10-16

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Park[62]

Ozawa[7]

Losev[72]

Boyd[69]

Bird[15]

Kang[66]

Rate Constant (m

Farbar[5]

3/molecule/s)

Temperature (K)

O + e- > O+ + e- + e10-23

10-21

10-19

10-17

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

NO + N+ > N + O + N+

10-30

10-28

10-26

10-24

10-22

10-20

10-18

10-16

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Park[62]

Ozawa[7]

Losev[72]

Boyd[69]

Bird[15]

Kang[66]

Rate Constant (m

Farbar[5]

3/molecule/s)

Temperature (K)

N + e- > N+ + e- + e10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Ozawa[9]

Bird[15]

Kang[66]

Gimelshein[67]

Rate Constant (m

Park[71]

3/molecule/s)

Temperature (K)

N2 + O > NO + N
118 

Figure C.35: Q-K reaction rates for 

reaction 35.

Figure C.36: Q-K reaction rates for 

reaction 37.

Figure C.37: Q-K reaction rates for 

reaction 39.

Figure C.38: Q-K reaction rates for 

reaction 36.

Figure C.39: Q-K reaction rates for 

reaction 38.

Figure C.40: Q-K reaction rates for 

reaction 40.

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Rate Constant (m

Gimelshein[67]

3/molecule/s)

Temperature (K)

NO + O > O2 + N

10-23

10-22

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Park[62]

Kang[66]

Park[71]

Rate Constant (m

Gupta[74]

3/molecule/s)

Temperature (K)

N + O > NO+ + e10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Kang[66]

Rate Constant (m

Phelps[75]

3/molecule/s)

Temperature (K)

N+ + N2 > N2

+

 + N

10-24

10-23

10-22

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]

Bird[15]

Park[62]

Rate Constant (m

Kang[66]

3/molecule/s)

Temperature (K)

N + N > N2

+

 + e10-28

10-26

10-24

10-22

10-20

10-18

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]

Bird[15]

Park[62]

Rate Constant (m

Kang[66]

3/molecule/s)

Temperature (K)

O + O > O2

+

 + e10-22

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]

Bird[15]

Kang[66]

Rate Constant (m

Losev[72]

3/molecule/s)

Temperature (K)

O+ + N2 > N2

+

 + O
119 

Figure C.41: Q-K reaction rates for 

reaction 41.

Figure C.42: Q-K reaction rates for 

reaction 43.

Figure C.43: Q-K reaction rates for 

reaction 45.

Figure C.44: Q-K reaction rates for 

reaction 42.

Figure C.45: Q-K reaction rates for 

reaction 44.

Figure C.46: Q-K reaction rates for 

reaction 46.

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Rate Constant (m

Park[76]

3/molecule/s)

Temperature (K)

O+ + NO > N+ + O2

10-22

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]

Rate Constant (m3/molecule/s)

Temperature (K)

NO+ + N > O+ + N2

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Rate Constant (m

Fujita[73]

3/molecule/s)

Temperature (K)

NO+ + O > O2

+

 + N

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]

Rate Constant (m

Bird[15]

3/molecule/s)

Temperature (K)

NO+ + N > N2

+

 + O

10-23

10-22

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]

Rate Constant (m

Kang[66]

3/molecule/s)

Temperature (K)

NO+ + O > N+ + O2

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]

Bird[15]

Kang[66]

Rate Constant (m

Losev[72]]

3/molecule/s)

Temperature (K)

NO+ + O2 > O2

+

 + NO
120 

Figure C.47: Q-K reaction rates for 

reaction 47.

Figure C.48: Q-K reaction rates for 

reaction 49.

Figure C.49: Q-K reaction rates for 

reaction 51.

Figure C.50: Q-K reaction rates for 

reaction 48.

Figure C.51: Q-K reaction rates for 

reaction 50.

Figure C.52: Q-K reaction rates for 

reaction 52.

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Rate Constant (m

Park[76]

3/molecule/s)

Temperature (K)

O2

+

 + N > N+ + O2

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[77]

Bird[15]

Rate Constant (m

Kang[66]

3/molecule/s)

Temperature (K)

O2

+

 + O > O+ + O2

10-48

10-47

10-46

10-45

10-44

0 5000 10000 15000 20000 25000 30000 35000

Q-K Analytic

Q-K Sampled

Park[65]*

Farbar[6]

Gupta[78]

Rate Constant (m3/molecule/s)

Temperature (K)

N + N + T > N2 + T

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Rate Constant (m

Park[76]

3/molecule/s)

Temperature (K)

O2

+

 + N2 > N2

+

 + O2

10-47

10-46

10-45

10-44

0 5000 10000 15000 20000 25000 30000 35000

Q-K Analytic

Q-K Sampled

Park[65]*

Farbar[6]

Gupta[78]

Rate Constant (m3/molecule/s)

Temperature (K)

O + O + T > O2 + T

10-47

10-46

10-45

10-44

10-43

0 5000 10000 15000 20000 25000 30000 35000

Q-K Analytic

Q-K Sampled

Fujita[73]*

Farbar[6]

Gupta[78]

Rate Constant (m3/molecule/s)

Temperature (K)

N + O + T > NO + T
121 

Figure C.53: Q-K reaction rates for 

reaction 53.

Figure C.54: Q-K reaction rates for 

reaction 55.

Figure C.55: Q-K reaction rates for 

reaction 57.

Figure C.56: Q-K reaction rates for 

reaction 54.

Figure C.57: Q-K reaction rates for 

reaction 56.

Figure C.58: Q-K reaction rates for 

reaction 58.

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]*

Bird[15]

Gimelshein[67]

Boyd[69]

Rate Constant (m

Carlson[80]

3/molecule/s)

Temperature (K)

NO + N > N2 + O

10-18

10-17

10-16

10-15

10-14

10-13

0 5000 10000 15000 20000 25000 30000 35000

Q-K Analytical

Q-K Sampled

Park[65]*

Park[62]

Bird[15]

Ozawa[7]

Boyd[79]

Boyd[69]

Kang[66]

Rate Constant (m3/molecule/s)

Temperature (K)

N2

+

 + e- > N + N

10-16

10-15

10-14

10-13

10-12

10-11

10-10

0 5000 10000 15000 20000 25000 30000 35000

Q-K Analytical

Q-K Sampled

Fujita[73]*

Park[62]

Bird[15]

Ozawa[7]

Boyd[79]

Boyd[69]

Gupta[74]

Rate Constant (m3/molecule/s)

Temperature (K)

O2

+

 + e- > O + O

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]*

Ozawa[9]

Bird[15]

Gimelshein[67]

Rate Constant (m

Boyd[69]

3/molecule/s)

Temperature (K)

O2 + N > NO + O

10-18

10-17

10-16

10-15

10-14

0 5000 10000 15000 20000 25000 30000 35000

Q-K Analytical

Q-K Sampled

Park[65]*

Park[62]

Bird[15]

Ozawa[7]

Boyd[79]

Boyd[69]

Kang[66]

Gupta[74]

Rate Constant (m3/molecule/s)

Temperature (K)

NO+ + e- > N + O

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[65]*

Bird[15]

Kang[66]

Rate Constant (m

Boyd[79]

3/molecule/s)

Temperature (K)

N2

+

 + N > N+ + N2
122 

Figure C.59: Q-K reaction rates for 

reaction 59.

Figure C.60: Q-K reaction rates for 

reaction 61.

Figure C.61: Q-K reaction rates for 

reaction 63.

Figure C.62: Q-K reaction rates for 

reaction 60.

Figure C.63: Q-K reaction rates for 

reaction 62.

Figure C.64: Q-K reaction rates for 

reaction 64.

10-23

10-22

10-21

10-20

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

rx_59

QK_59

Park[76]*

Bird[15]

Kang[66]

Rate Constant (m

Boyd[79]

3/molecule/s)

Temperature (K)

N2

+

 + O > O+ + N2

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]*

Bird[15]

Boyd[79]

Rate Constant (m3/molecule/s)

Temperature (K)

N2

+

 + O > NO+ + N

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]*

Boyd[69]

Kang[66]

Rate Constant (m3/molecule/s)

Temperature (K)

N+ + O2 > NO+ + O

10-19

10-18

10-17

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]*

Rate Constant (m

Boyd[79]

3/molecule/s)

Temperature (K)

N+ + O2 > O+ + NO

10-21

10-20

10-19

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]*

Boyd[79]

Rate Constant (m3/molecule/s)

Temperature (K)

O+ + N2 > NO+ + N

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Fujita[73]*

Boyd[79]

Rate Constant (m3/molecule/s)

Temperature (K)

O2

+

 + N > NO+ + O
123 

Figure C.65: Q-K reaction rates for 

reaction 65.

Figure C.66: Q-K reaction rates for 

reaction 67.

Figure C.67: Q-K reaction rates for 

reaction 66.

Figure C.68: Q-K reaction rates for 

reaction 68.

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]*

Ozawa[9]

Kang[66]

Boyd[79]

Rate Constant (m

Carlson[80]

3/molecule/s)

Temperature (K)

O2

+

 + NO > NO+ + O2

10-18

10-17

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]*

Boyd[79]

Rate Constant (m3/molecule/s)

Temperature (K)

N2

+

 + O2 > O2

+

 + N2

10-16

10-15

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[76]*

Rate Constant (m3/molecule/s)

Temperature (K)

N+ + O2 > O2

+

 + N

10-22

10-21

10-20

10-19

10-18

10-17

10-16

10-15

10-14

0 5x103 1x104 2x104 2x104 3x104 3x104 4x104

Q-K Analytical

Q-K Sampled

Park[77]*

Bird[15]

Kang[66]

Boyd[69]

Rate Constant (m

Boyd[79]

3/molecule/s)

Temperature (K)

O+ + O2 > O2

+

 + O
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Publication

 4. TITLE AND SUBTITLE

Extension of a Kinetic Approach to Chemical Reactions to Electronic 

Energy Levels and Reactions Involving Charged Species with Application 

to DSMC Simulations

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Liechty, Derek Shane

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20376

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 34

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

The ability to compute rarefied, ionized hypersonic flows is becoming more important as missions such as Earth reentry, landing high mass payloads on Mars, 

and the exploration of the outer planets and their satellites are being considered. Recently introduced molecular-level chemistry models that predict equilibrium 

and nonequilibrium reaction rates using only kinetic theory and fundamental molecular properties are extended in the current work to include electronic energy 

level transitions and reactions involving charged particles. These extensions are shown to agree favorably with reported transition and reaction rates from the 

literature for near-equilibrium conditions. Also, the extensions are applied to the second flight of the Project FIRE flight experiment at 1634 seconds with a 

Knudsen number of 0.001 at an altitude of 76.4 km. In order to accomplish this, NASA’s direct simulation Monte Carlo code DAC was rewritten to include the 

ability to simulate charge-neutral ionized flows, take advantage of the recently introduced chemistry model, and to include the extensions presented in this 

work. The 1634 second data point was chosen for comparisons to be made in order to include a CFD solution. The Knudsen number at this point in time is 

such that the DSMC simulations are still tractable and the CFD computations are at the edge of what is considered valid because, although near-transitional, 

the flow is still considered to be continuum. It is shown that the inclusion of electronic energy levels in the DSMC simulation is necessary for flows of this 

nature and is required for comparison to the CFD solution. 

15. SUBJECT TERMS

Computational fluids dynamics; Computerized simulations; Hypersonic flow; Ionized; Rarefied

18. NUMBER

 OF 

 PAGES

139

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

432938.11.01.07.43.40.05

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TP-2014-218254

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 04 - 2014