Hugh Douglas Perkins

Glenn Research Center, Cleveland, Ohio

Development and Demonstration of a Computational 

Tool for the Analysis of Particle Vitiation Effects in 

Hypersonic Propulsion Test Facilities

NASA/TM—2010-216765

October 2010
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c 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 Offi cer. 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 Reports 

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 signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings 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. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

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 via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
Hugh Douglas Perkins

Glenn Research Center, Cleveland, Ohio

Development and Demonstration of a Computational 

Tool for the Analysis of Particle Vitiation Effects in 

Hypersonic Propulsion Test Facilities

NASA/TM—2010-216765

October 2010

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Acknowledgments

I would like to express my appreciation for the ongoing support, counsel and patience of my graduate advisor, Dr. Chih-Jen Sung, 

in the preparation of this work. I would also like to thank my colleagues Dr. Dan Paxson and Mr. Scott Thomas for their detailed 

review of this dissertation. Lastly, I would like to thank the NASA Glenn Research Center for its ongoing support of this research 

and of my overall graduate program.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://gltrs.grc.nasa.gov

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.

Level of Review: This material has been technically reviewed by technical management. 
NASA/TM—2010-216765 iii

Contents

Abstract ......................................................................................................................................................... 1

1.0 Introduction .......................................................................................................................................... 1

2.0 Code Description .................................................................................................................................. 6

2.1 General Approach ....................................................................................................................... 6

2.2 Gas Phase Governing Equations ................................................................................................. 6

2.3 Thermodynamic Properties ......................................................................................................... 7

2.4 Transport Properties .................................................................................................................... 8

2.5 Finite-Rate Kinetics .................................................................................................................. 10

2.6 Solid Phase Governing Equations and Numerical Solution Method ........................................ 11

2.7 Gas Phase Numerical Solution Method .................................................................................... 20

3.0 Validation ........................................................................................................................................... 23

3.1 Gas Phase Flow and Chemical Kinetics Validation .................................................................. 23

3.2 Particle Drag and Heat Transfer Models and Gas-Solid Momentum and Heat Exchange 

Models Validation ..................................................................................................................... 30

3.3 Particle Internal Heat Transfer, Heat Release and Supersonic Drag/Convection Validation ... 34

4.0 Generation of Graphite Kinetic Constants ......................................................................................... 37

5.0 Model Facility and Scramjet Simulation Geometry and Flow Conditions ........................................ 39

5.1 Test Facility Components ......................................................................................................... 39

5.2 Scramjet Inlet and Isolator Sections.......................................................................................... 40

5.3 Modeling Approach .................................................................................................................. 44

5.4 Grid Descriptions ...................................................................................................................... 45

6.0 Results ................................................................................................................................................ 45

6.1 Inert Particle Calculations ......................................................................................................... 45

6.2 Single Graphite Particle Calculations ....................................................................................... 57

6.3 One Percent Mass Fraction Graphite Particle Simulations ....................................................... 67

6.4 Ignition Calculations ................................................................................................................. 78

7.0 Conclusions ........................................................................................................................................ 80

8.0 Recommendations for Future Work ................................................................................................... 81

8.1 Modeling Improvements ........................................................................................................... 81

8.2 Input Data Improvements ......................................................................................................... 82

8.3 Additional Modeling ................................................................................................................. 82

Appendix A.—Symbols List ....................................................................................................................... 83

Appendix B.—Gas Phase Kinetic Mechanism ........................................................................................... 87

References ................................................................................................................................................... 89

List of Tables

Table 1.—Mixer Section Inflow Total Pressure, Total Temperature and dynamic pressure for 

Mach 5, 6 and 7 Simulated Flight Conditions .................................................................................... 40

Table 2.—Summary of Model Scramjet Inlet/Isolator Geometry at Mach 7 Design Point ....................... 41

Table 3.—Summary of Flow Field Conditions for Model Scramjet .......................................................... 42

Table 4.—Test Article Configuration for Mach 5, 6 and 7 Test Conditions .............................................. 44

Table 5.—Flow Conditions Along the Pre-Compression Plate for Mach 5, 6, and 7 Test 

Conditions .......................................................................................................................................... 44

Table 6.—Gas Temperatures (K) at the Exit of the Isolator Section .......................................................... 78

Table 7.—Species Mass Fractions at the Exit of the Isolator Section ........................................................ 79

Table 8.—Particle Surface Temperature, Velocity, and Porosity for Surviving Graphite Particles 

at the Exit of the Isolator Section ....................................................................................................... 79

Table 9.—Ignition Delay Time in Microseconds with Vitiates .................................................................. 80
NASA/TM—2010-216765 iv

Table 10.—Ignition Delay Time in Microseconds without Vitiates ........................................................... 80

Table B.1.—Gas-Phase Kinetic Mechanismc

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

Table B.2.—Third Body Efficiencies ......................................................................................................... 88

List of Figures

Figure 1.—GASL Leg IV diagram. .............................................................................................................. 4

Figure 2.—HTF diagram. ............................................................................................................................. 5

Figure 3.—Particle heat transfer and reaction model.................................................................................. 13

Figure 4.—Effectiveness factor as a function of temperature. .................................................................... 17

Figure 5.—Detonation case geometry. ....................................................................................................... 24

Figure 6.—Detonation calculation reference frames. ................................................................................. 24

Figure 7.—Detonation wave speed convergence. ....................................................................................... 25

Figure 8.—Comparison of pressure profiles between Q1D and ZND codes. ............................................. 26

Figure 9.—Comparison of temperature profiles between Q1D and ZND codes. ....................................... 26

Figure 10.—Comparison of velocity profiles between Q1D and ZND codes. ........................................... 27

Figure 11.—Comparison of H2O mole fraction profiles between Q1D and ZND codes. ........................... 27

Figure 12.—Comparison of detonation OH mole fraction profiles between Q1D and ZND codes. .......... 28

Figure 13.—Comparison of detonation CO2 mole fraction profiles between Q1D and ZND codes. ......... 28

Figure 14.—Comparison of detonation NO mole fraction profiles between Q1D and ZND codes. .......... 29

Figure 15.—Comparison of detonation NO2 mole fraction profiles between Q1D and ZND codes. ......... 29

Figure 16.—Gas temperature profile from Q1D calculation from point of particle insertion 

compared to gas temperature profile from Maxwell/Seasholtz analysis. ........................................... 32

Figure 17.—Gas velocity profile from Q1D calculation from point of particle insertion compared 

to gas velocity profile from Maxwell/Seasholtz analysis. .................................................................. 32

Figure 18.—Particle temperature profile from Q1D calculation from point of particle insertion 

compared to particle temperature profile from Maxwell/Seasholtz analysis. .................................... 33

Figure 19.—Particle velocity profile from Q1D calculation from point of particle insertion 

compared to particle velocity profile from Maxwell/Seasholtz analysis. .......................................... 33

Figure 20.—A comparison of the post-shock gas-particle relative Mach number profiles between 

the Q1D and Baek codes. ................................................................................................................... 35

Figure 21.—Internal temperature profiles at 2 μsec intervals of particle heating up to ignition 

from Baek code. ................................................................................................................................. 35

Figure 22.—Internal temperature profiles at 2 μsec intervals of particle heating up to ignition 

from Q1D code. .................................................................................................................................. 36

Figure 23.—PGX graphite thermal conductivity. ....................................................................................... 38

Figure 24.—PGX graphite heat capacity. ................................................................................................... 38

Figure 25.—Q1D kinetic constant fit to Froberg data. ............................................................................... 39

Figure 26.—Scaled axisymmetric geometry for facility settling chamber/mixer and nozzle. .................... 40

Figure 27.—Scramjet geometric design parameters. .................................................................................. 41

Figure 28.—Scaled inlet/isolator at Mach 7 design point. .......................................................................... 41

Figure 29.—Scaled inlet/isolator with shock locations corresponding to Mach 6 operation. ..................... 42

Figure 30.—Scaled inlet/isolator with shock locations corresponding to Mach 5 operation. ..................... 42

Figure 31.—(a) Truncated configuration for wind tunnel testing with forebody simulation 

(precompression) plate and (b) truncated configuration with forebody simulation plate 

shown schematically with facility mixer and nozzle indicating computational domain. ................... 43

Figure 32.—Nozzle area distribution where X = 0 indicates the beginning of the convergent 

section of the nozzle. .......................................................................................................................... 45

Figure 33.—Tabular alumina heat capacity. ............................................................................................... 46

Figure 34.—Tabular alumina thermal conductivity. ................................................................................... 46
NASA/TM—2010-216765 v 

Figure 35.—Comparison of gas and single particle surface temperature profiles for various 

particle sizes for the Mach 5 total temperature flow condition. ......................................................... 48

Figure 36.—Comparison of gas and single particle velocity profiles for various particle sizes for 

the Mach 5 total temperature flow condition. .................................................................................... 48

Figure 37.—Comparison of gas and single particle surface temperature profiles for various 

particle sizes for the Mach 6 total temperature flow condition. ......................................................... 49

Figure 38.—Comparison of gas and single particle velocity profiles for various particle sizes for 

the Mach 6 total temperature flow condition. .................................................................................... 49

Figure 39.—Comparison of gas and 10 μm particle surface temperature profiles for Mach 5 and 6 

total temperature flow conditions. ...................................................................................................... 50

Figure 40.—Comparison of gas and 100 μm particle surface temperature profiles for Mach 5 and 

6 total temperature flow conditions. ................................................................................................... 50

Figure 41.—Comparison of gas and 1000 μm particle surface temperature profiles for Mach 5 

and 6 total temperature flow conditions. ............................................................................................ 51

Figure 42.—Comparison of gas and 10 μm particle velocity profiles for Mach 5 and 6 total 

temperature flow conditions. .............................................................................................................. 51

Figure 43.—Comparison of gas and 100 μm particle velocity profiles for Mach 5 and 6 total 

temperature flow conditions. .............................................................................................................. 52

Figure 44.—Comparison of gas and 1000 μm particle velocity profiles for Mach 5 and 6 total 

temperature flow conditions. .............................................................................................................. 52

Figure 45.—Radial temperature profiles of 100 μm particles for the Mach 5 and 6 conditions at 

the exit planes of the mixer/nozzle (MN) and the isolator (IS). ......................................................... 53

Figure 46.—Radial temperature profiles of 1000 μm particles for the Mach 5 and 6 conditions at 

the exit planes of the mixer/nozzle (MN) and the isolator (IS). ......................................................... 54

Figure 47.—Comparison of gas and particle temperature profiles for the Mach 6 total temperature 

flow condition for a single 10 μm particle and for a 1 percent 10 μm particle mass fraction. ........... 54

Figure 48.—Comparison of gas temperature profiles for the Mach 6 total temperature flow 

condition for a single 10 μm particle and for a 1 percent particle mass fraction of 10 μm 

particles over the last 1.0 m of the nozzle, the precompression plate, the final external 

compression ramp, and the isolator. ................................................................................................... 55

Figure 49.—Comparison of particle temperature profiles for the Mach 6 total temperature flow 

condition for a single 10 μm particle and for a 1 percent particle mass fraction of 10 μm 

particles over the last 1.0 m of the nozzle, the precompression plate, the final external 

compression ramp, and the isolator. ................................................................................................... 55

Figure 50.—Comparison of gas and particle velocity profiles for the Mach 6 condition for a 

single 10 μm particle and for a 1 percent 10 μm particle mass fraction. ........................................... 56

Figure 51.—Comparison of gas velocity profiles for the Mach 6 total temperature flow condition 

for a single 10 μm particle and for a 1 percent particle mass fraction of 10 μm particles over 

the last 1.0 m of the nozzle, the precompression plate, the final external compression ramp, 

and the isolator. .................................................................................................................................. 56

Figure 52.—Comparison of particle velocity profiles for the Mach 6 total temperature flow 

condition for a single 10 μm particle and for a 1 percent particle mass fraction of 10 μm 

particles over the last 1.0 m of the nozzle, the precompression plate, the final external 

compression ramp, and the isolator. ................................................................................................... 57

Figure 53.—Mass fractions of minor species from Mach 7 total temperature reacting flow 

simulation ........................................................................................................................................... 58

Figure 54.—Comparison of mass fractions of major species from Mach 7 total temperature 

simulations with and without finite-rate chemistry. ........................................................................... 58

Figure 55.—Comparison of static temperature profiles from Mach 7 total temperature simulations 

with and without finite-rate chemistry. .............................................................................................. 59
NASA/TM—2010-216765 vi

Figure 56.—Mach 5 graphite versus alumina surface temperature comparison. ........................................ 60

Figure 57.—Mach 6 graphite versus alumina surface temperature comparison. ........................................ 60

Figure 58.—Mach 7 surface temperature results. ....................................................................................... 61

Figure 59.—Radial temperature profiles of the 1000 μm particle for the Mach 5, 6 and 7 

conditions at the exit planes of the mixer/nozzle (MN) and the isolator (IS)..................................... 62

Figure 60.—Mach 5 graphite versus alumina velocity comparison............................................................ 62

Figure 61.—Mach 6 graphite versus alumina velocity comparison............................................................ 63

Figure 62.—Mach 7 velocity results. .......................................................................................................... 63

Figure 63.—Graphite particle porosity. ...................................................................................................... 64

Figure 64.—Mach 5 gas-particle relative mach number. ............................................................................ 65

Figure 65.—Heat flows for 10 μm graphite particle at Mach 5 total temperature condition...................... 66

Figure 66.—Heat flows for 100 μm graphite particle at Mach 5 total temperature condition.................... 66

Figure 67.—Heat flows for 1000 μm graphite particle outer shell at Mach 5 total temperature 

condition. ............................................................................................................................................ 67

Figure 68.—Comparison of particle temperature and velocity profiles between a single 10 μm 

graphite particle and a 1 percent mass fraction of 10 μm graphite particles at the Mach 5 

total temperature flow condition. ....................................................................................................... 68

Figure 69.—Comparison of particle temperature and velocity profiles between a single 10 μm 

graphite particle and a 1 percent mass fraction of 10 μm graphite particles at the Mach 6 

total temperature flow condition. ....................................................................................................... 69

Figure 70.—Comparison of particle temperature and velocity profiles between a single 10 μm 

graphite particle and a 1 percent mass fraction of 10 μm graphite particles at the Mach 7 

total temperature flow condition. ....................................................................................................... 69

Figure 71.—Comparison of particle temperature and velocity profiles between a single 100 μm 

graphite particle and a 1 percent mass fraction of 100 μm graphite particles at the Mach 5 

total temperature flow condition. ....................................................................................................... 70

Figure 72.—Comparison of particle temperature and velocity profiles between a single 100 μm 

graphite particle and a 1 percent mass fraction of 100 μm graphite particles at the Mach 6 

total temperature flow condition. ....................................................................................................... 71

Figure 73.—Comparison of particle temperature and velocity profiles between a single 100 μm 

graphite particle and a 1 percent mass fraction of 100 μm graphite particles at the Mach 7 

total temperature flow condition. ....................................................................................................... 71

Figure 74.—Comparison of particle temperature and velocity profiles between a single 1000 μm 

graphite particle and a 1 percent mass fraction of 1000 μm graphite particles at the Mach 5 

total temperature flow condition. ....................................................................................................... 72

Figure 75.—Comparison of particle temperature and velocity profiles between a single 1000 μm 

graphite particle and a 1 percent mass fraction of 1000 μm graphite particles at the Mach 6 

total temperature flow condition. ....................................................................................................... 72

Figure 76.—Comparison of particle temperature and velocity profiles between a single 1000 μm 

graphite particle and a 1 percent mass fraction of 1000 μm graphite particles at the Mach 7 

total temperature flow condition. ....................................................................................................... 73

Figure 77.—Comparison of particle porosity profiles between single graphite particles and a 

1 percent mass fraction of graphite particles at the Mach 5 total temperature flow condition. .......... 73

Figure 78.—Comparison of particle porosity profiles between single graphite particles and a 

1 percent mass fraction of graphite particles at the Mach 6 total temperature flow condition. .......... 74

Figure 79.—Comparison of particle porosity profiles between single graphite particles and a 

1 percent mass fraction of graphite particles at the Mach 7 total temperature flow condition. .......... 74

Figure 80.—Comparison of NO mass fraction profiles between no graphite particles and a 

1 percent mass fraction of graphite particles. ..................................................................................... 75
NASA/TM—2010-216765 vii

Figure 81.—Comparison of NO2 mass fraction profiles between no graphite particles and a 1 

percent mass fraction of graphite particles. ........................................................................................ 76

Figure 82.—Comparison of O mass fraction profiles between no graphite particles and a 1 

percent mass fraction of graphite particles. ........................................................................................ 76

Figure 83.—Comparison of CO2 mass fraction profiles between 1 percent graphite particle mass 

fraction cases. ..................................................................................................................................... 77

Figure 84.—Comparison of CO mass fraction profiles between 1 percent graphite particle mass 

fraction cases. ..................................................................................................................................... 78

Figure 85.—Example of ignition point for a H2/vitiated air/particle mixture from the Mach 6 total 

temperature simulation with a 1 percent mass fraction of 1000 μm graphite particles. ..................... 79

NASA/TM—2010-216765 1 

Development and Demonstration of a Computational Tool for the 

Analysis of Particle Vitiation Effects in Hypersonic 

Propulsion Test Facilities

Hugh Douglas Perkins 

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135 

Abstract

In order to improve the understanding of particle vitiation effects in hypersonic propulsion test 

facilities, a quasi-one dimensional numerical tool was developed to efficiently model reacting particle-gas 

flows over a wide range of conditions. Features of this code include gas-phase finite-rate kinetics, a global 

porous-particle combustion model, mass, momentum and energy interactions between phases, and 

subsonic and supersonic particle drag and heat transfer models. The basic capabilities of this tool were

validated against available data or other validated codes. 

To demonstrate the capabilities of the code, and to provide initial insight into the effects of various 

particle laden flows on ignition, a series of computations were performed for a model hypersonic 

propulsion test facility and scramjet. Parameters studied were simulated flight Mach number (Mach 5, 6, 

and 7), particle size (10, 100, and 1000 μm diameters), particle mass fraction (single particle and 1 

percent) and particle material (alumina and graphite). For the alumina particles, it was found that the 

presence of particles up to 1 percent mass fraction had very little effect on the gas phase, even though

only the 10 μm particles closely followed the gas flow velocity and temperature. With the graphite 

particles, the 10 μm particles were either quickly quenched, or were quickly consumed, depending on the 

gas temperature. As the particle size was increased to 100 μm, the particles did not quench, but were still 

typically consumed within the model test facility. For the 1000 μm particles, combustion was diffusion 

limited, so particle and gas temperature had little effect on the combustion rate. When the particle mass 

fraction was increased to 1 percent, the main change was the addition of significant heat release. In those 

cases where low graphite reaction rates were observed for single particles, the increase to 1 percent mass 

fraction had very little impact. 

Hydrogen/vitiated air ignition delay calculations for the 1 percent mass fraction of graphite particles 

cases showed significant decreases in ignition delay in cases where higher graphite combustion rates were 

observed. Further calculations showed that this was due primarily to increased combustor inlet 

temperature, not the gaseous or solid vitiate species present in the flow.

1.0 Introduction

The principle challenge associated with ground testing of hypersonic airbreathing propulsion systems 

is the requirement to heat the air supply to temperatures in excess of 1300 K prior to expansion through 

the facility nozzle. Over the last 40 years, several different methods have been employed to accomplish 

this heat addition. The most commonly used method is combustion heating, wherein fuel is burned in the 

air stream to bring the flow temperature up to the required stagnation temperature, and then makeup 

oxygen is added to bring the flow composition back to the normal 21 percent oxygen mole fraction. 

Several facilities have been built which utilize an electric arc to heat the air in a stilling chamber to very 

high temperatures. This flow is then mixed with additional high pressure air to create an air flow at the 

desired stagnation temperature which is then supplied to the facility nozzle. Lastly, heat sink storage 

heaters have also been built for several additional facilities. In these facilities, a large mass of a refractory 

material in the form of pebbles or drilled blocks is heated to a temperature significantly above the desired 
NASA/TM—2010-216765 2 

air stagnation temperature prior to the propulsion system test. The test air then flows through the 

refractory material which heats the air to a temperature slightly above the required stagnation 

temperature. Additional high pressure air is then added to the heated air to achieve the test stagnation 

temperature set point. It should be noted that heat exchangers and electrical resistance heaters that are 

typically used for lower Mach number propulsion test facilities cannot be used at the significantly higher 

temperatures required for hypersonic propulsion testing due to structural and heat flux limitations.

Detailed descriptions of the operational procedures and limitations of each type of hypersonic propulsion 

test facility can be found in the JANNAF Scramjet Testing Recommended Practices and Guidelines

(Weber, 2002). 

Each of the heating methods described above introduces contamination into the facility air flow, 

typically referred to as vitiation. This flow vitiation results in some level of departure of the ground test 

results from that which would be seen in free flight. For the combustion heated facility, the products of 

combustion remain in the air flow. If hydrogen is used as the heater fuel, the product is primarily water. If 

a hydrocarbon fuel is used, carbon dioxide would also be present. If the fuel/air mixing scheme is not 

adequate or if the heater combustion length is insufficient for the given test condition, unburned fuel may 

also be present in the flow. As the simulated Mach number is increased, the higher temperature of the 

mixture results in an increasing mole fraction of combustion radicals, such as H, O, and OH, as well as 

NO and NO2, exiting the heater and entering the nozzle, in addition to the increasing fraction of water 

vapor and carbon dioxide. If a hydrocarbon fuel is used, the fraction of carbon monoxide present in the air 

stream also increases with increasing simulated Mach number. The rapid expansion of the heater outflow 

through the facility nozzle freezes the composition at a point downstream of the nozzle throat, so that the 

mixture is not allowed to equilibrate. Under some circumstances, water droplets may form as the water 

vapor becomes supercooled due to the expansion to high altitude flight conditions. This condensation is 

rarely seen at the test point, in practice. 

As the combustion heated air enters the engine inlet, the changes in the ratio of specific heats, γ, and 

the gas molecular weight alter the shock structure from what would be seen in clean air. A variety of other 

physical characteristics of the flow are also altered due to the flow vitiation, such as heat transfer rates, 

vibrational energy relaxation rates, turbulence, and inlet flow distortion, among others. As the flow enters 

the combustor and fuel is added, the combustion kinetics can be significantly altered, effecting both 

ignition and flame holding in ways that can be hard to predict, due to a variety of competing mechanisms. 

Lastly, nozzle expansion is also effected by the flow composition. As can be seen by the list above, 

combustion heater vitiation adds a large amount of uncertainty to the test results. Additional detail on 

these effects can be found in the summary papers of Powell and Stallings (1998), Pellett, et al. (2002) and 

Fry (2004). Additional difficulty is incurred with combustion heated facilities in terms of determining 

which primary variables should be matched during the testing. The changes in gas properties makes it 

impossible to match both total temperature and total enthalpy, both Mach number and velocity, or both

total pressure and static pressure, for instance. The experimentalist must evaluate which properties are 

key(s) for each test, a subject of much debate in the hypersonic engine testing community through the 

years.

If an arc heater is used, high levels of NOx are formed in the vicinity of the arc, which can persist in 

the mixed flow entering the nozzle at mole fractions above 3 percent. These compounds have little effect 

on the shock structure or other aerodynamic features, but can significantly catalyze the ignition process

(Slack and Grillo, 1977). 

Lastly, a storage heater generally represents the lowest level of chemical vitiation. Like the arc heater, 

a storage heater can generate some NOx, but not to the same level. At most, an equilibrium concentration 

of NOx might be generated, depending on the residence time of the air mixture at high temperature. As 

the thermal generation of NOx is fairly slow, it is possible that NOx levels below equilibrium might be 

present in some facilities, depending on the specifics of the facility configuration. 

The above discussion gives a brief overview of the gas-phase vitiation issues for hypersonic 

airbreathing propulsion ground test facilities. However, each of these facility types can produce

particulate vitiation as well. For the combustion heated facility, the particulate source is most likely soot
NASA/TM—2010-216765 3 

formed during the combustion process due to the existence of fuel rich regions in the hydrocarbon fuel/air 

mixing zone. Additional particulates may be present due to piping corrosion caused by hot water vapor. In 

an arc heated facility, electrode erosion releases a small amount of metallic particulates into the air 

stream. In a storage heater facility, eroded particles from the refractory material, a ceramic or graphite, are

found in the air stream. While the quantities and size distributions of these particulates have generally not 

been characterized, it can be assumed that the level of particulate vitiation is highest with storage heaters

due to the extremely large surface area of non-surface abrasion resistant material exposed to the flow. 

While the most recent surveys of the state of the art regarding the effects of vitiation (Fry, 2004; Pellett, et 

al., 2002) acknowledge the importance of particulate vitiation, only a single study has ever been 

conducted of the potential effects of particulates (Mitani, 1995). In contrast, the surveys of Fry (2004) and 

of Pellett, et al. (2002), reference dozens of past studies addressing various aspects of gas-phase vitiation, 

reflected in the summary of effects described above. 

Mitani (1995) performed a series of theoretical analyses to determine the effect of low micron-sized 

inert particles, assumed to be in thermal equilibrium with the carrier gas, on scramjet combustor ignition 

and flameholding. Mitani concluded that a mass fraction above 0.1 percent can inhibit ignition through 

the mechanism of radical termination on the particle surface. This mechanism is only valid for particles 

smaller than 3 μm in diameter. For the particles to have a thermal heat sink effect, the study concluded 

that a mass fraction in excess of 10 percent would be required. 

In order to look at the effects of larger particles that are not in equilibrium with the carrier gas, as well 

as reacting particles, throughout the test facility flow path, a numerical analysis is necessary. However, up 

to this time, no numerical modeling has ever been done to examine the potential effects of particle 

vitiation in hypersonic propulsion ground test facilities. Reasons for this include: 

1. The effects of particles have been conveniently assumed to be of secondary importance compared to 

the effects of gas phase vitiation. 

2. The numerical analysis tools typically used for hypersonics flows do not contain particle flow models.

The existing tools are multi-dimensional codes with finite-rate chemistry capability characterized by 

long run times, and so the addition of particle flow models is not attractive as it would tend to require 

even longer computational turn-around times.

3. Particle-gas flow interaction models that cover the broad range of conditions experienced in a 

hypersonic propulsion test facility are not readily available.

4. The mass loading and size distribution of particle vitiates has not been determined, so key inputs for 

model formulation are not available. 

While the first assumption, that particle vitiation effects are of secondary importance compared to 

gas-phase vitiation effects, may well be valid for combustion heated facilities, and perhaps arc-heated 

facilities as well, it is not a good assumption for storage heater-based facilities. In storage heater-based 

facilities, the amount of particle vitiation is likely to be substantially higher than in other types of 

facilities, and the amount of gas-phase vitiation is significantly lower. The final three reasons given above 

can be addressed, all or in part, through the development of appropriate particle-gas flow models and 

applying them to a quick turn-around numerical tool. By integrating the typically computationally 

demanding particle flow analysis into a more rapid analysis tool, parametric studies of particles of various 

sizes and other pertinent characteristics can be performed in a timely manner to understand under what 

conditions particle vitiation may be of concern.

In this work, a quasi-one-dimensional (Q1D) numerical modeling tool and associated sub-models are 

developed that allow for a relatively rapid exploration of the effects of particle vitiation throughout the 

operating envelope of a storage-heated hypersonic propulsion test facility. A wide range of particle sizes 

are considered for both chemically inert particles (alumina) and reacting particles (graphite). In order to 

appropriately define the ranges of the parameters to be modeled, both a test facility and a scramjet engine 

model (up through the combustor) are laid out as part of this study. The test facility model is a 

compromise between the geometry and flow characteristics of the two large scale storage heater scramjet 
NASA/TM—2010-216765 4 

test facilities found in the United States. This same model facility and model scramjet will then be used to 

demonstrate the use of the Q1D code for simulating hypersonic propulsion test flows with both inert and 

reacting particles.

The first storage heater based hypersonic propulsion test facility is the GASL Leg IV facility. The 

facility consists of an alumina matrix storage heater followed by hydrogen fueled booster heater. The 

storage heater can heat dry air to a total temperature of 1390 K, simulating approximately Mach 5, with 

the booster heater then adding thermal energy up to a Mach 7 condition (Roffe, et al., 1997). The basic 

configuration is shown in Figure 1. 

The second configuration to be studied is the NASA Glenn Research Center Hypersonic Tunnel 

Facility (HTF) located at the NASA Plumbrook Station in Sandusky, Ohio. This facility consists of a 

drilled-core graphite block heater capable of heating gaseous nitrogen to 2500 K. The heated nitrogen is 

then mixed with ambient temperature gaseous oxygen to provide simulated air to the facility nozzle at a 

total temperature up to 2170 K (Perkins, et al., 1998). This basic configuration is shown in Figure 2. 

Figure 1.—GASL Leg IV diagram.

Settling

Chamber

Pebble Bed 

Heater

Air Supply

Booster

Heater

Facility

Nozzle

Air

Shutoff

Valve

Diluent

Air Supply
NASA/TM—2010-216765 5 

Figure 2.—HTF diagram. 

While the GASL Leg IV facility is vitiated above Mach 5, the model facility will assume non-vitiated 

flow for both the inert and the reacting particles at all test Mach numbers, as vitiated flow for the inert 

particle case would significantly affect the nozzle flow and inlet shock structure, making comparisons 

difficult between the different particle types.

Neither GASL Leg IV nor HTF have ever been characterized as far as particle size and concentration 

as a function of test condition. Therefore, these particle properties will be dealt with parametrically when 

the code is exercised through a range of test conditions. There is no reason to assume that particle size is 

limited to a few microns diameter, as assumed by Mitani (1995). In the case of the HTF facility, particles 

on the order of 1 mm have been observed adhering to the heater graphite blocks during facility 

maintenance. Therefore, particle diameters in the range of 10 μm to 1 mm will be considered. Neither 

inert nor reacting particles with diameters less than 10 μm were used in this analysis as such small 

particles would be expected to track the gas flow very closely (Maxwell and Seasholtz, 1974) and, as will 

be shown in the following analysis, the 10 μm graphite particles almost immediately quench or are 

consumed, so a smaller particle would not be expected to show any different behavior. A particle mass 

fraction of 1 percent of the total test flow will be used as an upper bound, since a higher mass fraction 

would represent such a significant mass loss from the storage heater that heater refurbishment would 

likely be required after each test program. Such heater degradation has never been reported. While the 

results of this parametric study will not conclusively determine whether or not particle vitiation is of 

concern at any particular test point for a particular test facility, it will provide insight into what particle 

concentrations and sizes may be of concern for future research test programs.

Graphite

Heater

Blocks

Diluent 

Oxygen & 

Nitrogen 

Injection 

Flange

Facility

Nozzle

Nitrogen

Inlet

Heater

Isolation

Valve

Flow

Mixing

Section
NASA/TM—2010-216765 6 

In the following sections, the Q1D code and its sub-models will be described and the results of a 

series of validation cases based on the above range of conditions will be given. The code is then exercised 

through a parametric study of particle and gas flow combinations appropriate to the model test facility and 

model scramjet. Results are presented and conclusions are drawn as to the potential effects of particle 

vitiates. Finally, recommendations for further study are presented. 

2.0 Code Description

2.1 General Approach

The need to model a large number of facility and test hardware components over a wide range of 

parameter space, as well as the need to consider two phase flow with mass, momentum and energy 

interactions between the phases and finite rate combustion kinetics, limits the number of dimensions that 

can be practically used for the analysis. For this reason, a quasi-one-dimensional (Q1D) analysis approach 

was selected for the underlying gas-phase code. No existing code was available with the requisite 

characteristics for this study. A custom Eulerian (fixed grid), gas-phase, multispecies, time-accurate, Q1D 

computational fluid dynamics (CFD) code with finite rate kinetics was therefore written and validated. To 

this code was then added a Lagrangian particle tracking scheme with a global particle reaction model.

Since larger particles are considered, an internal one-dimensional grid is generated for each particle being 

tracked which is solved at each time step for the particle internal temperature gradient and the heat 

transfer to the gas phase. If the particles are sufficiently small such that the Biot number is expected to 

remain less than 0.1 (based on the particle radius as the characteristic dimension, (Incropera and Dewitt, 

2002) throughout the simulation, the particle is treated as isothermal. This is easily accomplished within 

the code by setting the number of internal particle grid points to two, as the number of shells within the 

particle is equal to one less than the number of internal grid points. Mass, momentum, energy and species 

source terms are included in the gas-phase governing equations to provide feedback from the particle 

drag, heat transfer, and reaction models. For inert particles, the particle reaction model is simply turned 

off. When the particle mass fraction is high enough that the number of particles being tracked would be 

too high for efficient and timely computation, each particle tracked is defined as being a representative 

particle for a larger number of particles at that same location. The number of particles represented by each 

tracked particle is provided as a user input. The mass, momentum and energy contributions from each 

tracked particle are then multiplied by the number of particles in each set to get the gas-phase source 

terms. 

This modeling approach allows for a large number of computations to be performed, but with certain 

drawbacks. Boundary layer related effects cannot be captured, nor can oblique shocks. Perhaps the 

biggest drawback for this problem is that the lateral diffusion of the particles when passing through a 

divergent section at high flow velocity is not captured. Previous studies have shown that when two-phase 

flow passes through a convergent-divergent nozzle, the particle stream will “separate” from the wall, 

creating a particle-free zone near the wall (Chang, 1980; Ishii, et al., 1987). The concentration of particles 

therefore increases in the rest of the flow. This effect is more pronounced as particle size is increased.

However, since this study is parametric in nature, this effect can be accounted for when the results are 

applied to a specific configuration by correcting the particle concentration based on the nozzle flow 

conditions and the particle size under consideration. 

2.2 Gas Phase Governing Equations

The gas phase governing equations were adapted from Kuo (1986), re-derived to include area 

variation. Viscous dissipation and diffusion terms determined to be unimportant for these simulations are 

omitted. The φ term represents the local volume fraction of gas in the gas-solid mixture. The source 

terms, Sm, Sp, Se, come from the solid phase relations, while the species source terms, Sm,i, are a 

summation of both the gas phase and solid phase reactions. 
NASA/TM—2010-216765 7 

Continuity: 

( ) g ( ) gg ASAu m x

A

t

φ+φρ ∂

∂ −=φρ ∂∂ (1) 

Momentum: 

( ) ( ) x

A pASpAuA

x

uA

t gg gg p ∂

∂ φ+φ+φ+φρ ∂∂ −=φρ ∂∂ 2 (2) 

Energy: 

( ) ( ) φ+φ+φ+φρ π

∂

∂ −=φρ ∂∂ − AqASAqupAeuA

x

eA

t gg ggg cg we 4 (3) 

where 

x

T

kq g

gc ∂

∂ −= (4) 

and 

−= RThe ggg (5) 

Species Continuity: 

( ) ( ) ig igig ASAvuY im

x

AY

t , )( φ+φ+ρ ∂

∂ −=φρ ∂∂ (6) 

where 

x

Y

Y

D

v i

i

i

i ∂

∂ −= (7) 

The qw wall heat transfer term is currently a user input that can be used to “tune” the computational 

results to match available experimental data. If data are unavailable, an approach based on the ReynoldsColburn analogy relating wall friction to heat transfer has been shown to be effective for Q1D 

calculations (Paxson, 1993) and could be readily incorporated. 

2.3 Thermodynamic Properties

The basic thermodynamic properties for each species, enthalpy (hg,i), enthalpy of formation (hto,i), 

specific heat (cp,i), Gibbs energy of formation (gf,i

o

) and molecular weight (Mi), were all taken from the 

JANNAF Thermochemical Tables (Chase, et al., 1986). The temperature dependent properties were 

piece-wise curve fit over four temperature ranges using third-order polynomials. The four temperature 

ranges were 0 to 300 K, 300 to 600 K, 600 to 1100 K, and 1100 to 3000 K. The use of four temperature 

ranges provided high accuracy over the entire temperature range, typically within 0.1 percent of the 

original JANNAF data. From this data, along with the mass fractions of each species, the mixture 

properties can be defined, as follows. 
NASA/TM—2010-216765 8 

∑ ( ) =

= +

ns

i

igig hhYh to i

1

,, (8) 

∑

=

=

ns

i

CYC ipip

1

, (9) 

∑

=

=

ns

i

o

ifi of gYg

1

, (10) 

For some of the calculations, the mole fractions of the chemical species are needed instead of mass 

fractions, which can be easily calculated as follows.

∑

=

= ns

i i

i

i

i

i

M

Y

M

Y

X

1

 (11) 

Once the mole fractions are calculated, the mixture molecular weight can be obtained.

∑

=

=

ns

i

MXM ii

1

 (12) 

This leads to the simple relation for the mixture gas constant, given the universal gas constant, 

M

R R u = , (13) 

and the ratio of specific heats can then defined as

RC

C

p

p

− =γ . (14) 

2.4 Transport Properties

Species effective collision diameter (σi) and effective temperature (Tε,i) data, taken from Svehla 

(1962) and Hirschfelder, Curtiss, and Bird (1964), were combined with the Prandtl number (0.706) to 

generate the transport properties used in the calculations. Following the methodology found in White 

(1991), the binary diffusion coefficients are first calculated.

ji Dji

ji

ji

g

i j P

MM

MM

T

D

,

2 ,

2

1

2

3

,

001858.0

Ωσ















 +

= (15) 

Where the collision integral is defined as
NASA/TM—2010-216765 9 

( ) 2

,

*145.0

,

*

, 5.0

− − =Ω ++ D ji g ji g ji T T (16) 

and where 

ji T

T

T g

g ji ,

,

*

ε

= (17) 

( ) 2

1

, jiji = TTT εεε (18) 

and 

( ) σ+σ=σ jiji 2

1 , . (19) 

Once the binary diffusion coefficients have been calculated, an “average” diffusion coefficient for 

each species diffusing into the rest of the mixture can be calculated, using the following relation from 

Warnatz, et al. (1999). 

∑≠=

− = ns

ijj ji

j

i

i

D

X

Y D

,1 ,

1 (20) 

While the use of an average diffusion coefficient is somewhat less accurate than using a full multispecies 

diffusion model, it is much less computationally demanding and was deemed sufficiently accurate (~10

percent) for the purposes of this study. 

To calculate the individual species viscosity, the method described in Hirschfelder, Curtiss, and Bird 

(1964) is adopted. 

( )

2

2,2

2

1 993.26

i

gi

i

TM

σΩ

=μ (21) 

where the collision integral, Ω2,2, is represented by the curve fit 

( ) ( ) ( ) *1448.2 2,2 *3856.0 2,2 *0161.0 2,2 2,2 5715.25015.0 6391.0 3546.0 T T T e e e − − − +=Ω + + (22) 

and 

i

g

T

T

T ε

* 2,2 = (23) 

The mixture viscosity is then calculated using Wilke’s law, as given by Drummond (1990). 

∑

∑ 















+ φ

μ =μ

≠=

ns

ij

ns

ijj

j

i

i

X

X

1

,1

1 1

(24) 
NASA/TM—2010-216765 10

where

5.0

2 25.05.0

1

2

4

1

















+

















































ρ

ρ

μ

μ +

=φ

j

i

j

i

j

i

j

i

ij

M

M

M

M

(25) 

Lastly, the mixture thermal conductivity is simply calculated by the standard expression 

Pr

μ = p

g

C

k . (26) 

All of the above expressions assume laminar flow.

2.5 Finite-Rate Kinetics

Each of the N chemical reactions used in the gas-phase reaction mechanism, shown symbolically in 

Equation (27), is modeled using the standard Arrhenius expression for the forward reaction rate, 

Equation (28). 

i

N

i

i i

N

i

i S ∑∑ S

= =

′′ ν ′ ν⇔

1 1

(27) 

gu

a

TR

E

n

gff eTAk

−

= (28) 

The reverse reaction rate can be found from the relation

eq

f

b K

k

k = (29) 

where the equilibrium constant, Keq, is calculated from the Gibbs energy of formation for each species and 

the reaction coefficients from Equation (27), following the procedure found in Kuo (1986). 

For this study, it was necessary to assemble a chemical kinetic mechanism that included H, O, N, and 

C compounds in order to capture the basic hydrogen combustion process for a scramjet combustor as well 

as the effects of the flow vitiates of interest. Based on a review of the available literature on gas-phase 

vitiation, as outlined in the introduction, the H-O species H2, O2, H2O, OH, H, O, HO2, and H2O2 were 

selected. The nitrogen species N, HNO, NO and NO2 were added to the basic diluent N2 to capture 

dissociation and vitiation effects. Lastly, CO, CO2, and HCO were included as gas-phase products from 

the combustion of the carbon particles under study. The mechanism developed by Davis, et al. (2005), for 

H2-CO combustion was adopted for the bulk of the mechanism used for this study. The nitrogen 

chemistry was added from the GRI 3.0 mechanism (Smith, et al. 2005). The details of the combined 

mechanism are given in Appendix B. 

Once obtained, the forward and reverse reaction rate constants are used in the standard law of mass 

action expression to obtain the change in moles of each species for each reaction. These molar changes 

are then summed across all of the reactions, and then utilized as source terms within the governing 

equations. 
NASA/TM—2010-216765 11

2.6 Solid Phase Governing Equations and Numerical Solution Method

There are two fundamental approaches that can be taken when formulating the governing equations 

for the solid phase in a two-phase computation. The first is to assume that the particles form a continuum 

that can be treated in the same manner as the gas phase. A parallel set of conservation equations to the gas 

phase equations can then be set up for the solid phase, simply exchanging (1–φ) for φ as the volume 

fraction occupied by the flow, and making the source terms of opposite sign in Equations (1) through (3). 

This is the approach, referred to as Eulerian, taken by Kuo (1986), Ludwig and Roth (1997), and Di 

Giacinto et al. (1982), among many others, and is typical of computations of fluidized beds. 

Alternately, each particle can be tracked individually in a Lagrangian frame of reference. In this 

formulation, position, velocity, mass, and temperature of each particle are updated at each time step based 

on mass, momentum and energy interactions with the gas phase and any surface reactions that may occur.

This approach is typically used for very dilute particle concentrations with relatively small computational 

domains where the computational load is not unreasonable. This method was adopted by Egolfopoulos 

and Campbell (1999), and Zhou et al. (2002) in order to capture physical phenomena that could not be 

easily examined using the more traditional Eulerian approach. 

For this study, the Lagrangian formulation was adopted as the mass loading of particles was expected 

to be relatively low, and the consideration of larger particles made the continuum assumption required for 

the Eulerian approach appear questionable. Since these calculations are Quasi-1D, the calculations are 

additionally simplified by the lack of radial variation. This allows a single particle to be used 

computationally to represent all the particles at a given axial position, significantly reducing the number 

of particles that need to be analyzed at each time step. Particles are therefore inserted into the 

computational domain at the inflow boundary as sets of particles, with the number of particles in a set 

determined by the particle size and mass loading, for which only a single particle calculation must be 

performed per time step. The gas-particle interaction terms in the governing equations are then the results 

of the single particle calculation for a particle set at a given axial location multiplied by the number of 

particles in that particular set. This is particularly helpful for very small particles, where the number of 

particles can become very large even at small mass loadings. Once the gas-particle interaction terms have 

been calculated for a given set of particles, these values are linearly interpolated onto the gas-phase fixed 

grid for used in the gas-phase solution process. 

After a particle (or set of particles, hereafter referred to as though it were a single particle) is injected 

into the gas-phase flow with an initial mass, position, velocity, and temperature, a series of relations for 

spherical particles are used to determine the particle drag, heat transfer, and surface reaction rate that is to 

be applied to the particle. 

The first gas-particle interaction is aerodynamic drag. For laminar flow, White (1991) provides the 

following drag relation for a sphere. 

( ) 2

2

1

ppgpggDp π−−ρ= ruuuuCF (30) 

Many empirical relations have been proposed for the drag coefficient of a sphere. Carlson and 

Hoglund (1964) and Crowe (1967) developed correlations, based on the data available at the time, which 

have been used frequently for particle flow simulations. More current efforts, such as that of Igra and 

Takayama (1993) have provided more accurate expressions over a specific flow properties range based on 

detailed shock tube data. The drag coefficient used for this study is taken from Henderson (1976). This 

expression covers a very broad range of Reynolds number and Mach number with reasonable accuracy. 

As the relative velocity between the flow and the particles in this study could range from supersonic to 

near zero, and the static pressure from 80 atmospheres down to 0.1 atmospheres or less, it is important to 

cover the broadest possible range of the relative Mach and Reynolds numbers. The Henderson drag 

coefficient expression was developed from the extensive experimental data base compiled by Bailey and 

Hiatt (1971) through a series of ballistic range experiments and is accurate to within 16 percent across the 
NASA/TM—2010-216765 12

entire range considered. The expression is divided into three elements to cover the entire relative Mach 

number range. For a subsonic relative Mach number,

( )

e MSR

e M M

e

T

T

T

T

C MSR

rel

rel

rel

rel

rel

M

rel rel

rel rel

rel rel

M

MSR

g

p

g

p

D rel

1 6.0

2.01.0

03.01 Re 48.0 Re

03.038.05.4 Re 48.0 Re

353.01

165.3 53.

24 Re 33.4

Re

Re 2 8

5.0

1

Re .0 247





















−+

















+ +

+ +

+ +

+





















































































+

−

+= +















 −















 −

−





 





−

(31) 

where the molecular speed ratio is defined as

2

γ MSR = Mrel (32) 

For a supersonic relative Mach number greater than 1.75,

2

1

4

2

1

2

2

1

2

Re

86.11

058.12 1 2

Re

86.1 34.0 9.0

















+





















− 















++ 















++

=

rel

rel

g

p

rel

rel

rel

D

M

T MSR

T

MSR MSR

M

M

C (33) 

In both the subsonic and supersonic equations, the relative Mach and Reynolds numbers are evaluated 

at free stream conditions. If the relative Mach number lies between 1.0 and 1.75, then the drag coefficient 

is taken as a linear interpolation between the value calculated using the subsonic Equation (31) with a 

Mach number of 1.0 and the value calculated using the supersonic Equation (33) with a Mach number of 

1.75, expressed as,

( ) ( ) ( ) ( ) ( ) ( ) D MC relrel CD rel Mrel CD rel CD Rerel ,75.10.1 Re ,0.1 3

4 ,Re = ,0.1 Re −+ − (34) 

The relative Reynolds number is evaluated at the actual free stream conditions for use in Equation (34).

The gas-phase properties used in these equations are evaluated at the particle location, linearly 

interpolated between the nearest adjacent grid points. 

Once the force on the particle is calculated, the particle position and velocity can be updated by the 

simple relations,

tuxx i

p

i

p

i

p ∆+= +1 (35) 

and 
NASA/TM—2010-216765 13

t

m

F

uu

i

i

p

i p

p

i

p ∆+= +1 . (36) 

Performing the energy balance on a large reacting particle in a convective flow has a number of 

complexities that need to be considered. Since the particle can be relatively large, a uniform particle 

temperature cannot be assumed, and an internal grid must be used to capture the temperature gradient 

within the particle. Since significant velocity slip is expected between the particle and the gas, convection 

must be modeled. Radiation losses must also be modeled since the gas surrounding the very hot, high 

emissivity particles is expected to be optically thin, and the test facility walls are water cooled. And, of 

course, since the particles are reacting, a large amount of heat is input to this system via combustion. 

If a thin shell at the outer surface of the particle is considered, as shown in Figure 3, an energy 

conservation Equation (37) can be written for that shell that encompasses each of the heat flows described 

above, and which allows for the accumulation or depletion of heat energy in the outer shell. 

accum comb rad conv −−−= qqqqq cond (37) 

Heat energy will then be conducted into or out of the particle interior over time at a rate controlled by 

the thermal diffusivity of the particle material. For the conduction, convection, and radiation terms, 

simple textbook relations are used, such as those found in Incropera and DeWitt (2002). 

s

cond pS dr

dT rkq 



 





π−=

2 4 (38) 

Figure 3.—Particle heat transfer and reaction model. 

rp

qconv

qrad

qcond

qcomb

qaccum

1-D Radial Grid

∆r 2/

∆r

Outer Shell
NASA/TM—2010-216765 14

( ) 442

qrad − 4 ppBSrad −πσε−= TTr wall (39) 

conv ( ) −π−= TTrhq pgpp

2 4 (40) 

where

p

gd

p r

k

h

2

Nu = (41) 

As with the drag coefficient, a number of empirical expressions for the Nusselt number have been 

developed for spheres over a variety of flow conditions, characterized by the Mach number and the 

Reynolds number. Once again, an expression that covers a wide range of these parameters is necessary for 

this study, including the case of supersonic relative velocity. In the case of supersonic relative velocity, 

the effect of the shock wave in front of the particle needs to be considered. Drake and Backer (1952) 

determined that it is best to correlate the Nusselt number with the conditions immediately behind the 

shock wave in line with the direction of travel. Since this is essentially a standing normal shock, simple 

textbook relations, such as those found in Saad (1985), can be used to calculate the pressure, temperature, 

and density behind the shock wave. After adjusting the thermodynamic and transport properties to the 

post-shock conditions, the post-shock Mach and Reynolds numbers can be calculated, and an appropriate 

correlation, such as that proposed by Fox, et al.(1977), given below, can be used. 

























+

+

























+

=

−

−

5.1

5.01 459.0 RePr

Re

171

Nu 2

Re

17

33.0 55.0 rel

rel

rel

M

rel

rel

rel

M

d

e

M

e (42) 

For high speed flows, it is also necessary to calculate what is commonly referred to as the recovery 

temperature (White, 1988), given below for laminar flow. 







 





 



 



 − += 5.0 2

2

1 1 Pr rec TT gas Mrel

γ (43) 

This value is then substituted into the convection heat transfer Equation (40), in place of the free stream 

gas temperature, to determine the convective heat transfer term.

)(4 2 conv pp rec −π−= TTrhq p (44) 

The recovery temperature is a modified form of the total temperature which approximates the temperature 

in the boundary layer surrounding the particle. This temperature is based on the free stream temperature 

and free stream relative Mach number, even in the case of a supersonic relative Mach number. 

The calculation of the combustion heat addition term, qcomb, is not nearly so straight forward. A 

number of carbon particle combustion models have been developed, such as those of Blake (2002); Hurt, 

et al. (1998); Lee, et al. (1995); Makino and Law (1986); and Libby (1980). These models were typically 

developed to examine the oxidation of coal char. In addition to simplifications such as zero gas-particle 

velocity slip and uniform particle temperature, each of these studies was limited to combustion at lower 

pressures; typically atmospheric pressure. However, for the problem under consideration in this study, a 

large range of pressures must be considered. Only two studies were found that looked at carbon 
NASA/TM—2010-216765 15

combustion at pressures considerably above atmospheric pressure in a convective flow, those of Moors 

(1998) and Hong (2000). The methodology of Hong was adopted for this study. 

Calculation of the combustion rate for solid carbon is somewhat complicated by the fact that the 

products of reaction shift as a function of the reaction temperature. As the carbon surface temperature 

increases, the reaction shifts from 

C s O ⇒+ CO22)( (45) 

to 

O CO

2

1 C s 2)( ⇒+ (46) 

Hong provides the following Arrhenius-type relation, taken from Hurt and Mitchell (1992), for the mole 

ratio of CO to CO2. 

MR e 987.1 Ts

30000

CO

CO 40000

moles

moles

2

−

= = (47) 

For convenience, this relation can be recast to give the ratio of moles of CO2 produced by the surface 

reaction to the total moles produced. 

+ MR =ψ1

1 (48) 

Two other convenient parameters can be defined at this point that will be useful for simplifying the 

particle combustion equations. These are the stoichiometric coefficient of oxygen consumed for each 

mole of carbon consumed, ν0, and an intermediate variable, ξ. 

0 ( ) 15.0 ψ+=ν (49) 

1

1

+ψ

−ψ =ξ (50) 

Calculation of the combustion rate for a solid carbon particle is also complicated by the structure of 

the particle itself. For a typical carbon particle, the majority of the reactive surface area is made up of a 

labyrinth of micron to nanometer diameter pores internal to the particle. Oxygen must diffuse from the 

surface into the particle pore structure for significant reaction to take place. This process must be included 

in order to develop an accurate reaction rate model. This internal diffusion effect is captured through the 

use of an effectiveness factor, η, defined by Hong as the ratio of the actual observed reaction rate 

including intra-particle diffusion effects to what the reaction rate would be if the oxygen concentration at 

the particle surface persisted throughout the interior pore structure of the particle. The effectiveness factor 

is calculated as follows.

( ) 



 





=η − T TT c M MM f 3

1

tanh 3

11 (51) 

where
NASA/TM—2010-216765 16

( )2 1

2

1

2

2 2

2

1

2

1

1

mobs

T

T

c

M

M

f

−





























+

+= (52) 

The Thiele Modulus, MT, is defined as

2

1

,O

,O

10

2

2 1

1 2 3 























+

+





 





ν

=

sp

sp

e

g

up

p

T

PK

PK

D

T Rk

r

M (53) 

and the observed reaction order, mobs, is 

sp

obs PK

m

1 2 ,O

1

+ = (54) 

The above expressions contain two kinetic rates for which Hong gives the following Arrhenius-type 

relations for non-porous graphite flakes, based on the experimental data of Ranish and Walker (1993). 

TR ipu k p e ,

51000

8 1 1029.6

−

×= (55) 

TR ipu p eK ,

10100

4.13

−

= (56) 

The constants contained in these kinetic rate expressions must be re-derived from experimental data for 

each individual type of graphite or coal char, as they are dependent upon the physical structure of the 

particles as well as their chemical composition.

The Thiele Modulus expression also requires a value for the effective diffusion coefficient, De. The 

effective diffusion coefficient is a function of the particle porosity, ε (1.0 minus the ratio of the actual 

particle density to the density of a non-porous solid of the same material), the molecular diffusion 

coefficient of oxygen into the gas mixture at the local gas conditions, D 2 ,O M , and the Knudsen diffusion 

coefficient, Dk (for diffusion in small pores where wall collisions dominate over inter-molecular 

collisions). Hong gives the following relation for the combination, 

kM

e

DD

D 11

,O

2

2

+

ε = (57) 

where

2

1

O

,

2

9700 













 = M

T

rD ip

k p (58) 
NASA/TM—2010-216765 17

in units of cm2/s, as taken from Smith (1981). The value of D 2 ,O M is available from the gas-phase 

portion of the Q1D code. 

In the formulation of Hong, a uniform particle temperature was assumed, so a single value of the 

effectiveness factor was used. In this study, it was necessary to use an approach that accounted for an 

internal particle temperature distribution in what is otherwise a global reaction model. The use of a mass

averaged value of the effectiveness factor, equivalent to an internal surface area averaged value, was 

selected as an appropriate modeling approach. While no data was available to validate this modeling 

choice, this commonly used approach provides reasonable approximate values and reduces to the 

approach of Hong in the case of uniform particle temperature. The general behavior of the effectiveness 

factor as a function of temperature is shown in Figure 4, using the kinetic parameters for graphite flakes 

given above. 

In order to solve this complex system for the carbon burning rate, it is necessary to know both the 

overall particle oxygen consumption rate and the particle temperature profile. However, these values are 

inter-dependent, therefore requiring that the carbon burning rate be solved iteratively. Since the other heat 

flow terms from Equations (38), (39), and (44) depend on the particle surface temperature, and therefore 

the particle temperature profile, the energy balance given in Equation (37) must also be solved iteratively.

This is done by running an inner iterative loop that corrects the particle surface partial pressure of oxygen 

(at the particle surface temperature supplied from an outer iterative loop) by comparing the O2 diffusion 

rate from the free stream and the total carbon consumption rate of the particle based on the graphite 

kinetics, taking into account the temperature dependent stoichiometric coefficients described above 

(Eqs. (45) to (47)). An outer iterative loop is run that corrects the particle surface temperature by 

comparing the outer shell carbon consumption rate calculated by an energy balance on the outer shell to 

the same rate calculated using the carbon consumption kinetic expression using the updated particle 

surface partial pressure of oxygen from the inner loop. 

0.1

1

10

0 500 1000 1500 2000 2500

η

Temperature (K)

Figure 4.—Effectiveness factor as a function of temperature. 
NASA/TM—2010-216765 18

The procedure for this solution is as follows

1. Assume a value for the particle surface temperature. This is typically selected as the value from 

the previous time step.

2. Assume a value for the oxygen concentration at the particle surface. A value of 90 percent of the 

free stream oxygen concentration is used to start each iterative cycle.

3. Calculate the global kinetic carbon consumption rate, rk , for the given values of the particle 

surface temperature, the particle interior temperature distribution, and the surface oxygen 

concentration. An intrinsic Langmuir expression is used to calculate the local volumetric 

consumption rate within each shell of the particle. This expression, adapted from Hong, is valid 

over a wide range pressures, temperatures and porosities. 

Sip

cSip

tot

ext

avgik PK

MPk

S

S

S

S

r

,O,

int ,O,1

int

,

2

2

 1+











+η= (59) 

The natural units for this expression is g/cm3

-s, based on the particle volume, but it is 

typically given in terms of the particle external surface area, g/cm2

-s. This transformation is easily 

done by multiplying the rate by the particle volume and then dividing by the external surface 

area. In either case, the units of k1p and Kp are mole/cm3

-atm-s and 1/atm, respectively. The 

kinetic rates in Equation (59), k1p and Kp, are the same terms used in the Thiele Modulus 

expression. The interior, exterior, and total particle surface areas, Sint, Sext, and Stot, are used in the 

expression to proportion the amount of reaction that takes place on the particle exterior, and is 

therefore not affected by internal pore diffusion, versus that which occurs in the particle interior. 

It should be noted that the external surface area, Sext, is typically taken as the particle’s geometric 

external area multiplied by a roughness factor to account for surface irregularities. Hong suggests 

a value of 5 be used for the roughness factor. Using available data for porosity and pore diameter, 

and assuming a monodisperse pore size distribution within the particle, the particle interior 

surface area is readily calculated.

pore

p

d

r

S

3

int

16πε = (60) 

The local carbon consumption rate for each particle shell volume is calculated by multiplying rk,i

by the volume of the particle shell being evaluated. The global carbon consumption rate is then 

obtained by summing the individual shell rates. This value is then normalized by the particle’s 

geometric external surface area, so as to be comparable to the other rates used in the iterative 

procedure, which are expressed per unit external surface area.

4. Next, the rate of carbon consumption based on diffusion of O2 from the free stream to the particle 

surface must be calculated. The mass transfer coefficient for O2 diffusing to the surface, recast in 

terms of carbon consumption, can be expressed as, 

up film o

mC

D TRr

DM

K

ν = 2

,O Sh 2 (61) 

where

film 5.0 ( ) += TTT gp (62) 
NASA/TM—2010-216765 19

and where Sh is the Sherwood number, which is calculated from Equation (42) by substituting the 

Schmidt number for Prandtl number. A Schmidt number of 0.74 is used for this study. Then the 

rate of carbon consumption based on oxygen diffusion to the surface is given by, 





























ξ−

ξ−

ξ =

P

P

P

P

PK

r

fs

s

D

OX

,O

,O

2

2

1

1

ln ≠ξ )0( (63) 

or 

OX D ( fs PPKr ,O,O s ) 2 2 = − =ξ )0( (64) 

with the pressure terms in the above expressions given in atmospheres.

Assuming that the reaction and diffusion processes are in a quasi-steady state, the carbon 

consumption rate from Equation (63) or (64) and the rate of carbon consumption calculated from 

the kinetic expressions above should be equal. If they are not initially equal, the surface partial 

pressure of O2, P ,O s 2 , is iterated upon to bring these two rates into equilibrium. 

5. The carbon consumption-based O2 diffusion and O2 consumption rates are now in equilibrium at 

the particle surface temperature assumed in step 1. An energy balance of the particle shell using 

the terms shown in Figure 3 is then used to adjust the particle surface temperature as follows.

First, (37) is rearranged to solve for qcomb. 

comb accum rad conv +++= qqqqq cond (65) 

where 

)( 3 2

4 1

3

3 i

p

i

pp ps s accum TTC r

q rr − 



















 



 ∆ −− πρ = + . (66) 

Then we can write a carbon consumption rate expression for the outer shell based on the rate that

would be required to produce qcomb. 

( ) CO CO2 1 HH

q r comb e ∆ψ+∆ψ− = (67) 

If re is not equal to the kinetic rate for the outer shell calculated using Equation (59) within a 

given convergence criteria, Tp is adjusted upward or downward, as appropriate, and the procedure 

is repeated starting at step 2 until convergence is achieved. 

Now that the carbon consumption rate and particle outer shell energy balance have been calculated, 

the next step is to update the particle mass to account for the mass lost due to the interior and exterior 

surface reactions. Mass can be lost from both the exterior surface of the particle, reducing the particle 

diameter, and from the pore surfaces internal to the particle, increasing the diameter of the pores. To 

simplify the analysis for this study, it was decided to adjust the particle mass by increasing the pore 

diameter within the particle. This approach was selected because of the high porosity of the graphite 

material being studied and the difficulty in establishing a proper mass decrease to diameter decrease ratio.
NASA/TM—2010-216765 20

According to Hong, the relationship between mass decrease and particle diameter is generally determined 

experimentally, and this ratio changes during the combustion history of a particle and with ambient 

conditions. By changing the pore diameter only, a computational efficiency is obtained in that the particle 

grid does not have to be adjusted after each iteration. 

Once this has been done, the source terms for overall mass, momentum (including the momentum 

associated with the mass lost by the particle), energy, and species (CO, CO2, and O2) are calculated for 

use in Equations (1), (2), (3), and (6). These source terms are multiplied by the number of particles in 

each particle set, and then linearly interpolated onto the two adjacent grid points for use in the gas phase 

governing equations. 

Lastly, the internal particle temperature profile is updated using the transient, 1-D conduction 

equation in radial coordinates.

t

T C

r

T rK

rr

S SS ∂

∂  =



 





∂

∂

∂

∂

ρ 2

2

1 (68) 

The simple, 2nd-order accurate Crank-Nicolson semi-implicit method (Tannehill, et al., 1997) is used 

to solve the conduction equation. This method was selected because it is unconditionally stable and does 

not limit the time step used in the overall calculation. 

2.7 Gas Phase Numerical Solution Method 

The gas-phase governing equations (Eqs. (1), (2), (3), and (6)) are cast in the following form for 

solution

J

x

F

t

U

+

∂

∂ −= ∂∂ (69) 

where





























φρ

φρ

φρ

φρ

=

ig

gg

gg

g

YA

eA

uA

A

U (70) 

( ) 



























+φρ

φρ+φ+φ

φρ+φ

φρ

=

iigg

c gggg

gg

gg

YvuA

euAAPuAq

uAPA

uA

F

2

(71) 

and 





































φ

π−φ





 





∂

∂ φ+φ

φ

=

im

we

m

DA

AqDA

x

A A Dp P

DA

J

,

4

(72) 
NASA/TM—2010-216765 21

Two different schemes are available in the code for solving this set of equations. The first scheme is 

the standard MacCormack predictor-corrector scheme, as described in Anderson (1995). This solution 

method is 2nd order accurate in space and time, and follows the following steps. 

Predictor Step

1. The values in the F and J arrays are first calculated from the values contained in the U array from 

the previous time step.

2. Equation (69) is then used to obtain the time derivative of the U array using a forward difference 

in x. 

3. The U array is then updated to an interim value by the equation, 

 χ+∆

 





∂

∂

+= t

t

U UU j (73) 

where j refers to the time step and χ is a Total Variation Diminishing (TVD) artificial viscosity 

term that will be discussed later.

Corrector Step

1. The values in the F and J arrays are updated from values contained in the U array. 

2. Equation (69) is again used to obtain the time derivative of the U vector using a backward 

difference in x. 

3. An average of the two time derivatives is calculated as follows,





 











 







∂

∂  +



 





∂

∂  = 

 





∂

∂

t

U

t

U

t

U

avg 2

1 (74) 

4. Lastly, the U array is updated to the next time step by the following expression, 

 χ+∆ 

 





∂

∂

+= + t

t

U UU

avg

1 jj (75) 

This scheme can be easily improved to 4th order accuracy in space by making the following 

substitutions (Georgiadis 2000), with k representing the spatial grid coordinate.

1. Calculate the time derivative of U from predictor step 2 by the expression, 

( ) JFFF

xt

U

kkk +−+−

∆ −= ∂

∂ 87 ++ 21 6

1 (76) 

2. Calculate the time derivative of U from corrector step 2 by the expression, 

( ) JFFF

xt

U

kkk ++−

∆ −= ∂

∂ 87 −− 21 6

1 (77) 

The user can select either scheme within the code.
NASA/TM—2010-216765 22

Numerical viscosity, χ,which is generally required when using a MacCormack type scheme 

(Anderson 1995), is combined with a TVD routine in an easily implemented method developed by Davis 

(1987). The basic form of the TVD artificial viscosity term is 

( ) ( ) ( )−− 

  − + − + + + + =χ + + kk k kk k UUrKrK

1 1 2

1

2

1

[ ( ) ( )]( ) 2 1

1

2 1 1 − − − − − −− + ++

kk k k k k UUrKrK (78) 

where

( ) ( ){ })(15.0 ±± ± CrK θ−ν= r c (79) 

The leading variable in this expression is defined as 

( ) 



 ν−ν

=ν

25.0

)1( cc C c if

if

5.0

5.0

>ν

≤ν

c

c (80) 

where the Courant number, νc, is defined by 

x

t c k ∆∆

max λ=ν (81) 

and the eigenvalue at each point is defined as 

, +=λ au kkgk (82) 

The expression for the flux limiter )( ± θ r is given by

( ) ( )





 =θ

±

±

0

min r 1,2

r

if

if

0

0

<

>

±

±

r

r (83) 

where r

+

 and r

–

 are defined by the following inner products, 





 



 ∆∆





 



 ∆∆

=

+ +

− + +

2

1

2

1

2

1

2

1

,

,

k k

k k

i

UU

UU

r (84) 





 



 ∆∆





 



 ∆∆

=

− −

− + −

2

1

2

1

2

1

2

1

,

,

k k

k k

i

UU

UU

r (85) 

and ∆U represents a backward difference in x. 

The gas-phase solution method also requires boundary conditions. These boundary condition may 

need to represent subsonic, choked, or supersonic flows, depending on the simulation. For supersonic and 

choked inflow conditions, the velocity or Mach number is set along with either the static or total pressure 
NASA/TM—2010-216765 23

and temperature. For supersonic and choked outflow conditions, simple extrapolated boundary conditions 

are used. Subsonic boundary conditions are inherently more challenging, since information can travel 

both upstream and downstream. The methodology of characteristic boundary conditions of Poinsot and 

Lele (1992), valid for both the Euler equations as well as the Navier-Stokes equations, has been adopted 

for these conditions. The boundary condition can be pre-defined, or can be selected by the code as a 

function of the flow conditions as the calculation progresses. 

The computational grid used for these calculations is generally uniform in the axial coordinate x, with 

grid spacings ranging from approximately 4 mm down to 1 mm. Since the calculations are quasi-1D, the 

number of grid points required for any given calculation is typically reasonable, from 300 up to 4000. 

Departures from a uniform grid are pointed out where they occur. 

The time step used in these calculations was allowed to vary during the course of the computation to 

the largest value that did not exceed any of three limits put in place to ensure stability and accuracy of the 

solution. The first limit is the standard Courant, Friedrichs, and Lewy (CFL) stability limit for explicit 

numerical methods (Tannehill, et al., 1997). A time step limit of 50 percent of the maximum time step 

allowed by the CFL criterion was used to ensure stability. A gas-phase chemical reaction rate limit to the 

time step was also found to be necessary to prevent a chemical heat release “runaway” under some 

conditions. A limit of 1 percent change to any species mass fraction at any grid point was typically used.

Lastly, it was found during validation testing of the code that a time step limit based on the rate of change 

of the outer shell temperature of the particles was necessary to maintain accuracy. This limit was typically 

set to 0.02 K temperature change in the outer shell of any particle. The smallest of these three limits is 

selected at each time step as the value to be used for the succeeding time step. 

3.0 Validation

For the purposes of validating the above described Q1D 2-phase code, no single representative test 

case of known solution was found. Instead, different aspects of the code were validated individually using 

available solutions. Individual physical process validations were performed in addition to the higher-level 

validations included in these sections. These physical process validations included chemical equilibrium 

calculations, shock tube calculations, individual particle fixed outer temperature boundary condition 

calculations, individual particle constant heat flux calculations, rocket performance calculations, steady 

state nozzle flow calculations, and various carbon combustion calculations. Additionally, many test cases 

were run to determine maximum acceptable time steps and grid sizes and to determine appropriate 

convergence criteria for a wide range of problems. 

3.1 Gas Phase Flow and Chemical Kinetics Validation

As a first step, the gas-phase reacting flow capability was validated using a 1–D detonation 

simulation. The full kinetic mechanism was exercised by using a reactant mixture of 10 percent CO, 

3 percent H2, 32 percent O2 and 55 percent N2 by mass. As shown in Figure 5, the detonation was initiated 

by a small region of high pressure, high temperature H2 at the closed end of a tube filled with reactants at 

300 K and 100 kPa. The ignition region is initially at 3000 K and 4 MPa. The detonation wave stabilized 

after traversing the first 7 mm of the computational domain, so the calculation was terminated at that 

point. A time step of approximately 1.5×10–10 s was used.

These results were then compared to a calculation obtained from the 1-dimensional ZND code 

developed by Shepherd (1986). To run the ZND code, the desired kinetic mechanism is first input into the 

CHEMKIN code, described by Kee et al. (1993), that provides the kinetic and thermodynamic properties 

required. Next, the Gordon and McBride (1994) CEA code, or any other equilibrium detonation code, is 

used to determine the detonation (Chapman-Jouguet) velocity. With these inputs, ZND calculates the 1–D 

detonation structure by first assuming a normal shock based on the detonation velocity and quiescent 

reactant conditions and then reacts the mixture via the kinetic mechanism and gas properties from 

CHEMKIN using the post-shock conditions as the reaction initial conditions. ZND performs these 
NASA/TM—2010-216765 24

calculations in the reference frame fixed to the detonation wave, shown in Figure 6, so the velocity must 

be transformed back to the stationary reference frame by adding the detonation wave velocity. All of the 

static properties are unaffected by the change in reference frames. 

Figure 7 shows the convergence of the detonation wave speed upon the equilibrium value calculated 

using the CEA code. The wave speed at the termination of the calculation is approximately 1 percent

below the CEA value, which was judged to be sufficiently converged for the purposes of this validation. 

Figure 5.—Detonation case geometry. 

Figure 6.—Detonation calculation reference frames. 

Area = 0.025 m2 

Detonation Wave Position at end of Calculation

L = 10 cm

Initiation Region (0.2 mm)

Reactant Mixture

7.0 mm

Stationary Reference Frame 

(Q1D)

Reference Frame Fixed to

Detonation Wave

(ZND)

Detonation 

Front 

Gas

Stationary

Detonation 

Front 

Unburned Gas

V1=Vcj

Unburned Gas

 V1=0 2

Burned Gas

 V2

Vcj

Burned Gas

1 1 2 V2 
NASA/TM—2010-216765 25

Figure 7.—Detonation wave speed convergence. 

Figures 8 to 15 show the values of pressure, temperature, velocity, and H2O, OH, CO2, NO, and NO2

mass fractions calculated by Q1D at the end of the computation against those calculated by ZND using 

the same kinetic mechanism given in Appendix B. Since the results from ZND begin after the leading 

shock wave of the detonation front and therefore do not include part of the induction zone, leading to a 

difference in total induction time between the two codes, the two sets of results were aligned by matching 

the midpoints of the initial temperature rise region of the temperature profiles shown in Figure 9. Once 

the results were thus synchronized, the ZND and Q1D results generally agree quite well. The largest 

differences are observed in the region of the shock, particularly noticeable in the pressure and velocity 

profiles. This is caused by the chemical reactions and consequent heat release initiating within the shock 

wave in the Q1D model, whereas in the ZND model, it is assumed that the chemical induction time starts 

after the passing of the shock wave, as mentioned above. This difference in the time at which reactions 

start is evident in the H2O, OH, and CO2 profiles. The nitrogen reactions are somewhat slower, and 

therefore do not reflect this difference as significantly. The quicker rise in NO concentration shown in 

Figure 14 for the Q1D code is most likely due to differences in the thermodynamic and transport property 

data sets used by the two codes. Overall, the good agreement between these two results for the detonation 

case indicates that the gas phase fluid dynamic and kinetic models and solvers within the Q1D code are 

working properly. 
NASA/TM—2010-216765 26

0

5

10

15

20

25

30

0 0.02 0.04 0.06 0.08 0.1

Q1D

ZND

Pressure (atm)

Distance (cm)

Figure 8.—Comparison of pressure profiles between Q1D and ZND codes. 

0

500

1000

1500

2000

2500

3000

3500

0 0.02 0.04 0.06 0.08 0.1

Q1D

ZND

Temperature (k)

Distance (cm)

Figure 9.—Comparison of temperature profiles between Q1D and ZND codes. 
NASA/TM—2010-216765 27

0

500

1000

1500

2000

0 0.02 0.04 0.06 0.08 0.1

Q1D

ZND

Velocity (m/s)

Distance (cm)

Figure 10.—Comparison of velocity profiles between Q1D and ZND codes. 

0

0.05

0.1

0.15

0.2

0.25

0.3

0 0.02 0.04 0.06 0.08 0.1

Q1D

ZND

H O Mole Fraction

2

Distance (cm)

Figure 11.—Comparison of H2O mole fraction profiles between Q1D and ZND codes. 
NASA/TM—2010-216765 28

0

0.01

0.02

0.03

0.04

0.05

0.06

0 0.02 0.04 0.06 0.08 0.1

Q1D

ZND

OH Mole Fraction

Distance (cm)

Figure 12.—Comparison of detonation OH mole fraction profiles between Q1D and 

ZND codes. 

0

0.01

0.02

0.03

0.04

0.05

0 0.02 0.04 0.06 0.08 0.1

Q1D

ZND

CO2 Mole Fraction

Distance (cm)

Figure 13.—Comparison of detonation CO2 mole fraction profiles between Q1D

and ZND codes. 
NASA/TM—2010-216765 29

0

0.0001

0.0002

0.0003

0.0004

0.0005

0 0.02 0.04 0.06 0.08 0.1

Q1D

ZND

NO Mole Fraction

Distance (cm)

Figure 14.—Comparison of detonation NO mole fraction profiles between Q1D and 

ZND codes. 

0

5 10-8

1 10-7

1.5 10-7

2 10-7

2.5 10-7

3 10-7

3.5 10-7

4 10-7

0 0.02 0.04 0.06 0.08 0.1

Q1D

ZND

NO2 Mole Fraction

Distance (cm)

Figure 15.—Comparison of detonation NO2 mole fraction profiles between Q1D 

and ZND codes. 
NASA/TM—2010-216765 30

3.2 Particle Drag and Heat Transfer Models and Gas-Solid Momentum and Heat 

Exchange Models Validation 

The next step in validating the Q1D code is to verify that the particle momentum exchange and 

convective heat transfer models are working properly. This can be accomplished by studying a gas flow 

interaction with chemically inert particles. Maxwell and Seasholtz (1974) examined the behavior of 

micron-size inert particles passing through a Mach 1.6 normal shock wave for the purposes of 

determining the tracking behavior of these particles used in Laser Doppler Velocimetry (LDV). The 

physical situation modeled in the study was that of a constant area duct with a supersonic, 289 K static 

temperature, 1 atmosphere static pressure, particle laden inlet air flow and a standing normal shock wave 

interior to the duct. The particles were assumed to be at the gas static temperature and velocity up to the 

shock wave, and to occupy negligible volume. For the analysis, the post-shock conditions, as calculated 

using the standard Hugoniot relations, were used as the initial conditions for the gas, while the pre-shock 

temperature and velocity were used as the initial conditions for the particles. Since only small particles, 

less than 5 μm diameter, were considered for this study, uniform particle temperature was assumed 

throughout. Variations in particle diameter, particle density, gas and particle inflow velocity, particle 

specific heat and particle/gas mass loading ratio were examined. Constant thermodynamic and transport 

properties were assumed for the air and the particles, but the exact values were not given. Because of this, 

the calculations had to be replicated using the following equations adapted from Maxwell and Seasholtz 

(1974), as the values assigned for the thermodynamic and transport properties can significantly affect the 

results. 

A force balance on an individual particle yields the following relation, 

( ) pgpg

p

g

pp

p D uuuu

ur

C

x

u −− 















ρ

ρ















 = ∂

∂

8

3 (86) 

where the absolute value of the gas-particle velocity difference is used as shown to maintain the correct 

sign for the particle velocity spatial derivative. Likewise, an energy balance on an individual particle 

yields the following equation,

( ) pg

ppsp

p g TT

rCu

k

x

T − ρ = ∂∂

2

Nu

2

3 (87) 

Equating the momentum and energy transfers between the two phases, with α defined as the ratio of the 

particle to gas mass flows, yields the following two conservation equations. 

[ ( ) ] ( )

( ) g

g

gp

p s p

g p g

u

u

Tc

x

T

C

x

u

u u

x

u

−−γ

∂

∂

−γα− ∂

∂

−γα−αγ

= ∂

∂

1

1 1

(88) 

x

u

R

u

x

u

R

u

u

T

x

T pggg

g

gg

∂

∂α − ∂∂















 −= ∂

∂

(89) 
NASA/TM—2010-216765 31

In order to solve this system of differential equations, it is also necessary to use the following relation for 

the rate of change of pressure, generated by taking the spatial derivative of the ideal gas law and 

transforming the density term into a velocity term through the expression for continuity in a constant area 

duct. 

















∂

∂ − ∂∂ = ∂

∂

x

V

V

T

x

T

T

P

x

P g

g

gg

g

 (90) 

The particle drag coefficient and Nusselt number correlations used in the above equations are taken from 

Maxwell and Seasholtz, and differ only slightly in value from the expressions used in the Q1D code in 

this flow regime. These equations were then solved using a simple explicit space-marching routine, with a 

step size of 1×10–6 m. 

In order to fully validate the two-phase capability of the Q1D code, a particle/gas mass ratio, α, of 0.2 

was selected for comparison with a solution using the analysis method of Maxwell and Seasholtz. At this 

mass ratio, both the particles and the gas respond significantly to the presence of the other phase. A 

particle diameter of 2 μm and a particle/gas specific heat ratio of 1.125 was used in both analyses. The 

Q1D simulation was run by first establishing a steady, subsonic gas flow in the constant area duct that 

approximated the post-shock conditions, as determined from the standard Hugoniot relations. A total 

pressure, total temperature subsonic characteristic inflow boundary condition and a fixed static pressure 

subsonic characteristic outflow boundary condition, both based on the methodology of Poinsot and Lele 

(1992), were used to establish this initial solution. The particles were then introduced into the flow at the 

desired mass fraction from a fixed position in the computational domain as the solution continued to run. 

Each set of particles was inserted at the pre-shock static temperature and flow velocity. The solution was 

then allowed to run out to a new steady state solution. As the gas-solid interactions caused the gas mass 

flow to shift from the initial result, the particle mass flow rate was continuously adjusted to maintain a 

particle/gas mass ratio of 0.2. The gas phase inlet velocity, static temperature, and static pressure also 

shifted from the initial gas-phase solution due to the interaction with the particles. These new static/total 

ratios were then characteristic of a slightly different shock Mach number. Therefore, the analytical 

solution was run at the gas flow inlet conditions of the converged numerical solution for comparison. 

Figures 16 and 17 show excellent agreement between the Q1D simulation and the Maxwell/Seasholtz

analysis gas temperature and gas velocity from the point of particle introduction to the end of the 

computational domain. Figures 18 and 19 show the corresponding results for the particle “stream”, also 

with excellent agreement between the two methods. 

These results are consistent with the expected behavior of the two phases as given by Maxwell and 

Seasholtz. The slight differences between the results of the two methods observed in Figures 16 to 19 are 

due to a combination of the different particle drag coefficient and Nusselt number correlations used in the 

different calculations, and the use of varying thermodynamic and transport properties in the Q1D code 

versus the fixed values used in the Maxwell/Seasholtz analysis.
NASA/TM—2010-216765 32

380

390

400

410

420

0 0.005 0.01 0.015 0.02

Q1D

Maxwell

Temperature (K)

Distance (m)

Figure 16.—Gas temperature profile from Q1D calculation from point of particle 

insertion compared to gas temperature profile from Maxwell/Seasholtz

analysis. 

200

220

240

260

280

300

0 0.005 0.01 0.015 0.02

Q1D

Maxwell

Velocity (m/s)

Distance (m)

Figure 17.—Gas velocity profile from Q1D calculation from point of particle 

insertion compared to gas velocity profile from Maxwell/Seasholtz

analysis. 
NASA/TM—2010-216765 33

250

300

350

400

450

0 0.005 0.01 0.015 0.02

Q1D

Maxwell

Temperature (K)

Distance (m)

Figure 18.—Particle temperature profile from Q1D calculation from point of 

particle insertion compared to particle temperature profile from 

Maxwell/Seasholtz analysis. 

200

300

400

500

600

0 0.005 0.01 0.015 0.02

Q1D

Maxwell

Velocity (m/s)

Distance (m)

Figure 19.—Particle velocity profile from Q1D calculation from point of 

particle insertion compared to particle velocity profile from 

Maxwell/Seasholtz analysis. 
NASA/TM—2010-216765 34

3.3 Particle Internal Heat Transfer, Heat Release and Supersonic Drag/Convection

Validation 

Next, the results of Baek (1985), summarized also in Sichel, et al. (1985), were used to validate the 

particle internal heat transfer/heat release modeling in the Q1D code, as well as the supersonic relative Mach 

number drag and convective heat transfer correlations. Baek performed a series of experiments examining 

the shock ignition of various types and sizes of particles, which he then also examined using a 

computational model and asymptotic analysis. The most comprehensive set of results reported by Baek were 

obtained using 53 μm coal particles subjected to a Mach 4.8 traveling shock wave. The particles were 

injected into a shock tube containing pure oxygen at 1/3 of an atmosphere pressure and 295 K in the driven 

section. The mass of the particles was very small compared to the mass of gas present, so the particles were 

assumed to have no effect on the gas flow. The particles were also assumed to be of uniform size and 

spherical. The coal particles were characterized as having the following physical properties.

Density: 1200 Kg/m3

Thermal Conductivity: 0.887 J/m-s-K 

Heat Capacity: 987.4 J/Kg-K 

Heat of Combustion: 35812 KJ/Kg

Internal Surface Area: 426000 m2

/Kg

Under the conditions tested, the particle consumption rate, in g/cm2

-s (based on the particle external 

surface area), was assumed to follow a simple Arrhenius expression of the same form as shown in 

Equation (91), with an activation energy of 35,700 cal/g-mole and a pre-exponential factor of 

87,100 Kg/m2

-s-atm. Baek’s numerical single-particle combustion model was found to match the 

experimental data well over a range of shock Mach numbers and particle sizes. 

In the Q1D code, this scenario was modeled by looking at a single 53 μm coal particle, with the 

physical properties listed above, injected at zero velocity into a steady flow at the post-shock conditions 

as calculated from a standard ideal shock tube calculation, such as that found in Saad (1985). A particle 

porosity of 0.32 with 2 nm diameter pores was found to give the required internal surface area. The two 

Arrhenius rate expressions required for the Hong particle combustion kinetic rate, Equation (59) (restated 

below), were defined as shown in Equations (91) and (92) to match the behavior of the simple kinetic 

expression of Baek. The Q1D code was then run until ignition occurred. 

Sip

cSip

tot

ext

avgik PK

MPk

S

S

S

S

r

,O,

O,1 , int

int

,

2

2

 1+











+η= (59) 

TR ipu k p e ,

35700

9 1 1009.3

−

×= (91) 

TR ipu p eK ,

30300

4.13

−

= (92) 

The results from the Q1D code are compared to the results from Baek (1985) in Figures 20 to 22. The 

Baek results were recalculated using Baek’s original code, a listing of which was included in the report.

Figure 20 shows the gas-particle relative Mach number from just after particle insertion up until the time 

at which the particle ignited in the Baek calculation. The small difference in relative Mach number 

reflects differences in the thermodynamic properties between the two codes, resulting in a slightly 

different sound speed in the gas. Other than this slight offset, the Mach number decay is nearly identical, 

which is to be expected, as Baek also used the Henderson (1976) particle drag correlations, given in 

Equations (31) to (34).
NASA/TM—2010-216765 35

0

0.5

1

1.5

2

0 5 10 15

Q1D

Baek

Relative Mach Number

Time (microseconds)

Figure 20.—A comparison of the post-shock gas-particle relative 

Mach number profiles between the Q1D and Baek codes. 

200

400

600

800

1000

1200

0 0.05 0.1 0.15 0.2 0.25 0.3

Temperature (K)

Fraction of Particle Radius from Particle Surface

Figure 21.—Internal temperature profiles at 2 μsec intervals of particle 

heating up to ignition from Baek code. 
NASA/TM—2010-216765 36

Figure 21 shows a time series of particle internal temperature profiles as calculated using Baek’s 

code. The particle surface temperature is at the left hand side of the plot, with the temperature falling off 

with distance into the particle. The lowest temperature profile corresponds to 2 μsec after the passage of 

the shock wave, and each succeeding temperature profile represents the passage of another 2 μsec of time.

Ignition occurs just prior to the 14 μsec (uppermost) profile. 

Figure 22 shows a similar plot of results from the Q1D code. The overall ignition process took 

approximately 4 μsec longer in this simulation, due to a combination of small differences in the 

thermodynamic properties, transport properties, and heat transfer relations used in the different codes.

While this approximately 25 percent increase in ignition time seems significant, it actually represents a 

very small difference in particle heating rate. The continual decrease in the convective heating rate of the 

particle due to the decrease in the relative velocity between the gas and the particle and the decreasing 

temperature difference leads to the particle surface approaching the ignition temperature asymptotically.

This behavior translates to a high sensitivity of the ignition time to the heating rate, since it could take a 

significant portion of the ignition time to heat the last few percent of the required temperature rise.

Even with this longer ignition time, the Q1D result is within the scatter of the experimental data for 

ignition delay given by Baek for this condition. Taking this into consideration, the Q1D code results 

compare quite favorably with the analysis and experimental results of Baek.

200

400

600

800

1000

1200

0 0.05 0.1 0.15 0.2 0.25 0.3

Temperature (K)

Fraction of Particle Radius from Particle Surface

Figure 22.—Internal temperature profiles at 2 μsec intervals of particle heating up to 

ignition from Q1D code. 
NASA/TM—2010-216765 37

4.0 Generation of Graphite Kinetic Constants

The graphite combustion rate expression, Equation (59), repeated below, requires that the constants in 

the kinetic rate expressions, Equations (55) and (56), also repeated below, be derived from experimental 

data for the given material and particle shape. 

Sip

cSip

tot

ext

avgik PK

MPk

S

S

S

S

r

,O,

int ,O,1

int

,

2

2

 1+











+η= (59) 

TR ipu k p e ,

51000

8 1 1029.6

−

×= (55) 

TR ipu p eK ,

10100

4.13

−

= (56) 

The above expressions contain the two sets of kinetic rate constants given by Hong (2000) for nonporous graphite flakes, based on the experimental data of Ranish and Walker (1993). While the data used 

to generate this data covers a broad pressure range (from 0.1 MPa up to 6.4 MPa) over a solely kinetically 

controlled temperature range (733 to 842 K), the graphite flakes used were significantly different in 

composition, structure and shape than the graphite spheres of interest in this study. The Ranish and 

Walker graphite flakes averaged 30 μm in diameter and 0.5 μm thickness. They were also a high purity, 

non-porous natural crystalline graphite. In the HTF facility described earlier, the graphite used in the 

facility heater and the piping liners is Union Carbide PGX graphite (GrafTech, 2005). PGX is an artificial 

molded graphite material formed using lampblack particles as large as 0.8 mm held together with a coaltar pitch based binder. These materials are mixed together, poured into a mold and graphitized by baking. 

The resultant material has a density of 1730 Kg/m3

 (Graftech, 2005), which translates to a porosity of 

0.2345 (based on a non-porous graphite density of 2260 Kg/m3

) with an average pore diameter of 6 μm 

(Pierson, 1993). Observations made during test facility operations indicate that roughly spherical particles 

are to be expected. Because of these differences, an alternate data set from which to generate the kinetic 

constants was needed. However, no available data set was found in the literature for an appropriate 

material that covered the pressure and temperature range of interest. In particular, no data sets including 

data for pressures above 1 atmosphere were found. Because of this, the Hong kinetic constants were used 

as a starting point, since they included the effect of pressure, and then the constants in the k1p expression 

were modified to fit the available ambient pressure data. 

The frequently referenced graphite sphere combustion rate data of Tu, et al. (1934) was one possible 

data set, but upon examination, it was found that the graphite used was not sufficiently characterized so as 

to provide all the inputs needed for the combustion model of Hong. A number of other graphite oxidation 

data sets, such as that of Okada and Ikegawa (1953), were obtained for cylindrical graphite electrodes, but 

not for spheres. The only applicable data was that taken by Froberg (1967) in a set of experiments 

modeled after the work of Tu, et al. (1934), but for which significantly more graphite characterization 

information was available. The graphite used by Froberg had a density of 1650 Kg/m3 and an average 

pore diameter of 5 μm, both similar to the PGX graphite used at HTF. In the set of Froberg’s experiments 

that were used to generate the graphite rate constants for this study, large graphite spheres (~12.6 mm 

diameter) were exposed to heated air streams with air temperatures ranging from 813 to 1196 K. The air 

heater inflow velocity was held constant at 1.0 cm/s (at standard temperature and pressure), which 

resulted in an increasing air velocity as the air temperature was raised. The graphite sphere was allowed to 

come to equilibrium before combustion rate data was taken, which resulted in graphite surface 

temperatures as high as 1311 K for the highest temperature air flow. The combustion rate was directly 

measured via sphere weight loss.
NASA/TM—2010-216765 38

For the Q1D simulation of this configuration, thermal conductivity and heat capacity as a function of 

temperature were also needed for the graphite. As this data was not available from Froberg for the 

specific graphite studied, estimates were made based on literature information. The thermal conductivity 

of PGX is 120 W/m-K at room temperature (Graftech, 2005) and decays exponentially as temperature 

increases. It was assumed that the Froberg graphite had the same room temperature thermal conductivity 

as the PGX. The thermal conductivity decay rate, as well as graphite specific heat as a function of 

temperature, are given graphically by Pierson (1993), and were used to generate 3rd order polynomial 

curve fits of these properties for use in the Q1D code, as shown in Figures 23 and 24. For those cases 

where an internal particle temperature profile was used, mass-averaged values for heat capacity and 

thermal conductivity were used. 

Figure 23.—PGX graphite thermal conductivity. 

Figure 24.—PGX graphite heat capacity. 
NASA/TM—2010-216765 39

Figure 25.—Q1D kinetic constant fit to Froberg data. 

With this information, Q1D simulations of a number of data points from Froberg were made, and the 

kinetic constants adjusted until the simulations and the experimental results were in reasonable 

agreement. The resultant expression for k1p is given in Equation (93), which represents a factor of 11.3 

increase in the pre-exponential factor and a 24 percent increase in the activation energy over the 

expression from Hong (2000) for non-porous graphite flakes given in Equation (55). 

TR ipu k p e ,

63100

9 1 1011.7

−

×= (93) 

A graphical comparison between the Froberg data and the adjusted kinetic model is shown in Figure 25, 

based on the sphere surface temperature. Using this adjusted expression for k1p, PGX graphite can be 

simulated by adjusting the particle density, pore diameter, and particle size to the appropriate values.

5.0 Model Facility and Scramjet Simulation Geometry and Flow Conditions

5.1 Test Facility Components

In the introduction, two different hypersonic propulsion test facility storage-heater types were 

described. The first utilized ceramic pebbles heated by flowing combustion exhaust gases through the 

bed. The second consisted of inductively heated drilled graphite blocks. The first type is found in the 

GASL Leg IV facility (Roffe, et al., 1997), while the second type is used at the NASA Plumbrook HTF 

facility (Perkins, et al., 1998). Schematics of these facilities were presented in Figure 1 and Figure 2.

Rather than attempting to simulate both of these existing facilities for this code demonstration, a single 

facility geometry will be used to allow for a direct comparison of results for non-reacting and reacting 

particles. As the goal of this study is not to determine the state of particulates entering the combustor of a 

particular facility at a specific test condition, but instead to parametrically determine the range of 

particulate states that could be encountered, and the dominant physical processes that determine these 

states, the use of a generic facility geometry does not detract from the objectives of the study. 
NASA/TM—2010-216765 40

Figure 26.—Scaled axisymmetric geometry for facility settling chamber/mixer and nozzle. 

The selected generic facility geometry, shown in Figure 26, consists of a 3 m long low speed 

(~ Mach 0.08 gas flow velocity) mixing section followed by a 3 m long Mach 5 exit freejet nozzle. The 

scale of the mixer and nozzle sections is slightly less than that of the NASA HTF facility, and 

considerably larger than that of GASL Leg IV facility. The nozzle contour was generated using a spline 

fit to the inlet, throat, and exit axial locations and diameters, estimated from the facility descriptions given 

in the above cited facility references. The same nozzle geometry is used for testing at simulated flight 

Mach numbers of 5, 6, and 7 by adjusting the mixer section inflow total pressure and total temperature, 

and by adjusting the flow deflection angle and length of a pre-compression plate that will be described in 

detail in the next section. The inflow conditions for the Mach 5, 6 and 7 simulated flight conditions are 

given in Table I, along with the representative flight dynamic pressures used to calculate the facility 

conditions, and include corrections for real gas effects taken from NASA compressible flow tables (Ames 

Research Staff, 1953). 

TABLE 1.—MIXER SECTION INFLOW TOTAL PRESSURE, TOTAL TEMPERATURE AND 

DYNAMIC PRESSURE FOR MACH 5, 6 AND 7 SIMULATED FLIGHT CONDITIONS

Ptotal

(KPa)

Ttotal

(K) 

Dynamic pressure

(KPa)

Mach 5 3082 1219 89

Mach 6 4498 1652 54

Mach 7 8247 2174 47

5.2 Scramjet Inlet and Isolator Sections

After the facility components are defined, it is necessary to define the geometry of the scramjet model 

to be simulated. In practice, the entire scramjet inlet is not used in ground testing due to its length. The 

multiple external compression ramps are instead replaced by a single pre-compression plate that provides 

the same inflow condition at the engine inflow plane. However, it is still necessary to define the full inlet 

geometry in order to determine the appropriate pre-compression plate geometry. Figure 27 shows a 

generic, unscaled scramjet inlet and isolator, along with a short section of the combustor. Table 2 defines 

the geometric parameters shown in the figure for a Mach 7 design point inlet/isolator system. The values 

upstream of the engine throat were determined by first selecting the throat height, H1, to set an overall 

engine scale consistent with the generic facility size, and then iterating on the compression ramp angles 

and lengths such that all of the oblique shocks converge on the cowl lip, and the cowl lip oblique shock in 

turn impinges on the shoulder that defines the inlet throat, as shown in Figure 27. This design was further 

constrained by the typical design practices of keeping the compression ramp turning angles approximately 

equal to minimize overall total pressure loss, and of targeting a throat Mach number of approximately 

one-half of the free stream Mach number. The isolator section was sized using the methodology described 

by Heiser and Pratt (1994) for constant area isolators, with an isolator pressure ratio of 2. The slight 

18.517 cm

36.1 cm

49.4 cm

300 cm 300 cm

Mixer 99.7 cm

Nozzle
NASA/TM—2010-216765 41

divergence angle within the isolator section accounts for boundary layer growth that would otherwise 

restrict the flow in the duct. A series of linked spreadsheets were used to automate this design process.

Figure 27.—Scramjet geometric design parameters. 

Figure 28.—Scaled inlet/isolator at Mach 7 design point. 

Figure 28 shows a scaled representation of the inlet/isolator and the shock system at the Mach 7 

design point based on the geometric parameters given in Table 2. Figures 29 and 30 show the shock 

structure moving out from the cowl lip as the flight Mach number is decreased to Mach 6 and Mach 5, 

respectively. Table 3 then provides the static pressure, static temperature, and Mach number at key points 

through the geometry for the selected Mach 5, 6 and 7 flight conditions. It should be noted that the 

conditions at the isolator entrance are a result of the combination of the cowl lip oblique shock with an 

expansion fan emanating from the inflection point at the inlet throat , which is not shown in Figures 27 

to 30. 

TABLE 2.—SUMMARY OF MODEL SCRAMJET INLET/ISOLATOR 

GEOMETRY AT MACH 7 DESIGN POINT

L1 = 353.2 cm H1 = 6.0 cm β1 = 6.4° 

L2 = 120.2 cm α1 = 8° β2 = 8.2°

L3 = 99.6 cm α2 = 17° β3 = 10.4°

L4 = 33.2 cm α3 = 26° β4 = 13.2°

L5 = 547.6 cm α4 = 2° ------------

Cowl

Shock Wave

Isolator

Combustor
NASA/TM—2010-216765 42

Figure 29.—Scaled inlet/isolator with shock locations corresponding to Mach 6 operation. 

Figure 30.—Scaled inlet/isolator with shock locations corresponding to Mach 5 operation. 

TABLE 3.—SUMMARY OF FLOW FIELD CONDITIONS FOR MODEL SCRAMJET 

Mach 5 Mach 6 Mach 7

Pstatic

(Pa)

Tstatic

(K)

Mach 

no.

Pstatic

(Pa)

Tstatic

(K)

Mach 

no.

Pstatic

(Pa)

Tstatic

(K)

Mach 

no.

Free stream 5024 215 5 2110 227 6 1370 230 7

Along 1st ramp 12600 285 4.2 6200 319 5.0 4670 342 5.6

Along 2nd ramp 30100 372 3.6 16900 434 4.1 14200 483 4.6

Along 3rd ramp 63500 463 3.0 39300 555 3.5 35800 634 3.8

Isolator inlet 60400 512 2.8 37400 636 3.1 34300 751 3.4

Isolator exit 120800 798 1.8 74800 1000 2.0 68600 1237 2.2

Having established the scramjet inlet/isolator geometry and the flow conditions at each step through 

those components, the test configuration to be used can be defined. To shorten the scramjet geometry, the 

first two external compression ramps from Figure 27 are replaced by a single pre-compression plate, 

leaving the last external compression ramp, the cowl, and the isolator unchanged. In order to create the 

same conditions along the last external compression ramp before the throat as those achieved in the 

original configuration from Figure 27, the geometry has to be rotated slightly counter-clockwise, as 

shown in Figure 31(a). This rotation compensates for the change in the flow direction coming at the last 
NASA/TM—2010-216765 43

external compression ramp in the shortened configuration relative to that in the original configuration. 

Each of the geometric parameters shown in Figure 31(a) needs to be changed for each test Mach number 

to be able to match the flow conditions approaching the last external compression ramp calculated for the 

original configuration. The required values as a function of flight simulation Mach number are given in

Table 4. The flow conditions remain the same along the last external compression ramp and within the 

isolator, but a new set of flow conditions are created along the precompression plate, which are given in 

Table 5. Figure 31 (b) shows the test article from Figure 31(a) schematically “installed” in the model test 

facility, indicating the full computational domain.

(a) 

(b) 

Figure 31.—(a) Truncated configuration for wind tunnel testing with forebody simulation (precompression) 

plate and (b) truncated configuration with forebody simulation plate shown schematically with facility mixer 

and nozzle indicating computational domain. 

Mixing Length

N2

or 

Air

Refractory Material

(Pebbles or 

Drilled Blocks)

O2 & N2 Injection 

(As Needed)

Heater 

Shutoff 

Valve

Computational Domain for Code Demonstration

Facility Nozzle

Test Model

α5

α6

L6

β5

Precompression 

Plate 

Isolator

Final External

Compression Ramp

Cowl
NASA/TM—2010-216765 44

TABLE 4.—TEST ARTICLE CONFIGURATION FOR MACH 5, 6 AND 7 TEST CONDITIONS 

Mach 5 Mach 6 Mach 7

β5 9.3° 9.3° 10.1°

α5 15.0° 9.5° 4.5°

α6 2.0° 7.5° 12.5°

L6 100.6 cm 99.1 cm 84.0 cm

TABLE 5.—FLOW CONDITIONS ALONG THE PRE-COMPRESSION PLATE

FOR MACH 5, 6, AND 7 TEST CONDITIONS

Mach 5 Mach 6 Mach 7

Pstatic (Pa) 23700 21900 25200

Tstatic (K) 370 430 510

Mach no. 3.5 4.0 4.4

5.3 Modeling Approach

Simplifying assumptions must be made in order to model a gas-particle flow through the facility 

mixer, the facility nozzle, along two external compression ramps, through the inlet throat area and the 

isolator section using a quasi-1D code. Since the facility mixer and nozzle are assumed to be 

axisymmetric, it is straightforward to cast their geometry into a quasi-1D format. However, the 

precompression plate, final external inlet ramp and isolator geometries are inherently two-dimensional, as 

are the oblique shocks preceding each of those components and the expansion fan leading into the 

isolator. This difficulty was addressed by modeling each of these last three components separately from 

the mixer/nozzle section. The jump conditions across each of the oblique shocks were calculated from the 

exit conditions of the upstream component and the specified turning angle using classical methods, such 

as those found in Saad (1985). The calculated post-shock conditions were then used as the inflow 

boundary conditions for the next component. For the isolator, the post-shock conditions were used as the 

input to a classical expansion fan calculation to determine the inflow boundary conditions of that 

component. Each of these last components was then treated as a constant area section having a length 

equal to the corresponding body-side bounding surface, as shown in Figure 31. The isolator was 

considered as a constant area section, even though a slight area increase with length is shown in 

Figure 26. However, that area relief is designed to accommodate boundary layer growth so that the 

isolator can behave as a truly constant area section and not prematurely choke near the end of the duct. 

Since the 1-D analysis does not include boundary layer, a constant area duct was modeled. The crosssectional area for these constant area sections can be any arbitrary value, though a value of 0.1 m2 was 

generally used throughout. It should be noted that the length of each component is representative of the 

gas-particle flow near the bounding wall. Gas-particle flows further out from the wall will have a shorter 

residence time between shocks. 

One additional modeling step was required for the isolator. Isolator flow is generally characterized by 

a complex, three-dimensional series of reflected oblique shocks, possibly followed by a weakened normal 

shock in the ramjet operating mode. It is not realistic to try and model this flowfield using a 1-D code. 

Instead, a uniform negative momentum source term was added to the momentum Equation (2), at each 

grid point along the isolator length to reflect the total pressure loss generated by the internal shock 

structure. This approach yielded nearly linear changes in the primary flow variables between the isolator 

inlet and exit, ending at the desired isolator exit values given in Table 3. This behavior for the various 

flow parameters can be observed in the various plots in the later Results section. 

With these modeling assumptions, a 1 percent reacting particle mass fraction case had a typical total 

run time of approximately 80 hr on a dual-processor (2 GHz) PC workstation. This relatively long run 

time was driven primarily by the low convective velocity in the low-speed mixer section. A similar case 

with just a single reacting particle typically took approximately 40 CPU hours on the same machine. In 

both cases, a converged no-particle solution was used as the initial condition.
NASA/TM—2010-216765 45

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0 0.5 1 1.5 2 2.5 3

Area (m2

)

X (m)

Figure 32.—Nozzle area distribution where X = 0 indicates the beginning of 

the convergent section of the nozzle.

5.4 Grid Descriptions

Simple uniform 1-D grids were used for each component. For the mixer/nozzle simulation, an overall 

grid totaling 4000 points was split two sections. An 800 point uniform grid was used for the 300 cm 

constant-area mixer section and a 3200 point uniform grid was used for the 300 cm nozzle section. The 

exact area distribution used for all of the simulations for the 3 m long nozzle, as shown schematically in 

Figure 26, is given in Figure 32. Each of the remaining three components downstream of the nozzle was

modeled using uniform 300 point grids. 

6.0 Results

6.1 Inert Particle Calculations

A series of Q1D calculations were performed with inert alumina particles utilizing the test facility and 

scramjet geometries described in the previous section at the Mach 5 and 6 total temperature flow 

conditions defined therein. The Mach 7 total temperature condition was not run because this condition is 

near the melting point of alumina (2323 K). As mentioned in the introduction, particle diameters of 10 μm 

to 1 mm (1000 μm) were considered for this study. In addition to examining the two extreme particle 

diameters, an intermediate diameter of 100 μm was also included in the study to help define the effect of 

particle diameter on particle and gas flow conditions entering the model scramjet combustor. The 

majority of results shown below are for single particles traversing the facility mixer/nozzle, 

precompression plate, final inlet compression ramp and isolator components. Additionally, a single case 

with a 1 percent particle mass fraction was run to determine if the effect of increased particle mass 

fraction on the scramjet combustor inflow conditions warranted additional calculations at elevated inert 

particle mass fractions. 

The particles used in this portion of the study were assumed to be low porosity, high purity tabular 

alumina. All the required material properties, including temperature dependent properties, were available 
NASA/TM—2010-216765 46

from Munro (1997), with the exception of particle emissivity. However, Pluchino and Masturzo (1981) 

determined that the emissivity of small alumina particles approaches zero as the material purity 

approaches 100 percent. As the tabular alumina is a very high purity form of alumina, it was therefore 

assumed that the particle emissivity was in fact zero for these calculations, eliminating all radiative heat 

transfer from the particles to the water-cooled facility walls (500 K wall temperature). The published 

density for tabular alumina is 3984 Kg/m3

. Temperature dependant values of heat capacity and thermal 

conductivity are shown in Figures 33 and 34, respectively. 

Figure 33.—Tabular alumina heat capacity. 

Figure 34.—Tabular alumina thermal conductivity. 
NASA/TM—2010-216765 47

It should be noted that these property values had to be extrapolated above 1800 K. However, as each 

of these properties approach an asymptote by that temperature, no significant error is incurred. For 

particles with an internal temperature gradient, mass-averaged values of these temperature dependent 

properties were calculated and applied at each time step. 

Prior to running these inert particle cases, it was desirable to determine whether or not it would be 

necessary to include gas-phase chemistry in the calculations. Gas-only calculations with and without gasphase chemistry were run at the Mach 5 and 6 total temperature conditions. The Mach 5 and 6 

calculations showed no significant difference in the bulk fluid properties between the results with and 

without gas-phase chemistry. All minor species in the finite-rate chemistry calculations remained below 1 

part-per-million throughout the computational domain for both conditions. Based on these results, it was 

determined that the inert particle cases could reasonably be run without using finite rate chemistry. By 

running all of the inert particle cases in this manner, significant computational run-time improvements 

were achieved. 

A total of six cases were run using single 10-, 100- and 1000-μm alumina particles at the Mach 5 and 

6 total temperature flow conditions. Initial test runs showed that the Biot number for the 10 μm particles 

never exceeded 0.1 during the simulation, indicating that those particles should be modeled as isothermal.

The larger particles were modeled with internal grids of 51 grid points for the 100 μm particles and 101 

grid points for the 1000 μm particles. The number of internal grid points for each particle size was 

selected by running a series of test cases with increasingly refined grids until the solution became grid 

independent. Figures 35 to 38 show the gas and particle surface temperature and velocity results for these 

six cases, each plot at a single facility operating point total temperature. Only a single gas-phase solution 

is given in each of these figures as the gas-phase solutions at any given Mach number are not effected by 

the insertion of a single particle. The discontinuities in the gas flow results indicate the positions of the 

oblique shock waves indicated in Figure 31 at each flow turning point. The particle initial temperature 

was set 200 K higher than the initial gas temperature in the mixer section to simulate the overheating of 

the test gas typically required in storage heater facility operations, as described earlier. The particle initial 

velocity is set to the incoming gas flow velocity.

Examining these figures, several observations can be made. First, that while the 10 μm particles 

reasonably follow the gas flow temperature and velocity profiles, the 100 and 1000 μm particles 

significantly lag the gas flow due to their higher mass to surface area ratio that decreases the particle drag 

relative to the particle inertia. Second, that while the 10 μm particles respond to the step-change in gas 

flow properties at the oblique shock positions, albeit with significant delay, the 100 and 1000 μm particles 

show almost no response to the shocks for the same reason as the first observation above. And lastly, that 

increasing the facility total temperature from the Mach 5 condition to the Mach 6 condition does little to 

alter the response of the particles, other than to shift the temperature and velocity scales. This shift in 

scales is more evident in Figures 39 to 44, where the data is replotted keeping the particle size constant in 

each figure, while varying the simulated Mach number total temperature. 
NASA/TM—2010-216765 48

Figure 35.—Comparison of gas and single particle surface temperature profiles for 

various particle sizes for the Mach 5 total temperature flow condition. 

Figure 36.—Comparison of gas and single particle velocity profiles for 

various particle sizes for the Mach 5 total temperature flow condition. 
NASA/TM—2010-216765 49

Figure 37.—Comparison of gas and single particle surface temperature profiles 

for various particle sizes for the Mach 6 total temperature flow condition. 

Figure 38.—Comparison of gas and single particle velocity profiles for 

various particle sizes for the Mach 6 total temperature flow condition. 
NASA/TM—2010-216765 50

Figure 39.—Comparison of gas and 10 μm particle surface temperature 

profiles for Mach 5 and 6 total temperature flow conditions. 

Figure 40.—Comparison of gas and 100 μm particle surface temperature 

profiles for Mach 5 and 6 total temperature flow conditions. 
NASA/TM—2010-216765 51

Figure 41.—Comparison of gas and 1000 μm particle surface temperature 

profiles for Mach 5 and 6 total temperature flow conditions. 

Figure 42.—Comparison of gas and 10 μm particle velocity profiles for 

Mach 5 and 6 total temperature flow conditions. 
NASA/TM—2010-216765 52

Figure 43.—Comparison of gas and 100 μm particle velocity profiles for 

Mach 5 and 6 total temperature flow conditions. 

Figure 44.—Comparison of gas and 1000 μm particle velocity profiles for 

Mach 5 and 6 total temperature flow conditions. 
NASA/TM—2010-216765 53

Figures 45 and 46 show the internal particle temperature profiles for the 100 μm particles and 

1000 μm particles, respectively, at the exit plane of the mixer-nozzle section (MN) and the exit plan of 

the isolator (IS). The temperature values are normalized by temperature at which they are initially inserted 

in the flow stream, Tinsert, while the particle radius coordinate, r, is normalized by the particle outer radius, 

rp. It is of interest that the 100 μm particle internal temperature profiles shift somewhat with the change in 

simulated Mach number in Figure 45, indicative of significant particle temperature change, but that the 

particle internal temperature profiles remain almost constant with the change in simulated Mach number 

for the 1000 μm particles shown in Figure 46. It is also evident, as expected, that the larger particles, with 

their larger thermal mass, sustain a larger thermal gradient at the particle surface.

Lastly, the effect of a relatively high particle mass fraction is examined in Figures 47 to 52. In 

Figure 47, the gas and particle temperature profiles for a case using a single 10 μm particle are compared 

to the gas and particle temperature profiles for a case wherein the flow contained a 1 percent mass 

fraction of 10 μm particles inserted into the flow at the same conditions as the single particles. The 

profiles are nearly identical for both the gas and particles. Figures 48 and 49 show a closer view of the 

nozzle exit and the ramp and isolator sections where a small difference in gas and particle temperatures is 

evident, however. After the same manner, the overall gas and particle velocity profiles are shown in 

Figure 50, with more detailed views given in Figures 51 and 52. Looking at these results, one must 

conclude that while some difference does exist in gas and particle properties due to the presence of a 

1 percent mass fraction of particles, the differences are insignificant. Based on this result, no further 

1 percent inert particle mass fraction simulations were performed.

Figure 45.—Radial temperature profiles of 100 μm particles for the Mach 5 

and 6 conditions at the exit planes of the mixer/nozzle (MN) and the 

isolator (IS). 
NASA/TM—2010-216765 54

Figure 46.—Radial temperature profiles of 1000 μm particles for the Mach 5 

and 6 conditions at the exit planes of the mixer/nozzle (MN) and the 

isolator (IS). 

Figure 47.—Comparison of gas and particle temperature profiles for the 

Mach 6 total temperature flow condition for a single 10 μm particle and 

for a 1 percent 10 μm particle mass fraction. 
NASA/TM—2010-216765 55

Figure 48.—Comparison of gas temperature profiles for the Mach 6 total 

temperature flow condition for a single 10 μm particle and for a 1

percent particle mass fraction of 10 μm particles over the last 1.0 m of 

the nozzle, the precompression plate, the final external compression 

ramp, and the isolator. 

Figure 49.—Comparison of particle temperature profiles for the Mach 6 

total temperature flow condition for a single 10 μm particle and for a 

1 percent particle mass fraction of 10 μm particles over the last 1.0 m of 

the nozzle, the precompression plate, the final external compression 

ramp, and the isolator. 
NASA/TM—2010-216765 56

Figure 50.—Comparison of gas and particle velocity profiles for the Mach 

6 condition for a single 10 μm particle and for a 1 percent 10 μm particle 

mass fraction. 

Figure 51.—Comparison of gas velocity profiles for the Mach 6 total 

temperature flow condition for a single 10 μm particle and for a 

1 percent particle mass fraction of 10 μm particles over the last 1.0 m 

of the nozzle, the precompression plate, the final external compression 

ramp, and the isolator. 
NASA/TM—2010-216765 57

Figure 52.—Comparison of particle velocity profiles for the Mach 6 total 

temperature flow condition for a single 10 μm particle and for a 

1 percent particle mass fraction of 10 μm particles over the last 1.0 m 

of the nozzle, the precompression plate, the final external compression 

ramp, and the isolator. 

6.2 Single Graphite Particle Calculations

Using the kinetic and physical property information for graphite previously described, the Mach 5 and 

6 total temperature single inert particle simulations were repeated using single graphite particles of the 

same particle diameters. For graphite the radiative emissivity is non-zero, so radiation heat transfer from 

the hot particles to the facility wall is included in each graphite particle simulation. A value for the 

emissivity of graphite of 0.9 was used (Pierson, 1993). Also, since the allowable temperature for graphite 

is considerably higher than that for alumina, Mach 7 total temperature simulations were performed for all 

three particle sizes. As with the inert particle cases, gas-phase finite rate chemistry was not used for the 

Mach 5 and 6 calculations. For the Mach 7 total temperature simulations, examination of gas-phase only 

calculations with and without finite-rate chemistry showed that significant mass fractions of minor 

species, particularly NO, were generated, as shown in Figure 53. 

The formation of NO in turn had a measurable impact upon the O2 and N2 concentrations, as shown in 

Figure 54, which in turn had an effect on the static temperature, shown in Figure 55. Based on these results, 

the Mach 7 total temperature simulations were run with gas-phase finite-rate chemistry. In order to decrease 

computational time, those reactions given in Appendix B that include hydrogen-containing species were 

eliminated for these calculations, as there is no hydrogen in the system up through the isolator section.

One additional adjustment was made for the graphite particle calculations in comparison to the alumina 

particle calculations. For the alumina particles, an internal grid was used for both the 100 and 1000 μm 

diameter particles. However, the much higher thermal conductivity of graphite compared to alumina 

resulted in a Biot number for the 100 μm graphite particle below 0.1 during the entire particle trajectory. 

During the code validation portion of this study, it was shown that such low Biot numbers, particularly at 

high heat flux conditions, yielded lower accuracy heat transfer results. Therefore, the 100 μm graphite 

particle cases were run without an internal particle grid (i.e. with the whole particle at a single temperature). 
NASA/TM—2010-216765 58

Figure 53.—Mass fractions of minor species from Mach 7 total 

temperature reacting flow simulation

Figure 54.—Comparison of mass fractions of major species from Mach 7 total 

temperature simulations with and without finite-rate chemistry. 
NASA/TM—2010-216765 59

Figure 55.—Comparison of static temperature profiles from Mach 7 total 

temperature simulations with and without finite-rate chemistry. 

Figures 56 and 57 show the particle temperature results for the Mach 5 and Mach 6 total temperature 

conditions, respectively, along with the corresponding inert alumina particle results for comparison. 

Figure 58 shows the temperature results for the Mach 7 condition that was not run for the inert particles.

The 10 μm graphite particle is immediately quenched at the Mach 5 condition by the colder gas around it 

due to its high surface area to mass ratio leading to a high heat loss to the surroundings relative to the 

combustion heat addition at this condition, and stays quenched throughout the particle trajectory. The 

different thermal response compared to the 10 μm alumina particle in the nozzle is due to the differences 

in density, porosity, heat capacity, and thermal conductivity between the two materials. For the Mach 6 

and 7 simulations, the 10 μm particle is consumed within a very short distance as the combustion heat 

addition is now much greater than the heat loss to the environment. The code does not contain a specific 

particle burnout model, and so the final phases of particle combustion may be somewhat inaccurate, but 

considering how rapidly the particle burns out, it is unlikely that a more accurate modeling approach 

would yield a substantially different result. For all cases, burnout is assumed to occur when the particle 

density reaches 1 percent of the initial particle density, at which point the particle is removed from the 

simulation. 

The 100 μm graphite particles show similar behavior for the Mach 5, 6 and 7 total temperature 

conditions. In the constant-area low speed mixer, the first three meters, combustion is diffusion limited, as 

the particle and gas velocities are nearly identical. However, as the particle enters the convergent portion 

of the nozzle, a particle velocity lag is introduced, which greatly increases the transport of oxygen to the 

particle surface, and the temperature rapidly increases. The temperature rise then levels off, and for the 

Mach 6 case even drops off somewhat, based on the competing effects of lower free stream pressure in

the nozzle leading to lower oxygen concentration, decreased particle mass due to combustion, increased 

particle internal surface area, heat transfer and combustion kinetics. For the Mach 7 total temperature 

condition, the particle burned out just after entering the convergent portion of the nozzle. 
NASA/TM—2010-216765 60

Figure 56.—Mach 5 graphite versus alumina surface temperature comparison. 

Figure 57.—Mach 6 graphite versus alumina surface temperature comparison. 
NASA/TM—2010-216765 61

Figure 58.—Mach 7 surface temperature results. 

The 1000 μm particles show very little temperature change at any of the three simulated Mach 

numbers due to their size. The relatively low external surface area to particle mass ratio of the 1000 μm 

particles decreases the impact of radiation and convective heat transfer on particle temperature.

Simultaneously, intra-particle diffusion of oxygen and combustion products is greatly restricted, limiting 

combustion to the outer regions of the particle.

Figure 59 shows the internal particle temperature profiles for the 1000 μm particle at the exit plane of 

the mixer-nozzle section (MN) and the exit plan of the isolator (IS). The T/Tinsert values above 1.0 are 

indicative of ongoing particle combustion. 

Figures 60 to 62 show the velocity results for all three simulated Mach numbers and all three particle 

sizes, with the Mach 5 and Mach 6 results compared to their corresponding alumina particle cases. The 

Mach 6 and Mach 7 10 μm particle results are not readily discernable in these plots as the particles burn 

out so quickly that their trajectories do not differentiate themselves from the other cases. Since graphite is 

significantly less dense than alumina, it is not surprising to see that the graphite particles more closely 

follow the gas flow than the corresponding alumina particles do. This effect is accentuated due to mass 

loss from oxidation of the graphite particles, particularly for the 100 μm particles. In Figure 63, graphite 

particle porosity histories are given for each particle size and simulated Mach number. This plot shows 

that at each simulated Mach number, the 100 μm particles approach complete burnout, allowing the 

100 μm graphite particles to respond much more quickly to changes in flow velocity than the 

corresponding alumina particles. The 1000 μm graphite particles lose a relatively small amount of mass 

traversing the domain in each case, and so behave similarly to the alumina particles, but reach a higher 

velocity due to the density difference.
NASA/TM—2010-216765 62

Figure 59.—Radial temperature profiles of the 1000 μm particle for the 

Mach 5, 6 and 7 conditions at the exit planes of the mixer/nozzle (MN) 

and the isolator (IS). 

Figure 60.—Mach 5 graphite versus alumina velocity comparison. 
NASA/TM—2010-216765 63

Figure 61.—Mach 6 graphite versus alumina velocity comparison. 

Figure 62.—Mach 7 velocity results. 

It should be noted that for the 100 and 1000 μm graphite particles, the increase in porosity, 

corresponding to mass loss, does not increase as much as might be expected with simulated Mach 

number, despite the higher gas and particle temperatures. This is due to two causes. First, the particle 
NASA/TM—2010-216765 64

combustion is primarily diffusion limited throughout the particle trajectory. The higher temperature cases 

do not have significantly higher oxygen diffusion rates to and within the particle, so the mass loss cannot 

proceed at a much higher rate despite the increase in kinetic rate with temperature. Second, an increase in 

radiative heat transfer from the particle surface with increasing temperature tends to offset the internal 

heat generation from increased combustion rate.

To gain additional insight into the various factors effecting particle velocity and temperature profiles 

as a function of particle size, it is of interest to examine some additional parameters. Figure 64 give 

particle-gas relative Mach number (defined as Vgas–Vparticle divided by the local sound speed) data for all 

three particle sizes at the Mach 5 total temperature condition. The Mach 5 data was selected here because 

at the higher Mach numbers the 10 μm particle quickly burned out and therefore these values could not be 

compared throughout the particle trajectories. The relative Mach number traces shown in Figure 64 trend 

as expected, with the larger, heavier particles showing a larger velocity deficit. The 1000 μm particles 

have a supersonic relative Mach number through much of the nozzle, while the smaller particles never 

achieve a supersonic relative Mach number. While Mach number is generally defined as a scalar, for 

these Mach number calculations the relative Mach number is allowed to go negative when the particle 

velocity is greater than the gas flow velocity. Therefore, negative Mach numbers are shown in Figure 64, 

occurring in each case for the 10 μm particle after a shock jump-condition has been imposed on the gas 

flow and for the last two segments of the flowpath for the 100 μm particle.

Figures 65 to 67 show the different heat flows into, out of, and accumulating in the particle outer 

shell, as defined in Figure 3, for the three Mach 5 total temperature conditions. For the 10 and 100 μm 

particles, the “outer shell” actually encompasses the entire particle, since no internal grid was used. This 

means that the conduction term into or out of the outer shell is by definition zero, and is not shown. For 

the 1000 μm particle, a 101 point internal grid was used, so the outer shell has a thickness equal to 1 

percent of the particle radius (5 μm). It should be noted that in these figures, some of the heat flow terms 

have been multiplied by 10 to make them visible on the scale of the plot. 

Figure 63.—Graphite particle porosity.
NASA/TM—2010-216765 65

Figure 64.—Mach 5 gas-particle relative mach number. 

Looking at Figure 65, heat flows for the 10 μm particle, we see high initial heat flow rates as the 

particle cools from the particle insertion temperature (200 K higher than the gas temperature) down to the 

gas temperature. The particle quickly reaches steady state within the mixer section, up to the 3.0 m point, 

and the heat flows stay constant until the particle enters the nozzle section. As the gas temperature drops 

in the nozzle, combustion is quenched and the particle radiation drops to zero. At the same time, 

convection spikes, removing heat from the particle as shown by the negative heat accumulation term.

Finally, there is some particle heating across each of the oblique shocks and throughout the isolator 

section, located at approximately 6.0 m, 7.0 m, and 8.0 to 8.4 m, respectively. 

For the 100 μm particle, Figure 66, we see that the combustion term is always larger than the 

convection heat loss term (convection out of the particle is defined as being positive, as is radiation), 

resulting in continuous particle temperature rise, as the radiation heat loss remains small. In the area of 

the isolator, the sum of the convection and radiation heat flows do approach the value of the combustion 

term, so the temperature remains fairly constant, as shown in Figure 56. 
NASA/TM—2010-216765 66

Figure 65.—Heat flows for 10 μm graphite particle at Mach 5 total temperature condition. 

Figure 66.—Heat flows for 100 μm graphite particle at Mach 5 total 

temperature condition. 
NASA/TM—2010-216765 67

Figure 67.—Heat flows for 1000 μm graphite particle outer shell at Mach 5 

total temperature condition. 

Lastly, the heat flows (including internal conduction) for the 1000 μm graphite particle at the Mach 5 

condition are given in Figure 67. As this particle’s temperature changes very little throughout the 

simulation, we see that the radiative heat loss is essentially constant. As with the other particle sizes, the 

convection term increases rapidly entering the nozzle, but the actual change in particle temperature is 

small because of the large thermal mass represented by this relatively large particle. The ramp shocks and 

isolator section are discernable in Figure 67, but have only a very minor effect on the particle temperature 

profile. 

6.3 One Percent Mass Fraction Graphite Particle Simulations

The same simulations that were run for the inert particles and the single graphite particles were then 

run using a 1 percent mass fraction of graphite particles. All of the 1 percent mass fraction simulations 

were run using finite rate chemistry in the gas phase. As in the Mach 7 single graphite particle 

simulations, the hydrogen-containing species were eliminated from the kinetic mechanism to speed the 

calculations. Each 1 percent mass fraction simulation was run until a time invariant solution was obtained 

for both the particle and gas streams. 

Figures 68 to 76 show the temperature and velocity histories of 10, 100, and 1000 μm graphite 

particles and the corresponding gas stream solutions calculated using a 1 percent mass fraction of 

particles at each of the three simulated Mach numbers. These figures also contain the single graphite 

particle results for comparison. 

For the 10 μm particles at the Mach 5 total temperature flow condition, Figure 68, there is virtually no 

difference between the single particle case and the 1 percent mass fraction case. Since combustion of the 

particle is almost immediately quenched by the cooler surrounding gas, there is no impact on the gas 

stream that could cause a change in the overall solution. The insertion into the gas stream of a 1 percent
NASA/TM—2010-216765 68

ass fraction of particles at a temperature 200 K higher than the gas stream is insufficient to cause a 

significant change in the gas stream variables, as was observed with the inert alumina particles. At the 

Mach 6 and Mach 7 total temperature flow conditions, Figures 69 and 70, the 10 μm graphite particles 

undergo very rapid combustion. The heat release from the particles significantly raises the gas stream 

temperature in the mixer. This increase in mixer temperature results in an increased isolator exit 

temperature, as well as increased velocity in sections of the flow path. The increase in velocity is due to 

an increase in the sound speed due to the higher static temperature in regions of approximately constant 

Mach number. The somewhat lower relative temperature increase in the Mach 7 case, Figure 70, is due 

primarily to the increase in gas specific heat with increasing temperature. 

Figure 68.—Comparison of particle temperature and velocity profiles between a 

single 10 μm graphite particle and a 1 percent mass fraction of 10 μm graphite 

particles at the Mach 5 total temperature flow condition.
NASA/TM—2010-216765 69

Figure 69.—Comparison of particle temperature and velocity profiles 

between a single 10 μm graphite particle and a 1 percent mass fraction of 

10 μm graphite particles at the Mach 6 total temperature flow condition.

Figure 70.—Comparison of particle temperature and velocity profiles 

between a single 10 μm graphite particle and a 1 percent mass fraction of 

10 μm graphite particles at the Mach 7 total temperature flow condition. 
NASA/TM—2010-216765 70

Comparing the results for the Mach 5 and Mach 6 total temperature flow condition for the 100 μm 

particles, Figures 71 and 72, similar results are observed. While the single particles do not burn out within 

the computational domain, the higher gas temperature for the 1 percent mass fraction simulations from 

graphite combustion leads to more rapid particle consumption and eventual burnout. As would be 

expected, the particles burn out more rapidly at the higher Mach 6 total temperature condition. 

Some unique features are visible in the Mach 7 results for the 100 μm graphite particles, as shown in 

Figure 73. Here the single graphite particle and the 1 percent graphite particle stream burn out at 

approximately the same point. This occurs despite the increasing gas temperature found in the 1 percent 

mass fraction results. This occurs because at this very high temperature condition the combustion process 

is completely diffusion limited. Increasing the gas temperature does not translate into a higher particle 

combustion rate. 

Lastly, in Figures 74 to 76, little difference is observed between the single particle solutions and the 

1 percent mass fraction solutions for the 1000 μm particles at each of the three simulated Mach numbers. 

Since relatively little combustion is occurring in these cases, the gas streams are unaffected by the 

presence of the 1 percent mass fraction of particles.

Figures 77 to 79 show the particle porosity histories for the single graphite particles and the 1 percent 

mass fraction graphite particle stream as a function of particle size and simulated Mach number. 

Figure 71.—Comparison of particle temperature and velocity profiles between 

a single 100 μm graphite particle and a 1 percent mass fraction of 100 μm 

graphite particles at the Mach 5 total temperature flow condition. 
NASA/TM—2010-216765 71

Figure 72.—Comparison of particle temperature and velocity profiles between 

a single 100 μm graphite particle and a 1 percent mass fraction of 100 μm 

graphite particles at the Mach 6 total temperature flow condition.

Figure 73.—Comparison of particle temperature and velocity profiles between 

a single 100 μm graphite particle and a 1 percent mass fraction of 100 μm 

graphite particles at the Mach 7 total temperature flow condition.
NASA/TM—2010-216765 72

Figure 74.—Comparison of particle temperature and velocity profiles between 

a single 1000 μm graphite particle and a 1 percent mass fraction of 

1000 μm graphite particles at the Mach 5 total temperature flow condition.

Figure 75.—Comparison of particle temperature and velocity profiles between 

a single 1000 μm graphite particle and a 1 percent mass fraction of 

1000 μm graphite particles at the Mach 6 total temperature flow condition.
NASA/TM—2010-216765 73

Figure 76.—Comparison of particle temperature and velocity profiles between 

a single 1000 μm graphite particle and a 1 percent mass fraction of 

1000 μm graphite particles at the Mach 7 total temperature flow condition.

Figure 77.—Comparison of particle porosity profiles between single graphite 

particles and a 1 percent mass fraction of graphite particles at the Mach 5 

total temperature flow condition.
NASA/TM—2010-216765 74

Figure 78.—Comparison of particle porosity profiles between single graphite 

particles and a 1 percent mass fraction of graphite particles at the Mach 6 

total temperature flow condition.

Figure 79.—Comparison of particle porosity profiles between single graphite 

particles and a 1 percent mass fraction of graphite particles at the Mach 7 

total temperature flow condition. 
NASA/TM—2010-216765 75

Figure 80 shows the NO mass fraction profile for all the 1 percent graphite particle mass fraction 

cases where NO is produced above a mass fraction of 1×10–6 at any point in the computational domain. 

All of the Mach 5 total temperature flow condition cases and the Mach 6 total temperature flow condition 

case with no particles have NO levels below this threshold. The Mach 6 cases with particles all have very 

low NO mass fractions which only make it above the 1×10–6 mass fraction threshold because of particle 

combustion. However, all of the Mach 7 simulations have significant mass fractions of NO. The NO is 

generated by the high temperature in the mixing section, particularly when enhanced by particle 

combustion. The NO level decreases as the temperature drops with expansion, but the NO reactions are 

eventually “frozen” in the expansion section of the nozzle due to the rapid temperature drop along with 

relatively slow nitrogen chemistry. Looking at Figures 81 and 82, the NO2 and O mass fraction profiles, it 

is evident that a small mass fraction of O is generated by the high temperatures in the mixing section in 

addition to the NO, and that this atomic oxygen combines with the NO to produce NO2 during the 

expansion just upstream and downstream of the nozzle throat. These NO2 mass fractions are then also 

“frozen” as they progress downstream. Prominent in Figure 81 is a spike in the NO2 mass fraction early in 

the mixer for the 100 μm graphite particles at the Mach 7 total temperature flow condition. This spike 

represents a region wherein the kinetics favors the formation of NO2 from the NO and O formed at the 

initial high temperature. As the temperature increases due to graphite combustion, the kinetics transitions 

to favor the formation of NO and O instead of NO2, and the NO2 previously formed dissociates back to 

NO and O. In general, it can be concluded from these figures that the more rapid the particle combustion 

process, the higher the levels of NO and O that are generated, and therefore an increased mass fraction of 

NO2 is generated in the expansion.

Figure 80.—Comparison of NO mass fraction profiles between no graphite particles 

and a 1 percent mass fraction of graphite particles.

Mixer Nozzle

Precompression Plate

Final Inlet Ramp

Isolator
NASA/TM—2010-216765 76

Figure 81.—Comparison of NO2 mass fraction profiles between no graphite 

particles and a 1 percent mass fraction of graphite particles. 

Figure 82.—Comparison of O mass fraction profiles between no graphite 

particles and a 1 percent mass fraction of graphite particles.
NASA/TM—2010-216765 77

Figures 83 and 84 show the spatial distribution of CO2 and CO, respectively. As expected, the 

combinations of particle size and simulated Mach number that yield little or no particle combustion show 

low levels of CO2 and CO. This would include all the 1000 μm particle cases and the Mach 5, 10 μm 

particle case. Low temperature graphite combustion favors CO production at the particle surface and 

inhibits CO oxidation in the gas phase, as seen with the Mach 5 total temperature, 100 μm graphite 

particle case. As the temperature increases due to more rapid graphite combustion, such as with the Mach 

6 total temperature, 10 and 100 μm particle cases, CO is initially produced by the particle, but this 

transitions to CO2 production with rising particle temperature. With increasing gas phase temperature, the 

CO is oxidized to CO2. This gas phase CO oxidation takes much longer for the slower burning Mach 6, 

100 μm particle case, and is not completed before the nozzle expansion lowers the flow temperature. At 

that point, CO production is again favored over CO2 at the particle surface, and the CO level rises through 

the nozzle. This rise in the CO mass fraction for the Mach 6, 100 μm case in the nozzle is also contributed 

to by a small amount of CO2 dissociation, as indicated by the drop in CO2 concentration in Figure 83. 

Figure 83.—Comparison of CO2 mass fraction profiles between 1 percent 

graphite particle mass fraction cases.
NASA/TM—2010-216765 78

Figure 84.—Comparison of CO mass fraction profiles between 1 percent 

graphite particle mass fraction cases.

6.4 Ignition Calculations 

 As a final step in examining the effect of graphite particulate vitiation on hypersonic propulsion 

testing, a series of ignition delay calculations were performed based on the isolator section exit flow 

conditions for each 1 percent graphite particle mass fraction case. Hydrogen was used as the fuel at an 

equivalence ratio of 0.8 based upon the amount of molecular oxygen present at the exit of the isolator. 

The computational domain generally consisted of a constant area duct 0.5 meters long, open at each end. 

450 equally spaced grid points were used in each model. In some cases the length of the domain was 

extended to 1.0 m if ignition did not occur within the 0.5 m length. A supersonic inlet boundary condition 

was used with the pressure, temperature, velocity, and species mass fractions taken from the isolator exit 

conditions, as given in Tables 6 and 7. If particles were present at the isolator exit, they were included in 

the ignition delay calculation at their isolator exit condition, Table 8, and were allowed to continue to 

react based on the local flow conditions. The exit boundary conditions remained supersonic for each case.

The initial condition within the duct was 100 percent N2 at the temperature, pressure, and velocity of the 

isolator exit. Pure nitrogen gas was used as the initial condition for the calculation to help distinguish the 

location of the interface between the gas coming from the in-flow boundary versus and the gas initially in 

the duct for this plug-flow simulation. Ignition was defined as the point in time when a 200 K temperature 

rise above the inlet boundary condition was observed within the computational domain. Figure 85 shows 

a typical result at the point of ignition. 

 TABLE 6.—GAS TEMPERATURES (K) AT THE EXIT OF THE ISOLATOR SECTION 

Mach no. No particles 10 μm 100 μm 1000 μm

5 798 800 908 802

6 1008 1229 1201 1022

7 1264 1452 1441 1278
NASA/TM—2010-216765 79

TABLE 7.—SPECIES MASS FRACTIONS AT THE EXIT OF THE ISOLATOR SECTION

Mach no./

particle size

O2 N2 NO NO2 CO2 CO O

5 / None 0.233000 0.767000 ---------- ---------- ---------- ---------- ----------

6 / None 0.233000 0.767000 ---------- ---------- ---------- ---------- ----------

7 / None 0.229035 0.763448 0.007340 0.000177 ---------- ---------- ----------

5 / 10 0.232141 0.766549 ---------- ---------- 0.001175 0.000135 ----------

6 / 10 0.207284 0.757233 0.000347 0.000018 0.035111 0.000007 ----------

7 / 10 0.196415 0.743837 0.019493 0.000323 0.035401 0.000141 ----------

5 / 100 0.218665 0.756865 ---------- ---------- 0.002132 0.022322 0.000016

6 / 100 0.209855 0.756875 0.000084 0.000008 0.026443 0.006705 ----------

7 / 100 0.196023 0.744401 0.017977 0.000304 0.035747 0.000668 ----------

5 / 1000 0.232698 0.766849 ---------- ---------- 0.000271 0.000182 ----------

6 / 1000 0.232200 0.766511 0.000002 ---------- 0.000359 0.000923 0.000005

7 / 1000 0.227262 0.762568 0.007643 0.000183 0.001472 0.000872 ---------- 

TABLE 8.—PARTICLE SURFACE TEMPERATURE, VELOCITY, AND POROSITY 

FOR SURVIVING GRAPHITE PARTICLES AT THE EXIT OF THE ISOLATOR SECTION

Mach no./

particle size

Surface temperature 

(K)

Velocity 

(m/s)

Porosity

5 / 10 625 1115 0.264

5 / 1000 1387 785 0.247

6 / 1000 1911 941 0.273

7 / 1000 2421 1156 0.294

Figure 85.—Example of ignition point for a H2/vitiated air/particle mixture from the

Mach 6 total temperature simulation with a 1 percent mass fraction of 1000 μm 

graphite particles.
NASA/TM—2010-216765 80

The calculated ignition delay times for each condition with a 1 percent initial mass fraction of 

graphite particles, along with the ignition delay times for simulations done without any particles, are 

given in Table 9. In those cases where little graphite combustion occurred, only a slight reduction in 

ignition is observed, or even a slight increase in ignition delay (due to the heat-sink effect of the particles 

during the ignition delay calculation) for the Mach 5, 10 μm particle case. However, in the calculations 

where significant graphite combustion has occurred, the ignition delay time is greatly reduced from the 

cases with no particles. To examine this effect further, a second set of ignition delay calculations were

performed to determine if the decrease in ignition delay time was due to vitiate species being present 

(including particles) or graphite combustion heat release. These calculations were performed at the 

elevated temperatures given in Table 6, but using no gas phase or particulate vitiates, only regular air.

These results are given in Table 10. Comparing the results from Tables 9 and 10, it is evident that the 

majority of the graphite particle vitiation effect for these conditions is simply heat release. The presence 

of the vitiate species in the flow generally serves to slightly retard ignition. The relatively small impact of 

the combined vitiate species is consistent with other studies, such as Han, et al. (2004), where the 

significant effects of known ignition enhancing species such as NO are generally seen at lower 

temperatures and higher pressures than those studied here, and where the expected effects of a particular

vitiate is frequently cancelled out by the effects of other vitiate species when they are studied in 

combination. The purpose of these ignition delay calculations is only to look at the combined effect of all 

the vitiates. Further studies of individual vitiate species effects would be needed to determine more 

detailed information, but such studies are not warranted here. 

TABLE 9.—IGNITION DELAY TIME IN MICROSECONDS WITH VITIATES

Mach no. No particles 10 μm 100 μm 1000 μm

5 1105 1112 367 1055

6 183 82.7 91.2 169

7 56.1 43.1 41.9 52.1

TABLE 10.—IGNITION DELAY TIME IN MICROSECONDS WITHOUT VITIATES

Mach no. No particles 10 μm 100 μm 1000 μm

5 1105 1094 353 1055

6 183 75.8 84.4 169

7 56.0 39.1 37.8 51.7

7.0 Conclusions

In order to improve the understanding of particulate vitiation in hypersonic propulsion test facilities, a 

quasi-one dimensional numerical tool has been developed capable of efficiently modeling reacting 

particle-gas flows over a wide range of gas and particle conditions. Features of this code include finite 

rate kinetics for gas phase chemistry, a global porous-particle combustion model for the solid phase, 

mass, momentum and energy interactions between phases, subsonic and supersonic particle drag and heat 

transfer models, variable time step, and appropriate sets of subsonic and supersonic boundary conditions.

The basic capabilities of this tool have been validated against available literature data or other validated 

codes. 

This code has been demonstrated through a parametric study of particle vitiation effects in a model 

hypersonic propulsion test facility and scramjet. Parameters studied were simulated flight Mach number

(Mach 5, 6, and 7), particle size (10, 100, and 1000 μm diameters), particle mass fraction (single particles

and 0.01) and particle material (alumina and graphite). For the inert alumina particles, it was found that 

the presence of particles up to the 1 percent mass fraction had very little effect on the gas phase. As a 1

percent mass fraction would in reality be an extreme case, it is safe to conclude that hot wall ignition on 

the particle surface is likely to be the only potential inert particle vitiation effect, an effect that was 
NASA/TM—2010-216765 81

beyond the scope of this study. It was observed that only the 10 μm alumina particles tracked the gas flow 

reasonably well, and that when the particle size was increased to 1000 μm, the particle responded very 

little at all to changes in the gas flow, including shock waves, even though supersonic relative Mach 

numbers are observed. This is to be expected as the large particles have significant inertia against velocity 

changes, and have a relatively large volume-to-surface-area ratio that decreases heat transfer effects with 

the surroundings. 

With the graphite particles, it was observed that the smallest particles in the study will either quickly 

quench, or quickly be consumed, depending on the gas and particle temperatures. The combustion of 

these 10 μm particles is solely determined by the particle temperature-dependent kinetic rates and not by 

oxygen diffusion. As the particle size increases, the particle is less likely to quench, because it has a larger 

thermal inertia, and so will remain at a high temperature long enough for the particle to ignite. As the 

particle and surrounding gas temperatures increase, the particle consumption rate increases somewhat, but 

not as fast as might be expected. This is due to diffusion rate limitations, both to the particle surface and 

internal to the particle. For the largest particles, combustion is significantly diffusion limited, so particle 

and gas temperature have little effect on the combustion rate. As the particle mass fraction is stepped up 

to 1 percent, the main change is the addition of significant heat release. As the choice of 1 percent mass 

fraction was meant to be a limiting case, the amount of heat release is extreme and not expected to be 

observed in real facilities. However, it does serve to increase the particle combustion rate, and to increase 

the quantity of NOx generated in the flow. In those cases where low reaction rates of graphite were 

observed for single particles, the increase to 1 percent mass fraction had very little impact, similar to the 

inert particles. Ignition delay time calculations for the 1 percent mass fraction of graphite particles cases 

showed significant decreases in ignition delay in the cases where higher graphite combustion rates, and 

therefore higher gas temperatures, were observed, but further calculations showed that this was due 

primarily to the increased combustor inlet temperature, not the vitiate species present in the flow. What

chemical effects were present appeared to slightly increase the ignition delay. A review of several 

summary papers on gas-phase vitiate effects did little to clarify this result, as the findings are typically 

highly case-specific and not easily compared. 

While the results discussed above are by no means conclusive regarding the potential effects of 

particle vitiation, they do provide significant insight into the physical effects of gas flow conditions on the 

particles, and of the particles on the gas flow. Understanding of the relative effects of the particle drag, 

the various forms of heat transfer, and combustion throughout the test facility and test engine is an 

important first step. The tool developed as part of this effort can be useful in determining the likely effects 

of specific particle contamination situations at specific test conditions, when the required modeling inputs 

are made available. This numerical tool is also generally useful for modeling any gas phase or gas-solid 

two-phase fluid flow problem that can be modeled as quasi-one dimensional. 

8.0 Recommendations for Future Work

During the course of demonstrating the numerical tool developed for this study, a number of potential 

improvements were noted. These can be broadly categorized into three areas, as follows.

8.1 Modeling Improvements

1. Particle burnout model.—A particle burnout model should be added to the code if medium size 

particles are to be modeled. This is not important for very small particles which burn out quickly or for 

very large particles which suffer low mass loss.

2. Particle diameter shrinkage model.—This is another model that would be most critical for medium 

size particles. For smaller particles that burn out or quench quickly, the impact of incorrectly modeling 

the particle shrinkage is small, as there would not be expected to be a significant difference in the 

numerical solution at any significant distance downstream of the insertion point. For large particles, little 
NASA/TM—2010-216765 82

diameter change would be expected unless extremely long residence times were being considered. A 

particle diameter shrinkage model would need to be based on empirical data for the material, particle 

geometry, and test conditions under study. 

3. Multiple simultaneous particle sizes.—The code currently only models a single particle size for a 

given run. The capability for handling multiple particle sizes simultaneously would enhance the ability of 

the code to model real situations. 

4. Hot wall ignition model.—The ability to examine gaseous fuel-oxidizer ignition at the surface of an 

individual hot, reacting particle could be important. This type of modeling is most likely beyond the 

capability of a quasi-one dimensional code. A multi-dimensional code that addresses multi-species 

diffusion to and from the particle surface and particle surface chemistry would be needed to complement 

the current capability.

8.2 Input Data Improvements

1. Graphite Combustion Data.—For accurate results with reacting particles, it is necessary to have 

graphite combustion data, including combustion rate and particle shrinkage, across the full range of 

conditions to be experienced by the particles in the model. Such data could then be applied to the global 

reaction model in the code. 

2. Good particle characterization.—For modeling a real system, it is necessary to have characterized 

the particle contamination in the facility. Particle size distribution and particle mass fraction must be 

known. Also the particles themselves need to be fully characterized in terms of overall porosity and 

internal pore size.

3. Validation Data.—It was not possible to validate the quasi-one dimensional code as a whole 

because of a lack of relevant validation data. Each capability within the code was validated separately, to 

the greatest extent possible. Gathering applicable validation data would greatly enhance confidence in the 

accuracy of the numerical results.

8.3 Additional Modeling

1. Lower Particle Mass Fractions for Reacting Particles.—The upper limit of 1 percent particle mass 

fraction used for this study was unrealistically high, although it does provide an upper bound. Cases could 

be run at various particle mass fractions below this level to understand the effect of particle vitiation at 

more realistic levels for those facility test conditions at which particle vitiation effects were observed.

2. Additional Particle Sizes for Reacting Particles.—In this demonstration study, different modes of 

behavior were observed at the different particle sizes. Running more studies in the intermediate size range 

would be of significant interest, since it is in this size range that the particle behavior is the least 

predictable. 
NASA/TM—2010-216765 83

Appendix A.—Symbols List

cp Gas phase mixture specific heat

c ,ip Specific heat of gaseous species i

d pore Particle internal pore diameter

eg Gas phase specific internal energy

cf Effectiveness factor coefficient

o g f Gas phase Gibbs energy of formation

o g ,if Gibbs energy of formation of gaseous species i

hg Gas phase specific enthalpy

h ,ig Specific enthalpy of gaseous species i

hp Particle convection coefficient

hto Gas phase enthalpy of formation 

hto,i Enthalpy of formation of gaseous species i

kb Backward reaction rate

k f Forward reaction rate

kg Gas thermal conductivity

k1p 1st kinetic rate for carbon combustion 

mobs Observed carbon combustion reaction order 

n Arrhenius reaction temperature exponent

ns Number of chemical species

p Pressure

qaccum Heat accumulation within the particle outer shell 

qc Gas phase conductive heat flux 

qcond Conductive heat flow into the particle from the particle outer shell 

qconv Convective heat flow out of the particle outer shell

qcomb Combustion heat addition in particle outer shell

qrad Radiative heat flow out of the particle outer shell

qw Wall heat flux

r Particle radial coordinate

er Particle combustion rate based on particle outer shell energy balance 

kr Mass-averaged particle combustion rate based on carbon kinetics

ikr , Local internal particle combustion rate based on carbon kinetics

ox r Particle combustion rate based on oxygen diffusion 
NASA/TM—2010-216765 84

pr Particle radius

±r Inner product used to calculate artificial viscosity

t Time

ug Gas phase velocity 

up Particle velocity

x Axial distance coordinate

xp Particle axial position

A Cross sectional area

Af Arrhenius reaction pre-exponential factor

C Artificial viscosity interim variable

CD Particle drag coefficient

Cs Particle material heat capacity

De Particle internal diffusion coefficient

Di Diffusion coefficient of gaseous species i

D , ji Binary diffusion coefficient between gaseous species i and j

DK Knudsen diffusion coefficient

D 2 ,O M Molecular diffusion coefficient of oxygen into gas mixture 

Ea Arrhenius reaction activation energy

F Array of convective terms from the gas phase governing equations

Fp Drag force on a particle

H1 Scramjet inlet throat height

J Array of source terms from the gas phase governing equations

KD Mass transfer coefficient

Keq Equilibrium constant 

KP 2nd kinetic rate for carbon combustion

Ks Particle thermal conductivity

K± Artificial viscosity interim variable

Li Scramjet inlet and isolator ramp lengths

M Mixture average molecular weight

Mc Carbon atomic mass

Mi Species molecular weight

MR Ratio of moles of CO to moles of CO2 produced in carbon combustion 

Mrel Mach number based on particle/gas relative velocity

MSR Molecular speed ratio

MT Thiele modulus 
NASA/TM—2010-216765 85

N Number of chemical reactions 

Nud Nusselt number based on particle diameter

P ,O fs 2 Free stream partial pressure of oxygen

P ,O s 2 Partial pressure of oxygen at the particle surface

Pr Prandtl number 

R Mixture gas constant

Rerel Particle Reynolds number based on particle/gas relative velocity

Ru Universal gas constant

Se Energy equation energy source term

Sext Particle external surface area

Sh Sherwood number 

Si Reactant or product species

Sint Particle internal surface area

Sm Continuity equation mass source term 

S ,im Species continuity equation mass source term

Sp Momentum equation momentum source term 

Stot Particle total surface area

Tfilm Particle film temperature

Tg Gas temperature

T * , jig Reduced gas temperature for collision between species i and j

Tinsert The temperature at which a particle is inserted into the flow stream

Tp Particle surface temperature or bulk temperature for uniform temperature particles

T ,ip Particle temperature at internal grid point i

Trec Convective recovery temperature

Twall Facility mixer/nozzle wall temperature

Tε,i Effective temperature of gaseous species i

Tε ,, ji Average effective temperature between gaseous species i and j

T * 2,2 Reduced temperature for mixture viscosity

U Array of dependent variables from the gas phase governing equations 

Xi Mole fraction of gaseous species i

Yi Mass fraction of gaseous species i

α Ratio of particle mass flow to gas mass flow

αi Scramjet inlet/isolator ramp lengths 

βi Scramjet inlet/isolator shock angles
NASA/TM—2010-216765 86

χ Artificial viscosity

ε Particle porosity

εrad Particle emissivity

φ Gas phase volume fraction

φ , ji Averaging variable for mixture viscosity

γ Ratio of specific heats

η Particle internal diffusion effectiveness factor

ηavg Mass-averaged particle internal diffusion effectiveness factor

λ Artificial viscosity velocity eigenvalue

μ Gas phase viscosity

μi Viscosity of gaseous species i

νc Courant number 

νi Diffusion velocity of gaseous species i

νi

′ Stoichiometric coefficient of reactant species i

νi

′′ Stoichiometric coefficient of product species i

νo Stoichiometric coefficient of oxygen consumed per mole of carbon consumed 

θ Artificial viscosity flux limiter

ρg Gas phase density

ρs Particle material density

σi Effective collision diameter of gaseous species i

σ , ji Average effective collision diameter between species i and j

σ −BS Stefan-Boltzmann constant 

ζ Oxygen diffusion coefficient 

∆t Numerical time step

∆r Particle internal grid size

ΩD , ji Collision integral between species i and j

Ω 2,2 Collision integral for mixture viscosity 

ξ Intermediate variable for carbon combustion rate calculation 

Ψ Ratio of moles of CO2 to total moles produced in carbon combustion 
NASA/TM—2010-216765 87

Appendix B.—Gas Phase Kinetic Mechanism

The detailed gas phase kinetic mechanism used for this study is given in Table B.1. While the 

majority of the reactions follow the simple Arrhenius form, there are a number that have additional 

features. Ten of the reactions required the participation of a third body, for which the species third body 

efficiencies are included in Table B.2. Reactions 18, 19, 23, and 25 include a second duplicate reaction 

with different Arrhenius constants. The sum of the two reactions provides the correct temperature 

dependent behavior. Lastly, reactions 12, 14, and 24 are pressure dependent reactions. The first set of 

Arrhenius constants represents the high pressure limit, while the second set represents the low pressure 

limit. These two reaction rates are then combined via the Troe method, using the constants provided, for 

reactions 12 and 14, and via the Lindemann method for reaction 24. These methods are described in detail 

by Kee, et al (1993). 

TABLE B.1.—GAS-PHASE KINETIC MECHANISMc

Reaction 

ref. no.

Reaction Af n Ea Reaction

sourcea

1 2 +⇔+ OHOOH 2.644×1016 –0.6707 17041 05Davis

2 2 +⇔+ OHHHO 4.589×104 2.70 6260 05Davis

3 2 +⇔+ 2OHHHOH 1.734×108 1.51 3430 05Davis

4 +⇔ 2OHO2OH 3.973×104 2.40 –2110 05Davis

5 +⇔++ MH2MHH 1.78×1018 –1.0 0.0 05Davis

6 2 ⇔+ 2HH2H 2 9.0×1016 –0.600 0.0 05Davis

7 2 +⇔+ 22 OHHOHH2 5.624×1019 –1.25 0.0 05Davis

8 2 +⇔+ CO2HCOH2 2 5.50×1020 –2.0 0.0 05Davis

9 2 +⇔++ MOHMOHH 4.40×1022 –2.0 0.0 05Davis

10 +⇔++ MOHMHO 9.428×1018 –1.0 0.0 05Davis

11 2 +⇔++ MOMOO 1.20×1017 –1.0 0.0 05Davis

12 2 ( ) ( ) 2 +⇔++ MHOMOH 5.116×1012 0.44 0.0 05Davis

Low 6.328×1019 –1.4 0.0

TROE 0.5 1×10–30 1×1030

13 22 2 +⇔+ HHOOH 5.916×105 2.433 53502 05Davis

14b ( ) ( ) 22 +⇔+ MOHMOH2 1.11×1014 –0.37 0.0 05Davis

Low 2.01×1017 –0.584 –2293

TROE 0.0 345 10

15 2 +⇔+ 2OHOHHO 3.97×1012 0.0 671 05Davis

16 2 ⇔+ OH2HHO 7.485×1013 0.0 295 05Davis

17 2 +⇔+ OOHOHO 2 4×1013 0.0 0.0 05Davis

18a 2 +⇔+ 22 OHOOHHO 2.375×1013 0.0 –500 05Davis

18b 2 +⇔+ 22 OHOOHHO 1×1016 0.0 17330 05Davis

19a +⇔ OHOHO2 2222 1.3×1011 0.0 –1630 05Davis

19b +⇔ OHOHO2 2222 3.658×1014 0.0 12000 05Davis

20 22 +⇔+ HHOHOH 22 6.05×106 2.0 5200 05Davis

21 22 +⇔+ 2OHOHHOH 2.41×1013 0.0 3970 05Davis

22 22 +⇔+ HOOHOOH 2 9.63×106 2.0 3970 05Davis

23a 22 OHOH +⇔+ 22 OHHO 2×1012 0.0 427 05Davis

23b 22 OHOH +⇔+ 22 OHHO 2.67×1041 –7.0 37600 05Davis

24 ( ) ( ) 2 +⇔++ MCOMOCO 1.362×1010 0.0 2384 05Davis

Low 1.173×1024 –2.79 4191

25a +⇔+ H2COOHCO 8×1011 0.14 7352 05Davis

25b +⇔+ H2COOHCO 8.784×1010 0.03 –16 05Davis
NASA/TM—2010-216765 88

TABLE B.1.—GAS-PHASE KINETIC MECHANISMc

Reaction 

ref. no.

Reaction Af n Ea Reaction

sourcea

26 2 2 +⇔+ OCOOCO 1.119×1012 0.0 47700 05Davis

27 2 2 +⇔+ OHCOHOCO 3.01×1013 0.0 23000 05Davis

28 +⇔+ HCOHHCO 2 1.2×1014 0.0 0.0 05Davis

29 +⇔+ OHCOOHCO 3×1013 0 0 05Davis

30 2 +⇔+ HCOOHCO 3×1013 0 0 05Davis

31 +⇔+ 2OHCOOHHCO 3.02×1013 0 0 05Davis

32 ++⇔+ MHCOMHCO 1.87×1017 –1 17000 05Davis

33 2 ++⇔+ O2HHCOOHHCO 2.244×1018 –1 17000 05Davis

34 2 +⇔+ HOCOOHCO 2 1.204×1010 .807 –727 05Davis

35 2 +⇔+ ONNON 2.7×1013 0 355 GRI3.0 

36 2 +⇔+ ONOON 9×109 1 6500 GRI3.0 

37 +⇔+ HNOOHN 3.36×1013 0 385 GRI3.0 

38 2 2 +⇔+ OHNONOHO 2.11×1012 0 –480 GRI3.0 

39 2 +⇔++ MNOMONO 1.06×1020 –1.4 0 GRI3.0 

40 2 +⇔+ ONOONO 2 3.9×1012 0 –240 GRI3.0 

41 2 +⇔+ OHNOHNO 1.32×1014 0 360 GRI3.0 

42 +⇔++ MHNOMNOH 4.48×1019 –1.3 740 GRI3.0 

43 +⇔+ OHNOOHNO 2.5×1013 0 0 GRI3.0 

44 2 +⇔+ NOHHHNO 9×1011 .7 660 GRI3.0 

45 +⇔+ 2OHNOOHHNO 1.3×107 1.9 –950 GRI3.0 

46 2 2 +⇔+ NOHOOHNO 1×1013 0 13000 GRI3.0 

47 2 +⇔+ CONOCON 3×1012 0 11300 GRI3.0 a 05Davis refers to the mechanism developed by Davis, et al (2005), while GRI3.0 refers to the GRI 3.0 (2005) mechanism. b For this reaction, the optional fourth Troe constant is used with a value of 345.

c The units for the pre-exponential factor Af are Kg/m2

-s-atm, while the units for the activation energy Ea are cal/g-mole. The temperature 

exponent n is unit-less. 

TABLE B.2.—THIRD BODY EFFICIENCIES 

Species React

5

React

9

React

10

React

11

React

12

React

14

React

24

React

32

React

39

React

42

H2 0 2 2 2.4 .75 2 2 2 1 1

O2 1 1 1 1 .85 1 1 1 1 1

H2O 0 6.3 12 15.4 11.89 6 12 0 1 1

OH 1 1 1 1 1 1 1 1 1 1

H 1 1 1 1 1 1 1 1 1 1

O 1 1 1 1 1 1 1 1 1 1

HO2 1 1 1 1 1 1 1 1 1 1

H2O2 1 1 1 1 1 1 1 1 1 1

N 1 1 1 1 1 1 1 1 1 1

NO 1 1 1 1 1 1 1 1 1 1

NO2 1 1 1 1 1 1 1 1 1 1

HNO 1 1 1 1 1 1 1 1 1 1

CO2 0 3.6 3.6 3.6 2.18 3.6 3.6 3.6 1 1

CO 1 1.75 1.75 1.75 1.09 1.75 1.75 1.75 1 1

HCO 1 1 1 1 1 1 1 1 1 1

Ar .63 .38 .7 .83 .4 .7 .7 1 1 1 

N2 1 1 1 1 1 1 1 1 1 1
NASA/TM—2010-216765 89

References

Ames Research Staff “Equations, Tables, and Charts for Compressible Flow,” NASA Report 1135, 1953. 

Anderson, Jr., J.D. Computational Fluid Dynamics, McGraw-Hill, Inc., New York, 1995. 

Baek, S.W. “The Ignition of Dust Particles Behind Moving Shock Waves,” Ph.D. Dissertation, University 

of Michigan, Ann Arbor, MI, 1985. 

Bailey, A.B.; Hiatt, J. “Free-Flight Measurements of Sphere Drag at Subsonic, Transonic, Supersonic, 

and Hypersonic Speeds for Continuum, Transition, and Near-Free-Molecular Flow Conditions,” 

AEDC –TR–70–291, March 1971. 

Blake, T.R. “Low Reynolds Number Combustion of a Spherical Carbon Particle,” Combustion and 

Flame, Vol. 129, No. 1 / 2, pp. 87–111, April 2002. 

Carlson, D.J.; Hoglund, R.F. “Particle Drag and Heat Transfer in a Rocket Nozzle,” AIAA Journal, Vol. 

2, No. 11, pp. 1980–1984, November, 1964.

Chang, I.S. “One- and Two-Phase Nozzle Flows,” AIAA Journal, Vol. 18, No. 12, pp. 1455–1461; 

December 1980. 

Chase, M.W.; Davies, C.A.; Downey, J.R.; Frurip, D.J.; McDonald, R.A.; Syverud, A.N. JANAF 

Thermochemical Tables, Third Edition, American Institute of Physics, New York, 1986. 

Crowe, C.T. “Drag Coefficient of Particles in a Rocket Nozzle,” AIAA Journal, Vol. 5, No. 5, pp. 1021–

1022; May 1967.

Davis, S.F. “A Simplified TVD Finite Difference Scheme Via Artificial Viscosity,” SIAM Journal on 

Scientific and Statistical Computing, Vol. 8, No. 1, pp. 1–18, January 1987. 

Davis, S.G.; Joshi, A.V.; Wang, H.; Egolfopoulos, F. “An Optimized Kinetic Model of H2/CO 

Combustion,” Proceedings of the 30th International Symposium of The Combustion Institute, pp. 

1283–92, Chicago, IL, July 2002. 

Di Giacinto, M.; Sabetta F.; Piva, R. “Two-way Coupling Effects in Dilute Gas-Particle Flows,” 

Transactions of the ASME, Vol. 104, pp. 304–312, September 1982. 

Drake, R.M.; Backer, G.H. “Heat Transfer from Spheres to a Rarefied Gas in Supersonic Flow,”

Transactions of the ASME, Vol. 74, pp. 1241–1249; October 1952. 

Drummond, J.P. “Supersonic Reacting Internal Flow Fields,” NASA TM–103480, 1990. 

Egolfopoulos, F.N.; Campbell, C.S. “Dynamics and Structure of Dusty Reacting Flows: Inert Particles in 

Strained, Laminar, Premixed Flames,” Combustion and Flame, Vol. 117, No. 1 / 2, pp. 206–226, 

April 1999. 

Fox, T.W.; Rackett, C.W.; Nicholls, J.A. “Shock Wave Ignition of Magnesium Powders,” Proceedings of 

the 11th International Symposium on Shock Tube and Shock Wave Research, pp. 262–268, Seattle, 

WA, July 1977.

Froberg, R.W. “The Carbon-Oxygen Reaction: An Experimental study of the Oxidation of Suspended 

Carbon Spheres,” Ph.D. Dissertation, The Pennsylvania State University, State College, PA, 1967. 

Fry, R.S. “Fundamentals of Test Medium Effects in Ground Testing for Scramjet-Powered Flight,” CPIA 

CPTR 78, January 2004. 

Georgiadis, N.J. “A Hybrid Numerical Method for Turbulent Mixing Layers,” NASA/TM—2001-

210811, May 2001.

Gordon S.; McBride, B.J. “Computer Program for Calculation of Complex Chemical Equilibrium 

Compositions and Applications,” NASA Reference Publication 1311, 1994. 

GrafTech International, Ltd., “Technical Data Sheet 43: Grade PGX Molded Graphite,” 

www.graftech.com, 2005. 

Han, B.; Sung, C.J.; Nishioka, M. “Effects of Vitiated Air on Hydrogen Ignition in a High-Speed Laminar 

Mixing Layer,” Combustion Science and Technology, Vol. 176, No. 3, pp. 305–330, March 2004. 

Heiser, W.H.; Pratt, D.T. Hypersonic Airbreathing Propulsion, AIAA, Washington D.C., 1994. 

Henderson, C.B. “Drag Coefficients of Spheres in Continuum and Rarefied Flows,” AIAA Journal, Vol. 

14, No. 6, pp. 707–708; June 1976.
NASA/TM—2010-216765 90

Hirschfelder, J.O.; Curtiss, C.F.; Bird R.B. Molecular Theory of Gases and Liquids, John Wiley & Sons, 

New York, 1964. 

Hong, J. “Modeling Char Oxidation as a Function of Pressure using an Intrinsic Langmuir Rate 

Equation,” Ph.D. Dissertation, Brigham Young University, Salt Lake City, UT, 2000. 

Hurt, R.; Sun, J-K; Lunden, M. “A Kinetic Model of Carbon Burnout in Pulverized Coal Combustion,” 

Combustion and Flame, Vol. 113, No. 1/2, pp. 181–197, April 1998. 

Hurt, R.H.; Mitchell, R.E. “Unified High Temperature Char Combustion Kinetics for a Suite of Coals of 

Various Rank,” Proceedings of the 24th International Symposium of The Combustion Institute, pp. 

1243–1250, Sydney, Australia, July 1992. 

Igra, O; Takayama, K. “Shock Tube Study of the Drag Coefficient of a Sphere in a Nonstationary Flow,” 

Proceedings of the Royal Society-Series A, Vol. 442, No. 1915, pp. 1243–1250, August 1993.

Incropera, F.P.; DeWitt, D.P. Fundamentals of Heat and Mass Transfer, John Wiley & Sons, New York, 

2002. 

Ishii, R.; Umeda, Y.; Kawasaki, K. “Nozzle Flow of Gas-Particle Mixtures,” Physics of Fluids, Vol. 30, 

No. 3, pp. 752–760; March 1987. 

Kee, R.J.; Rupley, F.M.; Miller, J.A. “Chemkin-II: A Fortran Chemical Kinetics Package for the Analysis 

of Gas Phase Chemical Kinetics,” Sandia Report SAND89–8009B, September 1993. 

Kuo, K.K. Principles of Combustion, John Wiley & Sons, New York, 1986. 

Lee, J.C.; Yetter, R.A.; Dryer, F.L. “Transient Numerical Modeling of Carbon Particle Ignition and 

Oxidation,” Combustion and Flame, Vol. 101, No. 4, pp. 387–398, June 1995.

Libby, P.A. “Ignition, Combustion, and Extinction of Carbon Particles,” Combustion and Flame, Vol. 38, 

No. 3, pp. 285–300, August 1980. 

Ludwig, T.; Roth, P. “Modeling of Laminar Combustion Wave Propagation in Reactive Gas/Particle 

Mixtures,” International Journal of Multiphase Flow, Vol. 23, No. 1, pp. 93–111, February 1997. 

Makino, A.; Law, C.K. “Quasi-Steady and Transient Combustion of a Carbon Particle: Theory and 

Experimental Comparisons,” Proceedings of the 21st Symposium of the Combustion Institute, 

Munich, West Germany pp. 183–191, August 1986. 

Maxwell, B.R.; Seasholtz, R.G. “Velocity Lag of Solid Particles in Oscillating Gases and in Gases 

Passing Through Normal Shock Waves,” NASA Technical Note D–7490, March 1974. 

Mitani, T. “Ignition Problems in Scramjet Testing,” Combustion and Flame, Vol. 101, No. 3, pp. 347–

359, May 1995. 

Moors, J.H.J. “Pulverized Char Combustion and Gasification at High Temperatures and Pressures,” Ph.D. 

Dissertation, Technische Universiteit Eindhoven, Eindhoven, Netherlands, 1998. 

Munro, R.G. “Evaluated Material Properties for a Sintered Alpha-Alumina,” Journal of the American 

Ceramic Society, Vol. 80, No. 8, pp. 1919–1928, August 1997. 

Okada, J.; Ikegawa, T. “Combustion Rate of Artificial Graphites from 700 °C to 2000 °C in Air,” Journal 

of Applied Physics, Vol. 24, No. 9, pp. 1149–1150, September 1953. 

Paxson, D.E. “A Comparison Between Numerically Modeled and Experimentally Measured Loss 

Mechanisms in Wave Rotors,” AIAA Paper 93–2522, June 1993.

Pellett, G.L.; Bruno, C.; Chinitz, W. “Review of Air Vitiation Effects on Scramjet Ignition and 

Flameholding Combustion Processes,” AIAA Paper 2002–3880, July 2002.

Perkins, H.D.; Thomas, S.R.; DeBonis, J.R. “Rocket-Based Combined Cycle Propulsion System Testing,” 

Journal of Propulsion and Power, Vol. 14, No. 6, pp. 1065–1076, November 1998. 

Pierson, H.O. Handbook of Carbon, Graphite, Diamond, and Fullerenes, Noyes Publications, Park Ridge, 

NJ, 1993. 

Pluchino, A.B.; Masturzo, D.E. “Emissivity of Al2O3 Particles in a Rocket Plume,” AIAA Journal, Vol. 

19, No. 9, pp. 1234–1237, September 1981. 

Poinsot, T. J.; Lele, S. K. “Boundary Conditions for Direct Simulations of Compressible Viscous Flows,” 

Journal of Computational Physics, Vol. 101, No. 1, pp. 104–129, July 1992.

Powell, E.S.; Stallings, D.W. “A Review of Test Medium Contamination Effects on Test Article 

Combustion Processes,” AIAA Paper 98–0551, January 1998.
NASA/TM—2010-216765 91

Ranish, J.M.; Walker, P.L “High Pressure Studies of the Carbon-Oxygen Reaction,” Carbon, Vol. 31, No. 

1, pp. 135–141, 1993. 

Roffe, G.; Bakos, R.; Erdos, J. “The Propulsion Test Complex at GASL,” ISABE Paper 97–7096, 

September 1997. 

Saad, M.A. Compressible Fluid Flow, Prentice-Hall, Englewood Cliffs, NJ, 1985. 

Shepherd, J.E. “Chemical Kinetics of Hydrogen-Air-Diluent Detonation,” in Progress in Astronautics 

and Aeronautics, Vol. 106, p. 263–293, AIAA, New York, 1986. 

Sichel, M.; Baek, S.W.; Kauffman, C.W.; Maker, B.; Nicholls, J.A. “The Shock Wave Ignition of Dusts,” 

AIAA Journal, Vol. 23, No. 9, pp. 1374–1380, September 1985. 

Slack, M.; Grillo, A. “Investigation of Hydrogen-Air Ignition Sensitized by Nitric Oxide and by Nitrogen 

Dioxide,” NASA CR–2896, October 1977. 

Smith, G.P.; Golden, D.M.; Frenklach, M.; Moriarty, N.W.; Eiteneer, M.G.; Bowman, C.T.; Hanson, 

R.K.; Song, S.; Gardiner, W.C.; Lissianski, V.V.; Qin, Z. “GRI–Mech 3.0,” 

http://www.me.berkeley.edu/gri_mech/ , March 2006. 

Smith, J.M. Chemical Engineering Kinetics, McGraw-Hill, New York, 1981. 

Svehla, R.A. “Estimated Viscosities and Thermal Conductivities of Gases at High Temperatures,” NASA 

TR R–132, 1962. 

Tannehill, J.C.; Anderson, D.A.; Pletcher, R.H. Computational Fluid Mechanics and Heat Transfer,2nd

Edition, Taylor & Francis, Washington D.C., 1997. 

Tu, C.M.; Davis, H.; Hottel, H.C. “Combustion Rate of Carbon: Combustion of Spheres in Flowing Gas 

Streams,” Industrial and Engineering Chemistry, Vol. 26, No. 7, pp. 749–757, July 1934.

Warnatz, J.; Maas, U.; Dibble, R.W. Combustion, Springer-Verlag, Berlin, 1999.

Weber, J., Editor, “Scramjet Propulsion Testing Recommended Practices and Guidelines,” CPIA 

Publication 710, March 2002. 

White, F.M. Viscous Fluid Flow, McGraw-Hill, New York, 1991. 

White, F.M. Heat and Mass Transfer, Addison-Wesley, Reading, MA, 1988. 

Zhou, H.; Flamant, G.; Gauthier, D.; Lu, J. “Lagrangian Approach for Simulating the Gas-Particle Flow 

Structure in a Circulating Fluidized Bed Riser,” International Journal of Multiphase Flow, Vol. 28, 

No. 11, pp. 1801–1821, November 2002. 
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

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-10-2010 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Development and Demonstration of a Computational Tool for the Analysis of Particle 

Vitiation Effects in Hypersonic Propulsion Test Facilities 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Perkins, Hugh, Douglas 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 599489.02.07.03.03.03.03 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17390 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2010-216765 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 09, 34, and 07 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

In order to improve the understanding of particle vitiation effects in hypersonic propulsion test facilities, a quasi-one dimensional numerical 

tool was developed to efficiently model reacting particle-gas flows over a wide range of conditions. Features of this code include gas-phase 

finite-rate kinetics, a global porous-particle combustion model, mass, momentum and energy interactions between phases, and subsonic and 

supersonic particle drag and heat transfer models. The basic capabilities of this tool were validated against available data or other validated 

codes. To demonstrate the capabilities of the code a series of computations were performed for a model hypersonic propulsion test facility 

and scramjet. Parameters studied were simulated flight Mach number, particle size, particle mass fraction and particle material. 

15. SUBJECT TERMS

Combustion; Particulates; Hypersonic wind tunnels; Supersonic combustion ramjet 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

102 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18