Computation of high-speed reacting flows


Material Information

Computation of high-speed reacting flows
Physical Description:
viii, 185 leaves : ill. ; 29 cm.
Clutter, James Keith, 1966-
Publication Date:


Subjects / Keywords:
Combustion -- Measurement   ( lcsh )
Weapons systems -- Data processing   ( lcsh )
Aerospace Engineering, Mechanics and Engineering Science thesis, Ph. D   ( lcsh )
Dissertations, Academic -- Aerospace Engineering, Mechanics and Engineering Science -- UF   ( lcsh )
bibliography   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1997.
Includes bibliographical references (leaves 177-184).
Statement of Responsibility:
by James Keith Clutter.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 027588928
oclc - 37153585
System ID:

This item is only available as the following downloads:

Full Text







To my wife, Paula


I would like to express my gratitude to those individuals who have helped me

accomplish this goal. Dr. Bruce Simpson, Russ Klugg and everyone at Eglin AFB made it

possible for me to start this effort while on active duty. Also, the continued financial support

from the Armament Directorate of Wright Laboratory, Eglin AFB, has been essential as well

as the computer resources they have put at my disposal. The collaboration and interaction

with Gregg Abate and Rudy Johnson have both helped steer my research and provided key

technical input. The guidance I have received from Dr. Dave Belk has been essential and I

can not thank him enough for his encouragement and for serving on my committee.

I would like to thank Dr. Hartmuth Lehr for providing his experimental data and for

his inputs on the topic of shock-induced combustion. Also, Dr. Greg Wilson has generously

provided both reference material and insights that are greatly appreciated. Dr. Mark Janus

has been a good friend and provided both encouragement and education into the finer arts

of motorsports and statistics.

Several individuals at the University of Florida have helped me see this effort

through. First and foremost is Dr. Wei Shyy who has been both a mentor and friend.

Completion of this effort would not have been possible without his guidance,

encouragement, and tolerance. I thank Dr. Ed Milton for helping me begin this endeavor and

Dr. Renwei Mei, Dr. Corin Segal, and Dr. James Klausner for serving on my committee as

well as Dr. David Mikolaitis for his inputs and suggestions in the area of kinetics. I would

also like to thank my colleague, Dr. Venkata Krishnamurty, who has provided valuable

insight into the modeling of turbulence.

My family has provided endless support whether it be emotional or dropping in to

help with the kids during times of peak studying and work. I can not thank my mom and my

mother-in-law enough as well as my father-in-law, Richard Cooper, who recently passed

away. I must thank my father, the eternal student and consummate mechanic, for instilling

in me a desire to learn how things work and the moxie to take them apart. My children,

Courtney and Cooper, are without a doubt the finest children in the world and have been

immensely understanding and patient. They have also proven to be delightful distractions

and have helped keep me focused on what is truly important.

Finally, there is my wife, Paula, whose strength and encouragement are

immeasurable. She has never ceased to support me whether it meant her understanding in

me quitting my job to return to school or balancing her career and graduate school along with

mine. I will never be able to express all she means to me.


ACKNOWLEDGEMENT ............................................ iii

ABSTRACT ....................................................... vii

1 INTRODUCTION ............................... ............ 1
1.1 Reactive M unition ................................. ......... 1
1.2 Gun Muzzle Blast .......................................... 6
1.3 Phenomenological Scales .................................... 11
1.4 Summary ....................................... ........ 14
2 SHOCK-INDUCED COMBUSTION ............................. 17
2.1 Reaction Mechanisms for Hydrogen and Oxygen .................. 21
2.2 Previous Computational Studies of Lehr's Experiments ............. 22
2.3 Summ ary ................... .................... ........ 33
3.1 Governing Equations ....................... ............ .... 37
3.2 Equation of State and Temperature ............................ 41
3.3 Thermal Conductivity and Viscosity ................... ....... 43
3.4 Flux Vector Splitting Scheme ........................ ....... 44
3.5 Chemical Reaction Source Term ............................... 47
3.6 Special Treatment For Chemical Length Scales ................... 49
3.6.1 Transformed Governing Equations ........................ 50
3.6.2 Logarithmic Interpolation ............................... 51
3.6.3 Error Estimation Based on One-Dimensional Model Equations .. 53 Global Order of Accuracy Analysis ................. 57 Current Application of Scaling Factor Scheme ......... 61
3.6.4 Nonphysical Wave Speeds .............................. 63
3.7 Numerical Solution Methods ................................ 64
3.7.1 Fully-Implicit Formulation ................ .............. 64
3.7.2 Explicit-Implicit Splitting Scheme ........................ 67
3.7.3 ODE Solution Issues ................................... 69
3.8 Summary .............................................. 71
4 TURBULENCE MODEL ....................................... 72
4.1 Averaging Techniques and Averaged Variables .................... 75

4.2 Averaged Equations ........................................ 77
4.3 Turbulence Closure ................................ ..... .. 79
4.4 k-e Model ................................................ 81
4.5 Wall Function ............................................. 83
4.6 Turbulence Effects On The Reaction Source Term ................. 85
4.7 Summary ................................................ 86

5.1 Comparison of Solution Methods .................... ........ 94
5.2 Source Term Treatment ...................................... 96
5.3 Effect of Reaction Mechanism on the Computed Flow Field ......... 104
5.3.1 Pressure Effects on Reaction Rates ........................ 123
5.3.2 Summary of Reaction Mechanism Effects .................. 128
5.4 Computations Using Turbulence Model ......................... 129
5.5 Summary ...................................... ......... 136

6.1 Gun Muzzle Blast ........................................... 138
6.2 Gun Muzzle Flash ......................................... 142
6.3 Summary ................................................. 143

7.1 Large Caliber Howitzer Simulations ............................ 145
7.1.1 Flow Field Model and Boundary Conditions ................ 146
7.1.2 Experimental Overview ................................ 150
7.1.3 Computational Results ................................. 151
7.2 Small Caliber Sound Suppressor Simulations .................... .160
7.3 Summary ................................................ 165

8 CONCLUSIONS AND FUTURE WORK ......................... 168
8.1 Shock-Induced Combustion .................................. 168
8.2 Gun M uzzle Blast ................................ ...... 171
8.3 Summary ................................................ 172


REFERENCES ................................................... 177

BIOGRAPHICAL SKETCH ........ ............................. 185

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy




May 1997

Chairperson: Dr. Wei Shyy
Major Department: Aerospace Engineering, Mechanics and Engineering Science

A computational study has been conducted for high-speed reacting flows relevant to

munition problems, including shock-induced combustion and gun muzzle blast. The

theoretical model considers inviscid and viscous flows, multi-species, finite rate chemical

reaction schemes, and turbulence. Both the physical and numerical aspects are investigated

to determine their impact on simulation accuracy.

A range of hydrogen and oxygen reaction mechanisms are evaluated for the

shock-induced combustion flow scenario. Characteristics of the mechanisms such as the

induction time, heat release rate, and second explosion limit are found to impact the accuracy

of the computation. On the numerical side, reaction source term treatments, including

logarithmic weighting and scaling modifications, are investigated to determine their

effectiveness in addressing numerical errors caused by disparate length scales between

chemical reactions and fluid dynamics. It is demonstrated that these techniques can enhance

solution accuracy.

Computations of shock-induced combustion have also been performed using a k-E

model to account for the turbulent transport of species and heat. An algebraic model of the

temperature fluctuations has been used to estimate the impact of the turbulent effect on the

chemical reaction source terms. The turbulence effects when represented with the current

models are found to be minimal in the shock-induced combustion flow investigated in the

present work.

For the gun system simulations, computations for both a large caliber howitzer and

small caliber firearms are carried out. A reduced kinetic scheme and an algebraic turbulence

model are employed. The present approach, which accounts for the chemical reaction

aspects of the gun muzzle blast problem, is found to improve the prediction of peak

overpressures and can capture the effects produced by small caliber firearm sound


The present study has established the numerical and physical requirements for

simulating high-speed reacting flows. Also, key parameters useful in quantifying the roles

of these aspects have been assessed.


In the present work, a computational capability is developed to analyze high speed

reacting flows arising in munition problems. The particular problems focused on are a

munition concept and the occurrence of gun muzzle blast. These problems are difficult to

analyze given the fact many physical and chemical processes are involved. First a description

of each problem is presented to highlight these various processes and frame the motivation

of the current research. During these descriptions the phenomenological issues such as the

fluid dynamics, chemical reactions, and turbulence are highlighted. This is followed by a

description of how the time and length scales of these various aspects impact computational

models. Finally, the specific scope of the present thesis is defined and the approach adopted

to solve these problems is presented.

1.1 Reactive Munition

The first munition problem of interest is a weapon concept dependent on the mixture

of fuel and oxidizer and shock initiation and is similar to the concepts discussed by Gore et al.

(1993). This class of weapons differs from conventional munitions in the manner by which

they derive their destructive force. This destructive force has to be generated by some energy

source and in typical munitions, such as general purpose bombs, the source is a condensed

phase explosive. This explosive contains the chemical compounds needed to generate the

force once initiated by a fuse. Alternate munition concepts are based on the idea of

generating the destructive force from alternative means. This may be done to enhance the

energy source beyond that which is carried to the target but another focus is to use materials

on target as energy sources. This is highly attractive for munitions to be launched from both

ships and aircraft given the fact that the less material that must be carried to the target, the

lighter and smaller the munitions system must be. This directly impacts logistic concerns

such as storage as well as factors such as the fuel required to deliver the munitions.

This current study is part of an effort to study the various physical and chemical

phenomena which directly impact the effectiveness of new munition concepts. Another facet

of the effort is to quantify the modeling tools needed to simulate this class of weapons and

to aid in the development of engineering codes useful in design and evaluation of these

various concepts. The range of concepts and targets which fall in this category is quite

extensive and a complete review is beyond the scope of this document. Therefore, one target

and one munition concept are highlighted here. However, it should be noted that the

underlying physical and chemical processes to be presented in the current context are

common to all targets and munition concepts.

As mentioned, one of the premises of various new munitions is to take advantage of

materials on target which can be compromised and used to defeat the target. One component

common to all air vehicles whether missiles or aircraft is on-board fuel. Given the particular

vehicle, the fuel can be stored in various states. For brevity, here we consider an air vehicle

carrying fuel in a fluid state. For example, most aircraft contain numerous fuel cells in the

wings and fuselages. Depending on the stage in the mission, the fuel cells contain some level

of liquid fuel as well as fuel vapor. Therefore, the munition concept used to help define the

important physical and chemical phenomena to be studied centers around using the on-board

fuel in a combustion process as the source of energy to generate the destructive force.

Even though there can be an ample supply of fuel on target and an ample supply of

oxidizer, the ambient air, if the required conditions are not achieved, the munition-target

system can remain in a fuel or oxidizer rich state and no effective combustion is initiated. One

plausible system concept to attack the fuel cells is a weapon that generates a "cloud" of

objects designed to penetrate and destroy the fuel cells. The focus of the current project is to

study various key phenomena to aid in the development of an engineering level model with

which to quantify the destruction output once this "cloud" intercepts the fuel cell.

Continuing with the concept of attacking the fuel cell, the target can be thought of as a

control volume bounded by a solid wall with specific characteristics such as material and

thickness. In the context of the larger issue of weapon effectiveness, criteria which need to be

quantified and are dependent on output from any model of the munition-target interaction are

factors such as whether the walls are damaged such that this damage is rapidly transferred

throughout the entire air vehicle. Two tragic examples of such a process are the recent Valujet

and TWA Flight 800 accidents. In the first case, the predominant theory is that the damage

was initiated by the failure of one or several oxygen generators (McKenna, 1996). In the

TWA case, as of this writing, it is still unknown whether the damage to the aircraft had been

initiated by an explosive device or a system malfunction. Regardless of the initiation process,

damage in both cases had been tremendous due to the fact that conditions were such that the

destruction process was sustainable. In any combustion process, which both these cases

entailed, sustainment is dependent on the correct supply of fuel and oxidizer as well as the

required system settings. The first goal of this effort is to better understand the relationship

between such system settings and the physical and chemical phenomena which conspire

together to generate global effects such as the destructive force output. The second goal of the

present work is to quantify the attributes of the computational model required to simulate

these phenomena.

Our fuel cell scenario is depicted in Figures 1.1 and 1.2 and the steps leading to

destruction are as follows. The fuel cell contains some amount of liquid fuel (F1) as well as

fuel vapor (Fv). At stage 1, the projectile intercepts the fuel cell with a shock front (Sf) due to

the projectile flight through the ambient region which is composed of air, an oxidizer (0).

After the projectile enters the fuel cell, the structure of the established shock (Sf) changes

somewhat due to the properties of the fuel. There may also be a precursor compression (Sp)

wave initiated by the impact of the projectile with the fuel cell wall. The strengths of the

established shock in stage 2 and precursor compression wave depend on the characteristics of

the wall. If the wall is relatively thin, then the momentum transfer from the projectile to the

wall is negligible before the projectile penetrates. On the other hand, for a stronger wall, the

momentum transfer will be larger and the velocity of the projectile will decrease

substantially during penetration. This would result in a stronger precursor compression wave

and dramatically weaken the shock already established by the projectile.

Fuel Cell
Projectile Cloud

Figure 1.1. Fuel cell scenario showing projectile cloud. Single projectile
highlighted in Figure 1.2.

After penetration, as the projectile and shock proceed through the cell, the heat from

the shock can vaporizes any liquid fuel. This fuel vapor mixes with the oxidizer (0) which

has penetrated into the fuel cell and the mixing is enhanced by the wake of the projectile. The

projectile may be constructed out of an ablative material formulated to enhance the reaction

process. In that case, the projectile changes phase from a solid (Ms) to eventually a gas (Mg).

This compound reacts with the fuel vapor (Fv) and oxidizer (0) to produce a final product

(P). The final step leading to destruction is for the reaction to reach a state such that it will

progress through the fuel cell and combine with the same process initiated by the other


The overall effectiveness of this or any weapon concept is a function of the target and

weapon system parameters. The target parameters here include primarily the fuel chemical

composition, phase, and state which directly impact key aspects such as how much total

energy is needed to compromise the fuel volume. These aspects are controlled by a

combination of detailed phenomena such as chemical time scales. The second of the target

I 0i I




(1) Projectile intercepts fuel cell.
Shock (Sf) due to flight through
ambient air.
Fuel cell contains some level of
liquid fuel (FI) and fuel vapor (Fv)

Fv O


L -- **- jI

(2) Precursor (Sp) generated during
Heat release behind shock (Sf)
produces fuel vapor (Fv)
Oxidizer (0) penetrates cell


(3) Solid projectile (Ms) ablates (4) Reaction front (Rf) proceeds
producing gas (Mg) through cell
Fuel vapor (Fv), catalysis (Mg)
and oxidizer (0) mix and react
producing product (P)

Figure 1.2. Sequence of events surrounding an individual projectile.

parameters is the container characteristics such as wall material and thickness. In this

example, these parameters enter into the analysis process somewhere between the modeling

of the munition-target interaction assessment and a larger lethality model. Essentially, some

sort of effective destruction power is generated and the effectiveness of this power in

damaging the fuel cell wall and surrounding structures has to be quantified and this damage

output is fed into a higher level lethality model.

The second class of parameters is those associated with the munition. The particular

munition concept used in the current discussion is a scenario where the fuel cell is infiltrated

by a certain number of projectiles. Therefore, the munition parameters are first the number of

projectiles and their dispersion pattern. Next are the parameters of the individual projectiles

such as their size, shape, material composition and velocity. The importance of these

parameters on the system effectiveness is manifested in how they drive the important

detailed physical and chemical phenomena. From an application context, the quantities of

interest are the global effects such as the destructive force generated during the

munition-target interaction. In the fuel cell scenario one key global effect is whether

combustion is achieved and sustained. The correlation between system parameters and these

global effects is controlled by the detailed physical and chemical phenomena.

It is clear from the description of the fuel cell problem that among the key physical

phenomena is the mixing of the fuel and oxidizer. Given the flow scenario, turbulence plays

a key role in the mixing as well as the transport of heat. Among the processes of interest are

phase change and the reactions between the fuel and oxidizer. The modeling of these various

phenomena is the focus of the present work.

1.2 Gun Muzzle Blast

The second munition problem used in the current study is the simulation of gun

muzzle blast. This problem contains many of the same physical and chemical processes

discussed in the description of the reactive munition problem and requires much of the same

modeling tools. However, here computational efficiency is paramount given the desire to

develop an engineering tool useful in the design of muzzle devices. The primary concern in

this problem is the reduction of the pressure levels generated when a conventional gun is

fired. In the case of large caliber systems, the reduction is sought to alleviate loads the support

vehicle must sustain. In small caliber guns, the concern is the associated hearing loss from

repeated use and the ability to go undetected for certain applications.

This problem has been studied by several investigators. The primary tool of

investigation has been experiments with a limited computational effort. A literature review

is provided in a later chapter. Here the problem is described to highlight the need for further

research. The ballistic process determining blast is depicted in Figure 1.3 and can be divided

into two domains. The first is the interior ballistic cycle where the solid propellant in the

rounds is ignited and burned. This results in a build-up of chamber pressure which

accelerates the projectile down the barrel. Following the projectile is a volume of gun gas

which is created from the spent propellant. This gas is characterized by both high pressures

and temperatures as well as some level of particle loading. The point in the barrel where the

particle loading is not changing and the solid propellant has been combusted to essentially

an equilibrium state can be considered as the boundary between the interior ballistic cycle

and the launch cycle (point A in Figure 1.3). This is the point where the current code imposes

specific boundary conditions.

The launch cycle encompasses various physical occurrences such as the projectile

exiting the barrel and the progression of the precursor and main blast wave from the barrel

into the ambient environment. The precursor blast wave denotes the shock created by the

fly-out of the projectile while the main blast wave denotes the high pressure gun gas which

exits the muzzle with the projectile. The level of blast experienced is also directly related to

parameters such as the particular propellant used in the rounds.

Another physical phenomenon which often occurs when a conventional gun is fired

is muzzle flash. The occurrence of muzzle flash is a result of further combustion of the gun

gas with the ambient air. This chemical reaction process not only generates flash but also

Solid Propellant Projectile

Large Grains S
I..... V

(1) (3)

mall Particles Gun Gas / Air Interface

Precursor Blast Wave

Main Blast Wave

Precursor Blast Wave

(2) (4) 1

Figure 1.3. Interior ballistic and launch cycle of a typical conventional gun.

alters the overpressure signature. Even though the focus here is on the overpressure, it is

worthwhile describing the mechanics of the occurrence of flash to obtain a full appreciation

of the problem.

The gun gas which is discharged during the launch cycle is always high in

temperature, always contains some level of solid particle loading, and is typically rich in

hydrocarbons which can react with the oxygen in the ambient air. Flash associated with

conventional guns can be classified in three main regions; primary, intermediate, and

secondary flash. These regions not only occur at different stages in time after the gun firing,

but also in different spatial regions of the field downstream of the gun muzzle. These regions

are represented in Figure 1.4. The primary flash is located at the exit of the muzzle and

encompasses a small spatial region and is low in intensity. However, the intermediate flash

is a more extensive region of greater intensity and occurs just downstream of the shock disk

which forms some distance from the muzzle exit. The location of the intermediate flash

region is a direct result of the fact that it is the heating of the gun gas by the shock disk which

initiates the combustion producing the flash. Another influence of the shock disk is the

1 AZ


reduction in the velocity of the gas which increases the residence time of the gas and

promotes the initiation of the chemical reactions.

In addition to the heating from the shock disk and the reduction in the gun gas

velocity as it crosses the shock, a certain level of oxygen must be mixed with the gun gas

before the flash occurs. This is why there is typically a delay between the firing of the gun

and the occurrence of the intermediate flash. This delay is directly dependent on the time

it takes to entrain oxygen from the ambient air and to mix it through the volume of gun gas

toward the line of fire. The entrainment process is highly dependent on the turbulent

processes which operate along the interface between the gun gas and ambient air. Therefore,

any alteration to the gun system such as the addition of a muzzle device will impact the

mixing process and in turn the intermediate flash.

The level of flash generated from the intermediate region is typically repeatable from

firing to firing for a given gun system, i.e. gun caliber, ammunition, and muzzle device. The

caliber of the gun and the ammunition fired determine the volume of the gas discharged as

well as the gas properties when it exits the muzzle. These gas properties set the characteristics

of the shock disk as well as the conditions upstream of the disk such as the fuel richness.

Along with the influence of any muzzle device mentioned earlier, these parameters

determine the level of flash in the intermediate region.

Secondary flash is the final and typically the largest and most intense region of the

three. Like the intermediate flash, the secondary flash is initiated by the mixing of the oxygen

rich air and gun gas and also requires a source of energy to initiate the reaction. And as before,

the residence time of the mixture must be such to allow combustion to initiate and progress.

In the case of the intermediate flash, the shock disk provides the energy source as well as

reduced the velocity of the gas to ensure a high enough residence time. In the case of

secondary flash, it has been found (Klingenberg and Heimerl, 1992) that if the temperature

in the intermediate flash region does not reach a certain level, dependent on the ammunition,

then secondary flash does not occur. This supports the theory that the heat released in the




Intermediate Flash

Primary Flash

Interface of Gun
and Air

Main Blast Wave

Figure 1.4. Depiction of flash regions and the development of muzzle blast and flash.



intermediate flash is the energy source for initiating the secondary flash which occurs some

time after the intermediate flash and further downstream. The residence time criterion is met

by the ever increasing volume of gun gas. As the cloud of gas grows larger, the velocity

throughout the cloud is reduced. Therefore, as depicted in Figure 1.4, the secondary flash

covers a larger area of the field. Secondary flash from shot to shot is not as repeatable as

intermediate flash. This is due to the difference in the initiation process. In the case of

intermediate flash, the initiation is driven by the shock disk which is less influenced by small

perturbations in both the ambient field and the particular round fired. The initiation of

secondary flash is much more dependent on the transport of both species and heat in a much

larger region, further from the gun muzzle. This makes it more accessible to ambient

influences such as winds.

To simulate all the intricacies of flash requires a detailed model of the chemical

composition of the gun gas. Some modeling efforts have been carried out and these are

reviewed later. Since the focus of the current study is to predict the overpressure in regions

somewhat removed from the gun system, the key phenomenon associated with the reactions

which cause flash that must be reproduced is the energy released from the combustion. This

is due to the fact that the combustion occurs behind the main blast wave. This deposition of

energy behind this wave helps sustain the magnitude of the main blast wave as it progresses

out into the field. This is key since the peak pressure at any location in the field is a function

of the main blast wave magnitude. Therefore, the current model employs a reduced reaction

mechanism in an attempt to better reproduce the blast levels without greatly increasing the

computational requirements. However, the characteristics of the flash processes are to prove

advantageous in evaluating the model.

1.3 Phenomenological Scales

The description of the problems above makes it clear that there are a multitude of

physical and chemical processes key in high-speed, chemically reacting flows. Associated

with this myriad of mechanisms are multiple time and length scales. The time and length

scales refer to the characteristic time and length over which effects due to certain physical

and chemical mechanisms occur. Consequently, it is difficult to investigate this class of

problems both theoretically and experimentally. Therefore, computational models are

routinely used to study reacting flows and this is the methodology of the current work.

Any models used to simulate the various detailed phenomena involved must

correctly reproduce all phenomena and their corresponding scales, or at the least the role the

phenomena play in the larger problem. As for the numerical aspects of such models, the

requirement to resolve a wide range of scales determines simulation parameters such as the

computational time and computer memory needed to carry out a useful simulation. This is

due to the fact that when solving a problem using computational models, the governing

equations are usually solved at discrete points in time and space. Therefore, the temporal and

spatial discretization must be done with both the time and length scales of interest in mind.

A key parameter of chemically reacting flows is the flow speed or the fluid dynamic

time scale relative to the chemical processes. A nondimensional parameter useful in

classifying problems is the Damk6hler number defined as the ratio between the fluid

dynamic and chemical time scales, ) = Tf/Tc. For low-speed flows, the chemical processes

typically occur orders of magnitude faster than the fluid dynamic processes. This implies

) > 1 and adds validity to modeling assumptions such as local chemical equilibrium or thin

flame zones. From a modeling standpoint, such assumptions allow the use of a simplified

chemical reaction model. For instance, as a result of the equilibrium assumption, the finite

rate chemical processes can be ignored and the presence of intermediate species may be

neglected. Even if the flow is not in chemical equilibrium but the reactions can be considered

fast, either a partial equilibrium (Correa and Shyy, 1987) or a progress variable (Libby and

Williams, 1994) concept can be adopted. For the case of high-speed flows as in the case of

shock-induced combustion, such assumptions are generally not valid and a model of the

finite rate chemical reactions accounting for the detailed distribution of intermediate species

is needed. This has direct implications as to the complexity of the computational model.

Many combustion problems have to address the issue of turbulence. This is due to

the fact that the mixing characteristics of turbulence are routinely present with high Reynolds

number flows. Also, in many cases turbulence is used in reactive systems to enhance the

mixing and subsequently the combustion process. As in all turbulent flows, the aspects due

to turbulence must first be compared to the mean flow and with the addition of reactions,

the relative scales of the turbulence to the chemical processes are important. Turbulence is

known to possess a multitude of scales, particularly the large scale eddies which contain what

is labeled the turbulent kinetic energy (k) and the small scale eddies which work to dissipate

this energy. The dynamics of turbulence consist of various processes by which the large

eddies break down into small eddies and thereby transfer energy between the large and small


Depending on the relationship between the large and small scales of turbulence and

those of the chemical reaction some conclusions can be drawn as to the impact of turbulence

on the combustion process. A parameter useful in the analysis process is the ratio of the

large-eddy time scale to the chemical time scale, the large-eddy Damk6hler number

(I) = tj/Tc). The large-eddy time scale, t, = I/ k, is defined using the characteristic

length scale (1) and velocity (-k) of the eddy. Using this parameter, the regime of interest

runs from slow chemistry, small dt, to fast chemistry, large D/. For the case of slow

chemistry, the turbulent processes are working on a scale smaller than the chemical process.

This produces enhanced mixing of reactants or even intermediate species resulting in the

reaction zone being distributed in the flow field. In this case, the impact of turbulence on

mass and thermal diffusion must be considered. For the case of fast chemistry the reactions

take place on a scale much smaller than the large eddies. This results in distinct thin flame

fronts which lie within the large eddies and are transported and strained by the turbulence.

Such cases can be modeled using the laminar flamelet concept which is reviewed by Libby

and Williams (1994). For the high-speed flows to be studied here, it is probable that the

turbulent and chemical time scales are comparable. Therefore, one goal of the current study

is to incorporate the effects of turbulence on the chemical processes in the computation.

1.4 Summary

Here the problems to be studied have been described to identify the important

physical and chemical processes which must be modeled. Given the high-speed nature of the

flows, many assumptions used to reduce the complexity of models employed in low-speed

combustion are not applicable. Some of the length and time scales have been identified and

how these impact the numerical aspects of the model is discussed later. The first primary

focus of the current effort is to quantify the complexity of the models needed to accurately

simulate the munition problems of interest. The second focus is to identify the impact both

the phenomenological and numerical aspects of these models have on the accuracy of these

simulations. The remainder of this document describes the specific models and

computations to be carried out.

Chapter 2 describes in detail a shock-induced combustion scenario to be used as a

representation of the munition problem consisting of the fuel cell scenario. This is done in

part due to the availability of experimental data associated with the particular scenario. Also,

this scenario has been studied by previous investigators using a wide variety of physical,

chemical, and numerical models. A literature review of these earlier works is given to help

draw clear comparisons between these efforts and the current study. Using the code

developed here, many of these various models can be tested within a single computational

framework to ascertain the correlation between these various models and simulation

accuracy. This not only allows judgments to be made as to the validity of the various models

but also aids in quantifying the fidelity needed in a computational tool for the munition

concepts of interest.

The governing equations of the current class of munition problems are presented in

Chapter 3 along with a detailed description of the computational code developed here. It is

to become evident that there exist numerous modeling options for the various facets of the

code from both a phenomenological and numerical point of view. Where appropriate these

options are discussed. The assumptions which drive the model selection are detailed as well

as compelling reasons for the current model selection based on both the physical and

chemical aspects of the problems to be solved and the computational requirements the

solution process must meet. Also discussed in Chapter 3 are the various numerical

techniques used in the present computational code.

Chapter 4 describes the turbulence model used in the current effort as well as the

various terms which are impacted by the turbulence aspects of the flow scenarios. The model

selection for the representation of turbulence is driven in a large part by the computational

requirements of the code developed here. Various model modifications for addressing

attributes of the high-speed flows to be solved are also discussed.

Results for the shock-induced combustion scenario are presented in Chapter 5.

Among the data reviewed is analysis as to how the chemical model and numerical aspects

of the code impact accuracy. Computations using both inviscid and turbulence fluid dynamic

models are presented to both test the current code but also the ascertain to what extent the

turbulence aspects of the flow impact the combustion scenario. The current results are

reviewed in the context of previous findings both to support the conclusion of the present

effort and to fill in those aspects of solving this class of flows not directly addressed in earlier


A literature review of pertinent work in the area of gun muzzle blast is provided in

Chapter 6. Many of the conclusions and findings of these earlier efforts have helped steer

the development of the current code used to simulate this munition problem. It is to be made

evident that the models to date have ignored some of the physical and chemical processes

and it is the focus of the present effort to quantify the impact these processes have on the

generation of muzzle blast. This is followed in Chapter 7 by a description of how the model

discussed in Chapter 3 is used to simulate the muzzle blast scenario. Chapter 7 also presents

results using a computational model analogous to those used in earlier efforts and compares

these simulations to results from the current model.

Finally, Chapter 8 summarizes the findings of the present study into both the

shock-induced combustion and gun muzzle blast problems. The accuracy of the code

developed here is discussed as well as the computational requirements. Areas into which

more focus is needed are identified in the context of the munition problems of interest. Also,

in a more general sense, those issues in need of more attention to increase the capability of

simulating high-speed reacting flows for full scale engineering problems are identified and

some suggestions for the future direction are made.


The fuel tank scenario described earlier contains a multitude of physical and

chemical processes. The predominate processes are the shock dynamics, mixing, phase

change, and chemical reactions between the fuel and oxidizer. For the current look into this

class of problems, a premixed combustion scenario is studied. The shock dynamics are

explicitly present and the mixing aspects are manifested in the importance of the distribution

of intermediate species and heat. The chemical processes are represented by the chemical

reactions between the fuel and oxidizer. Aspects also important to the phase change process

are present. These include the competition between the fluid dynamics and chemical time

scales and the dependence of the process on local field properties.

The combustion scenario to be used in the present research for studying the reactive

munition problem is the shock-induced combustion scenario studied experimentally by Lehr

(1972). In his experiments, Lehr shot projectiles into stoichiometric hydrogen and air

mixtures as well as hydrogen and oxygen mixtures. Projectiles of various geometries have

been used but here only data from shots using hemispherically nosed projectiles are studied.

Lehr's primary findings are that the flow field features generated by the shock-induced

combustion can either be one of steady combustion or a situation in which unsteady

combustion occurs. The case of unsteady combustion can be further subdivided into low and

high frequency regimes. The regime achieved is found to be directly related to the ratio of

the projectile velocity to the detonation velocity of the mixture. The data available from

Lehr's study are qualitative in nature and are in the form of shadowgraphs. A representative

image from the three regimes has been provided by Lehr (personal communication, 1995)

and is shown in Figure 2.1. The shock in all three images is visible as is a second front due

to the chemical reactions and labeled the energy release, or reaction front. All three of these

cases are for a stoichiometric mixture of hydrogen and air at a temperature of 293 K and a

pressure of 320 Torr.

As is clear by the images in the figure, for those cases where the velocity of the

projectile is lower than that of the detonation velocity of the mixture (VD), an instability is

seen in the reaction front. As the velocity approaches the detonation velocity, the frequency

of the instability increases until the combustion becomes stable at velocities greater than VD.

It should be noted that at the other end of the velocity spectrum, for low Mach numbers, no

Vp = 1685 m/s
Vp/VD = .82
M= 4.18
freq = 150 kHz

Vp = 1931 m/s
Vp/VD = .94
M = 4.79
freq = 720 kHz

Vp = 2605 m/s
VPVD = 1.27
M = 6.46
freq = steady

Figure 2.1. Example shadowgraphs from Lehr's Experiments.

combustion is initiated. Several investigators have proposed theories as to the important

phenomena driving the instability and have supported their theories with computations.

Among these is the work of McVey and Toong (1971), Alpert and Toong (1971), Matsuo et

al. (1993), and Matsuo and Fujiwara (1993). There is also some question among the

investigators as to when the actual transition from unsteady to steady occurs. Sussman

(1994) presents arguments that the transition does not occur exactly at VpVD = 1.

Many investigators have simulated the cases presented in Figure 2.1 and a review is

provided in the next section. However, here some general comments are needed on the

experimental data themselves as well as any computational effort to simulate these

shock-induced combustion cases. The shadowgraph from Lehr for M=6.46 is shown in

Figure 2.2 with points used to denote the shock and reaction fronts. It has been verified

(Lehr, personal communication, 1995) that hemispherically nosed projectiles with

cylindrical bodies have been used in the experiments. However, in the image the nose of the

projectile does not appear to be a true hemisphere. One possible source for this alteration

is that the projectile is simply yawed. Such a yawing would also be evident in a

corresponding change to the projectile's baseline seen in the images and this is not the case.

There is a slight alteration to the baseline in the images but not nearly on the order of the

alteration of the nose. After consultation with Lehr it has been verified that the alteration is

due to the optical system (Lehr and Maurer, 1969) and he cautions against using the data for

precise measurements. However, the data can be used qualitatively to benchmark the

computational models to see if the flow field features are being captured. And given the fact

the current study is focused on the interplay of the fluid dynamics and chemical processes,

Lehr's experiments offer a solid set of data in which the scales of each of these aspects are


As is evident by now, the successful simulation of a flow scenario such as Lehr's

experiments (Lehr, 1972) or the reactive munition problem requires a multitude of models

for the various phenomenological aspects. Each model not only embodies the particular

Figure 2.2. Shadowgraph for M=6.46 with the shock and reaction front denoted.

phenomenology framework used to represent either the physics or chemistry but also their

numerical traits. For instance, the reaction process between hydrogen and oxygen in reality

occurs through countless intermediate processes. However, to represent this phenomenon

with a numerical model, a distinct number of reaction steps must be selected and reaction

rates for these steps must be specified. The selection of a particular model, i.e. the number

of steps and rates, inherently specifies the chemical scales produced by the model. These

scales may or may not match those which exist in reality. If a discrepancy exists between the

computation and experimental data, the model can be altered until a good match is made.

This can increase the confidence in the ability of the model to reproduce key

phenomenological characteristics such as the chemical time scales in the case of the reaction

model. However, such alteration should be done with caution. For example, many

computations have been carried out by previous investigators for the shock-induced

combustion scenario captured by Lehr (1972). In these earlier studies, several reaction

mechanisms have been used as well as a variety of grids. In most cases, an individual study

uses a single reaction mechanism and the grid is adapted or other numerical aspects are

altered until good agreement with experimental results is achieved. The accuracy of the

computation is dependent not on a single model but the combination of all models. It is to

become clear that for this class of problems, the phenomenological and numerical aspects

must be considered in conjunction.

2.1 Reaction Mechanisms for Hydrogen and Oxygen

One of the predominate differences among the earlier computational studies into the

shock-induced combustion scenario is the representation of the chemical reactions.

Therefore, here a description of various reaction mechanisms used for hydrogen and oxygen

is provided. As detailed earlier, the chemical processes in the problems of interest have time

scales such that a finite rate chemistry model is needed. The chemical reaction processes

involved are incorporated into the governing equations by way of the source terms to be

defined later.

There have been many reaction mechanisms proposed for the reaction between

hydrogen and air. Here, just those used in the previous computational studies to be reviewed

and those used in the current study are presented. They range in complexity from a 2-step

global model proposed by Rogers and Chinitz (1983) to a 32-step model proposed by

Jachimowski which is based on his original 33-step model (Jachimowski, 1988) but

incorporates a slight modification he and others suggest (Oldenborg et al., 1988). These

mechanisms are presented in Tables 2.1 through 2.7. The constants are used to calculate the

rate of the individual reaction using the Arrhenius expression
P-- Ean
kn = AnTb exp[ E (2.1)

which is discussed in more detail later. In the review, a 19-step mechanism is mentioned. This

mechanism is composed of the first 19 steps of the Jachimowski. Note, the two 8-step

mechanisms presented use the same steps but different rate constants.

Table 2.1. Two-step hydrogen-oxygen reaction mechanism from Rogers and Chintz

n Reaction An bn Ean / Ru

1 H2 + 02 = 20H (8.917) + 31.433/(-28.950) x 1047 -10 2449.77
2 H2 + 20H = 2H20 (2.0 + 1.333/(P-0.833() x 1064 -13 21400.87

( = equivalence ration, set to 1
in current computations

Units are moles, seconds, centimeters, and Kelvin

Table 2.2. Seven-step hydrogen-oxygen reaction mechanism used by Shang et al.
(1995) based on eighteen-step mechanisms from Drummond et. al. (1986).

n Reaction An bn Ean / Ru

1 H2 + 02 = 20H 1.7 x 1013 0 24233.
2 H + 02 = OH + O 1.42 x 1014 0 8254.
3 H2 + OH = H + H20 3.16 x 107 1.8 1525.
4 O + H2 = OH + H 2.07 x 1014 0 6920.
5 20H = H20 + O 5.5 x 1013 0 3523.
6 H + OH + M = H20 + M 2.21 x 1022 -2 0.
7 2H + M = H2 + M 6.53 x 1017 -1 0.

Units are moles, seconds, centimeters, and Kelvin

2.2 Previous Computational Studies of Lehr's Experiments

Both steady and unsteady combustion cases captured in shadowgraphs by Lehr

(1972) have been studied by several investigators using numerical techniques. The actual

Table 2.3. Eight-step hydrogen-oxygen reaction mechanism from Moretti (1965).

n Reaction An bn Ean / Ru

1 H + 02 = OH + O 3.0 x 1014 0 8810.
2 O + H2 = OH + H 3.0 x 1014 0 4030.
3 H2 + OH = H + H20 3.0 x 1014 0 3020.
4 20H = O + H20 3.0 x 1014 0 3020.
5 H2 + M = 2H + M 1.85 x 1017 -1 54000.
6 H20 + M = OH + H + M 9.66 x 1018 -1 62200.
7 OH + M = O + H + M 8.0 x 1016 -1 52200.
8 02 + M = 20 5.8 x 1016 -1 60600.

Units are moles, seconds, centimeters, and Kelvin

Table 2.4. Eight-step hydrogen-oxygen reaction mechanism from Moretti (1965)
with rate constants proposed by Evans and Schexnayder (1980).

n Reaction An bn Ean / Ru

1 H + 02 = OH + O 1.2 x 1017 -.91 8315.
2 O + H2 = OH + H 7.5 x 1013 0 5586.
3 H2 + OH = H + H20 2.0 x 1013 0 2600.
4 20H = O + H20 5.3 x 1012 0 503.
5 H2+M=2H+M 5.5 x 1018 -1 51987.
6 H20 + M = OH + H + M 5.5 x 1018 -1.5 59386.
7 OH + M =O + H +M 8.5 x 1018 -1 50830.
8 02 + M = 20 7.2 x 1018 -1 59340.

Units are moles, seconds, centimeters, and Kelvin
parameter range studied by Lehr exceeds that presented by Lehr (1972) and these additional

data have been provided by Lehr (personal communication, 1995). Here a brief review of

previous numerical studies is provided. In addition to the impact the reaction mechanism has

on the computed results, the numerical aspects including solution method and grid resolution

strongly affect the computed flow fields. Therefore, a description of these facets of each

study is included in the review along with key results.

Yungster et al. (1989) have computed the M=6.46, hydrogen-air case captured by

Lehr (1972) as well as the shots into hydrogen and oxygen presented by Lehr (1972). The

8-step reaction mechanism for hydrogen and air presented in Table 2.4 by Moretti (1965)

is used but the rate constants are replaced by values given by Evans and Schexnayder (1980).

Table 2.5. Eighteen-step hydrogen-oxygen reaction mechanism from Drummond
et al. (1986) used by Ahuja and Tiwari (1993).

n Reaction An bn Ean / Ru

1 02 + H2 = 20H 1.70 x 1013 0 24245.93
2 02 + H = OH + O 1.42 x 1014 0 8258.22
3 H2 + OH = H20 + H 3.16 x 107 1.8 1525.76
4 H2 + O = OH + H 2.07 x 1014 0 6923.81
5 20H = H20 + O 5.50 x 1013 0 3524.85
6 OH + H + M = H20 + M 2.21 x 1022 -2 0.
7 2H + M = H2 + M 6.53 x 1017 -1 0.
8 H + 02 + M = HO2 + M 3.20 x 1018 -1 0.
9 HO2 + OH = H20 + 02 5.00 x 1013 0 503.55
10 HO2 + H = H2 + 02 2.53 x 1013 0 352.48
11 H02 + H = 2 OH 1.99 x 1014 0 906.39
12 H02 + O = OH + 02 5.00 x 1013 0 503.55
13 2HO2 = H202 + 02 1.99 x 1012 0 0.
14 HO2 + H2 = H202 +H 3.01 x 1011 0 9416.39
15 H202 + OH = H20 +HO2 1.02 x 1013 0 956.75
16 H202 + H = H20 + OH 5.00 x 1014 0 5035.50
17 H202 + O = H2 + OH 1.99 x 1013 0 2970.95
18 H202 +M = 20H 1.21 x 1017 0 22911.53

Units are moles, seconds, centimeters, and Kelvin

The governing equations are solved using a fully coupled equation set and the total variation

diminishing (TVD) algorithm by Yee and Shinn (1989), also known as the "point implicit

TVD MacCormack." Grids of resolution 42x44 and 32x32 are used with the only difference

reported between the solutions being the "thickness of the captured shocks." Good

agreement between the computed shock and reaction front location and the experimental

value is reported. However, only a limited region encompassing the nose is modeled.

Lee and Deiwert (1989) have computed the M=6.46 case using the F3D code by Ying

(1986) which is based on the implicit flux vector splitting method. The governing equations

are solved in a loosely-coupled approach which assumes that in the Jacobians of the fluxes,

the effective specific heat ratio and pressure are not affected by the changes in species mass

fractions. The reaction mechanism used is the same as that of Yungster et al. (1989). Some

of the calculations are carried out with the rate coefficients reduced one and two orders of

magnitude in an effort to quantify the impact of reaction rates on the flow field. They do find

Table 2.6. Nineteen-step hydrogen-oxygen reaction mechanism used by Oldenborg
et al. (1990).

n Reaction An bn Ean / Ru

1 02 + H2 = H + H20 2.16 x 108 1.51 1726.
2 H + 02 = O + OH 1.91 x 1014 0. 8273.
3 O + H2 = H +OH 5.06 x 104 2.67 3166.
4 H + HO2 = H2 + 02 2.50 x 1013 0. 349.
5 H + HO2 = 20H 1.50 x 1014 0. 505.
6 O + HO2= OH +02 2.00 x 1013 0. 0.
7 OH + HO2 = H20 + 02 2.00 x 1013 0. 0.
8 H + 02 + M = HO2 + M 8.00 x 1017 -.8 0.
9 H + OH + M = H20 + M 8.62 x 1021 -2. 0.
10 2H + M = H2 + M 7.30 x 1017 -1. 0.
11 H+O+M=OH+M 2.60 x 1016 -.6 0.
12 20 + M = 02 + M 1.14 x 1017 -1. 0.
13 20H = O + H20 1.50 x 109 1.14 0.
14 20H + M = H202 +M 4.73 x 1011 1. -3206.
15 OH + H202 = H20 + H02 7.00 x 1012 0. 722.
16 O + H202 = OH + HO2 2.80 x 1013 0 3220.
17 H + H202= H2 + HO2 1.70 x 1012 0 1900.
18 H + H202= H20+ 02 1.00 x 1013 0 1800.
19 2HO2 = H202 + 02 2.00 x 1012 0 0.
Units are moles, seconds, centimeters, and Kelvin

that the slower rates do impact both the computed shock location and the temperature field

behind the shock. Their results are discussed further in the context of the results from the

current study. A grid of resolution 57 x 41 over the nose region is used and the results are

reported to be very similar to those by Yungster et al. (1989).

Wilson and MacCormack (1990) compute the steady combustion cases at M=5.11

and 6.46 using the 32-step reaction mechanism given in Table 2.7. The mechanism, as

presented by Wilson and MacCormack (1990), is in a slightly different order and there is a

typographical error in the pre-exponential constant for reaction 21. The constant An for

reaction 10 is different than that used by Jachimowski (1988). Wilson (personal

communication, 1997) reports this alteration has been made based on the rate information

available in other sources such as the work of Warnatz (1984). In the current study, the rate

for step 10 as defined by Wilson and MacCormack (1990) is used to facilitate comparisons.

A fully-implicit flux-splitting technique based on the work by MacCormack (1985) and

Table 2.7. Hydrogen-oxygen reaction mechanism from Jachimowski (1988).

n Reaction =An_ bn Ean / Ru

Units are moles, seconds, centimeters, and Kelvin
Third-body efficiencies: (6) H20 =6.0, (7) H20= 6.0, H2

: 2.0, (8) H20= 5.0, (9) H20 = 16.0, H2 = 2.0, (19) H20 = 15.0

Candler (1988) is used along with Steger-Warming flux splitting (Steger and Warming,

1981). The governing equations also include vibrational nonequilibrium effects; however,

Wilson (personal communication, 1996) reports these effects to be negligible. The study

evaluates increasing the induction delay by changing the reaction rate for H + 02 OH +

O, step 2 in Table 2.7. The rates

(a) k2 = 2.60 x 1014 exp(- 8459.64/T)

(b) k2 = 1.91 x 1014 exp(- 8277.35/T) (2.2)

(c) kf = 1.2 x 107T--91 exp(- 8315.12/T)

H2 + 02 = HO2 + H
H + 02= OH + O
O + H2 = OH + H
OH+ H2 = H + H20
20H = O + H20
H + OH + M = H20 + M
2H + M = H2 + M
H + 02+ M = HO2 + M
20+ M = 02+ M
HO2 + H = 2 OH
HO2 + H = H20 + O
HO2 + 0 = 02 + OH
HO2 + OH = H20 +02
2HO2 = H202 +02
H + H202 = H2 + HO2
O + H202 = OH + HO2
OH + H202 = H20 + HO2
H202 + M =20H + M
2N + M = N2 + M
N + 02 + NO + O
N + NO = N2 + O
N + OH = NO + H
H + NO + M = HNO + M
H + HNO = NO + H2
O + HNO = NO + OH
OH + HNO = NO + H20
HO2 + HNO = NO + H202
HO2 + NO = NO2 + OH
H + NO2 = NO + OH
O + NO2 = NO + 02
M + NO2 = NO +O + M

L ____________________ I

1.0 x 1014
2.6 x 1014
1.8 x 1010
2.2 x 1013
6.3 x 1012
2.2 x 1022
6.4 x 1017
6.0 x 1016
2.1 x 1015
6.0 x 1013
1.4 x 1014
1.0 x 1013
1.5 x 1013
8.0 x 1012
2.0 x 1012
1.4 x 1012
1.4 x 1013
6.1 x 1012
1.2 x 1017
2.8 x 1017
6.4 x 109
1.6 x 1013
6.3 x 1011
5.4 x 1015
4.8 x 1012
5.0 x 10"11
3.6 x 1013
2.0 x 1012
3.4 x 1012
3.5 x 1014
1.0 x 1013
1.16 x 1016



are evaluated and it is determined choice (c), the rate suggested by Warnatz (1984) gives

the best match for Lehr's experiments at M=5.11. This same rate is used to compute the

M=6.46 case resulting in a good match of both computed shock and reaction front to the

experiments. The M=6.46 case is computed on a 321 x 65 adapted mesh. In the current study,

the 32-step mechanism with the rates per Warnatz is labeled the modified 32-step model.

Wilson and MacCormack also evaluate altering the time of reaction by changing the rate

constants for H + OH + M H20 + M, step 6 in Table 2.7. This alteration affects the heat

release rate and the time required to reach the peak temperature. Reducing the original rate

per Jachimowski (1988) by a factor of 4 and using choice (a) above for H + 02 ; OH + O

produces no significant change in shock and reaction front location.

Ahuja and Tiwari (1993) use an 18-step mechanism, given in Table 2.5, to compute

the M=6.46 case. The mechanism is reported to be based on the work of Jachimowski (1988)

and matches the mechanism used in Drummond's work (1988) except for inconsistencies in

steps 9 and 10 given by Ahuja and Tiwari (1993) which are believed to be simply

typographical errors. The rate constants are not provided by Ahuja and Tiwari (1993) but

are presented here based on the data available from Drummond (1988). The equations are

solved using MacCormack's method and a Lax-Wendroff type scheme is used for the fluxes.

The grid resolution used is 197x152 over a domain comprising a quarter sphere. The

computed flowfield is reported to show a shock and reaction front but it is difficult to make

a detailed comparison to experiments given the limited computational domain.

Matsuo et al. (1993) use a 19-step reaction mechanism based on 32-step used by

Wilson and MacCormack (1990). In the 19-step mechanism the nitrogen reactions are

neglected and as mentioned, the specific reactions used are steps 1 through 19 of Table 2.7.

The governing equations are integrated explicitly, and the chemical reaction source term is

treated in a linearly point-implicit manner. Yee's non-MUSCL TVD upwind explicit scheme

is used and computations have been carried out for the M=4.18 and 4.79 cases captured by

Lehr (1972; personal communication, 1995). A grid of resolution 161x321 over a quarter

sphere domain is used. For the M=4.79 case, computations show the frequency of the

unsteadiness to be 725 kHZ where the reported experimental frequency is 712 kHZ. For the

M=4.18 case, a frequency of 159.8 kHz is computed where the experiments of Lehr (1972)

report a frequency of 150kHz. Also important in this study is their characterization of the

chemical time scales which is discussed later.

Sussman (1993) develops a source term modification scheme for use in combustion

computations. Using a model problem, the error associated with the representation of the

reaction source term is estimated. A scaling factor based on this error estimation is used in

computing both a quasi-one dimensional shock-induced combustion scenario as well as the

M=6.46 case captured by Lehr. This source term modification is discussed in detail in a later

chapter and is evaluated in the current study. The 19-step reaction mechanism presented in

Table 2.6 (Oldenborg et al., 1990) is used. The fluid dynamic terms are solved explicitly

using the Harten-Yee (Yee, 1989) TVD algorithm and the source terms are integrated using

a point implicit method. A grid resolution of 44x49 is used. In the quasi-one dimensional

case, the source term modification results in a reduction in the number of grid cells needed

to resolve the induction zone. This same scheme when used in Lehr's problem does not give

improved results without the addition of a grid adaptation technique.

Wilson and Sussman (1993) use the 32-step mechanism as in Wilson and

MacCormack (1990) with the rate suggested by Warnatz (1984). A logarithmic form of the

species conservation equations is used for which a description is provided later. A

comparison between Steger-Warming flux-vector splitting (Steger and Warming, 1981)

based on work by MacCormack (1985) and Candler (1988) and Harten-Yee (Yee, 1989)

upwind TVD scheme is made. The convective terms have been treated explicitly and source

terms implicitly. Both steady and unsteady cases have been computed. For the steady case

of M=6.46, a 52x52 grid is used with the logarithmic form of the equations. The computed

shock and reaction fronts from both flux schemes compare well with experiments. For the

unsteady case of M=4.79, a 375x161 grid is used but the standard form of the species

conservation equations is solved. The frequency of the unsteadiness is calculated to be 530

kHz where experiments are reported to be 720 kHz. The reported experimental value is

slightly different from the value of 712kHz reported by Matsuo et al. (1993). The

computations are repeated using the original rate for H + 02 OH + O per Jachimowski and

the frequency is computed to be 820 kHz. Matsuo et al. (1993) cite a private communication

with Wilson and Sussman that after running the computation for a longer time, the results

match the experimental frequency better than the 820 kHz cited. The findings that a different

rate for the H + 02 OH + O reaction are needed for the unsteady and earlier reported steady

cases is one area of concern regarding the available reaction mechanisms.

Matsuo and Fujiwara (1993) compute the unsteady shock-induced combustion flow

using the conditions from Lehr's experiments for M=4.18. The diameter of the nose is

changed to generate different cases. A two-step mechanism using progression variables is

used and the parameters are based on shock-tube data for 2H2 + 02 + 7Ar. No direct

comparison to Lehr's data is made but the computations are used to help explain the

phenomena causing the instabilities. The solution scheme used is the same as in Matsuo et

al. (1993). No details on the grid resolution are given.

Sussman (1994) presents results of computations for both the low and high frequency

cases studied by Lehr. The equations are solved using the technique of Wilson and Sussman

(1993) including the logarithmic form of the species governing equations. All computations

are performed on a 192x252 grid over a quarter sphere domain and the 19-step mechanism

of Table 2.6 (Oldenborg, et al., 1990) is used. Excellent agreement between the computed

and experimentally measured frequencies is reported for all cases except at the lowest and

highest projectile velocities. The discrepancy at the lowest frequency case is attributed to the

fact that the energy release front is closest to the body at this velocity and either a finer grid

in that region is needed or the neglect of the viscous effects is most prevalent here. The

discrepancies at the high frequency end of the spectrum are suspected to be due to the

difficulty in measuring the scales from the experimental results given their small amplitude.

Also of interest in Sussman's study is the discussion of the mechanism determining

the instabilities. Sussman reports finding an instability for projectile velocities of 2119 m/s

where Lehr (1972) found no evidence of unsteadiness. Sussman finds that this transition

point from an unsteady to a steady combustion scenario it directly related to the point at

which the induction time of the reaction becomes smaller than the energy release time. This

criterion is different from the projectile-to-detonation velocity criterion mentioned earlier.

The frequency of the instabilities is directly related to the induction time and given the fact

that the computed values compare well with the experiments over the majority of the velocity

regime, Sussman deduces that this aspect of the chemical process is being adequately

modeled by the reaction mechanism being used. However, since the transition point from

unsteady to steady combustion is over predicted, the energy release rate produced by the

model is greater than that in the actual experiments. This conclusion that the ratio of energy

release time, or reaction time, to the induction time is a key parameter is in agreement to the

findings of Matsuo et al. (1993) and those of McVey and Toong (1971). Sussman also

presents detailed computations of a one-dimensional shock-induced combustion scenario

and demonstrates the advantages of solving the logarithmic form of the species equations.

Yungster and Radhakrishanan (1994) computed the unsteady cases captured by Lehr

at M=4.18, 4.48, and 4.79. The 19-step mechanism of Matsuo et al. (1993) is used with the

original rates per Jachimowski. The governing equations are solved fully coupled using

Yee's 2nd order TVD scheme (Yee, 1989) and a grid resolution of 220x220 over half the

hemispherical nose region. Good agreement between computed frequencies of instabilities

and those measured by Lehr are achieved. No computation of the steady case is performed.

Shang et al. (1995) compute Lehr's M=6.46 case using the FDNS code (Chen et al.,

1993), a pressure based finite difference code using the TVD scheme of Chakravarthy and

Osher (1985). An operator-splitting scheme is used where the chemical kinetics source terms

are integrated implicitly with an ODE solver in a predictor step and then coupled with the

fluid dynamics in a corrector step (Rhie et al., 1993). The M=6.46 case is computed on a

uniform (87x97) and adapted grid (44x49) using a 7-step mechanism, presented in Table 2.2,

based on the 18-step mechanism by Drummond et al. (1986) presented in Table 2.5. This

mechanism is also documented in other work by Drummond (1988, 1990). Various

interpolation techniques for defining the species concentrations in the predictor step are

evaluated. Good comparisons with experimental results are achieved when either an adapted

grid or a high order interpolation is used. These interpolation techniques are discussed in

detail later and are evaluated in the current study.

Matsuo and Fujii (1995) compute unsteady cases using the conditions from

experiments of Lehr for M=4.79. The same solution scheme of Matsuo et al. (1993) and a

2-step reaction model similar to that of Matsuo and Fujiwara (1993) are used. The 19-step

mechanism (Matsuo et al., 1993) is used to perform a time integration of a reduced equation

for the species conservation with no convection or diffusion. Different activation energies

for the exothermic reaction of the 2-step mechanism are evaluated to generate different heat

release rates or time of reaction. Results show the reaction time directly impacts the

flowfield. Of the three activation energies evaluated, the one with a slow and the one with

a medium heat release rate generate periodic unsteady results. A rapid release (Ea=0)

produces a larger, more random disturbance.

Matsuo and Fujii (1996a) perform computations using the experimental conditions

of Lehr (1972) and Ruegg and Dorsey (1962). A grid resolution of 401x401 is used to

represent a quarter sphere domain and the 19-step mechanism used by Matsuo et al. (1993)

is used. The pressure, temperature and velocities of the experiments are used but the diameter

of the projectile is altered to generate steady, low frequency unsteady, and high frequency

unsteady cases. A zero-dimension time integration of the governing equations is used as in

Matsuo and Fujii (1995) to determine correlation between a defined Damkohler parameter

and the resulting combustion situation. The integration uses the conditions along the

stagnation line and assumes internal energy and total density are constant. This results in a

temporal profile depicted in Figure 2.3. The definition used for the Damk6hler parameter

S= tf/t
= [d/a2]/[T2/(dT/dt)max (2.3)

where d is the projectile diameter and the quantities with subscript 2 denote values behind

the shock. Based on their results, a critical value of T = 80 is defined below which the
unsteadiness is of a high frequency and above which it is a low frequency instability. No
universal parameter to classify steady as well as unsteady cases is found. Such classification
may require using the parameter defined in Matsuo and Fujii (1996a) as well as the critical
parameter of Vp/VD as determined by Lehr (1972) where Vp is the projectile speed and VD
is the detonation velocity of the mixture.
Matsuo and Fujii (1996b) compute unsteady cases from experiments of Ruegg and
Dorsey (1962). The 19-step mechanism used by Matsuo et al. (1993) is employed and the
governing equations are solved explicitly with the chemical reaction source term treated in

by shock dT/dt


Figure 2.3. Depiction of the temporal profile from a time integration of the
governing equations using the conditions behind the shock on the stagnation line.

a linearly point-implicit manner. Yee's non-MUSCL TVD upwind explicit scheme is used

on a 601x601 grid over a quarter-sphere domain. Results from grids of 201x301, 401x401,

and 801x801 showed similar flow features. No direct comparison to experimental results is

presented but conclusions as to the mechanisms causing the instabilities are made.

As is clear by the review here, the steady and unsteady shock-induced combustion

scenarios captured by Lehr have been used in many computational studies. The current study

computes the case at M=6.46 in a stoichiometric mixture of hydrogen and air. Therefore,

those previous studies which have computed this same case are summarized in Table 2.8.

Included in this summary are the particular fluid dynamics model and reaction model used

as well as the grid resolution employed. Also, a depiction of the flow domain computed and

resulting fronts are given. These images have been generated using the computed density

field presented in each reference and orientating a line depicting the shock and energy release

fronts based on this data. It is clear that those computations which do not cover a substantial

portion of the projectile flow domain do not provide enough information to judge the

accuracy of the computation. Therefore, the current study computes the entire length of the

projectile as depicted in Figure 2.4 which is a larger domain than any of the previous

computations. Also denoted in Figure 2.4 are the domains used in earlier computational

studies in which the M=6.46 steady combustion case has been simulated. The utility of the

larger domain is to be evident later.

2.3 Summary

As detailed here, the shock-induced combustion problem experimentally studied by

Lehr (1972) and computationally studied by several investigators is to be used as an initial

problem in developing the needed models for simulating the munition concepts discussed

earlier. As is evident in the review, both the phenomenology of the various models and the

numerical aspects are key in determining the accuracy of a simulation for this class of

problems. For example, the chemical time scales of induction and reaction time are found

to be key. However, this trend in the physical and chemical processes is only detectable if

Table 2.8. Summary of previous computational studies into steady shock-induced

Lehr (1972)

Yungster et al. (1989)
Mechanism: 8-step, Table 2.4
Grid : 42x44
Scheme : Fully implicit
Lee and Deiwert (1989)
Mechanism: 8-step, Table 2.4
Grid: 57x44
Scheme : Fully implicit

Wilson and MacCormack (1990)
Mechanism : 32-step, Table 2.7
Grid : 321x65 adapted
Scheme : Fully implicit
Ahuja and Tiwari (1993)
Mechanism: 18-step, Table 2.5
Grid: 197x152
Scheme: MacCormack

Sussman (1993)
Mechanism: 19-step,Table 2.6
Grid : 44x49
Scheme : Explicit-Implicit
Notes : Comparison improves with adapted grid.

Wilson and Sussman (1993)
Mechanism: 32-step, Table 2.7
Grid : 52x52
Scheme : Explicit-Implicit
Notes : Solve logarithmic form of equations.

Shang et al. (1995)
Mechanism: 7-step, Table 2.2
Grid: 87x97
Scheme : Explicit-Implicit
Notes : Comparison improves with extapolation scheme.

(1) Extent of domain used by
Yunster et al. (1989)
Lee and Deiwert (1989)
Ahuja and Tiwari (1993)
(2) Extent of domain used by
Wilson and MacCormack (1990)
Sussman (1993)
Wilson and Sussman (1993)
Shang et al. (1995)
(3) Extent of domain used
in current study

Figure 2.4. Computational domain used in previous investigations and
throughout current study.

an adequate spatial discretization and reaction model is used. In the current study, not only

are the computations compared to experiments, but the change in the dependent variables

as a result of alternative models is used to help quantify the roles of the individual models

on the accuracy of the complete computational tool. This not only aids in determining the

impact these aspects have on the computation but also helps quantify the phenomenological

and numerical aspects involved in solving a class of problems such as the fuel cell scenario.

This study uses the steady combustion case shown in Figure 2.2 to assess the aspects

of modeling high-speed reacting flows. If the conditions for M=6.46 are used in the

Damkihler parameter definition of Matsuo and Fujii (1996a), Eqn. (2.3), a value of 18 is

calculated. This supports the use of this flow scenario to study the computation of reacting


flows where 2) = 0(1). As mentioned, the optical distortion in the image raises some

concerns in a one-to-one comparison between the computations and image. However,

instead of re-interpreting the image, the locations denoted in Figure 2.2 are used as in the

case of previous investigations. It should also be noted that the experimental data is not

symmetric about the stagnation line which suggests the projectile had some slight angle of

attack during its flight. Results for the steady shock-induced combustion case are discussed

following the description of the computational models developed here.


In this chapter the governing equations and various numerical techniques needed to

solve the current class of problems are detailed. After the presentation of the governing
equations, the individual approaches are discussed. For each model there are various options

as to the particular phenomenological and numerical form. Where appropriate these options

are discussed.

3.1 Governing Equations
The governing equations for the current problem are the axisymmetric form of the

Navier-Stokes equations with appropriate source terms. The integral form of the governing

equations is

f Q dT + F.dS = dT (3.1)

where Q is the dependent variables in volume Y, F is the flux of those variables across the

surface (S) which defines the volume T and Q is any source term per unit volume for the

dependent variables. The flux, F, is a non-linear function of the dependent variables Q. The

governing equations once discretized for a general curvilinear coordinate system are

a J Q a J (F, Fy) 9 J (GI Gy)
at + + J (H Hy) = J 2Q (3.2)

where the dependent variables are

Q = [Q, Qu, v,E,QXi ....QNs-1] (3.3)

The variables Q, u, v, and E are respectively the density, x component of velocity, y

component of velocity, and energy. Also, ai is the mass fraction of ith species with the fluid

being defined by NS total species. Note that the mass fraction of the NSth species is not
explicitly modeled since the total density is included and the relationship p = Qpti holds.

Another modeling option is to not solve an equation for total mass but to solve an additional

species conservation equation for the NSth species and calculate the total mass using the

summation given above. With the addition of the relationship between the individual species

densities and the total mass, both implementations require attention to numerical round-off

issues. This is evident in previous investigations such as the work of Sussman (1993) and

Wilson and Sussman (1993) in which it is noted that species densities can sometimes become

negative due to numerical errors and that limiters must be used.

In the curvilinear directions and r, the inviscid flux vectors are

puU + SxP puV + YTxP
QvU + yP QvV + ryP
FI = U(QE + P) GI = V(QE + P) (3.4)
0aU galV


where P is pressure, and the viscous flux vectors are

0 0
Jxtxx + lytxy Txtxx + TlyXxy
Jxtxy + yTyy x'Txy + 1y'tyy
FV = xx + MyIy G = lxnx + ryy (3.5)
QalUl QaiVi

eNS- i NS-l QaNs-IVNs-

The third dimension is incorporated using the axisymmetric terms

Ogv 0
Quv txy
pv2 Tyy Too
H = v(E + P) H (3.6)
H y
Qalv eQiv1


where y is the distance from the line of symmetry and the vector of chemical reaction source

terms is

Q = [0,0,0,0,60,,...,oNS-_ (3.7)

The grid Jacobian J, an indication of the cell volume, and the contravariant velocities are

defined as

J = xyr) xy
U = xu + yV (3.8)
V = lTxu + T]yV


x -y x
X J y T)

y, x (3.9)
fix j 11y = -J-

The components of the viscous flux vectors are defined as

Iix = uTxx + vTxy qx
Iy = utxy + VTyy qy

Ty = 2 + + + (3.10)
aay (ay x Y

ee = 2 + (ay ax Y)

where the bulk viscosity is represented using Stokes' hypothesis, = The heat flux

components are

q= -k + h(aG^), qy= k T + _hi(aivi) (3.11)
qx ax + ay
i=l i=l

where k is the thermal conductivity and ai i and aiVij are the diffusion velocities defined in

the next section.

3.1.1 Molecular Diffusion
The vector components Ui iV are the contravariant diffusion velocities for species

i and are represented using Fick's law as described by Anderson (1989), and Libby and

Williams (1994). The exact form of the modeled terms are
Ui = xUi + yVi Vi = TIxUi + TlyVi
A Oa aO (3.12)
aui = D- aiiv = Di ay

The use of Fick's law is made possible by some simplifying assumptions as to the molecular

diffusivity of individual species. A more exact representation of the diffusion of a single

species i in a mixture composed of NS species is to calculate the diffusivity of the species

i as a function of its diffusivity in each of the other species. One level from this in

simplification is to define a background gas or the gas which comprises the majority of the

mixture. This allows the diffusivity of each species in the predominant gas to be calculated

and used for Di. For example, combustion problems involving hydrogen and air typically

contain a large amount of nitrogen which can be defined as the background gas. The next

level in simplification is to assume a single diffusion coefficient such that Di=D for all

species, and that the binary diffusion can be related to the viscosity using the Schmidt


Sc (3.13)

This simplification, used in the current study, has proven reliable in many computational

studies and is considered acceptable when all the species involved have approximately the

same molecular weight. The only concern in the context of the present study is the hydrogen

involved given its molecular weight is one order of magnitude less than the other species.

However, this simplification is employed given the predominance of the other transport

mechanisms for the species in the current flow scenario and considering the reduction in the

computational requirement. Also, in the context of diffusion itself, the molecular diffusivity

is secondary to the turbulent enhancement which is discussed later.

3.2 Equation of State and Temperature

The equation of state used in the current computations is derived by assuming the

ideal gas equation is valid for each individual species and has the form (Anderson, 1989)

P = RuT ZMi (3.14)
where Ru is the Universal Gas Constant and Mi is the molecular weight of species i. There

are flows for which an alternate equation of state is more applicable and from a

computational model standpoint, such alternate equations simply require an alternate

definition of terms for which the derivative of pressure with respect to the conserved

variables are needed. These modeling aspects are presented in a later section.

The temperature during the calculations must be extracted from the internal energy

using the relationship

e=E -u = ah
T (3.15)
hi = hf0 + CpidT

where TR is the reference temperature and hf the heat of formation of species i at TR. The

specific heat, Cpi, of each species is a known function of temperature and is available from

various references such as the JANNAF data (Stull and Prophet, 1971). The dependence of

Cpi on temperature can be represented in various forms ranging from a constant up to a high

order polynomial. In the shock-induced combustion computations presented here, a third

order polynomial of the form

Cpi = c1 + c2T + c3T2 + C4T3 (3.16)

is used. The coefficients have been determined by fitting curves to the JANNAF data up to

6000K and are given in Table 3.1. Higher order polynomials have been used by previous

Table 3.1. Coefficients used to calculate specific heats for the various species.

n c1 C2 c3 C4 hf

02 6.72433 1.75132 x 10-3 -3.35085 x 10-7 2.40485 x 10-11 0.
H2 6.18833 1.28193 x 10- -1.59252 x 10-7 9.05493 x 10-12 0.
OH 6.96378 5.90649 x 10-4 2.75234 x 10-8 -1.03921 x 10-11 9.432
H20 6.9803 3.6993 x 10-3 -6.55318 x 10-7 4.07621 x 10-" -57.798
O 5.25211 -2.0219 x 10-4 3.79786 x 10-8 0 59.559
H 4.968 0. 0. 0. 52.1
HO2 7.44933 7.3918 x 10-3 -2.8673 x 10-6 3.955 x 10-'1 5.
H202 6.6606 1.49214 x 10-2 -7.7563 x 10- 1.3478 x 10-9 -32.53
N2 6.51307 1.65586 x 10- -3.74607 x 10-7 2.8942 x 10-" 0.
Cpn in cal / (mole K) hf in kcal / mole TR = 298 K

investigators but here, the third-order has been found to provide a good match using a single

set of constants across the temperature range. Given this representation of specific heats, an

iterative Newton-Raphson procedure must be used to extract the temperature in each cell at

each time step as detailed by Molvik and Merkle (1989). The procedure solves the equation

em+ = em + e mAT (3.17)

where m denotes the iteration step. The initial guess for the temperature in each cell at a given

time step is simply the temperature at the previous time step and convergence is typically

achieved in only a few iterations if an over-relaxation scheme is used.

Even though convergence is typically rather fast, it is clear to see that for a

computation with a larger number of computational cells, this added aspect for reacting

flows and for that matter high temperature flows can become costly. Therefore, for some

applications a lower form for Cpi is worth using. In fact, if a linear representation for

dependence of specific heat on temperature is used, the temperature can be calculated by

solving a simple algebraic equation with no iteration. This representation enables the aspect

of temperature dependence to be included with no added computational time. However,

before using such a scheme, attention should be paid to the temperature range over which

the computation is expected to vary to ensure the constants used best match this regime.

3.3 Thermal Conductivity and Viscosity

For the current study, it is assumed that the thermal conductivity can be determined

by specifying a Prandtl number where

Pr = p (3.18)
This as the case with the diffusivity is a simplification which has proven applicable to the

class of flows modeled here. Note, the Schmidt and Prandtl number are related through the

Lewis number (Le = Sc/Pr) which is often near unity. In the current study Le= 1 is assumed

and the Prandtl number is set to .72. We now have the binary diffusion coefficient and the

thermal conductivity related to the viscosity of the mixture by a constant of order one.

The viscosity of the mixture is modeled as a function of temperature and the species

composition using Sutherland's law as described by White (1974),

Vi T 5 Ti + Si
oi T + Si (3.19)

and Wilke's formula

X j(ij j]
= Nsij = + 8( 1)2()][( + 2 (3.20)


where Xi is the molecular concentration of species i and Toi and Si are constants. The

constants for use in Eqn. (3.19) are given in Table 3.2. One of the requirements for the fluid

mechanics model is to limit the numerical errors in the models of quantities such as the

Table 3.2. Coefficients used to calculate viscosity for the various species.

diffusion of species and heat. Given the fact that both these quantities are related to the

molecular viscosity with an order one constant, it is reasonable to conclude that the

numerical errors in these two quantities are of the same order as numerical errors in

molecular viscosity. The apparent errors in viscosity for typical computational models are

not a function of the phenomenological model such as Eqn. (3.20) but are a result of

numerical viscosity. Numerical viscosity is a direct function of the numerical scheme used

in representing the governing equations. The particular scheme used in the current study is

described next and the issue of numerical viscosity is addressed in that context. However,

it is safe to assume that if the numerical viscosity is minimized, so will be the numerical

errors in diffusion and conductivity.

3.4 Flux Vector Splitting Scheme

For the solution of compressible flows, the most popular schemes for representing

the inviscid fluxes are flux difference and flux vector splitting. These upwind schemes are

attractive due to their ability to minimize numerical viscosity. A review of these various

schemes for multi-species flows is provided in the work of Montagne et al. (1987), Glaister

(1988), Cinnella (1989), Grossman and Cinnella (1990), Shuen et al. (1990) and Liou et al.

(1990). Here the scheme used to define the inviscid numerical fluxes is the Steger-Warming

flux vector splitting (Steger and Warming, 1981) which has been extended to model

multi-species flows by Liou et al. (1990). The flux vector splitting algorithm decomposes

n -o (Ns/m2) To (K) S (K)

02 1.919x 10-5 273 139
H2 8.411 x 10-6 273 97
HzO 1.120 x 10- 350 1064
N2 1.663 x 10-5 273 107

the inviscid fluxes into components based on the eigenvalues of the Jacobian AI = and

likewise for GI. The split fluxes take the form

K* = kXK, + X:K2 + O:K3 (3.21)


Qu (u kxc)
K _- 1 ( c2 1 (v kyc)
yK h y 1 K2,3 y Q(ht ) (k22)
Qal QaI

L QaNS-1 QaNS-1

with the eigenvalues being

Xk1 = (k Ikl)
S= 13k c(3.23)
2 = Pk + clVkI
k3 = Pk clVkl


Ok = kxu + kyV Ok = kxu + kyv
k- kY 1- -- (3.24)
kx =kx ky k IVI =/k2x + k2
IVkl IVkl IVkl k- +k(

The above formulation gives K=FI when k=E and K=GI when k=ri. For the multi-species

chemically reacting flow, c is the frozen speed of sound and y is the effective specific heat

ratio. The calculation of c is detailed in the Appendix.

In the finite volume formulation, the fluxes are evaluated at the cell faces and are

either first or second-order spatial representations. For the class of flows to be solved here,

a second-order scheme is desirable to limit the numerical viscosity. The flux at the face is

a function of the states in the neighboring cells,

(F),i+,j = F+(QL+Q + F-(QR+ (3.25)

If a first-order spatial representation is used, then Q = Qi, Q ,i = Qi+,j .To

achieve second-order accuracy, a MUSCL approach described by Shuen (1992) is used in
which cell-center values are extrapolated to the interfaces. Also, to guard against the

interpolation introducing any nonphysical extrema into the field in the region of large

gradients, a limiter must be used. The formula for the neighboring states takes the form
QL Qij + Q i+,j

iQR, Q J (3.26)
S',i+j = ij i+,J

where the limiting function is

S j = ,j minmod[AQ+i+I j,AQ-i+1,j
2 .( 3 .2 7 )
+j =i 2 minmod[AQ AQ ,j]


2(Q i+ ,j i- Q,j
AQ+ij i+ j + ij(

2(Qi, Q(3.28)
AQ-'J fi,j + i-lj

Various limiters can be used with this implementation but here the minmod limiter is used


minmod [X, Y] = sign(X) max[0., min(IXI, Ysign(X))] (3.29)

Note e ij is the cell-length to provide weighting for nonuniform grid spacing. The same

extrapolation procedure is carried out for the fluxes in r~ and can be performed on either the

dependent or primitive variables. Previous investigations by Shuen (1992) have shown that
using primitive variables gives better performance for flows with strong shocks and this is
the method used here.

Note, the second-order accuracy is achieved by the assumption that the field
quantities vary linearly across a cell. Higher order schemes are achievable if a higher order

variation is assumed. However, in all cases, limiters must be used in conjunction with the

higher order schemes to ensure against the introduction of any nonphysical extrema in

regions of large gradients. The viscous fluxes are evaluated at the cell faces using central

differencing and second order representations of the derivatives.

3.5 Chemical Reaction Source Term

The particular reaction mechanisms used in the current study have been discussed

in the earlier literature review of shock-induced combustion simulations. The production of

each species is incorporated into the governing equations by way of the source terms in Eqn.

(3.7). This source term is calculated using the finite rate chemistry model described here. The

reaction between hydrogen and oxygen which drives the shock-induced combustion must

be modeled by a specific number of reactions (NR) using the form

NS kfi NS
vi Xj <, j1j I i = NR (3.30)
j=l kbi j=l

where vi,' and vii" are the stoichiometric coefficients for species j in reaction i, and Xj is

the molecular concentration of species j. The parameters kf, and kb, are, respectively, the

forward and backward reaction rate constants for the ith reaction and can be expressed as a

function of temperature with the Arrhenius expression

ki = ATbexp Ea (3.31)

where Ru is the universal gas constant, Ea, is the activation energy and Ai and bi constants

for the ith reaction. The production rate of an individual species is calculated by summing

the contribution due to each reaction

j = M i i Vij ')x k XVn' kb, XV i (3.32)
i=l n=l n=l
The forward rate constants are available in various sources for most reactions of

interest. The backward rate for a reaction is sometimes given explicitly using Arrhenius rates

but is usually obtained using the equilibrium constant and the relationship

kei = -. (3.33)

The equilibrium constant can be expressed as a function of the stoichiometric coefficients

and the Gibbs free energies (gn) of the species in the reaction through the formula
ke,= [RuT] ,= exp RT (3.34)


Ag = (ni" Vni') gn (3.35)
The Gibbs free energy (gn) associated with each species is given as a function of temperature,

gn = gn(T). Using the JANNAF data (Stull and Prophet, 1971), it is clear that a linear

function of the form

gn = gn,, + gn,T (3.36)

is adequate. The values of the constants for the species to be used throughout the current work

are provided in Table 3.3 as well as their molecular weights. From a solution perspective,

the chemical reaction source terms are numerically stiff and require special attention which

is discussed later.

One note from the perspective of implementing a finite-rate reaction model as

described here is needed. The necessary routines can be coded up to handle a generic reaction

mechanism. This requires multiple do loops to calculate the source term as well as the

Jacobian of the source term (aQ2/Q) which is discussed later. This presents a problem when

running such a code on a vector machine, or any machine for that matter, and literally orders

of magnitude in computational time can be saved if the routine is coded up for a specific

reaction mechanism. This is easily done if an analysis tool such as Mathmatica (Wolfram,

1988) is used to preprocess the needed equations and is well worth the time when large

mechanisms or large grids are used.

Table 3.3. Molecular weights and coefficients used to calculate the Gibbs free energy for
various species.

n Mn go gl

02 32. 0. 0.
H2 2.016 0. 0
OH 17.08 9.1646 -3.4 x 10-3
H20 18.016 -59.29 1.35 x 10-3
O 16. 60.464 -1.574 x 10-2
H 1.008 53.3078 -1.396 x 10-2
HO2 33.008 4.195 1.194 x 10-2
H202 34.016 -33.251 2.646 x 10-2
N2 28.016 0. 0.
N 14.008 114.407 -1.608 x 10-2
NO 30.008 21.468 -2.94 x 10-3
HNO 31.016 23.2391 1.173 x 10-2
NO2 46.008 7.738 1.514 x 10-2

gn in kcal/mole

3.6 Special Treatment For Chemical Length Scales

In the discussion of the flux-vector splitting scheme used in the current code, a

detailed description is provided as to how to represent the numerical flux with a second-order

spatial model. Such a model reduces the numerical dissipation of field quantities and

improves the accuracy of the scheme. The second-order model is approximating the

variation in variables across a cell as a linear function. As noted earlier, higher order schemes

are achievable if a higher order variation is assumed. Second order schemes have proven to

provide acceptable accuracy for simulations of compressible turbulent flows with little

increase in computational requirements. However, this is due to the fact that the

computational grid used in the simulations employ discretization of the appropriate size to

capture the scales of the mean flow. With the addition of chemical reactions, the length scales

which must be resolved can be orders of magnitude smaller. Of course one option is to

construct the computational grid using a finer discretization. However, unless the location

of key events such as the reaction front are known a priori this finer discretization must be

used throughout a large portion of the flow field. An alternative is to employ a grid adaptation

technique which is only advantageous if it is closely coupled with the fluid dynamics model.

Both these options present an increase in the computational requirements.

Other options of addressing the error associated with the small chemical length scales

have been evaluated. Wilson and Sussman (1993) have used a transformed governing

equation while Sussman (1993) and Shang et al. (1995) have tested source term modification

techniques. These techniques take advantage of the known profiles of those species created

in the chain-branching reactions. The variation of these species, such as OH, H, and 0, is

more logarithmic in nature. The transformed governing equation used by Wilson and

Sussman (1993) is discussed briefly but attention is paid to the schemes which are an

alteration only to the source term. Typically the source term is calculated using the cell

averaged value for temperature and mass fractions, a piecewise representation. The modified

treatment calculates the source term using an altered value for the mass fractions,

mij = 0(as)ij. The scheme used by Shang et al. (1995) differs from that of Sussman (1993)

so both are described here. All the source term treatments are an effort to address the fact

that the cell averaged source term is not necessarily equal to the source term calculated based

on the cell averaged dependent variables. This is due to the highly nonlinear nature of the

Arrhenius source term.

3.6.1 Transformed Governing Equations

The modification option evaluated by Wilson and Sussman (1993) is a logarithmic

transformation of the entire species conservation equation minus the diffusion term. The


nn -= ln() Qln(an) (3.37)

is used with the total mass conservation equation to specify the transformed equation for

species n,

An+ a anuj= n (3.38)

In practice these transformed equations are solved in conjunction with conservation

equations for the individual elements to ensure the conservation of mass. Solutions by

Wilson and Sussman (1993) for the M=6.46 case of Lehr using this set of transformed

equations produce good results on somewhat coarse grids. Wilson and Sussman state this

formulation is no less robust than the conventional form of the equations but do warn against

the added distortion at contact discontinuities between different gases. This could be key in

non-premixed combustion scenarios. The computation time is cited to be 10-20 % more per

time step above the standard scheme. Also, of interest is the error trend when this scheme

is used on coarse grids. If the shock-induced combustion flow scenario is considered, the

flow structure is such that behind the shock, a small induction zone can exist which lies

between the temperature rise due to the shock and the temperature rise resulting from the

reaction front. If the standard form of the governing equations are solved, a coarse grid

results in a computed flow field such that the shock and reaction fronts are merged and the

induction zone is not resolved. Sussman and Wilson (1991) have shown that a coarse grid

used in conjunction with the transformed governing equations results in an overprediction

of the induction zone. In short, the computed reaction front is positioned too far down stream.

Therefore, even though the transformed governing equation method has shown promise

when applied to the shock-induced combustion scenario (Wilson and Sussman, 1993), this

grid dependence should be considered when analyzing the results.

3.6.2 Logarithmic Interpolation

The scheme of Shang et al. (1995) uses a logarithmic interpolation weighted by the

mass flux into the cell and can be explained using the image of Figure 3.1. In this case the

mass flux (rm) into cell (i,j) is from cell (i-l,j) and (i,j+l). The scheme first estimates the

mass fraction at the cell faces using

Log (i-,j = Log a-1,'j) + Log(ai)]

Logi+) = Log ) + Log ) (3.39)
L W = ILogij) + Lg(aij)]

Figure 3.1. Schematic of mass fraction interpolation scheme per Shang et al. (1995).

and then calculates the mass fraction to be used in the source term evaluation by

m. i_,Log (i + mi. ,Log(t +
2 J, ij +
cs = Exp +. (3.40)
m. .+ mn.
-? j '+

where a represents each species mass fractions. The interpolations are a sum of all cell faces

for which the mass flux is into the cell. In actual implementations, Log functions are singular

when the mass fractions are zero. Therefore, after simple arithmetic, the interpolated value

for the example case is simply

r/ \m \ 2!m 1 +m 1 \
as = 2" i-lj)' T,+l) \ ++ (m i j) (3.41)

Shang et al. (1995) present both solutions of a one-dimensional auto-ignition case and the
M=6.46 case captured by Lehr (1972). In the one-dimensional case, the modified source

term treatment produced "exact" solutions using only about a fourth of the grid points

required to achieve the same accuracy with the piecewise representation. This interpolation

scheme also improves the comparison between the experimental data and the computed flow



field for the M=6.46 case. This is accomplished with no addition of grid cells or adaptation.

These results are discussed further in the context of the results of the current study.

3.6.3 Error Estimation Based on One-Dimensional Model Equations

The source term modification proposed by Sussman (1993) is based on the analysis

of a one-dimensional model problem of the governing equations. The analysis is

summarized here for completeness. Consider the model problem

aQ aQ
+ u = kQ
at ax (3.42)
Q(, t) = Qo

and when discretized using a first-order upwind scheme is

Q+l = Q f At Qn-) + Atn (3.43)

where the source term is

Qn = Q(Q) = kQ (3.44)

Note, in Sussman's (1993) analysis, u is set to 1. The steady-state solution can be

approximated by

Q(x) = Qoekx/u (3.45)
and the source term by

2(Q) = kQ (3.46)
where k=kf. The values of the discretized function and source term, using a constant spatial

discretization, are

Qi = Qoek(i-1)Ax/u
gi = kQek(i-1)Ax/u (3.47)

Therefore, the approximation of the spatial derivative of Q at point i, using first-order

upwinding, is

Di- = Q = Qoek( 1)Ax/u[ -e-kAx/u] 1 (3.48)

The error present in a steady-state solution of the model problem can be estimated by

analyzing the ratio of this operator to the source term,

D- [1 e-kAx/u]
QTi kAx/u (3.49)
DZi kAx/u

Recall that k represents the reaction rate and 1/k is the characteristic time scale of the

reaction. Therefore, kAx/u is the ratio of the actual chemical time scale to the fluid dynamic

time scale, associated with the grid spacing, used in the computation. This is a cell

Damk6hler number, D)cfd. The variation of the ratio in Eqn. (3.49) with kAx/u is given in

Figure 3.2. Both the figure and Eqn. (3.49) demonstrate that as kAx/u oo the ratio

0.9 -

0.8 -

DI/W 0.7-

0.6 -

0.5 -

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Figure 3.2. Discretization error for Sussman's (1993) model problem.

Di-/Qi is asymptotically proportional to [kAx/u]-1. This relationship has direct

implications to the solution of chemically reacting flows and demonstrates the expected

trend. In a given flow scenario, the velocity is set by the solution as is the rate k, through the

field variables and the reaction mechanism used. The advantage of adapting the grid to

reduce Ax in regions of large k is obvious given the fact that the error changes exponentially

with kAx/u. However, the effectiveness of such adaptation is dependent on the magnitude

of k and it is clear for large k, Ax must be very close to zero to reduce the error. We next

consider alternative ways of reducing the error introduced through the discretization of the

exponentially varying function.

Sussman (1993) uses the information from the error analysis to propose a

modification to the source term evaluation using a scaling factor f such that a, = f ij. In

the context of the model problem,

O(Qi) = Q(fQi) = kfQi (3.50)

where f is a scaling factor which gives Di-/Wi = 1 if

f [ e-kAx/u]
f = (3.51)
In a general problem f is not known a priori because k is not a constant. Sussman (1993)

estimates k using the exact solution of the model problem,

k= u ln( (3.52)

Using this estimate, which is exact if Q is truly an exponential function in space, the scaling

factor when using a first-order upwind scheme is

f = n(r) r = Qn /Q (3.53)

and for a second-order upwind scheme, is
3 2r + ir2
f = 2 )2 r /Qp (3.54)
f ln(r) r I (3.54)

As presented by Sussman (1993), to achieve the exact solution of the model problem with

this modified source term treatment and a first-order scheme requires only 1/8th the number

of grid points if a piecewise treatment of the term is used with a second-order upwind


In applying this modification to a quasi-one-dimensional shock-induced combustion

case for hydrogen-air, the modification is applied to only those species produced by the

chain-branching reactions. For application to two-dimensional flows, it is assumed the
exponential variation occurs in the direction of the flow and k is estimated by

k = -- In(ri) ln(r) (3.55)
Ax Ay J

cos 0 = u sin = (3.56)
(u2 + v2)2 (u2 + v2)

This results in a modification factor for use with second-order upwind schemes of

Cos1 2ri + r2? + in ( 2r + r2?
f Ay (3.57)
csn(ri) ln(rj)

where the parameters ri and rj are determined using upwinding,

( Ji-lj/Qacij, u > 0; fQi,j- /Qi,j, v > 0;
ri = [Oi+l,j/ai,j, u < 0 =i lai,j+l/eai,j' v < 0 (3.8)

A separate factor is calculated for each species a produced by the chain-branching reactions.
To ensure the modification factor is bounded, the limits ri < 1 rj < 1 are imposed.

Sussman (1993) demonstrated this technique for the M=6.46 case of Lehr but found the
scheme did not produce results comparable to the experimental data unless an adapted grid
is also used. These results are discussed in more detail later in conjunction with the current
Given the above discussion, it is of interest to compare the modification of Shang et
al. (1995) in the same context. Beginning with Eqn. (3.41), and assuming the mass fractions
in the neighboring cells enter the formula through a similar representation of k in the two
dimensions as that used in Eqn. (3.52), we have

e-kAx/u Ci- ,j e-kAy/v + (3.59)
"ij ,i,j

and after simple arithmetic, Eqn. (3.41) becomes

[-m kAx/u -m ikAy/v 2m ( .,+m) ) ,
as = e -TJ e -+ 2

where v is used to represent the velocity in the second direction. For comparison purposes,

consider the one-dimensional case as done by Sussman, we have Ay=0 and mn = 0.

Therefore, the logarithmic extrapolation results in a scaling factor of

f = e-2Ax/u (3.61)

which at first looks quite different from that of Eqn (3.51) but as evident in Figure 3.3, when


2 2.5 3 3.5 4

Figure 3.3. Comparison of the scaling factors from the logarithmic extrapolation
and error estimation.

plotted for a range of kAx/u is comparable to the scaling factor by Sussman (1993) for

first-order upwind. Also presented in Figure 3.3 is the scaling factor for the second-order

upwind scheme of Eqn. (3.54). Global Order of Accuracy Analysis

The scaling factor proposed by Sussman (1993) has been derived by analyzing the

local error introduced by the piecewise representation of the source term. It is also

enlightening to consider the global order of accuracy and the impact of this scaling on the


magnitude of the global error. The global order of accuracy is dependent on the accumulation
of error at the locations upstream of the local computational stencil (Shyy, 1994) as depicted
in Figure3.4. This is particularly of interest in the current reactive flow problems since one



Figure 3.4. Depiction of region over which error accumulation affects the global
error of accuracy.

of the key flow field features is the variation in mass fractions.
Consider the model problem, Eqn. (3.42), used by Sussman (1993) which has an
exact solution for the steady-state case of the form

Q(x) = Qoea (3.62)
where a = k/u. Employing a first-order upwind difference operator with the piecewise
source term representation, the difference equation used to represent the model differential
equation is

q(x + h) q(x)
h a q(x) (3.63)
where h is the grid spacing and for the current analysis constant spacing throughout the
domain is assumed. The exact solution of the difference equation is

qn = (1 + ah) qn-1
= (1 + ah)[(l + ah)qn_2]
= (1 + ah)n qo (3.64)
qn = (1 + ah)x"/h oq

and to estimate the global error the solution of the difference equation is subtracted from the
exact solution of the differential equation
E(xn) = Q(xn) qn
= Qoeax (1 + ah)X"/hqo (3.65)
Since the initial conditions for the difference equation are specified, qo = Q0 and the

leading order of the global error is dependent of eaxn (1 + ah)x"/h. The second term can
be expanded as

(1 + ah)xn/h = exp I--Log(l + ah)

= exp [(ah 2 + (h3))
a2X h )1 (3.66)
= exp axn a2h + O(h2)
= exp(axn) exp(a2h) exp(O(h2))

where it has been assumed ah < 1 in the expansion of the logarithmic function. In
implementation this may not be true but here the global order of accuracy, and specifically
how the error decreases as h is reduced, is being assessed so the expansion is applicable.
Using the expansion

ex = 1 + x + x2/2! + x3/3! +...
Eqn. (3.66) becomes

(liahx)n/h i1- an2x~h
(1 + ah) /h = exp 1 a2 O h2)
ea2 (3.67)
= e2 h (eax" + O(h2)

and the global order of accuracy is

E(xn) = Qo h eax" O(h2)] (3.68)

This results in a normalized error, when the piecewise source term is used, of

Enorm(xn) = ( 0O(h) (3.69)
which is an indication of the error at locations far from the boundary. This error, an

accumulation of errors across the region from where the boundary conditions are imposed

to the point of interest, can be key in problems such as the current shock-induced combustion

flows where features such as induction zones and reaction fronts are of interest. This is due

to the fact that regions such as induction zones are areas over which quantities such as

radicals are growing. Since factors such as the mass fraction of radicals directly impact the

resulting chemical process, any error in their growth pattern introduce inconsistences

between the simulated field and what is actually experienced.

When the scaling factor proposed by Sussman (1993) is used, the actual difference

equation being solved is somewhat different from that of Eqn. (3.63) since the scaling factor

is a function of the dependent variables. The exact form, using a first-order upwind scheme,


q(x + h) q(x) 1 q(x+h)
h = a q(x) qx) (3.70)
Log(x + h)

which can be simplified to

qnLog( qn-1Log ( ) + ah(q- qn-1) = 0 (3.71)

and further reduced to

(qn qn-) Log( + ah = 0 (3.72)

which has two roots,

qn = qn-1 qn = eah+Log(qn) (3.73)
Dismissing the first root due to the nature of the problem being solved, the exact solution

of the difference equation is

qn = eah qn-1
= 1eah[e n-2

(eah) nq (3.74)

qn = (eah)n/h q

and when used to calculate the global error results in

E(xn) = Q(xn) qn = Qoeax, q0eax = 0 (3.75)

given the fact the initial conditions for the difference equation are specified so qo = Q0.

Therefore, the exact solution of the difference equation being solved when the scaling factor

is introduced matches the exact solution of the original differential equation. This infers the

accumulation error when the scaling factor is included is minimized and is further reason to

use the scaling factor in the shock-induced combustion computations to be performed here. Current Application of Scaling Factor Scheme

In the current study the two scaling methods described above are compared by

applying each to the solution of the shock-induced combustion flow. It is obvious by this

point that the error is introduced in the solution of the problem through the representation

of the spatial derivative. This is due to the fact that the mass fractions of the species in many

cases vary exponentially over a region much smaller than the typical grid spacing. Therefore,

the source term for a given cell calculated based on the cell average of the dependent

variables is not equal to the cell average of the true source term. Given the use of the scaling

schemes, it is worth assessing how close to the governing equations is the model problem

used by Sussman (1993).

The species conservation equation can be reduced to

__ + U = c (3.76)
at ax
using the continuity equation. In the actual governing equations the source term 6) is a

function of all species involved in the reactions. Therefore the factor k in Sussman's model

problem (Sussman, 1993) is an effective rate of production of a per concentration of a, or

the derivative of the source term with respect to the species ac. This information is available

in the form of the source term Jacobian, given in the Appendix. Alternate estimations for k

may prove useful but here we use the form as proposed by Sussman (1993) to facilitate

comparison to these earlier results.

In the documented scheme by Sussman (1993), the criterion for selecting the upwind

direction in specifying the ratio of mass fractions, Eqn. (3.58), uses the Cartesian velocity

components. It is not clear if this has been done simply for notation but would be very

erroneous if Cartesian velocities are used in conjunction with a curvilinear grid. Therefore,

here the contravariant velocities are used to determine the upwind directions. Also, the

particular form of Eqn. (3.57) used by Sussman (1993) in the two-dimensional computations

determines the error estimate using a weighted sum from the two curvilinear directions.

However, the Cartesian velocities are used in the weighting which may introduce some

inconsistencies. Therefore, the current study evaluates the scaling by Shang et al. (1995) and

Sussman (1993), with curvilinear upwinding, and also evaluates the scaling factor in the


mw 2ri + r) + r n 2rj + r2)

rTw ln(ri) finnln(rj)
2-,j i~j +
mw = mn =
S.2 2 2 2
m_,j +Im m + m
+ ij +- ij+

where the components are as defined earlier and the flow scenario is that of Figure 3.1. The

schemes of Shang et al. (1995), Sussman (1993) and the current scaling factor, Eqn. (3.77),

based on Sussman's work (1993) are evaluated here and the results are presented later. When

these scaling are used, the effective reaction rate is not just the change in the mass fraction

over the distance between neighboring cells but a velocity scaling is included as detailed in

the derivation using the model problem.

As is evident by the attention to both adaptive gridding techniques, transformed

equations, and modified source term treatments, the additional scales introduced by the

chemical reaction processes in the current problem set can not be ignored. Various studies

have evaluated these three options for treating the problem and here the modified source term

treatment is tested. This is done in part because if such a treatment proves worthwhile, it is

easily integrated into the solution process. The scaling factor offers an attractive option to

grid adaptation but the error estimation required is only as good as the estimate of the

effective rate k. Given the behavior of the error as kAx/U increases, any errors in this estimate

are heightened for large k.

3.6.4 Nonphysical Wave Speeds

LeVeque and Yee (1990) also document the problems associated with highly

nonlinear source terms. In their study a model problem is used to study the behavior of the

linear advection equation with a source term. Nonphysical wave speeds are noted for large

At. The source of the nonphysical speeds is found to be solely a function of the presence of

the source term since the conservative numerical methods used correctly reproduce the wave

speed if no source term is present. LeVeque and Yee's (1990) analysis is summarized here.

Consider the governing equation of the current reactive flow problems. If no source term is

present, the numerical scheme is conservative if the summation over the domain results in

the cancellation of the defined fluxes. In a word, for a finite volume scheme, the defined

fluxes at the cell faces must balance. The addition of a source term results in a quantity that

does not cancel when summed over the domain and therefore, if this term is incorrectly

modeled, the accuracy in characteristics such as wave speed are reduced.

In the numerical model of both the model problem studied by LeVeque and Yee

(1990) and the full set of governing equations to be solved for reacting flows, the source term

is typically represented as a function of the cell averaged a, w(Ui). In the context of the

integral form of the governing equations which is used to validate conservatism, this replaces

the average value of the source term, ii, with a source term which is calculated base on an

average a. The discrepancy between these two is small for a smoothly varying a but in the

case of reacting flows where sharp variations exist, the difference can be large. The

conclusions of LeVeque and Yee (1990) on how to treat the occurrence are much the same

as those already presented. First is the option of refining the grid which is not appealing as

mentioned earlier. A second option proposed by LeVeque and Yee (1990) is to better

represent the integral of the source term across a given cell. This requires a better model for

the variation across the cell which is in the same line as the logarithmic variation

approximations evaluated here.

3.7 Numerical Solution Methods

Many options exist for the integration and numerical solution of the governing

equations. These options typically fall in the categories of explicit or implicit schemes.

Problems including chemical reaction source terms as here normally require an implicit

method for the treatment of the stiff source terms. This can be achieved in two ways. The

first option is to solve the governing equations in a fully coupled manner which requires the

derivatives of both the fluid dynamic and source terms with respect to the dependent variable

vector, Q. These derivatives, Jacobian matrices, can be quite large in the case of a

multi-species flow. The second option is a splitting method in which the fluid dynamic terms

are treated explicitly and only the source terms are treated in an implicit manner. Both of

these schemes are to be used in the current study so a description of each is provided here.

3.7.1 Fully-Implicit Formulation

As detailed in the description of the flux scheme, a second-order spatial accurate

model is used. Likewise, it is advantageous to use a second-order temporal scheme as well.

A detailed review of the fundamentals of achieving such a scheme is available in Hirsch

(1990a, 1990b) and the current implementation follows that of Shuen (1992). A concise

review is provided here.

The governing equation for the current problem, Eqn. (3.2), is of the form

d = M(Q, t) (3.78)

where M is a non-linear function of Q and represents the differential space operators and
source terms of Eqn. (3.2). A one-step scheme for solving the system of equations is

Qn+l Qn = 8t [{Mn+' + (1 )Mn] (3.79)

and when 1 = 1/2 the scheme is second order in time. Since the whole purpose of solving

the set of equations is to determine Qn+ the question is first how to represent

Mn +1 = M(Qn +1). One option is to use a predictor and corrector approach which results

in an explicit scheme. The second option is to linearize the non-linear term M +1 which
results in a fully implicit scheme. The term is linearized as

Mn+ = Mn + ) t + O(6t2)
QJ\ n 6t
= Mn + at) at + 0(8t2) (3.80)

= M + Nn 6Q" + O(t2)

where 6Qn = Qn+1 Qn and the Jacobian, N = aM/aQ, has been introduced.
Introducing the linearized terms the equation becomes
6Q = 6t [~(Mn + Nn"Q) + (1 )M"]
= 6t [(Nn8Q) + Mn] (3.81)

Using At = 6t/J and the definitions

F = J(F, + Fv)
G = J(G + G) (3.82)
H = J(H, + Hv)
S = JQ
the governing equation, Eqn. (3.2) becomes

I + '-At(8A + B + H' S')1 AQ= -At (tF + G+ H-S) (3.83)

where 6N and 68 represent the spatial discretization in the two curvilinear directions,

A = aF/aQ, B = 8G/cQ, H' = aH/aQ, and S' = aS/aQ. Note, has been set to 1/2

for second order accuracy. When implemented, the Jacobians of the viscous terms are

neglected which has proven to save computational time with no negative impact on the

solution process.

Introducing the flux vector splitting described previously, Eqn. (3.83) can be

symbolically expressed as

[D +L + U] AQ = At R (3.84)

where the components of the implicit operator are
D1 + (A+t -A- +B+ -B +H' -)
D = I + IAt A j A j + B j B-i,j + Hi,j S i,j

L At (Ai-1j+ Bij-1) (3.85)

U = At A-i+l.j + B-i,j+l)

and the explicit operator is

R = At (Fi+ F.,. + G. ,L G.. + Hi )Si (3.86)

Due to the linearization, all the components of the implicit and explicit operator are

evaluated at time level n. As detailed by Janus (1989), when solving unsteady flows, many

times sub-iterations are used in the solution process to minimize the linearization error.

The Jacobian matrices of the split fluxes are exactly A = 8F: /Q but are

approximated here as in previous studies (Shuen, 1992) using

A: = TAT-' (3.87)

The matrix A* is a diagonal matrix composed of the eigenvalues X* given in Eqn. (3.23).

The matrix T is composed of the right eigenvectors of the matrix A, in columns, and T-

the corresponding left eignevectors in rows. The eigenvectors are provided in the Appendix.

The exact Jacobians are available in Janus et al. (1992).

Here the system of equations represented in Eqn. (3.84) is solved using a lower-upper

(LU) approximate factorization scheme described by Shuen (1992) and Shuen and Yoon

(1989) and based on the scheme developed by Jameson and Yoon (1987). The left-hand side

of Eqn. (3.84) is approximately factored into the product of two operators

[D + L]D-'[D + U] AQ = At R (3.88)

and when implemented can be solved in the sequence

[D + L]AQ* = At R,
[D + L]AQ = DAQ* (3.89)
Qn+1 = Q" + AQ

Various stability advantages of this solution scheme are detailed by Shuen (1992).

3.7.2 Explicit-Implicit Splitting Scheme

As mentioned, options for the solution method for the current problem range from

fully implicit to fully explicit schemes. Implicit schemes are essential when solving

problems with strong source terms as is the case here. However, explicit methods are

attractive due to their reduced computational requirements. Also, there is some question as

to the advantage of fully implicit schemes for this class of problems. For high-speed reacting

flows, small temporal steps are used to either capture the unsteady nature of the flow or in

light of the small time scales of the chemical processes.

The need for an implicit scheme for the treatment of the source terms and the

efficiency of an explicit formulation can be met using the Strang splitting method described

by Le Veque and Yee (1990). The scheme computes the dependent variables (Q) at the new

time level (n+l) by performing various operations on the variables at the current time level

(n). The splitting can be represented by

Qn+I = SQ(At/2) SF(At) Su(At/2) Qn (3.90)

where SQ represents the numerical solution operator for the source term and SF the

numerical solution operator for the fluid dynamics conservation equations. The Strang

splitting has two features which make it very attractive for the current problem. First the

splitting maintains second-order accuracy which is needed to reduce the numerical

dissipation of the important phenomena. The second feature is the fact that the splitting

allows proven solution schemes for each aspect, the chemistry and fluid dynamics, to be

easily integrated. Such a splitting may raise some concerns about capturing the coupling
between the fluid dynamics and the chemical processes. However, as demonstrated in the

combustion instability simulation work of Shyy et al (1993), this approach yields a very

satisfactory treatment for the source terms.

As mentioned, in the splitting technique the fluid dynamics aspects of the problem

are modeled using explicit schemes. To maintain second-order accuracy, the fluid dynamics

operator must be second-order and here a prediction-correction scheme is used in the form

given by Hirsch (1990b). To clarify the steps in the solution process, the sequence of

equations are as follows,

S,(At/2) : I- Q' AQ* =-At

Q* = Q + AQ*

Sp(At) : AQ" = 8[aF + aG)' + H]

Q** = Q* + AQ*
AQ*** = At[aF(2) + G(2) + H]

Q*** = Q* + AQ***

Sn(At/2) : I -Q' AQ****= At***
Qn+l = Q*** +AQ**

where the superscript *'s denote what value of Q is used to calculate the fluxes, source terms,

and Jacobian of the source term, Q'(Q"i,j) = aQ(Qnij)/Qni, j. The superscripts (1) and (2)

denote the spatial order of the numerical fluxes. The term a F = Fi+. F. and
S 1+ .1 1 2-'j
likewise for 9,G. The scheme used to define the inviscid numerical fluxes is the same as that

used in the fully-implicit method. However, here the Jacobians of the fluid dynamic terms

do not have to be calculated which reduces computational time as does the fact that a banded
system of equations does not have to be solved.

3.7.3 ODE Solution Issues

Even when implicit integration is used difficulties can arise when implementing the
solution process for problems with large reaction rates. Here some of these issues are
discussed in the context of solving the ODE associated with the source term operator in the
splitting solution method. One difficulty is due to the fact that the Jacobian matrix Q' can
become ill-conditioned which can introduce errors when the equations are solved. The
Jacobian matrix can be changed somewhat to help alleviate this problem. Consider the two
step reaction
H2 + O2 20H (1)
H2 + 20H 3 2H20 (2) (3.92)

which is used later in this study and has been presented in Table 2.1. The source term for each
species, when expanded, is

)H = MHj kf,(XH)(Xo + kb(XoH)2 kXH)(XOH)2 + kb XHo)2

o, = Mo{- kf, (X,)(Xo) + kb(XoH)2

OH = 2MoH kf,(XH2)(Xo) k OH kb(XH(XoOH)) + kb(XHo)} (3.93)

HO = MH20 2k,(2XH)(XoH)2 2kb2(XHo)2

and after some simple arithmetic those for H2 and OH can be written in the form
MH6. MH [b0

2MOH [o2] OH (3.94)
2"OH = M 2[ H20
02 HO
The Jacobian matrix of the source term can be altered to reflect the relationships of Eqn.
(3.94). This implementation has been used for reaction mechanism composed of two and

eight steps and has proven to increase the robustness of the computational code in the case

of the two-step mechanism. This is due to the fact that the reaction rate constants for the

two-step mechanism are orders of magnitude larger than those for the higher order

mechanisms and result in a Jacobian matrix with a much larger condition number. For the

higher order mechanisms, the Jacobian matrix using the standard form is well conditioned.

Even if the Jacobian matrix is well conditioned and a solution is obtained, the

resulting state of Q can be nonphysical. This is typically manifested in the mass fractions of

various species being less than zero. These nonphysical solutions are a result, once again,

of the fact that quantities such as mass fraction vary in a logarithmic fashion. Because the

solution process utilizes a linearization for the derivative of Q, these overshoots are

produced. In the solution of problems such as the reacting flows currently studied, such

non-physical solutions must be guarded against. This can be done by either taking smaller

time steps, which is not desirable, or by employing a scaling technique such as the method

demonstrated by Ramshaw and Chang (1995). Here a technique is used which scales AQ if

the solution of the ODE in Eqn (3.91) would result in a non-physical Q. In the current code,

the solutions from the ODEs of Eqn. (3.91) are evaluated to see if Q* or Qn+ 1 contain species

mass fractions such that ai < 0 for all i. If so, a single appropriate scaling factor is calculated

and applied to all AQ's such that 0 < ai < 1 for all i. The change in Q due to the fluid dynamic

operator presents no problem. This is due to the fact that the At is calculated based on the

CFL < 1 condition necessary for explicit schemes.

When applying such a scaling, caution must be used to ensure the actual equation

being solved is not altered. In the case of steady-state solutions, typically the scaling is

needed only in the early stages of the evolution of the flow field. If the scaling factor

converges to 1 long before the convergence of the governing equations to the steady-state

solution, then there are no concerns that the equations have been altered. In the case of

unsteady problems, an iterative process can be used to solve the ODE to ensure non-physical

solutions are not encountered. However, typically when solving unsteady problems, the time

step is set small to resolve the evolution of the flow field and such problems are not

encountered. It should be noted that these overshoots are different than the round-off errors

which can appear as mentioned earlier and as addressed in previous investigations (Sussman,

1993 and Wilson and Sussman, 1993).

3.8 Summary

Here the various models needed to successfully simulate the problems of interest

have been described. As is clear, there are numerous phenomenological and numerical

models available. Those selected to be used in the current development offer the latest in

methodology and have a wide range of applicability. This is important in the context of the

current study given the fact that a multitude of scales are intertwined. Those simplification

that have been made have been done so after analyzing the phenomenology of the process

and based on previous investigations. From the numerical aspect, one of the key items to be

highlighted is the source term treatment. This is due to the fact that previous investigations

have shown the accuracy is highly dependent on the grid resolution. The source term

treatments offer a viable alternative to the computationally intensive grid adaptation


The primary chemical aspect to be focused on is the impact of the reaction

mechanism. The model described here is used to compute the shock-induced combustion

problem with various mechanisms in a effort to correlate mechanism characteristics to

simulation accuracy. The impact of turbulence in the reactive flow scenarios presented here

is also quantified. Both the modeling of the turbulence and the coupling of these effects with

the chemical processes are detailed in the next chapter. The turbulence effects are known to

be manifested in the governing equations through the viscous aspects. It is clear from the

description of the governing equations that these aspects impact field quantities but also the

chemical process through the dispersion of species and heat.


The turbulence aspects of the problems of interest like the other physical and

chemical aspects require a phenomenological model. As eluded to earlier, turbulence

impacts many of the aspects already discussed as well as the particular models that have been

described. Therefore, this aspect of the problem is dedicated its own chapter to emphasize

both the phenomenological importance of turbulence but also the attention which must be

paid to its modeling.

A turbulent flow, like all continuum fluid problems, is governed by the

Navier-Stokes equations. But to model these equations on the length and time scales to

resolve all pertinent physical phenomena is not a practical method given today's

computational platforms. Simulations which do resolve all the scales ranging from the small

to the large eddies are conducted and are labeled Direct Numerical Simulations (DNS). A

review of DNS is provided by Galperin and Orszag (1993). To highlight the computational

resources required to conduct a DNS let us estimate the number of grid points that would

be required for the fuel cell scenario.

As mentioned, the phenomena of turbulence is known to be driven by a multitude

of scales, particularly the large scale eddies which contain what is call the turbulent kinetic

energy and the small scale eddies which work to dissipate this energy. The dynamics of

turbulence consist of various transformation processes which transfer this energy between

the eddies. To capture all the important phenomena, the governing equations must be

discretized spatially such that all scales are resolved. This infers that the computational grid

used must contain points spaced on the order of the size of the small eddies. The

characteristic length scale of the small eddies is the Kolmogoroff scale (1) and the ratio of

this scale to the characteristic length scale of the large eddies (le) is related to the large eddy

Reynolds number by Yl/le = (Re ) where Rel = -. Therefore, to capture all scales of

interest, the grid spacing must be on the order of I] = le(v3/u31e3 ). Given the fact that these

small eddies can have an arbitrary orientation in the flow, this same spacing must be used

in all three physical dimensions. An estimate of le is needed and for the fuel cell scenario,

the large scale eddies are on the order of the thickness of the various shear layers which mix

the fuel and oxidizer. These shear layers are on the order of the size of the projectile so this

length scale can be estimated as le = O(d) where d is the projectile diameter. Using the

Reynolds number of the mean flow (Re = ud/v) we can estimate the Kolmogoroff scale to

be 1 = 0 d(Re) ]. Recalling that the discretization must be on this order in all directions

and considering the computational domain, at the least, must be on the order of the projectile

diameter, the number of points to carry out a DNS study of this problem is on the order of
(Re)4. Using quantities for the viscosity, length, and velocities comparable to those in the

problems to be solved here, we have Re = O(1.0x107) which would require on the order

of 5x1015 spatial points on which the governing equations must be solved. This is clearly

not practical for typical engineering analysis or much less easy given the current

computational platforms. However, DNS is pursued for moderate Reynold's numbers and

simple geometries in hopes of better understanding the dynamics of turbulence thereby

advancing optional modeling techniques. These optional techniques are required to

complement the solved set of governing equations in an attempt to reproduce the dynamics

of those scales known to not be adequately resolved.

The temporal scale of the small eddies must also be considered if a DNS is to be

performed. This time scale can be thought of as the turnover time of the small eddies. Like

the spatial requirements impact the number of points at which the governing equations must

be solved and in turn the amount of computer memory needed, this requirement on the

temporal discretization effects the required computational time as does the small time scales

of the chemical reactions.

The next level from DNS in the context of grid requirements is Large Eddy

Simulation (LES) in which grid refinements are used such that the large scales of turbulence

are resolved but a model is needed for the small scale phenomena. LES is also characterized

by the use of spatial filters on the governing equations. This terminology is used because the

computational grid has to be constructed such that the spatial discretization is on the order

of the inertial subrange, or that region which lies between the large and small eddies and over

which the primary physical process is the transfer of energy. Since the primary role of the

small eddies is known to be one of energy dissipation, this is the phenomena which must be

modeled. A review of LES is provided by Erlebacher et al. (1992) and Germano (1992).

In the context of solving a full scale engineering problem, the term large eddy is

somewhat misleading in that the size of the large scales in turbulence can still be orders of

magnitude smaller than the domain over which the solution must be carried out. Therefore,

the LES requirement to resolve only the large scales is still computationally restrictive.

However, LES is currently being pursued aggressively and has been used for many problems

of interest. LES will continue to become more viable as computational platforms become

faster and cheaper.

Based on the description of DNS and LES, it is clear that an optional modeling

technique is needed for problems composed of large Reynolds numbers and large

computational domains. The problem in modeling turbulent flows is that the highly

nonlinear dynamics, represented by the Navier-Stokes equations, occurring at these

unresolved scales produce phenomena when observed from a larger frame of reference

appear random and non-systematic. The model of choice for problems such as the reactive

munition and shock-induced combustion scenarios remains a class of engineering models

labeled moment methods which are solved in conjunction with a set of averaged governing

equations. There are many options as to both how to average the equations as well as the form

of the engineering models. More details on various engineering models is provided by

Wilcox (1992) and a review in the context of reacting flows is provided by Libby and

Williams (1994). Here the work of Krishnamurty (1996) is followed to describe the

implementation of one particular model in an attempt to quantify the role of turbulence in

the reacting flow problems of interest.

4.1 Averaging Techniques and Averaged Variables

As mentioned, when a turbulent flow is observed, an apparent randomness exists.

Conceding the fact that all scales are not to be resolved, it is assumed that a solution on the

resolvable scales can provide a statistical average of the exact flow field. In LES a spatial

averaging, or filtering, procedure is used. Here, however, the averaging is cast as either a

time or ensemble average. A time average of a quantity with instantaneous value (d is

to +T

= lim j dt (4.1)
T-.00 T

where T is the period over which the instantaneous values are sampled. This results in a

formula for the instantaneous value at a given time of the form

+ = + (4.2)

where (' is the deviation from the mean and by definition 4' = 0.

When solving unsteady flows, it is more applicable to employ a methodology of

ensemble averaging. The ensemble average of a quantity c( is an average over a number of

realizations of the problem and mathematically is

= lim Z n (4.3)
N--co N

where again (' is the deviation from the mean and by definition 0' = 0. A realization is

simply the same problem restarted from the same initial conditions with the quantity 4 being

measured at the same point in space and at the same time after the start of the problem.

The mental picture of this type of averaging is an excellent way to frame the problem

of turbulence. Consider the task of collecting an ensemble average of say the velocity at a

specific point in the flow field about a model in a wind tunnel at a specific time after the

tunnel is started. Assuming the measurement can be made precisely at the same point in time

and space, the variation from measurement to measurement is a function of the nonlinearity

of the governing equations. Small perturbations can be generated in any quantity whether

it is caused by the startup process, a slight movement of the model, a small contaminant on

the tunnel wall, or a miniscule change in the environment. These small changes can grow

within the framework of the Navier-Stokes equations and effect the field quantities at various

points in space and time. And as mentioned several times now, to capture this evolution with

a computational model is no easy task.

Continuing with the premise that when the Navier-Stokes equations are solved on a

grid known to introduce some statistical averaging, the equations which are actually solved

need to be cast into a form accounting for this averaging process. The first step in doing so

is defining the averaged variables based on the averaging techniques discussed earlier. When

solving turbulent flows, averaged variables in the form of Eqn. (4.2) are typically defined

as Reynolds averaged variables. With the added desire to model compressible flows, a

second averaged variable definition is introduced. This is in light of the fact that unlike

incompressible flows, the added aspect of compressibility can introduce deviations, or

fluctuations, in density. By introducing a mass weighted average, denoted as a Favre

average, the resulting equations are somewhat simpler. Using this form, the instantaneous

value is given by

4)= 4) + )" (4.4)


which results in the relationships

= = 0. (4.6)

Note, the tilde and double primes are used to distinguish Favre averaged variables from

Reynolds averaged ones. One reason this approach is taken in modeling compressible

turbulent flows is that with the density weighting introduced, the final form of the governing

equations is closely analogous to the incompressible form and it is in the context of

incompressible turbulence that most models have been developed.

Following the typical convention, a mix of Reynolds and Favre averaged variables

are used. Specifically, Reynolds averaging is used for the density and pressure and Favre

averaging for the velocity components and scalars such as temperature and species mass

fractions. The forms are

Q = +Q'
p P + p'
ui = uf + u," (4.7)
T = T+ T"
ai = di + a"ill

The next section presents the form of the governing equations once these averaged variables

are introduced.

4.2 Averaged Equations

A review of the averaging process is provided by Wilcox (1992). Here the final form

of the equations are presented. This final form is obtained by introducing the defined average

variables of Eqn. (4.7) into the governing equations of Chapter 3 and making use of the

relationships given above. Note, for convenience, the Cartesian form of the equations with

tensor notation is used in the discussion. The continuity equation takes the form

i() + )= 0 (4.8)

and the momentum

(ugiJ) + (ax u) = + x i u"u ") (4.9)
t i -xi axj ij


i v vi' + i I +X-1 X 6
axx axi axk (4.10)

[au." u." [ "1 (4.10)
Tj"=v a- + +- + x- 6,uk"
j -ax -ax- "- ixk

and the term u, "uj" is yet to be defined. The energy equation, with the introduction of total

enthalpy, H = E + and specific enthalpy, h = e + becomes

)+ a-x( ) = "i + Tij Qu"uj") + ui"ti + ui Tij
qj quj"h" Quj'(,u,"u,")

Recall for reacting flows, the heat flux is also a function of species diffusion velocities as
given in Eqn. (3.11). Also, as a result of the averaging process, the the total energy, E, now

includes the turbulent kinetic energy which is defined later. The final equation is the species

conservation equation and when averaged for a single species n is

(n) + nj c ax ean"uj = (4.12)

where the source term is denoted as an average and the particular representation is discussed

The first terms to be addressed are those that are dependent on the product of the

molecular viscosity and the fluctuating terms. This includes Tij" in the momentum equation,

ui "ij and qj" in the energy equation, and -Sc I in the species conservation equation.

These terms are small in magnitude when compared to the terms that are a function of the
product between molecular viscosity and the mean field quantities. Therefore, these terms
are ignored in the modeling. The remaining terms are significant and must be modeled.

4.3 Turbulence Closure

The first term in need of modeling, or closure, is Qu "uj" which is a symmetric tensor

with six independent components. Analysis shows that this term behaves much like the

molecular stress term and is labeled the Reynolds-stress tensor. Boussinesq first proposed

modeling this term in an analogous fashion, dependent on the gradients of the mean field.

This requires an effective viscosity dependent on the turbulence, and the estimation of this

turbulent viscosity varies based on the model employed. The first level in sophistication is

the mixing-length model proposed by Prandtl (Tennekes and Lumley, 1972) which uses the

gradient of the mean velocity field and a yet undetermined mixing length to define the

viscosity. In a shear layer flow, the form for the turbulent viscosity is

tt = Qmix (4.13)

where u is the streamwise component of velocity and y is normal to the flow direction. This

model performs well for some flows in which the mixing length, Imix, has been empirically

determined. However, this model offers limited potential for complex flow fields in which

a single length scale does not hold. A more appropriate model for such flows is one which

incorporates the dynamics of turbulence in the estimation of turbulent viscosity. This

requires the definition of a velocity scale and length scale.

Various models have been used and a review is provided by Wilcox (1992). Here, one

particular model is discussed and integrated with the other models required for the

high-speed reactive flow problems. The model used here is the k-e model, a two-equation

model with its foundation in the work of Jones and Launder (1973). The exact equations used

to govern these quantities are presented next but first the closure of the various terms in the

averaged equations are discussed. The quantity k is the turbulent kinetic energy mentioned

earlier and e is the dissipation of k and both are defined as

k ui"u," 3ui," u,"'
k = E = v aX (4.14)
Q axi ax

Using k and E to define the velocity and length scale of turbulence, the turbulent viscosity

kt = Cg- (4.15)
where C. is a constant found to be .09. The Reynolds stress term is modeled as

r[aiu. 8a'xi [ik 1
Quu = t + axL k ij 2 (4.16)

with the last term introduced to maintain the proper trace of the Reynolds stress tensor. The

turbulent bulk viscosity is represented as is the molecular, Xt = -- t.

The next term in need of closure is Quj"h" which is modeled, using the Reynolds'

analogy (Hinze, 1959), as
CpIt aT(
Quj"h" Prt axj (4.17)

where Prt is the turbulent Prandtl number. Using this closure, a single turbulent Prandtl
number, assumes the turbulent transport of heat and momentum are dependent on a single
time scale, i.e. the ratio between the turbulent heat and momentum flux are the same
throughout the flow field. This assumption is known to not hold in many flow scenarios and
an optional modeling concept such as the work of Nagano and Shimada (1996) is to define
a second velocity and length scale associated with the heat flux. This allows the thermal eddy
diffusivity to be calculated explicitly throughout the flowfield. For the current study which
is an initial investigation into the coupling of the turbulence and chemical reacting aspects
of the flow, the closure of Eqn. (4.17) is deemed acceptable. In most flows, the turbulent
Prandtl number is close to unity and that is the value used in the current study. The remaining

term in the energy equation, Eqn. (4.11), is ui'"Tij" Quj"(ui"'ui") which represents the

diffusion of energy due to turbulent fluctuations and is modeled as
U7 -t ak
"" "(2 u,"u") = k (4.18)
Ui Ci -Qj u kaxj

where ok is a Prandtl number associated with the turbulent kinetic energy and is set to unity.

The one term in the species conservation equation needing closure is QCn"uj".

Using the gradient treatment already introduced for the turbulent heat flux, the term is

modeled as

Qcg"u" -- It aan (4.19)
Sct axj

where the turbulent Schmidt number, Sct, has been introduced. As in the case of the turbulent

Prandtl number, assuming a single turbulent Schmidt number prescribes that the scales

associated with the turbulent transport of species are the same throughout the flow relative

to those of the momentum. Also, similar to the modeling option for the turbulent heat flux,

additional equations can be solved for the correlations of Eqn. (4.19) but an additional

equation is needed for each species which can be computationally expensive. Here, Sct=1

which also implies the turbulent Lewis number is 1. All the modeled terms are now expressed

as a function of k and E so the next issue is the modeling of these two quantities.

4.4 k-E Model

The equation for k is derived by manipulating the momentum equation. A description

is given in Krishnamurty (1996) but here only the final form of the equation is given,

=- p --u1 --
(-xk) + -\li)jk -. Quif"uj" ax uap aui
t axI ax i ax axu "
_,, (4.20)
+ a ui "ij Qu tu"ui"i) p'j] + p aX

where some terms needing closure have already been defined. The definitions for the

production of turbulent kinetic energy and the rate of dissipation of this same energy are such


Pk = Qui uj a

uax (2
^E = ~-axp"

Note, the term p'uj" is not explicitly modeled and is included in the closure of Eqn. (4.18).

Bu." ap
The terms p' and u" are, respectively, the pressure dilatation term and the

enthalpic production term. A model has been proposed by Krishnamurty (1996) for these

terms and has been evaluated for nonreacting blunt body flows with shocks. The work of

Krishnamurty (1996) shows these terms may prove important in flows where

compressibility effects are significant. The shock-induced combustion flow studied here is

a high-speed, compressible flow and it may prove necessary to model these terms. However,

for this initial investigation, the standard form of the k-E model is used. This is due in part

to establish a benchmark using the baseline model but also in light of the fact that the

derivation of the models proposed by Krishnamurty (1996) need to be investigated in the

context of multi-species reacting flows. There have also been several modifications

proposed for addressing known compressibility effects and how they impact factors such as

the dissipation rate of k. These modifications are also discussed in detail by Krishnamurty

(1996). Here, the final form of the k equation used is

S(k)t +a (ijk) = Pk +E + V ( + ~) (4.22)

The governing equation for E is not as cleanly derived from the Navier-Stokes

equations but is based on the assumption the time rate of change and the convection of the

dissipation (E) is balanced by a production, dissipation and diffusion term in the form

a-- + ijE) = PE E + Te + B (4.23)

Unlike the k equation which is derived from manipulation of the Navier-Stokes equations,

the terms in the E equation have been derived based on dimensional analysis and physical

reasoning. The terms in the standard model are

PE = E, = E2 Te t + (4.24)

where Oj is a Prandtl number associated with the E and is set to 1.3. The coefficients in the
production and dissipation terms can be modified to account for the imbalance between the
production and dissipation of turbulent kinetic energy. These modeling option are discussed
by Chen and Kim (1987) for C.1 and Thakur et al. (1996) for CE2. Here, the standard values
of CE1 = 1.44 C,2 = 1.92 are used.

The term BE represents the baroclinic torque term which arises due to gradients in

density and pressure, () x (Vp) which do not appear in incompressible problems. A

model for this term, which may prove important to flows such as the current shock-induced
combustion scenario, is proposed by Krishnamurty (1996). Here, though the standard form
of the equation,

a + a ) = (Ck p -CE2) + J + (4.25)

is solved. Now that the equations for k and E have been presented, the next section specifies
the needed boundary conditions.
4.5 Wall Function
The k-E model used here requires a wall treatment due to the low Reynolds number
associated with this region. Here, the standard wall function (Wilcox, 1992) is used. This
function is based on the known correlation between the nondimensional velocity near the
wall and the nondimensional distance from the wall. This velocity is specified as

U + = Up/UT where UT = wal/ and the nondimensional height is y + = (QUT yp)/t

where yp is the distance normal to the wall and [t is the molecular viscosity. In the viscous
sublayer, the region nearest the wall, the relationship

U+ = y+ (4.26)
holds. Outside this region, the relationship between the two nondimensional quantities is

U = ln(y+) + B (4.27)

where x is the von Kdrmdn constant set to .41 and B=5.1. Therefore, this region is labeled

the log layer. Where the curves defined by Eqn. (4.26) and Eqn. (4.27) meet is defined yc,

the critical value of the nondimensional height.

In a finite volume scheme, such as the one developed here, the flow properties are

defined at the cell centers. Therefore, when implementing the wall function, the y +

corresponding to each cell which bounds the wall must be calculated based on the half-height

of the cell. Using this value along with velocity and viscosity, the boundary values for k and

E are calculated, if y + < y+, using (Krishnamurty, 1996)

U2 2 1 3
wall = wall wa (4.28)
c+ w


l iY- cp (C R Q wail YP
S' ci = x Re = (4.29)
1 + I, It

If y + > yc+ then the boundary conditions are specified as

U2 U3
wall Ewall =- (4.30)

Now that the equations for k and E have been specified as well as the boundary conditions

for both quantities at the wall, the final information needed is the boundary condition for k

and E at the inflow. This is typically specified in terms of a turbulent intensity, the magnitude

of the velocity fluctuations relative to the mean value. Here a value of 3% is used base on

data for similar flow scenarios used in earlier studies (Krishnamurty, 1996).

The governing equations for k and E are solved in conjunction with the Navier-Stokes

equations. These equations provide the necessary information with which to define pt. With

the addition of the turbulent viscosity, an estimate of the enhanced diffusion of species and

heat is obtained. This alone can impact the mean flow field and the impact this has on the

shock-induced combustion scenario is studied. In addition to these effects, the turbulence

aspects of the flow can also directly impact the reaction source term. This is the focus of the

next section.

4.6 Turbulence Effects On The Reaction Source Term

The chemical reaction source term, On, in Eqn. (4.12) is a function of temperature

and species mass fractions. The chemical source term is a function of the mean values of

these quantities and any deviations. The effect of turbulence on the source term is first

represented by modeling the impact turbulence has on the mean temperature and mass

fraction field. This is done in the framework of the two-equation type model by applying the

averaging and closures described above. The next step in modeling the impact of turbulence

on the source term is to account for the effect of the fluctuations. This is needed given the

highly nonlinear dependence of the source term on the dependent variables. The analysis of

Sagara and Tsug6 (1970) shows that the fluctuation in velocity does not effect the mean

reaction rate, the effect of density fluctuations is relatively minor, and the primary effect of

turbulence on the source term is manifested through the temperature fluctuations. They find

the mean reaction rate in a turbulent field to be

o = (1 + A)~exp + (1 A)Texp f l t) (4.31)

where A = 6T/T and F = E/RT. The term 6T = T-'T" is the root-mean-square of the

temperature fluctuation. This effective source term can be written in the form

S= [(1 + A)w(T + 6T) + (1 A) (T 6T) (4.32)

after simple arithmetic. Given the highly nonlinear nature of o, the resulting reaction source

term can be quite larger than simply the order of magnitude of the temperature fluctuation,

especially for large activation energies. This model has been used by Mizobuchi et al. (1997)

to simulate a scramjet combustor scenario. In their computations the factors (1 A)2 are

neglected but here, these are included. Their results show an improvement in the computed

combustion efficiency when compared with experimental data.

To this point in the modeling of the various terms in the averaged equations, the

presence of temperature fluctuations has been solely in the context of the correlation of these

fluctuations with those of velocity. As mentioned in the context of the closure of Eqn. (4.17),

one modeling option is to solve additional equations for the temperature fluctuations scales.

If such a model is used then the estimate of the root-mean-square of the temperature

deviation is explicitly provided. Since we have chosen to specify a single turbulent Prandtl

number, these deviations are to be estimated based on the deviations in the velocity field.

This is considered appropriate given the goal of the current study is to determine if any

estimate for the effect of turbulence on the temperature field impacts the computed flow field

in the problems of interest.

Using the strong Reynolds analogy (Gaviglio et al., 1977), the ratio of the

root-mean-square of the temperature fluctuations to the mean value is related to the velocity

field as

8T 2
S= (y 1)M2 (4.33)

where U denotes the velocity magnitude and k is used to estimate the root-mean-square of

the velocity fluctuations.

4.7 Summary

The model employed to represent the turbulent aspects of the flows has been

described. The two-equation model has been selected due to its utility in solving full-scale

engineering problems. The various modifications to the k-E for compressibility effects and

the imbalance between the production and dissipation of the turbulent kinetic energy have

been presented. In the initial computations, the standard form of the equations is used to

produce a benchmark to which further computations using these modifications can be

compared. Turbulence is known to enhance the dispersion of species and heat and these

aspects in turn affect the chemical processes. An additional effect explicitly accounted for

in the present effort is the impact of temperature fluctuations. This is incorporated through


a modified source term dependent on the root-mean-square of the temperature fluctuations.

The impact turbulence has on the problems of interest is detailed in the context of results

using these models presented in the next chapters.


Through the review of the earlier computational studies into the shock-induced

combustion scenario and model development discussion, it is clear that many issues related

to both the numerical and phenomenological aspects of the models used can affect the

outcome of the simulation. For example, it has been quantified in the discussion of the source

term treatments that the numerical error due to the presence of the reaction source term is

dependent on the cell Damkohler number, D)cfd = kAx/u. Using this parameter to frame the

computational investigation into high-speed reacting flows, the phenomenological aspects

of importance is the reaction rate, k, and the local flow speed, u, while the numerical aspect

is the mesh spacing, Ax. To date, the numerical aspect has been the primary focus of

investigations into the shock-induced combustion flow scenario.

In the context of the shock-induced combustion problem, various techniques for

addressing this numerical aspect have been evaluated including grid adaptation (Wilson and

MacCormack, 1990, Sussman, 1993, Shang et al., 1995), source term treatments (Sussman,

1993, Shang et al., 1995), and solving a logarithmic form of the governing equations

(Sussman and Wilson, 1991, Wilson and Sussman, 1993). In these earlier studies a single

reaction mechanism is selected and used to simulate the M=6.46 shock-induced combustion

case captured by Lehr (1972). A main focus of these earlier computations has been to

quantify the effectiveness of particular numerical treatments. Wilson and MacCormack

(1990) do compare the effect of altering the rate of the chain-branching reaction H + 02

OH + O and note the impact the alteration has on the induction time. However, none of the

studies perform a comprehensive comparison over a wide range of reaction mechanisms

which is one focus of the present investigation. A subset of those numerical techniques

developed in these earlier studies is evaluated, particularly the source term treatments

including the logarithmic interpolation proposed by Shang et al. (1995) and the scaling

factor scheme proposed by Sussman (1993). Also evaluated is the explicit-implicit splitting

scheme used to integrate the governing equations.

Given the variety of reaction mechanisms that has been proposed for the

hydrogen-oxygen reaction process, the models used here have been selected based on their

use in earlier computations of the shock-induced combustion scenario. These models are

used to derive key parameters useful in quantifying reaction mechanisms and shown to be

directly related to the computed flow fields. This establishes comparison parameters useful

in future studies into the aspects of computing high-speed reacting flows and the coupling

of the numerical and physical aspects.

To assess the impact the numerical and physical aspects have on the accuracy of such

computations, a starting point must be selected. The steady combustion scenario at M=6.46,

captured by Lehr (1972) is chosen as the case to study. Based on the success of earlier

computations, neither wall nor turbulence seems to have order one effects in the overall flow

field and a logical starting point is to perform the computations with an inviscid model using

piecewise representation of the source term. This is done in part to compare the current

computations with those of earlier investigations.

Based on the review of previous computational studies of Lehr's experiments (1972)

it has been decided to first consider the 7-step mechanism on a moderate grid and with no

special treatment of the source term. The present computational grid is composed of 80x 120

cells with 120 being in the direction normal to the body. The current grid is shown in Figure

5.1. To simplify the issue related to the grid distribution, uniform spacing is used in the

normal direction and the tangential spacing along the body is approximately uniform.

The 7-step hydrogen-oxygen mechanism used in the baseline computation is given

in Table 2.2 and is from the work of Drummond et al. (1986). The fully implicit solution

Figure 5.1. Baseline 80 x 120 grid used for the M=6.46 case.

scheme is used for the computations presented in this section. Also, for all shock-induced

combustion computations, a third-order representation for the dependence of specific heats

on temperature is used. The density contour from this baseline computation is shown in

Figure 5.2 along with the reference points from the experimental shadowgraph presented

earlier. The computation is carried out for only half the field but a mirror image is also

presented to facilitate comparison to the experiments. The first noticeable feature of the

comparison is the asymmetry in the experimental data, especially in the axial locations near

the baseline of the projectile. It should be reiterated that there is also some optical distortion

in the shadowgraph and the reference points have been altered in no way. They have simply

been placed such that the experimental shock location on the stagnation line corresponds to

* Experimental Shock Front
o Experimental Energy-Release Front

Figure 5.2. Solution on the 80 x 120 grid for M=6.46 using the 7-step mechanism
and the piecewise constant representation of the reaction source term.

that of the computations. Also, the data has been scaled 1:1 in length and height until the

diameter of the experimental projectile matches that of the computational geometry. As is

evident, this results in the reference points not extending as far to the base as in the

experimental image of Figure 2.2. The data from the experimental shadowgraph can be

altered to better match the length of the projectile in the computations but this requires

additional assumptions on the optical distortion which are not made here.

The comparison between the computed shock and reaction fronts are quite good in

the nose region but in the region further downstream the reaction front is somewhat closer

to the shock than the experimental data infers. This is one reason, unlike many of the previous

studies, the current study models the complete length of the projectile to provide a more

complete comparison. This thinning trend is also seen in the computations of Shang et al.

(1995) in which the same reaction mechanism is used on a comparable grid. It is noted that

all computations presented here have been converged to the point where the average density

residual is reduced six orders of magnitude.

Before studying other modeling aspects, the dependence of the solution on the grid

resolution used in the computation must be measured. The 80x120 grid has been chosen

based on the results of Shang et al. (1995) and the results from this grid as well as three

alternate grids are shown in Figure 5.3. For comparison purposes, only the experimental data

from the upper half of the shadowgraph in Figure 2.2 is used. To make a more precise

comparison, in Figure 5.4 the variation in density along a line normal to the cylindrical body

is presented. The x location of the lateral variation is one radius downstream of the

hemispherical nose-cylindrical body junction. All lengths are nondimensionalized by the

projectile diameter. A similar comparison along the stagnation line results in lines that are

essentially indistinguishable. The solutions on the coarser grids of 40x90 and 40x120 do

alter the shock and reaction front location slightly, but even then the relative location of the

reaction front to the shock are similar. Therefore, the thinning at the larger x must be due to

factors other than the grid.

Although Figure 5.4 demonstrates that as the grid is refined, the solution seems to

asymptotically converge, caution must be exercised when determining grid dependency for

high-speed reacting flows. This is due to the role of the chemical reaction source term which

is dependent on the distribution of species and temperature. In the context of solving the

multi-dimensional problem, this source term is calculated in each computational cell, or

control volume. The discussion on source term treatments has shown that error introduced

through the representation of the source term in each cell is proportional to kAx/u, the cell

Damk6hler number, where k denotes the reaction rate and Ax the grid spacing. Figure 3.3

Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd