Automated Nozzle Design through Axis-Symmetric Method of Characteristics
Coupled with Chemical Kinetics
by
Reid Britton Young
A thesis submitted to the Graduate Faculty of
Auburn University
in partial ful llment of the
requirements for the Degree of
Master of Science
Auburn, Alabama
August 4, 2012
Keywords: rocket, propulsion, nozzle, design
Copyright 2012 by Reid Britton Young
Approved by
Roy J. Hart eld, Jr., Chair, Woltosz Professor of Aerospace Engineering
Andrew Shelton, Assistant Professor of Aerospace Engineering
Brian Thurow, Reed Associate Professor of Aerospace Engineering
Abstract
A method of designing axis-symmetric nozzles with included  nite-rate chemical kinetics
has been developed. The use of the axis-symmetric method of characteristics coupled with
 nite-rate chemistry allows for the design and performance prediction of nozzles for liquid
rocket engines. The use of a particle swarm / pattern search optimizer enables viable pressure
distributions to initiate the axis-symmetric method of characteristics. The method allows
for the nozzle boundary to be formed through a full method of characteristics solution. The
contour formed allows for a smooth throat, shock free nozzle.
This development allows the process for designing rocket nozzles to be automated, fa-
cilitating the design of axis-symmetric nozzles to meet a speci c performance goal, while
simplifying the e ort of the designer. This thesis will describe the chemical kinetic cou-
pled axis-symmetric method of characteristics, and the associated automated nozzle design
methodology.
ii
Acknowledgments
I would like to especially thank Dr. Roy Hart eld for his support and guidance through-
out my graduate work at Auburn University.
I am also grateful for all of the support of my parents, they have been there through all
of my endeavors.
I would like to dedicate this work to my grandfather, Donald R. Young. \Illigitimi Non
Carborundum," these words will never let me down.
iii
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 A Brief History of Rockets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 The Liquid Rocket Engine . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Thrust Chamber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.2 Thrust Chamber Performance . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Liquid Propellants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Combustion of Liquid Propellants . . . . . . . . . . . . . . . . . . . . . . . . 10
2 Review of Literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1 Contoured Nozzle Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Transonic Flow Zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Method of Characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.1 Axis-symmetric Method of Characteristics . . . . . . . . . . . . . . . 24
2.3.2 Valid Pressure Distribution Selection . . . . . . . . . . . . . . . . . . 29
2.4 Chemical Kinetics with Finite Reaction Rates . . . . . . . . . . . . . . . . . 32
2.4.1 One Dimensional Chemical Kinetics . . . . . . . . . . . . . . . . . . . 35
3 Axis-Symmetric Method of Characteristics Nozzle Design . . . . . . . . . . . . . 39
3.1 Nozzle Contour Design Process . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.1 Comparison of Processes . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.1.2 Converging Section Geometry . . . . . . . . . . . . . . . . . . . . . . 44
iv
3.1.3 Analysis of Transonic Region . . . . . . . . . . . . . . . . . . . . . . 44
3.1.4 Select Number of Charcteristic Points . . . . . . . . . . . . . . . . . . 48
3.1.5 Select Pressure Distribution And Expansion Length . . . . . . . . . . 48
3.1.6 Perform Method of Characteristics . . . . . . . . . . . . . . . . . . . 49
3.1.7 Check Solution for Viability . . . . . . . . . . . . . . . . . . . . . . . 49
3.1.8 Analysis Near Throat . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.1.9 Determining The Boundary . . . . . . . . . . . . . . . . . . . . . . . 55
3.2 Pressure Distribution Optimization . . . . . . . . . . . . . . . . . . . . . . . 58
3.2.1 Fitness Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.3 Uniform Flow Nozzles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.4 Nozzle Truncation and Performance Evaluation . . . . . . . . . . . . . . . . 70
3.4.1 Optimization of Truncated Nozzles . . . . . . . . . . . . . . . . . . . 71
4 AXSMOC Nozzle Design with Chemical Kinetics . . . . . . . . . . . . . . . . . 74
4.1 Combining Chemical Kinetics with AXSMOC . . . . . . . . . . . . . . . . . 74
4.2 Nozzle Contour Design Process with Chemistry Included . . . . . . . . . . . 85
4.2.1 Nozzle Contour Design Process . . . . . . . . . . . . . . . . . . . . . 86
4.2.2 Converging Section Analysis . . . . . . . . . . . . . . . . . . . . . . . 86
4.2.3 Isentropic Nozzle Formation . . . . . . . . . . . . . . . . . . . . . . . 92
4.2.4 Perform Chemically Coupled MOC . . . . . . . . . . . . . . . . . . . 92
4.3 E ects of Finite Rate Chemistry on Nozzle Contour . . . . . . . . . . . . . . 95
5 Conclusions and Suggestions for Future Work . . . . . . . . . . . . . . . . . . . 98
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
A Sauer?s Solution to Transonic Region . . . . . . . . . . . . . . . . . . . . . . . . 104
B Chemical Equilibrium with Applications . . . . . . . . . . . . . . . . . . . . . . 108
C Thermodynamic Properties Database . . . . . . . . . . . . . . . . . . . . . . . . 109
D Reaction Mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
v
E One-dimensional Finite Rate Chemistry Data . . . . . . . . . . . . . . . . . . . 114
E.1 LOX-Hydrogen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
vi
List of Figures
1.1 Flow Zones in a Typical Thrust Chamber[1] . . . . . . . . . . . . . . . . . . . . 6
2.1 Typical Outer Mold Line of a Conical Nozzle . . . . . . . . . . . . . . . . . . . 14
2.2 Typical Outer Mold Line of a Contoured Nozzle . . . . . . . . . . . . . . . . . . 15
2.3 Typical characteristics net in a Rao nozzle with illustration of control surface . 16
2.4 Formation of internal shocks in conical and contoured nozzles[2] . . . . . . . . . 18
2.5 Pressure Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.6 Flow in the throat region of the nozzle with Sauer?s solution to the sonic line. . 23
2.7 Characteristic Point and its Associated Angles[3] . . . . . . . . . . . . . . . . . 26
2.8 Characteristic Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1 Method of characteristics solution zones . . . . . . . . . . . . . . . . . . . . . . 41
3.2 Initial step for both nozzle design approaches . . . . . . . . . . . . . . . . . . . 41
3.3 Establishing the boundary conditions for both approaches . . . . . . . . . . . . 42
3.4 Kernel solution developed for both design approaches . . . . . . . . . . . . . . . 42
3.5 Form method of characteristic solution to  ll the non-simple region . . . . . . . 43
3.6 Form boundary region solution with the method of characteristics . . . . . . . . 43
vii
3.7 Formation of nozzle contour for both approaches . . . . . . . . . . . . . . . . . 44
3.8 Transonic Region of Nozzle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.9 Unde ned Region[3] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.10 Viable Kernel solution, MOC mesh contains no crossing characteristics . . . . . 50
3.11 Non-viable method of characteristics solution, crossing like characteristics in
MOC mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.12 Close inspection of a non-viable method of characteristics solution, crossing like
characteristics in MOC mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.13 Speci cation of the unde ned region . . . . . . . . . . . . . . . . . . . . . . . . 53
3.14 Launching characteristic o of the constant Mach number line . . . . . . . . . . 54
3.15 Solution to the non-simple region . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.16 Illustration of boundary formation through extended characteristic mesh . . . . 56
3.17 Determination of the boundary segments . . . . . . . . . . . . . . . . . . . . . . 57
3.18 Normalized pressure distribution selected from optimizer . . . . . . . . . . . . . 61
3.19 Resulting static pressure distribution . . . . . . . . . . . . . . . . . . . . . . . . 62
3.20 Evaluation of  tness through intersection of characteristics . . . . . . . . . . . . 63
3.21 Evaluation of  tness for non-simple region and boundary . . . . . . . . . . . . . 64
3.22 Mach 2 Nozzle Contour solved with 20 characteristics,  = 1:23 . . . . . . . . . 67
3.23 Mach 3 Nozzle Contour solved with 20 characteristics,  = 1:23 . . . . . . . . . 67
viii
3.24 Mach 4 Nozzle Contour solved with 20 characteristics,  = 1:23 . . . . . . . . . 68
3.25 Mach 5 Nozzle Contour solved with 20 characteristics,  = 1:23 . . . . . . . . . 68
3.26 Comparison of the above contours . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.27 Example of contours that can be used for mapping thrust coe cient at a given
truncation length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.28 Coe cient of thrust versus truncation length for shock free contours . . . . . . 73
4.1 Solution can be developed through streamlines[4] . . . . . . . . . . . . . . . . . 84
4.2 Converging section geometry with dimensions normalized by throat radius . . . 87
4.3 Converging section pressure distribution for LOX H2 with Pc = 27:186 Bar . . . 90
4.4 Procedure of  nite rate chemistry analysis . . . . . . . . . . . . . . . . . . . . . 90
4.5 Initial Point Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.6 Points determined along the centerline through symmetry . . . . . . . . . . . . 93
4.7 Update the boundary region characteristic for chemistry . . . . . . . . . . . . . 94
4.8 Chemically reacting solution mesh for Mach 4 nozzle . . . . . . . . . . . . . . . 95
4.9 Comparison of isentropic kernel to kernel including  nite rate chemistry . . . . 96
4.10 Speci c heat ratio across expansion length . . . . . . . . . . . . . . . . . . . . . 96
4.11 Comparison of isentropic Mach 4 nozzle to nozzle including  nite rate chemistry 97
A.1 Transonic Flow Zone[5] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
E.1 Mass fraction of LOX Hydrogen system in converging section . . . . . . . . . . 114
ix
E.2 Mass fraction of H in converging section . . . . . . . . . . . . . . . . . . . . . . 115
E.3 Mass fraction of H2 in converging section . . . . . . . . . . . . . . . . . . . . . 115
E.4 Mass fraction of H2O in converging section . . . . . . . . . . . . . . . . . . . . 116
E.5 Mole fraction of O in converging section . . . . . . . . . . . . . . . . . . . . . . 116
E.6 Mole fraction of O2 in converging section . . . . . . . . . . . . . . . . . . . . . 117
E.7 Mole fraction of OH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
x
List of Tables
1.1 Performance factors of di erent hydrocarbon fuel compared to liquid hydrogen[6] 9
3.1 Pressure distribution optimization parameters for a 20 characteristic solution . . 60
3.2 Fitness weights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.3 Directly prescribed normalized pressure points . . . . . . . . . . . . . . . . . . . 66
4.1 Method of Characteristics Equations . . . . . . . . . . . . . . . . . . . . . . . . 83
4.2 Elemental breakdown of propellant fuel oxidizer pairs . . . . . . . . . . . . . . . 91
D.1 Simple reaction mechanism for Hydrogen/Oxygen combustion[7] . . . . . . . . . 110
D.2 Third-body e ciency ratios for Hydrogen/Oxygen . . . . . . . . . . . . . . . . . 110
D.3 Simple reaction mechanism for CHON combustion[7] . . . . . . . . . . . . . . . 111
D.4 Third-body e ciency ratios for CHON[7] . . . . . . . . . . . . . . . . . . . . . . 111
D.5 Reaction mechanism for C,Cl,F,H,N,O system combustion[7] . . . . . . . . . . . 112
D.6 Third-body e ciency ratios for C,Cl,F,H,N,O system[7] . . . . . . . . . . . . . . 113
xi
Chapter 1
Introduction
The performance of a rocket engine is dependent upon the momentum imparted to the
gaseous products of combustion by ejecting them through a nozzle. As these exhaust gases
pass through the nozzle they are continually accelerated. The speed of the gases begins in
the low subsonic range, in the head end of the thrust chamber, and increases to a highly
supersonic  ow. The thrust produced by a rocket engine accrues losses due to, among other
things,  ow divergence from the axial direction of the nozzle. Thus a nozzle that produces
maximum thrust for a given exit condition is one with a uniform  ow  eld across the exit
that is directed axially. To achieve this maximum thrust condition, a nozzle contour must,
also correctly expand the  ow. However, contoured nozzles required to achieve such  ows
are generally extremely long, heavy, and undesirable for rocket nozzle applications.
The contour of a converging diverging nozzle is comprised of three sections: convergent
section, throat section, and divergent section. Each section a ects the  ow di erently and
must be considered to produce the desired nozzle performance. Analysis of nozzle  ows
can be performed assuming that the exhaust gases expand adiabatically as an ideal gas
with constant speci c heat ratio. In addition, consideration of chemical reactions due to
incomplete combustion, chemical properties in the exhaust gases, viscous boundary layer
e ects and radiative heat loss can be considered.
The design of a rocket engine nozzle can be developed using the method of characteris-
tics. The method allows for physical boundaries to be found directly making the MOC an
attractive design tool for  nding the nozzle contours. The method can also be coupled with
 nite rate chemistry, allowing for the nozzle contour to be designed with chemical reactions
1
accounted for, including consideration of incomplete combustion, where the chemical prop-
erties of the exhaust gases are captured by analyzing the reactions and recombination of the
combustion products.
The use of the axis-symmetric method of characteristics coupled with chemical kinetics
of the combustion products has been developed by combining the method with the use of
an optimization routine, which allows for an automated design approach to liquid rocket
engines to meet desired performance goal.
1.1 A Brief History of Rockets
As early as 600 AD the Chinese had found that a mixture of charcoal, sulfur, and
saltpeter could be used as a propellant allowing the creation of their rocket-propelled  re
arrows[8]. These were the  rst true rockets. The use of such weapons slowly spread to western
civilization. Despite the rockets early advent, technical understanding and development of
the rocket did not begin until the end of the 19th century[9].
The Russian school teacher Konstantin Eduardovitch Tsiolkovsky is credited to be the
 rst true rocket scientist[8]. As a student he was infatuated with space travel, he carried this
with him into his career, and slowly worked on rocket theory without institutional support.
His most known contribution is the Tsiolkovsky rocket equation.
 v = Ispgolnmim
t
(1.1)
He also worked on a design for a liquid hydrogen, liquid oxygen rocket engine. He
published his work on the design and the rocket equation in 1903, the same year as the  rst
 ight of a powered airplane[8]. However, due to his lack of funds he performed no practical
experiments and generated no data on his design.
Convinced that rocketry was the key to space  ight, Dr. Robert H. Goddard was the
 rst to apply theory and develop practical working designs. As with Tsiolkovsky, Goddard
2
also determined that liquid hydrogen, liquid oxygen would be an e cient rocket propellant.
In 1914, Goddard was granted patents on his designs for propellant feed systems, combustion
chambers, nozzles, and multistage rockets. Obtaining a small grant from the Smithsonian
Institution in 1917, he was able to secure his career in rocketry, allowing him to publish
his work \Reaching Extreme Altitudes" in 1919[8]. Goddard went on to show how de Laval
nozzles could be applied to the combustion exhaust gases to accelerate the  ow. Most modern
day rocket engines use some form of a Laval nozzle.
During the early 1920?s, Goddard spent much of his time performing activities in his
laboratory. This time was spent with testing and development of what would become the
 rst liquid propellant rocket engine. Through these experiments, Goddard proved that a
rocket can work in a vacuum, by showing that the rocket needs no air to push against[10].
On March 26, 1926 Goddard launched his 10ft long rocket that used a liquid oxygen and
gasoline propellant. The rocket  ew to an altitude of 184 ft and reached a maximum speed
of 60 mph[8]. Seeming to be only a modest feat, this  rst  ight is equivalent to that of the
 ight of the Wright brothers in importance to rocketry.
Larger and larger rockets were designed and tested during the 1930?s, however the de-
velopments were slow and steady. The outbreak of World War II greatly accelerated the
advancement of the rocket engine. This especially occurred in Germany with the develop-
ment of the A-4 rocket, commonly known as the V-2 [11]. Much of this development was
due to the work of Herman Oberth, who was considered the \founding father" of German
rocketry. On October 3, 1942 the V-2 took its  rst  ight, the signi cance of this  ight was
the demonstration that a rocket had the potential to obtain orbit[11].
Following World War II intense research lead to the evolution of liquid rocket propulsion
from its rudimentary basis to a sophisticated science that has produced the engineering feats
of modern day[12]. Arguably the most important development in liquid-propellant rocket
engines came with the development of the Pratt and Whitney Aircraft RL10 powerplant.
The RL10 was the  rst rocket engine to use liquid hydrogen as a fuel[13]. This technological
3
breakthrough ultimately allowed the  rst lunar landing and lead to the development of the
Viking Spacecraft to Mars, and Voyager (to Jupiter and beyond[12]).
Today the state-of-the-art in high performance liquid-propellant rocket engines is rep-
resented by the Space Shuttle Main Engine. A single SSME delivers a vacuum thrust of
512,300 lb while the propulsion system weighs about 7775 lb [12, 14]. The SSME has a
vacuum speci c impulse of 452 seconds. This staged-combustion cycle engine is impressive
machine from head to toe. The hydrogen fuel is fed using a 71,140 hp turbopump. The
SSME has achieved a 100%  ight success and demonstrated reliability of 0.9995[14]
1.2 The Liquid Rocket Engine
Liquid rocket engines are chemical rockets that use liquid propellants. A complete
rocket engine is comprised of many highly integrated systems that allow for the operation
and control of the engine. There are four major subsystems that a liquid-propellant rocket
engine consists of:
 Thrust chamber
 Propellant feed system
 Propellant control system
 Thrust vector control system
The thrust chamber is used to convert stored chemical energy to directed kinetic energy
through the combustion of propellants and the acceleration of the exhaust gases. The pro-
pellants are supplied to the thrust-chamber at the correct pressure by the propellant feed
system.
The propellant feed system is comprised of storage tanks for the propellant, feed lines,and
turbopumps or pressure regulators. There are two basic types of propellant feed systems,
4
direct pressurized gas fed and turbopump fed. Direct pressurized gas feed systems store pro-
pellants in pressurized tanks, limiting the chamber pressure to the design of the tank. High
pressure tanks become heavy, due to tank thickness required to withstand such pressures[12].
Turbopump feeds use turbopump that are fed propellant at relatively low pressures.
The pumps can raise the pressure of the propellant to high pressures. This allows for high
pressures at much lower weight cost[12].
The turbine-drive system directs the appropriate gases through turbines in order to
power the turbopumps used in the propellant feed system. The gases used to drive the
turbans are directed from gas generators, preburners, heat exchangers, or tapped from the
main combustion chamber.
The propellant control system provides regulation of the propellants in order to control
the propellant mixture ratio and the total amount of propellant. This allows thrust control
and proper operation. The propellant control system is typically a system of valves.
Other control systems allow for proper and desired operations of the rocket. The focus
of the present development lies within the thrust chamber, which is be discussed in Sec 1.2.1
1.2.1 Thrust Chamber
The thrust chamber is the single most important system of a rocket engine as it is
the system of the rocket designed for the generation of thrust. The thrust chamber for
a liquid engine is comprised of three main sections: injector, combustion chamber, and
supersonic nozzle. It is here that the liquid propellants are injected, atomized, vaporized,
mixed, combusted, accelerated and ejected. Figure (1.1) displays the di erent zones of  ow
within a typical thrust chamber.
The injector injects both the fuel and oxidizer into the combustion chamber in such away
that the liquids are separated into small droplets, this process is referred to atomization. The
injectors are designed in such away that the injection streams of the fuel and oxidizer impinge
on one another to form a well mixed solution in the correct mixture ratio for combustion.
5
Supersonic
Expansion
Zone
Transonic 
Flow Zone
Rapid 
Combustion 
Zone
Streamtube 
Combustion 
Zone
Sonic Line
Injection &
Atomization
Zone
Combustion Chamber
With Subsonic Flow
Figure 1.1: Flow Zones in a Typical Thrust Chamber[1]
This impingement in many cases aides to the atomization, through molecular interaction[1].
The design of these injectors is a complex subject beyond the scope of this investigation. It
will be assumed henceforth that the injectors are designed in such away that proper mixing
occurs at any speci ed mixture ratio.
The combustion chamber is the section of the thrust chamber where the burning of
the propellants takes place. The combustion chamber is designed in such a way that it
retains the mixture of propellants for enough time that it ensures complete mixing and
combustion. The time that is required for this to happen is referred to as the combustion
residence time and is e ected by many parameters. Combustor geometry, injection condition,
propellant combination are a few such parameters that must be considered in combustion
chamber design[12]. Modern design of the combustion chamber typically calls for a cylindrical
chamber. For this development, the combustion chamber will be assumed to be cylindrical
and of the proper volume to ensure complete combustion to equilibrium conditions. Also,
ignition processes are beyond the scope of this investigation therefore it will be assumed that
the combustion process is taking place during operating conditions.
6
Following the combustion chamber is a converging diverging nozzle designed to produce
supersonic exit conditions. This section of the thrust chamber is used to convert the chemical
energy released in the combustion process to directed kinetic energy. This conversion of
energy ultimately results in the production of thrust. De Laval, or converging diverging
nozzles, are the focus of this investigation, and although other nozzle types exist they will
not be discussed. A good discussion of annular nozzles can be found in Ref. [1],and Ref. [12].
The development of supersonic nozzle design will be discussed in great detail in Chapter 2.
1.2.2 Thrust Chamber Performance
The performance of a liquid propellant rocket engine is directly dependent on the com-
bination of propellants used, the combustion e ciency of the propellants, and the expansion
of the gases formed from the combustion process. For the most part the performance of a
rocket engine can be expressed through the performance parameter called speci c impulse,
ISP. The speci c impulse is de ned as the thrust divided by the weight  ow rate measured
at sea level gravity[1]. This de nition provides that the units of speci c impulse are in time,
typically in seconds.
Isp =
Rt
0 Fdt
goRt0 _mdt (1.2)
The thrust force generated by the rocket can be expressed as
F = _mue +Ae(pe pa) (1.3)
Where the exit velocity can be found through
ue =
s
2 RuMWTo
  1
 
1 pep
o
  1
 
 
(1.4)
7
From these relationships it is apparent that the speci c impulse is dependent on the
speci c heat ratio, molecular weight and combustion temperature of the propellant. This
can easily be shown by rewriting Isp as
Isp =
vu
ut 2 
  1
RuTo
MW
"
1 
 p
po
 (  1)
 
#
(1.5)
Maximizing speci c impulse has an important impact on the design of the total vehicle.
From the de nition of speci c impulse above and the rocket equation, Eq. (1.1), the rela-
tionship between speci c impulse and payload is apparent. Large values of Isp result in an
increase in allowable payload. This is achieved by the minimization of the speci c heat ratio
and the molecular weight of the working  uid.
1.3 Liquid Propellants
As with all rockets, liquid propellant rocket engines are supplied internally with both
the fuel and oxidizer. Unlike an air breathing propulsion device, which supplies the oxidizer
in the form of air through some means of intake, rockets internal supply of both the fuel
and oxidizer allows for the combustion of the propellant in a vacuum. Furthermore, the
rocket is not limited by aerodynamics when supplying the oxidizer, allowing for a rocket
to operate at very high speeds while in the atmosphere. There are three basic types of
liquid rocket propellants: monopropellant, bipropellant, and tripropellant. The selection
of the propellant is typically based on the application of the rocket. The main focus of
this development is bipropellants, with liquid oxygen, LOX, as the oxidizer. However the
development is designed so that other propellants can be used. A list of propellants that are
readily available for analysis will be provided in a section 4.2 .
Clearly, from Eq 1.5, the performance of a rocket is dependent on the propellants used.
The exit velocity, speci c impulse, and thrust of a rocket are dependent on the speci c heat
ratio, molecular weight, and combustion temperature. This can easily be seen in Eqn. (1.2)
8
and Eqn. (1.4) and the thrust equation. From these relationships, it is clear that careful
consideration of the chemical kinetics of the propellants will lead to a high  delity solution
to the design problem.
Selection of the propellants to be used is based on trade studies or optimization of
mission design. This allows for the relationship between rocket performance, total weight
of the vehicle, and other design considerations to be investigated. One characteristic of the
propellant that lends to the selection is the density speci c impulse. De ned as
Id = Isp o f(1 +r)r 
f + o
(1.6)
Where  is the speci c gravity of the propellant component, and r is the mixture ratio.
Performance of LOX Oxidizer Propellants
Fuel Mixture Bulk Density Combustion Temp Vacuum Isp Density Speci c Impulse
Ratio kg=m3 K sec kg s=L
Hydrogen 6.0 358 3610 455.9 124
Methane 3.0 801 3589 368.3 235
Propane** 2.7 1014 3734 361.9 367.1
RP-1 2.5 1026 3803 354.6 294
*Isp values based on 100:1 expansion into vacuum. Assumed chamber pressure 20 MPa.
**Propane data based o of sub-cooled propane.
Table 1.1: Performance factors of di erent hydrocarbon fuel compared to liquid hydrogen[6]
This characteristic of the propellant, can be useful in the overall design of a vehicle
where any reduction in system mass can translate into an increase in payload mass for a
given thrust capability. A comparison of a few hydrocarbon fuels and hydrogen is given in
Table 1.1. It is clear hydrogen as a fuel yields the highest performance. However, the low
density allows for the hydrocarbons to be a better choice for situations when mass savings
outweighs overall performance. Ultimately the selection of the propellant is chosen as a
vehicle design parameter and a rocket engine is designed for the speci c propellant types[6].
9
1.4 Combustion of Liquid Propellants
For a chemical rocket, the working  uid of exhaust gases is produced by active combus-
tion of propellant. Combustion allows the release of stored chemical energy in the propellant.
This transfers the energy into the form of thermal, and random kinetic energy that can be
directed by the use of a nozzle to impart momentum on the system for the generation of
thrust[1]. The behavior of the combustion is dependent on the propellants being used. Inves-
tigation into the combustion process of the propellants results in more accurate predictions
of the performance parameters of the rocket engine as well as the thrust produced[15, 1].
Typically, combustion of liquid propellants in thrust chambers is very e cient. E -
ciencies in well-designed combustion chambers for liquid rocket engines can be found in the
range of 95 to 99.5%[1]. The ability of the combustion chamber to provide thorough mixing
through correct injection and atomization of the propellants along with the high reaction
rates at the high chamber temperatures allow for such e ciencies. The nature of combustion
provides that the exhaust gases are chemically active. The chemically active exhaust gases
can be modeled by treating the gases as equilibrium, frozen, or non-equilibrium chemically
reacting. Solutions to nozzle performance for equilibrium  ow and frozen  ow can readily be
obtained using generalized one-dimensional  ow analysis. However, these are the extreme
cases. Equilibrium  ow corresponds to the best case where the reaction rates of the  ow are
in nite. The lower limit corresponds to frozen  ow where the reaction rates are zero. Both
the equilibrium  ow and frozen  ow cases can be shown to be isentropic processes. For frozen
 ow the composition of the gas mixture has no change resulting in the entropy change in the
system to be zero, while with equilibrium the sum of products of the chemical potential and
changes in mass fraction of the reactants is equal to the sum of the changes in the products
of the chemical potential and changes in mass fraction of the products resulting in no change
in entropy in the system[4]. The two conditions can be shown through the e ect of chemical
reactions on the entropy of a mixture and is discussed in detail by Zucrow[4]
10
Tds = 
NX
i
 idCi (1.7)
In reality, the chemically reacting  ows have  nite rates of reaction somewhere between
these two limiting cases, and because the chemical reactions occur at a  nite non-zero rate
the entropy of the system will change.
The equilibrium composition and  nal equilibrium thermodynamic state can be found
through the minimization of Gibbs free energy of the system, the minimization of Helmholtz
free energy or the maximization of entropy of the system[16]. Characterizing the thermody-
namic state by temperature and pressure lends itself to using the minimization of Gibbs free
energy for equilibrium calculation due to temperature and pressure being natural variables
of Gibbs free energy[16]. Gibbs free energy represents the amount of energy available to do
work on an open system in which mass is exchanged within the surroundings. For a general
reaction the Gibbs free energy change can be written as
 G =  H T S (1.8)
which can become
 G =
X
[ j( fG )j]products 
X
[ j( fG )j]reactants (1.9)
Gordon and McBride[16] used a Lagrangian multiplier technique for the minimization of
Gibbs free energy to  nd the equilibrium composition of a reaction. Gordon and McBride?s
technique is used for the equilibrium calculation needed in this development. Their method
and the program which they developed for a generic equilibrium composition calculation is
discussed in Appendix B and a complete explanation of the their process can be found in
[16][17].
When the exhaust gas velocities become primarily axial the gases are passing through
the streamtube combustion zone. Although the gas composition is not truly uniform when
11
traversed radially across the combustion chamber, by time the exhaust gases become striated
they have been well mixed and can be considered uniform composition for following the
reactions through the subsonic portion of the nozzle. This allows to treat the problem as
quasi-one-dimensional  ow until the Mach number of the  ow reaches just beyond unity,
based on a prescribed pressure distribution through the combustion chamber, converging
section and throat. As these streamlines of chemically reacting gases  ow through the
converging section of the nozzle they are accelerated more. Depending on the propellant this
zone may become nearly chemically frozen if the axial velocity of the  ow is at a rate that
is faster than that at which the products can react.
12
Chapter 2
Review of Literature
2.1 Contoured Nozzle Design
It can easily be shown that in order to expand  ow through a duct from subsonic  ow
to supersonic  ow the area of the passage that the  uid is passing must  rst decrease in
area and then increase in area. This area relationship is the basis in nozzle design given
in Eq. (2.1). The relationship between local Mach number and the local area ratio was
found through the study of quasi-one-dimensional  ow. Although this relationship provides
no information for the contour of such a duct[18] or the losses that are associated with a
multidimensional  ow  eld.
A
A =
1
M
 2
 + 1
 
1 +   12 M2
   +1
2(  1)
(2.1)
In early rocket engine designs, conical nozzles were used almost exclusively[12]. The
geometry of a conical nozzle, depicted in Fig. 2.1, allows for ease of manufacturing. The
simple design also allowed for any given expansion ratio to be created with simply extending
the length of the nozzle.The  ow in a conical nozzle is not turned to become axial and thus
no strong shocks will form in the nozzled during operating conditions.
However, the conical nozzle has associated performance losses because the  ow has a
non-axial component[19]. In order to save weight on a potentially very long nozzle, large
cone angles must to be used. As the cone angle is increased the loss of thrust,due to non-
axial  ow conditions at the exit, increases. The momentum losses are typically incorporated
through the theoretical correction factor below[1]. This yields the percentage of the ideal
13
a
R
t
R
e
C
L
L
n
Figure 2.1: Typical Outer Mold Line of a Conical Nozzle
exhaust velocity that is axial. The in uence as the cone half angle increases can be easily
determined through [12, 1]
 = 12(1 +cos ) (2.2)
These losses and the extreme weights associated with small divergence angles make the
conical nozzle a poor choice of geometry. In order to impart the momentum axially at the exit
the contour nozzle or bell nozzle was developed, as depicted in Fig. 4.1[12]. The contoured
nozzle increases performance the  ow must be turned so that the exit  ow of the nozzle
is parallel to the centerline of the nozzle at the exit. This imparts all of the momentum
axially. To develop this condition a contoured nozzle can be designed so that immediately
behind the the throat the  ow undergoes a pronounced expansion. As the  ow continues
throughout the nozzle the expansion continues to occur but the convex shape allows for
the rate of expansion to diminish. Expansion waves in the  uid allow for the direction of
the supersonic  ow to change. As the gas  ows through these expansion waves the  uid
is allowed to change its direction and increase its velocity with little loss of energy, and a
contour can be formed to follow these changes. The physicality of the  ow lends itself to
the method of characteristics for analysis, allowing the method to model these expansion
14
R
sub
R
sup
R
e
R
t
L
N
Figure 2.2: Typical Outer Mold Line of a Contoured Nozzle
waves as mathematical characteristics[1][20]. A detailed description of the method will be
presented in a later section.
The  rst successful implementation method of characteristics for nozzle design was per-
formed by Ludwig Prandtl and Adolf Busemann[18] in 1929. Since the implementation by
Prandtl and Buseman, the method of characteristics has become a fundamental basis in noz-
zle design. This is because the method allows for physical boundaries to be located. Prandtl
and Busemann implemented the method graphically to solve two-dimensional nozzle prob-
lems. A comprehensive presentation on both the graphical and the analytical approach to
two-dimensional method of pharmacogenetics was completed by Shapiro and Edelman in
1947[21]. In 1949, Foelsch proposed a method for developing solutions to axis-symmetric su-
personic streams using a method of characteristics approach[21]. Antonio Ferri extended
the approach to axis-symmetric  ows in 1954, through a theoretical adaptation to the
mathematics[22, 21].
In 1959, Guerntert and Netmann implemented the analytical approach for the devel-
opment of supersonic wind tunnels with desired mass  ows[3, 21]. This implementation de-
veloped a solution based o of initial conditions along the nozzle centerline. The Guerntert
15
and Netmann consideration of the approach for wind tunnel design had no length require-
ment but required uniform exit  ow. Their solution resulted in di culties designing short
length and large expansion ratio nozzles[3]. Their work also showed that truncation of nozzle
lengths resulted in only a small reduction of vacuum speci c impulse from the uniform  ow
case[3].
The work of G.V.R. Rao used the method of characteristics as part of his solution in
developing contour nozzle designs. Rao developed a method of designing a contoured exhaust
nozzle for optimum thrust of a  xed length nozzle. Rao?s solution used a combination of
Lagrangian multipliers and method of characteristics[19, 23].
Rao?s solution used a control volume that does not span across the nozzle?s exit. The
prescribed control volume follows a characteristic. This creates a boundary on the control
volume be at an oblique angle to the axis of the nozzle[19, 23], as illustrated in Fig (2.3).
The control volume being set to an angle allowed Rao to integrate along a characteristic line
to obtain mass  ow rate and thrust.
L
N
Control
Surface
Figure 2.3: Typical characteristics net in a Rao nozzle with illustration of control surface
Rao?s method for obtaining a nozzle contour, begins by prescribing the geometry of the
throat. This sets the subsonic and supersonic radius as two unique values. The subsonic
16
radius of curvature is set to allow for substantially supersonic  ow at the throat of the nozzle
and was obtained on the basis of Sauer?s solution[19], which will be explained in depth later.
The supersonic radius of curvature is prescribed to a value that is large enough to ensure that
 ow separation at the throat is not a concern, yet kept short to allow for the minimization
of total nozzle length[19, 23]. Rao?s objective was to form nozzle contours that allowed for
the maximization of thrust for a given nozzle length.
When examining a nozzle without uniform exit  ow, following a line perpendicular to the
centerline that intersects the nozzle exit, the Mach number will not be found to be constant.
Starting at the exit, the Mach number along such a line will increase as the centerline is
neared. With the exit Mach number selected, Rao  nds the the wall slope at the exit under
these  ow conditions. From this point Rao?s method constructs a characteristic net in the
region between the exit, a calculated point on the control volume and the transition point
on the nozzle wall where the contour will transform from convex to concave geometry with
respect to the  ow. The contour of the nozzle is then found following a stream line between
the transition point and the exit point.
The method Rao presented in 1958 was the basis for the Rao parabolic approximation
to bell nozzle contours. The parabolic nozzles are presented in detail in both Sutton[1] and
Huzzel and Huang[12]. The Rao approximation method prescribes a set of parabolas that
allow for a contour that is near the optimum thrust contours developed by his method.
These nozzle contours are chosen by specifying the fractional nozzle length of a 15 degree
half-angle conical nozzle, and the expansion area ratio. The appropriate parabola is chosen
from these parameters. The parabola speci es the exit  ow angle and the initial turning
angle. The throat geometry is then found using a subsonic radius of curvature of 1:5rt
and a supersonic radius of curvature of  0:4rt[1, 12]. Neither Sutton or Huzel and Huang
explain the selection of the supersonic radius of curvature. Sutton goes on to note that small
changes in the supersonic radius have only small e ect on the overall approximation. The
method presented by Rao was widely used in the 1960?s. However it has been found that
17
Figure 2.4: Formation of internal shocks in conical and contoured nozzles[2]
the contours formed by the Rao approach and the Rao parabolic approximation produce
shocks. Shocks occur due to a pressure rise in the supersonic  ow  eld. In nozzles, the
formation of shocks can be caused by the presence of a downstream high pressure condition,
as in ignition transients, or by the how the contour is constructed to control the expansion
of the  ow[2]. The formation of shocks due to the wall construction occurs as a result of
turning the wall into the  ow, allowing a compressive wave to form. Interaction between the
compressive waves and expansion waves can cause a pressure rise. Figure (2.4) shows the
formation of shocks in di erent contours[2]. The TIC nozzle is a truncated ideal contour,
which is the focus of this thesis. TOC is a thrust optimized contour nozzle or Rao nozzle
as described above. TOP is a parabolic bell nozzle contour. The nozzle wall can be formed
such that a condition can be developed in which the compressive turn is exactly canceled
18
by an expansion wave. A viable solution can be obtained by speci cation of the centerline
pressure distribution so that this condition for compressive cancelation occurs.
Obtaining viable axis-symmetric method of characteristics solutions can be a challenging
exercise for large expansion ratio nozzles[21]. The solution is heavily dependent on a speci ed
boundary condition along the nozzle centerline. This condition must be assumed prior
to the implementation of the method of characteristics, as a result small deviations from
viable conditions can cause characteristic lines in the solution to cross. Crossing of like
characteristics may produce a solution however this solution cannot be considered valid, for
a shock free contour.
In order to form a solution with axis-symmetric method of characteristics using the
centerline as a boundary condition, a center line pressure distribution must be assumed.
Three basic pressure distributions that have been used with success [21]. However, not every
pressure distribution will yield a valid solution for any given nozzle. The pressure distribution
must be  tted to the problem at hand. Three distributions that have previously been used
with success are: [3]
Cosine: ln p = 12(ln pt ln pe)cos xl + 12(ln pt + ln pe) (2.3)
Parabolic: ln p = l 2(x l)2(ln pt ln pe) + ln pe (2.4)
Cubic: ln p =
 ql + 2(ln p
t ln pe)
l3
 
x3 
 2ql + 3(ln p
t ln pe)
l2
 
x2 +qx+ ln pt (2.5)
The cubic pressure distribution has an additional factor that allows it to be manipulated.
This can be used to create a valid pressure distribution. Examples of pressure distributions
are illustrated in Figure (2.5) for a Mach 4 Nozzle with an expansion length that is 5:75rt.
Chamber pressure is 4400 KPa and chamber temperature is 3000 K. As can be seen the
cubic pressure distribution can be manipulated with the changes in q giving the cubic more
 exibility
19
Distance From Throat
normalized by r
t
P
r
essu
r
e, 
P
a
Figure 2.5: Pressure Distributions
Further examination into these pressure distribution shows that they are only valid for
certain cases. The cosine pressure distribution is the least useful of all. For most expansion
ratios the cosine distribution does not yield a realistic pressure distribution after the throat.
The parabolic pressure distribution is valid for small expansion ratios. The cubic pressure
distribution must be manipulated to the expansion ratios. For small expansion ratios the
value of q must be close to zero to provide valid solutions. As the expansion ratio increases
the value of q must decrease to provide valid pressure distributions. However as the expan-
sion ratio increases further all of these distributions tend to form solutions with crossing
characteristics. In order to produce a pressure distribution that allows for valid solutions
through out a wide range of expansion ratios a di erent approach must be used to specify
the centerline pressure distribution.
20
Haddad and Mosst[24] used a method of characteristics approach to axis-symmetric
supersonic nozzle design, this was met with limited success when solving the pure axis-
symmetric design requirements[21]. Wang and Ho man[25] made use of the two-step pre-
dictor corrector algorithm on pre-speci ed  ow  eld data. This has regularly been used
since[21]. The work of Wayne and Mueller applied the rotational axis-symmetric method of
characteristics employing a technique that allows the down stream characteristic point loca-
tions to be selected to  t the  ow- eld as it is developed. This technique was developed by
Hartree. The technique is still subjected to the problem created by the centerline boundary
condition. The problem can be overcome but the process to overcome the centerline bound-
ary condition has become case dependent[21]. These methods result in heavy user interface
and hit-and-trial solutions to  nd viable solutions to the case at hand[21].
In 1978 Allman and Ho man presented a procedure for the design of a maximum thrust
contours by a direct optimization method. The contour used was a second-degree polynomial
 tted to a prescribed initial expansion contour. The contours produced were similar to that
of a Rao nozzle. Allman and Ho man showed that a polynomial could be used to develop the
nozzle boundary with comparable results to a Rao nozzle, however the  ow  eld was solved
in much the same way, and the solution di ered only slightly in the formation of the nozzle
boundary. Essentially, their solution was a Rao nozzle with a polynomial  tted boundary
instead of a boundary determined by the solution[26].
In 2011 Hart eld and Ahuja[21] put forth an e ort to automatically determine a viable
centerline pressure distribution when the centerline pressure distribution is chosen as the
governing boundary condition. The technique used by Hart eld and Ahuja was to remove
the ability of the characteristics to cross despite using a non-viable axial centerline pressure
distribution. The method is a non-physical solution to the internal  ow, however the resulting
contour of the nozzle allowed for the appropriate expansion lending the contour to a physical
solution.
21
2.2 Transonic Flow Zone
Transonic  ow is the  ow regime where the  uid transitions from subsonic to supersonic
velocities[5]. The transonic  ow regime has been intensely studied. The work of Sauer[5]
detailed the complexities and the mathematical treatment of such  ows, especially as applied
to the passage of  ow through Laval nozzles. The  ow in the throat region of a converging-
diverging nozzle under choked  ow conditions is transonic[4].
The work of Sauer has been the primary basis for the treatment of the transonic  ow
zone in supersonic nozzle design, because Sauer?s method is a closed form solution for the  ow
 eld in the nozzle throat[4], and it can produce excellent approximate solutions for nozzles
with a large subsonic radius of curvature relative to the throat radius. In 1958, Rao[19]
used Sauer?s approach in his treament of transonic  ow region. The work of Nickerson et
al. [7] based their work in the transonic analysis o of Sauer, applying the analysis to zoned
expansions. Nickerson et al. applied this transonic model to striated  ow for  ows with
variable mixture ratios[7]. Because of the importance of Sauer?s work on the treatment of
the transonic  ow zone his work is outlined and described in Appendix A.
The transonic solution is important to the method of characteristics solution. It allows
for the determination of a subsonic radius of curvature that allows for substantially supersonic
 ow at the nozzle wall at the minimum area point and also locates the position where like
 ow is on the axis of symmetry. In determining this line of constant substantially supersonic
Mach number, an initial value line can be formed so that the  ow satis es the wall boundary
condition at the throat exactly. Hence, the  ow at the wall boundary of the throat will have
a radial component of velocity exactly equal to zero.
22
r
x
e
h
Fl ow  D irect ion
 ?
sub
 ?
sup
S up er so ni c F lo wS ub son ic Fl ow
S on ic Li ne
Figure 2.6: Flow in the throat region of the nozzle with Sauer?s solution to the sonic line.
2.3 Method of Characteristics
The method of characteristics is a technique for solving partial di erential equations
based on absolute di erential calculus. The method is a numerical approach, but the nature
of the approach di ers from  nite di erencing schemes in that the solution provides the
information necessary for the location of a physical boundary in which produces shock free
 ow within the nozzle. The method was developed by the mathematicians Jaques Saloman
Hadamard in 1903 and by Tullio Levi-Civita in 1932[21]. The method of characteristics uses
a technique of following propagation paths in order to  nd a solution to partial di erential
equations. This is done through the reduction of the partial di erential equation to ordinary
23
di erential equations or algebraic expressions if appropriately restricted. These propaga-
tion paths are called the characteristics. These characteristics can represent such things as
geometrical form, discontinuities, or physical disturbances[20]. For the case of irrotational
supersonic  ow the characteristics follow the propagation of Mach waves. Because the  ow
is supersonic and these Mach waves have no in uence of the  ow upstream of themselves
and the method of characteristics may be used.
In 1949, it was proposed that an analytical solution to an all axial supersonic stream
could be developed using the axis-symmetric method of characteristics[21]. Working at
NASA Lewis, Guentert and Neumann implemented the analytical approach for designing
axis-symmetric supersonic windtunnels[3]. The approached used may be extended to form
the contour of a propulsive nozzle as well.
The applications of the method of characteristics for nozzle  ows are not limited to
the design of contours. The method may also be used to analyze the  ow  eld inside a
known contour as well[4]. The method is also not limited to the  ow within the nozzle.
The approach can be extended to analyze the exhaust plume for both underexpanded and
overexpanded nozzle  ow using the free pressure boundary of the exhaust plume.
2.3.1 Axis-symmetric Method of Characteristics
Assuming that the  ow through the nozzle is comprised of an ideal gas that is both
steady and maintains constant total enthalpy, the second law of thermodynamics, the con-
tinuity equation and the momentum equation can be used to model the  ow. With the
assumption that there are no shocks or viscous e ects present in the  ow the second law
of thermodynamics can be reduced to the condition for irrotationality. Thus the governing
equations are given as[21]:
~r ~V = 0 (2.6)
~r ( ~V) = 0 (2.7)
24
 D
~V
Dt  
~r= 0 (2.8)
It can be shown that combining these equations and solving using Cramer?s rule for
the indeterminate case yields the characteristic equations[21]. They can also be found using
a linear algebra approach as shown by Ferri[22]. These equations are devided into two
families, one that follows the left running characteristics and one that follows the right
running characteristics. The characteristic equations are[3]:
dy
dx = tan(   ) (2.9)
From the governing equations the compatability equations are also derived. These two
equations are grouped into the two families. The compatability equations are:
dV
V  tan d  
tan sin sin 
cos(   )
dx
y = 0 (2.10)
A complete derivation of the characteristic and compatibility equations can be found in
Ferri[22] and Zucrow[27, 4]. Both present a similar approach to the derivation.
The angle  is the Mach angle. Because  is dependent on the speci c heat ratio,  and
Mach number, M, it is a variable at every point in the  ow  eld. The characteristic points
will be described by the Mach angle, the  ow angle, and the velocity at each point in the
characteristics web. Although the characteristics are actually curves, they can be represented
as a line segment between characteristic points when distance between characteristic points
is kept small[3]. Figure (2.7) shows how each characteristic point is described.
Practical implementation of the above equations requires that they be discretized. The
di erential equations are discretized by rewriting them into a  nite-di erence solution. This
allows for point to point calculation of the  ow properties throughout the given  ow  eld.
Figure (2.8), shows a typical characteristic mesh. In order to form the solution information
must be known at points A and B in order for the information at point C to be found[3].
25
Radial Axis
Nozzle
Centerline
m
m
L
1
I
1
L
2
I
2
Q
V
u
v
Figure 2.7: Characteristic Point and its Associated Angles[3]
For a  rst approximation it is common to treat the location of the characteristics as the
characteristics themselves. That is the location of C may be calculated using the  rst order
 nite-di erence form of the characteristic equations if the location of points A and B are
known[3].
yC yA
xC xA = tan( A + A) (2.11)
yC yB
xC xB = tan( B  B) (2.12)
Solving these equations simultaneously results in the location of C. This essentially is
using the slope of the left running characteristic from A and the slope of the right running
26
A
B
Q
Q
A
B
C C'
m
m
B
A
Nozzle Centerline
V
B
V
A
Q?m
Q+m
Figure 2.8: Characteristic Mesh
characteristic fromB and their associated angles to  nd the intersection of the characteristics.
It is most practical to solve for the x location of C and then the y location. This is simple
algebra and results in the following[3]
xC = yA yB +xB tan( B  B) xA tan( A + A)tan( 
B  B) tan( A + A)
(2.13)
yC = yA + (xC xA)tan( A + A) (2.14)
The compatibility equations may be discretized and combined to form a  nite di erence
form of the di erential equations. The result of the discretization is then solved for the
velocity of the point C. Because the velocity at the point C is dependent on the  ow angle
at the point, velocity is determined through elimination of  c resulting in the following
equations[3].
VC = 1cot A
VA +
cot B
VB
fcot A(1 + $(xC xA)) + cot B(1 + M(xC xB)) + B  Ag
(2.15)
27
where
$ = tan A sin A sin Ay
A cos( A + A)
(2.16)
and
M = tan B sin B sin By
B cos( B  B)
(2.17)
In turn the  ow angle of the characteristic may be obtained by solving the compatibility
equations for  c yeilding:
 C =  A + cot A
 V
C VA
VA
 
 $(xC xA)) (2.18)
The solution at point C with the equations in presented thus far yields a  rst order
solution. However, it may be necessary to use a second order approach to form a valid
solution. The second order solution can be developed by solving Eqs 4.50, 4.53 for velocity
and  ow angle and then performing a second iteration upon the point C producing point C?.
In turn the  ow velocity of the characteristic may be obtained by solving the compati-
bility equations for  c yeilding[3, 22]:
x0c = yA yB +
1
2(tan( B  B) + tan( C  C))xb 
1
2(tan( A + A) + tan( C + C)xA
1
2(tan( B  B) tan( A + A) +tan( C  C) tan( C + C)
(2.19)
y0C = yA + (xC xA)12(tan( A + A) + tan( C + C)) (2.20)
V0C = F + G + B  A4((V
A +VC) 1(tan A + tan C) 1 + (VB +VC) 1(tan B + tan C) 1)
(2.21)
where
F = 2tan 
A +tan C
 2V
A
VA +VC +
$A + $C
2 (xC xA)
 
(2.22)
G = 2tan 
B + tan c
 2V
B
VB +VC +
MB + MC
2 (xC xB)
 
(2.23)
28
and
 0C =  A + 2tan 
A + tan C
 2(V0
C VA)
VC +VA  
$A + $C
2 (xC xA)
 
(2.24)
Using these equations the  ow properties are determined at C?, which becomes the
location of the characteristic and the starting point for future calculations. Because the
characteristic points are described by velocity and Mach angle the Mach angle and thermo-
dynamic conditions are determined through the Mach number at this point. This can easily
be calculated by using the de nition of enthalpy
To = T + V
2
2cp (2.25)
From Eq. (2.25) here the Mach number, pressure, and density of the  ow may be found
through the isentropic process and de nition of Mach number.
2.3.2 Valid Pressure Distribution Selection
In order to produce a pressure distribution that allows for valid solutions throughout
a wide range of expansion ratios a direct approach must be used to specify the centerline
pressure distribution. The functions presented in section 2.1 are di cult to implement while
obtaining solutions without crossing characteristics, and have the potential to produce non-
physical pressure distribution. A direct selection method together with a simple particle
swarm/pattern search optimizer[28] can produce valid pressure distributions for a much
larger range of expansion ratios than with the three provided in section 2.1. Two direct
methods are presented here.
Lagrange Polynomial Pressure Distribution
A method for developing pressure distributions is to form the pressure distribution with a
Lagrange polynomial, through the direct selection of pressure points. The associated method
for developing the pressure distribution is similar to that of polynomial interpolation. With
29
this method it is essential to use an optimization process so that the polynomial is generated
to the case and known to be valid.
Much like polynomial interpolation, the process in developing a pressure distribution
using Lagrange polynomials use points spread across an interval that describes the function
 tting a curve that passes through each point[29]. Although the exact pressure distribution
is not known, points along the interval are given a value that is assumed to be the known
pressure at that point and the distribution is determined by the curve that best  ts these
points. The location of the points along the interval is held constant yet the amplitude of the
points is allowed to change. Small changes are made until the pressure distribution allows
for a viable solution.
To select a pressure distribution using a lagrange polynomial a design interval [0;l] is
used. Along this interval there are set n+ 1 distinct nodes x0;:::;xn and the corresponding
values y0;:::;yn. It can be shown that with these distinct nodes and their corresponding
values a unique polynomial exists which exactly intersects each node at its corresponding
value[29]. This polynomial is expressed in the form
L(x) =
nX
i=0
yili(x) (2.26)
where
li(x) =
nY
0 j k
j6=i
x xj
xi xj (2.27)
If the distinct nodes x0;:::;xn are distributed uniformly across [0;l] then for functions
of the type needed to model a viable pressure distribution
limn!1jf(x) Lf(x)j6= 0 (2.28)
30
This is the phenomenon that occurs in Runge?s counterexample[29]. This phenomenon
occurs due to the equally spaced nodes. In order for the interpolating polynomial to have
uniform convergence an appropriate choice of the distribution of the nodes must be made.
Use of Gaussian quadratures considered with respect to the Chebyshev weight can pro-
duce a distribution of the nodes that will allow for uniform convergence. With this consid-
eration the nodes on the interval [ 1;1] are distributed by[29]
xi = cos (2i 1) 2 n (2.29)
This distribution is not de ned on the appropriate interval. The distribution of the
nodes over the desired interval, [0;l], can be found by an a ne transformation. The result
of the transformation gives for an arbitrary interval [a;b]
xi = 12 (a+b) + 12 (b a) cos
 (2i 1) 
2n
 
(2.30)
With the distribution of the nodes selected. The corresponding values of the nodes are
arbitrarily selected while following the scheme
y0 >y1 >y2 >:::>yn (2.31)
while
y0 = 1andyn = 0 (2.32)
Once a set of values is selected the corresponding pressures are set when both the exit
pressure and the initial method of characteristics pressure is known. Piecewise multiplication
and scaling as follows will provide the correct pressure values for the polynomial produced
above
p = (pt pe) [y] +pe (2.33)
31
Direct Selection Method
A very simple approach that can easily be implemented for solutions using a moderate
number of characteristics is direct selection of the pressures at each point. In order to produce
a valid pressure distribution a optimizer should be used.
Choosing the expansion length and the number of characteristics, the center line points
are distributed evenly across the expansion length. For each of these points a unique pressure
distribution is chosen, and the method of characteristics mesh can be generated and checked
for viability
This is a simple brute force approach that can be implemented quite successfully for
moderate numbers of characteristics using an optimizer.
2.4 Chemical Kinetics with Finite Reaction Rates
When considering a chemically reacting gas mixture with  nite reaction rates, a system
of reactions that model the mixture must be developed. This system of elementary reactions
is referred to as the reaction mechanism[4]. The reaction mechanism can be written as a
system of arbitrary chemical reactions with Eq. (2.34), where  and  0 are the stoichiometric
coe cients, A is the molecular formula, i is the ith species of the mechanism and j is the
jth reaction.
NX
i=1
 ijAi,
NX
i=1
 0ijAi (2.34)
(note:  and  0 are the stoichiometric coe cients of the reactants and products respectively)
For example a LOX Hydrogen system the reaction mechanism is described by the fol-
lowing reactions, where M denotes a third body [30]:
2H +M  H2 +M (2.35)
32
H +OH +M  H2O +M (2.36)
O +H +M  OH +M (2.37)
2O +M  O2 +M (2.38)
O2 +H O +OH (2.39)
H2 +O H +OH (2.40)
H2 +OH H2O +H (2.41)
OH +OH H2O +O (2.42)
It is clear that this reaction mechanism has a total of eight reactions and the total num-
ber of species for the mechanism is six. For a given reaction mechanism the production rate
 can be calculated for each species of the mechanism. The production rate of the species of
a reaction rate is called the species source function and is a function of local thermodynamic
properties as well as properties that are unique to the individual reaction, the equilibrium
constant, and the reaction rate. From the Law of Mass Action it can be shown that the
reaction rate is proportional to the collision rate between appropriate molecules and the
collision rate is proportional to the product of concentrations, thus the source function is
also proportional to product of concentrations of the species in the gas mixture[15].
 i = MWi
kX
j=1
  0
ij  ij
 "K
fj
NY
i=1
  ci
MW
  ij
 Kbj
NY
i=1
  ci
MW
  0ij#
Mj (2.43)
Where Kf and Kb are the forward and reverse reactions rates respectively and Mj is
the third body in uence of the reaction. The rate constants are equated by the equilibrium
constant which is the ratio of reverse to forward reaction rates
Kc = KfK
b
(2.44)
33
Typically either the forward or the backward rate for a given reaction is known, or available.
From the relation above the unknown rate can be calculated once the equilibrium constant
is calculated. There are many forms of the equilibrium constant and selection of the correct
de nition is essential to a correct solution. The form to be used to  nd the species production
rate must be expressed in terms of concentration. Data available for reaction rates typically
results in corresponding units of cm, mol, K, s for units. Thus the equilibrium constant must
be found in the form
Kc = 101:325R
uT
  
e S
 
Ru e
  H 
RuT (2.45)
where the universal gas constant Ru = 8:31447Jmol 1K 1 the exponent   is the change
in total moles of the reaction[31]. For reaction in Eq.2.42 this value would be 1. The value
101:325 comes from the value of 1 standard atmosphere in units of N=cm2 or kPa. The
enthalpy and entropy exponentials are equivalent to
 G f =
NSX
i=1
( 0i  i)g f (2.46)
where  G f is the standard state change in Gibbs free energy for the reaction at temperature[32].
The gibbs free energy of formation for a individual species g f can be found simply through
the enthalpy and entropy for the species at the speci ed temperature through
g f =  h f;298 +hT h298 TsT (2.47)
The values above can be found in the JANAF database or through the CAP program[33].
See Appendix C.
The reaction rate can be expressed as a function of temperature,  t to the reaction data
found through experiments. Data is usually provided as forward rates taking the same form
as the equilibrium constant equation above. The rate equation is
34
Kj = ajT nje( bj=RuT) (2.48)
where aj is the pre-exponential coe cient, nj is the temperature dependence of the pre-
exponential factor, and bj is the activation energy, these coe cients are tabulated for each
elementary reaction for many reaction mechanisms. They can readily be found in Kushida[30]
and Nickerson[7] for many oxidizer/propellant combinations.
The third body in uence of a reaction determines how the rates are changed with the
interaction of a third body to the reaction. Nominally if a M is introduced in the reaction
all other bodies besides the species of the reaction itself act as third bodies. However, some
species have more of an e ect on the reaction rates than others. The ratio that expresses this
in uence is given for each reaction that is in uenced by a third body and the the in uence
is calculated through
Mj =
NSX
i=1
mthr  ciMW
i
(2.49)
2.4.1 One Dimensional Chemical Kinetics
The development of the one-dimensional governing equations for non-equilibrium chem-
ically reacting gas is presented. Considering steady  ow that progresses in only a single
positively increasing direction the following equations are obtained.
Momentum:
 VdVdx + dpdx = 0 (2.50)
Energy:
h+ V
2
2 = constant (2.51)
35
In the above equation h is the total enthalpy which comprises of the sum of the sensible
enthalpy plus the energy of formation. Eq (2.51) can be rewritten on a mass basis with the
following equation.
h =
NX
i=1
cihi =
NX
i=1
ci
 Z T
To
cpidT +hoi
 
(2.52)
Continuity Equation:
Because of the chemical reactions the composition of the gas varies along the  ow
path. Thus conservation of mass must be applied to individual chemical species in the
gas. Assuming no mass di usion or thermal non-equilibrium occur the continuity equation
becomes
d _m = d( iAV) =  i(Adx) (2.53)
Where  is the source function as shown in Eq (2.43).
Substituting  i = ci into Eq (2.53) yields
d(ci AV) = ( AV)dci +cid( AV) =  i(Adx) (2.54)
Combining Eq (2.54) with the global continuity equation _m =  AV yeilds the di erential
equation
dci
dx =
 i
 A (2.55)
Equation of State:
36
The equation of state for a gaseous mixture is written as the sum of the product of the
species massfraction and the species spi c gas constant yeilding the equation:
p =  T
NX
i=1
ciRi =  RT (2.56)
The equations may be solved numerically for a non-equilibrium chemically reacting gas
mixture. In order to avoid singularities near the sonic point it is recommended to use a
speci ed pressure method to solve the system of equation simultaneously. This results in the
following equations and requires a known pressure distribution to solve the system[4]
dV
dx =
 1
 V
dp
dx (2.57)
d 
dx =
 1
 p
dp
dx 
 A
 
 (2.58)
dT
dx =
   1
 
1
p
dp
dx 
 B
 
T (2.59)
dci
dx =
 i
 V (2.60)
Area = _m V (2.61)
37
Where
 A = 1
PV
" NX
i=1
 iRiT   1 
NX
i=1
 ihi
#
(2.62)
and
 B =   1
 
1
pV
NX
i=1
 ihi (2.63)
It is worth noting that for non-equilibrium chemically reacting gases the constant pres-
sure speci c heat and hence the speci c heat ratio of the gaseous mixture cannot be assumed
constant[34]. To calculate the constant pressure speci c heat the CAP program developed
by Zehe, Gordon, and McBride[33] was found to be useful. This program allows for the
calculation not only of constant pressure speci c heat of numerous constituents but many
other thermodynamic properties as well.
38
Chapter 3
Axis-Symmetric Method of Characteristics Nozzle Design
The general approach in the design of an axis-symmetric nozzle based on a prescribed
centerline pressure distribution uses a three zone method of characteristics solution. This ap-
proach produces a \smooth throat" design, where the entire expansion contour is determined
through the method of characteristics. This di ers from Rao?s approach and the approaches
outlined by Zucrow, in that the initial expansion section of the nozzle was prescribed for
these approaches. If this initial expansion section is not prescribed correctly, crossing of like
characteristics, will occur. To avoid crossing characteristics with a prescribed initial turn,
small radius turns near the throat are required. The coalescence of like characteristics repre-
sents the formation of an oblique shock wave. If the oblique shock waves are relatively weak
they can be ignored without serious error in the solution[4]. However by using the method
of characteristics solution that contains no crossing of like characteristics a shock free nozzle
contour can be formed.
3.1 Nozzle Contour Design Process
The following describes the process for axis-symmetric method of characteristics nozzle
design, assuming isentropic  ow. Providing the chamber conditions and propellant prop-
erties, a contour will be formed to produce a nozzle of speci ed exit Mach number, with
uniform axial exit  ow.
1. Select Converging Section Geometry
2. Analysis of Transonic Region
3. Select Pressure Distribution
39
4. Select Number of Characteristic Points
5. Select Expansion Length
6. Determine Centerline Distribution
7. Perform Method of Characteristics
8. Check Solution for Viability
9. Perform Method of Characteristics Near the Throat
10. Form Nozzle Boundary
In order to design the smooth throat nozzle using the method of characteristics transonic
 ow at the throat must be analyzed. The smooth throat design takes inconsideration the
geometry in the throat region bring importance into how the geometry e ects the  ow both
in the subsonic region and the supersonic region of the nozzle. The sonic line cannot be
assumed to be a vertical line spanning the throat. The geometry of the throat will be curved
both in the subsonic and supersonic regions of the nozzle, thus the geometry e ect on the
transition from subsonic to supersonic  ow must be considered.
This entails a three fold analysis of the method of characteristics. The  rst being the
method of characteristics performed along the prescribed pressure distribution, called the
\Kernel". The second method of characteristics analysis is performed between the line of
substantial supersonic  ow from the transonic solution and the  rst left running characteristic
from the  rst analysis, this in the \Non-simple region". The third is the extension of the
characteristic mesh in the rest of the  ow  eld, called the \Boundary Region". These zones
are depicted in Fig 3.1.
40
r
x
uniform flow
Boundary Region
Kernel
Non-Simple Region
Sauer's Solution
Figure 3.1: Method of characteristics solution zones
3.1.1 Comparison of Processes
The process described in Sec. 3.1 is not fundamentally di erent than the process used
by Rao, Zucrow, and others [19][4][12][1]. However the approach to the solution has some
key di erences. A comparison to the method at hand and the method used by Rao to design
nozzle contours shows the key di erences between the two approaches.
Both approaches begin with the selection of the converging section geometry. With the
geometry known, an analysis of the transonic region is performed to produce a set of initial
conditions for the solution. Figure 3.2
The next step in the process is to set the initial conditions. In an approach such as
Roa?s approach, the geometry of the initial expansion region downstream of the throat is
(a) Rao Method (b) Proposed Method
Figure 3.2: Initial step for both nozzle design approaches
41
(a) Rao Method
Expansion Length
S
t
a
t
i
c 
P
r
e
ssur
e
,
 B
ar
(b) Proposed Method
Figure 3.3: Establishing the boundary conditions for both approaches
(a) Rao Method (b) Proposed Method
Figure 3.4: Kernel solution developed for both design approaches
prescribed. In the proposed method, the pressure distribution along the centerline of the
nozzle is prescribed, as can be seen in Fig. 3.3
With the initial conditions and the boundary conditions for both approaches deter-
mined, the next step is to form and analyze the \kernel" solution through the method of
characteristics. As shown in Fig. 3.4, an approach such as Rao?s, the \kernel" will contain
the solution from the transonic solution through the expansion length. For the proposed
approach the \kernel" solution will contain only the solution produced along the expansion
length.
For the proposed approach the method of characteristics is then used to form a solution
to the non-simple region, as depicted in Fig. 3.5. The solution is allowed to extend beyond
where the physical boundary will be, allowing for the physical geometry near the throat to
be determined. This is a main di erence in the result between the two approaches at hand.
42
(a) Rao Method (b) Proposed Method
Figure 3.5: Form method of characteristic solution to  ll the non-simple region
(a) Rao Method (b) Proposed Method
Figure 3.6: Form boundary region solution with the method of characteristics
An approach like Rao?s, the geometry just down stream of the throat was prescribed. With
the proposed approach and the extension of a method of characteristics solution near the
geometric throat, the geometry near the throat can be found.
From here, the solutions to each approach are carried on in the same manor. The
boundary region solution is formed using the method of characteristics, as shown in Fig. 3.6.
This provides for a solution space that allows for the boundary to be located.
With the boundary region solution formed, the location of the physical boundary can
be found, as depicted in Fig. 3.7. The location of the boundary allows for the contour of
the nozzle to be described. Either  ow tangency or characteristic continuity can be used
to obtain the location of the boundary along the characteristic lines. The location of the
physical boundary completes the solution to the nozzle design.
43
(a) Rao Method (b) Proposed Method
Figure 3.7: Formation of nozzle contour for both approaches
3.1.2 Converging Section Geometry
In order to design the smooth throat nozzle using the method of characteristics the
geometry throughout the throat region is considered. The subsonic geometry has a large
role in the transonic region of the nozzle and must be selected so that the location of where
the  ow becomes substantially supersonic, where the  ow downstream has does not a ect
the current position in the  ow. If the  ow is not substantially supersonic the method of
characteristics solution can not produce a viable solution. A rule of thumb for the converging
geometry that allows for substantially supersonic  ow near the throat is a radius of curvature
that is 2:0rt. This geometry allows for Sauer?s solution to hold, and is very similar to that of
the RL-10. This is not necessarily the optimal geometry for the subsonic region but clearly
suits its purpose of generating substantially supersonic  ow near the throat.
3.1.3 Analysis of Transonic Region
An analysis of the transonic region must be completed to provide a boundary condition
that is is substantially supersonic near the throat and to provide an axial location to start
the method of characteristics. For the transonic region analysis will be performed by using
the technique presented by Sauer[5]. The assumption of the sonic line being a vertical line
spanning the minimum area point is no longer considered. The sonic line will be a curve
that starts at some point upstream of the throat, crosses through the minimum area point,
44
and lies up-stream of the throat at the centerline of the nozzle.
r
x
e
h
Fl ow  D ir e ct ion
Figure 3.8: Transonic Region of Nozzle
Given a subsonic radius of curvature, the shape of the sonic line can be found. With
this information The Mach number at the boundary of the throat will be calculated and de-
termined if substantially supersonic. If so, the constant Mach number line will be developed
and followed to the centerline. This will be the initial point for the method of characteristics
solution as well as the starting point of the pressure distribution. Note that the coordinate
system for Sauer?s solution has an origin located where the sonic line intersects the nozzles
axis of symmetry[5].
From Sauer[5] the following equations are needed
 = r
 
8
s
2( + 1)r 
 sub (3.1)
45
Equation [3.1] is the axial displacement from the throat in which the sonic line inter-
sects the nozzle wall. For an axis-symmetric nozzle this distance is equal to that of the
displacement of the sonic line upstream of the throat. Resulting in
 =  (3.2)
The location of the throat is
x =  =  (3.3)
The Mach number in the transonic region is then calculated by  nding transonic components
u =  x +  + 14  2r 2 (3.4)
v =  + 12  2x r + 3r 3 ( + 1)
2
16 (3.5)
where
 =
s
2
( + 1) subr (3.6)
The Mach number can now be found from
M =
p
(1 +u)2 +v2 (3.7)
Using the above equations the position where the same Mach number, as found above, is
located on axis of symmetry can be found. Recall that the coordinate system for the transonic
analysis is based o of the location that the sonic line crosses the axis of symmetry. Thus
 must be added to the supersonic location found, representing this point in throat  xed
coordinates.
Solving the above equations for the sonic condition yields the sonic line which can be
projected as
46
xK =  + 14  r2k (3.8)
At every point along this line the Mach number is unity. The solution method of
characteristic solution cannot begin where the sonic line intersects the centerline at this
point. The  ow along this line is at Mach 1 which causes singularities in the Prandtl Myer
Expansion equation as well as with the slope of the characteristic because the  ow angle
along the centerline is  = 0. When the solution is initiated at a point where the  ow
is substantially supersonic the right running characteristic quickly intersects the sonic line.
Because the method of characteristics is only valid for supersonic  ow, the solution does
not encompass the entire supersonic portion of nozzle[3]. This zone has been referred to
in literature as the unde ned region or the non-simple region of the  ow, as illustrated in
Fig. 3.9.
Undefined 
Region
C
L
Sonic
Line
No
zzle
 
Wal
l
Flow Direction
Q+m
Characteristic
Figure 3.9: Unde ned Region[3]
To determine the location of the substantially supersonic point on the centerline the
Mach number is determined at the boundary of the geometric throat by Sauer?s solution.
The geometric throat lying upstream of the sonic line will have a  ow that is substantially
supersonic. Using Sauer?s solution a line of constant Mach number eminating from the
physical minimum area point along the wall. can be found and its intersection with the
47
centerline provides the location and Mach number of the initial centerline characteristic
point.
3.1.4 Select Number of Charcteristic Points
The number of characteristic points is chosen arbitrarily. The use of a small number
of allows for greater spacing in the characteristic mesh allowing for greater accumulation
of error. However, when a large number of characteristics is used the computational time
increases. When developing a method of characteristics code, the use of small numbers of
characteristics can be helpful for debugging purposes but will not result in highly accurate
solutions.
3.1.5 Select Pressure Distribution And Expansion Length
Through the isentropic process the pressure at the initial characteristic point can be
found. The location of the initial point will be the location that the constant Mach number
line emanated from Sauer?s solution. Starting the pressure distribution at this known location
avoids the singularities. The even distribution of expansion points can be
dx = ln
c
(3.9)
Where nc is the number of characteristics used in the solution and l is the expansion
length.
The expansion length l is the distance along the centerline that it takes for the  ow
to expand from the Mach number of the initial expansion point to the desired exit Mach
number along the centerline of the nozzle. It is important to understand that the expansion
length is not the length of the physical nozzle. This is because the  ow will expand to
the desired Mach number at the centerline upstream of where it becomes fully expanded
at the boundary. Thus, the most downstream point in the initial characteristic mesh is
located at l + 2 along the centerline and the pressure distribution chosen describes the
48
change of pressure between the constant Mach line and this point. At  rst this length can
arbitrarily be selected, but for for a given exit condition, pressure distribution and number
of characteristics there exists a unique expansion length that provides the correct expansion
through the nozzle, while minimizing the nozzle length. This can be found using an iterative
technique or an optimization process.
With the pressures selected by the chosen pressure distribution for each characteristic
point along the centerline the problem has been initialized for the the method of character-
istics
3.1.6 Perform Method of Characteristics
Following the scheme presented above the method of characteristics solution along the
centerline can be determined. This  rst method of characteristics solution will be referred
to as the \kernel".
3.1.7 Check Solution for Viability
The viability of the  rst method of characteristics solution is checked before additional
analysis is to be performed. Checking kernel for viability allows for the solution to be
discarded before trying to perform the second method of characteristics solution. Validity
of the solution as a whole is only possible if both solutions are viable. When the solution is
found to be nonviable the pressure distribution prescribed for the analysis is incorrect.
The contour being designed is to be a nozzle contour that produces shock free  ow.
Due to this goal a viable solution can only occur if none of the like characteristics intersect.
The viability of the solution can be performed easily. The characteristic lines between the
points of the characteristic mesh are treated as line segments. By calculating the location
where corresponding lines with equivalent slopes intersects, it can be determined if the two
characteristic segments cross. The segments cross if the location of the intersection is between
the points in question.
49
Figure 3.10: Viable Kernel solution, MOC mesh contains no crossing characteristics
As the solution unfolds, the characteristics may not initially cross but the grid points
are aligned to result appropriate spacing of characteristic lines, resulting in coalescence. The
result shown in Fig. 3.11 shows a MOC mesh with crossing characteristics. This mesh was
created using the same pressure distribution as that in Fig. 3.10, but with a di erent exit
Mach numbers.
Crossing Characteristics
Figure 3.11: Non-viable method of characteristics solution, crossing like characteristics in
MOC mesh
50
Figure 3.12: Close inspection of a non-viable method of characteristics solution, crossing like
characteristics in MOC mesh
It is important to note that the viability of the entire solution is not solely dependent
on the kernel The prescribed pressure distribution may produce a initial analysis that allows
for none of the characteristics to coalesce. However, this pressure distribution may allow
for the characteristics near the throat to cross. This requires a check on the viability of the
solution after the non-simple region has been formed as well.
To count the number of characteristics that cross the following checks are made:
1. Let x and r be the positions of the characteristic points of the MOC solution. The x
location of the intersection of the line segments conjoining between four characteristic
points can be found by
x =
  
  
  
  
  
  
  
  
  
  
 
xi;j ri;j
xi;j+1 ri;j+1
  
  
  
 
xi;j xi;j+1
  
  
  
 
xi+1;j ri+1;j
xi+1;j+1 ri+1;j+1
  
  
  
 
xi+1;j xi+1;j+1
  
  
  
  
  
  
  
=
  
  
  
 
xi;j xi;j+1 ri;j ri;j+1
xi+1;j xi+1;j+1 ri+1;j ri+1;j+1
  
  
  
 
(3.10)
51
The characteristic lines between these four points intersect in the solution if xi+1;j <
x xi+1;j+1
2. If ri+1;j <ri;j!n = n+ 1
3. If ri+1;j <ri;j and ri+1;j+1 <ri;j + 1!n = n+ 1
4. If xi+1;j <xi;j!n = n+ 1
3.1.8 Analysis Near Throat
With a viable method of characteristics solution carried out along the prescribed pressure
distribution. The analysis of the  ow between the kernel and the line of substantial supersonic
 ow can be analyzed.
Using Sauer?s method a line of constant Mach number is formed, beginning at the wall
of the geometric throat and extending to the centerline of the nozzle. This is the same line
that provided for the location of the  rst characteristic point. The zone between this line and
the  rst characteristic line running from the centerline is considered part of the non-simple
region, see Fig. (3.9). The  ow along the Sauer constant Mach number line is found to be
substantially supersonic, allowing for the method of characteristics solution to be applied to
the region. The line of constant Mach number can be produced using Sauer?s solution to the
transonic region. And can be well modeled by
x =  ( + 1)4  r2 + (3.11)
Through Sauer?s solution the Mach number components U and V are known anywhere
along this line. From these components the  ow angle at any point along the line can be
found. All of the other  ow variables can be found through the isentropic relations from
the Mach number. Although the  ow is known at every location along this line, where
characteristic lines can be launched so that the solution is in uenced by the prescribed
52
Constant 
Mach Number
Viable MOC
Solution
C
L
Figure 3.13: Speci cation of the unde ned region
centerline pressure distribution has not been determined. Having the  ow known along this
line allows for a solution to be found for the non-simple region. Although the line of constant
Mach number is not a characteristic line as has been previously de ned, it will be treated
as one for the solution of the unde ned region. This can be done because the  ow variables
are known along the line. When  nding the solution of the unde ned zone, when a right
running characteristic line intersects the constant Mach number the point of intersection
becomes the characteristic point and a left running characteristic line is launched. The
solution of the unde ned region begins by using the  rst left running characteristic line of
the viable MOC solution as the initial conditions. From the solution the slopes of the right
running characteristics are known at each characteristic point along this line. This allows
the points along this characteristic to be characteristic points in both the initial solution and
the unde ned solution.
53
C
L
New Characteritic
Launched
Figure 3.14: Launching characteristic o of the constant Mach number line
The process for the solution in the unde ned region is the same as the method of
characteristics solution presented with one distinct di erence. When the solution nears the
constant Mach line, a characteristic points must lie along this line. These points can be
found through the intersection of the right running characteristics from the initial solution
and the constant Mach line. They are found through
x =
 
( +1)
2   (  xA +rA) + 1
 
 
r 
 ( +1)
2   ( xA +rA) 1
 2
+ (( + 1) )2 2(  xA +rA)2
 ( +1)
2   
2
(3.12)
where
 = tan(   ) (3.13)
54
At this point the  ow properties are known from the conditions of the constant Mach
number line and a characteristic line can be launched back into the solution space, allowing
the method of characteristics is to be carried out in the usual way.
Figure 3.15: Solution to the non-simple region
3.1.9 Determining The Boundary
In the design of a two-dimensional planar nozzle the region between the kernel and
where the boundary is treated as a simple wave region. The left running characteristics
continue outward and a stream line that matches tangency at each of these characteristics
can be treated as the boundary. This is not the case for axis-symmetric  ow, and the
axis-symmetric method of characteristics must be applied to this region. Once applied the
boundary can then be determined by  ow tangency through the solution or by applying
continuity along left running characteristics as developed by Rao[19].
To form the extended characteristic mesh the initial conditions are based o of the last
right and left running characteristics. Downstream of the last left running characteristic the
 ow is become uniform and axial. Along this characteristic the  ow angle and the Mach
55
number are constant. To initiate the characteristic mesh a number of points are distributed
along the characteristic line and form the initial value points. The other set of initial value
points lies along the last right running characteristic and all needed information is known.
r
x
uniform flow
Figure 3.16: Illustration of boundary formation through extended characteristic mesh
Following the axis-symmetric method of characteristics, the left and right running char-
acteristic lines are propagated from each initial value point. Forming a mesh similar to that
in Fig. (3.17). The  nal initial point placed on the last left running characteristic should be
set to achieve the correct area ratio for the chosen exit Mach number. This is found through
r =
r
A
A (3.14)
x =
q
A
A 
tan( + ) +xc (3.15)
This point is also the point at which the contour will end. Insuring that the exit area
is correct. As with the other portions of the solution the prescribed pressure distribution
must be such that there are no crossing characteristics in the solution. The same process as
above can be used to determine if any crossing occurs.
In order to form the boundary a streamline is traced through this mesh. This can be done
by  nding the position and  ow angle at each intersection point between the boundary and
56
the characteristic mesh. The boundary is allowed to be found through  ow tangency because
the nature of the method of characteristics solution satis es continuity. Thus continuity is
satis ed for any stream line in the  ow  eld provided the correct pressure distribution is
provided. When the expansion ratio on the boundary is met than continuity has been
satis ed for the given contour.
x
1
,r
1
,?
1
x
2
,r
2
,?
2
x
A
,r
A
,?
A
B
B'
L
B
L
Figure 3.17: Determination of the boundary segments
The point B can be found through
xB = rA r1 tan AxA +tan( 1  1)x1tan( 
1  ) tan A
(3.16)
rB = rA + (xB xA)tan A (3.17)
Note that point 1 is always on the inside of the boundary and point 2 is always in
the outside of the boundary. This allows the same scheme to be used when the boundary
intersects a right running characteristic with only changing the sign of  in the characteristic
57
slope term in the above equations. The  ow angle at point B is then found by linearly
interpolating along the length L between 1 and 2. The  ow angle at point B is found
through
 B =  2  1L LB + 1 (3.18)
To meet  ow tangency the segment of the boundary must contain all  ow angles between
B and A. The boundary must then have a slope that is the average slope of B and A. This
changes the slope of line  BA which is currently tan A. The new location of B is then found.
xB0 = rA r1 
1
2(tan A +tan B)xA +tan( 1  1)x1
tan( 1  1) 12(tan A +tan B) (3.19)
rB0 = rA + (xB0 xA)12(tan A +tan B0) (3.20)
These steps can be repeated until there is no change between point B and B0. The
process is then repeated for the next intersection with a characteristic until the boundary is
formed at the throat.
3.2 Pressure Distribution Optimization
A viable axis-symmetric method of characteristics solution can be found using an opti-
mizer to  nd a viable centerline pressure distribution that allows for both the viable solution
to the initial MOC mesh and the MOC mesh in the unde ned region. The use of Genetic
Algorithm/Pattern Search/Particle Swarm Optimization developed by Hart eld, Jenkins
and Albarado[28] and the Evolving Swarm Optimizer developed by Jenkins have been used
e ectively to produce results.
The optimization method of Hart eld, Jenkins and Albarado combines the attributes of
particle swarm strategy of Kennedy and Eberhart with the direct search method technique
of Hooke and Jeeves[28]. The technique has been proven to be successful in matching solid
58
rocket motor performances. The method was implemented because of the complexity and
interdependence of variables within the design space.
The evolving swarm optimizer of Jenkins is a hybrid particle swarm pattern search
optimizer with the addition of the survival of the  ttest ideology used in genetic algorithms.
There are two fundamental di erences between the particle swarm/pattern search optimizers
and genetic algorithms. The particle swarm / pattern search optimizer is computationally
more e cient with the use of smaller generations. Jenkins recommends 20 - 50 members
for a population, whereas binary GA?s are often more e cient with large populations of
members. When smaller populations are used with a GA the results tend to the same value.
The second di erence is that the particle swarm pattern search optimizer may take any
value for the parameters it is allowed to manipulate. Where a GA can only take discrete
values. These values are determined by the a set resolution and the minimum and maximum
value limits coded into binary form. Jenkins showed that the %RMS error can be reduced
drastically using both a particle swarm / pattern search optimizer and the evolving swarm
optimizer. The evolving swarm optimizer has a substantially lower number of calls to the
objective function which allow a solution to be found much quicker.
The direct selection method was chosen for  nding a viable pressure distribution. Given
the chamber conditions, desired exit Mach number, and the number of characteristics that
the solution is to be developed over the optimizer selects a normalized pressure for each
of the centerline points. The solution is then carried out regardless of its viability as the
solution progresses. Once the solution is completed the  tness of the solution is calculated
to determine the viability of the pressure distribution.
For both optimizers the same parameters were used to determine a pressure distribution
associated with a valid MOC solution. The use of these optimizers is well suited to the direct
selection method for determining the pressure distribution. The parameters chosen by the
optimizer are the expansion length and the pressure values at each characteristic point that
lies on the centerline.
59
Parameter Maximum minimum
Length 1 25
P1 0.1 .99
P2 0.1 .99
P3 0.1 .99
P4 0.1 .99
P5 0.1 .99
P6 0.1 .99
P7 0.1 .99
P8 0.1 .99
P9 0.1 .99
P10 0.1 .99
P11 0.1 .99
P12 0.1 .99
P13 0.1 .99
P14 0.1 .99
P15 0.1 .99
P16 0.1 .99
P17 0.1 .99
P18 0.1 .99
Table 3.1: Pressure distribution optimization parameters for a 20 characteristic solution
The parameters shown above are for a 20 characteristic solution. The  rst pressure
point is always selected as 1 and the last point is always selected as 0. The parameters P1
through P18 are selected at random. Excluding the length selection, the physical meaning
to the parameters as they are selected is the percentage of the last parameter selected. this
ensures that the pressure distribution selected is always decreasing in pressure. The pressure
values are then determined by multiplying the succeeding parameter by the current value.
This value is then multiplied by the its succeeding parameter and so on. The result is a
percentage of pressure range from the Sauer?s solution pressure to the pressure that results
in the desired Mach number. With a pressure distribution selected by the optimizer the
method can be performed and the solution evaluated. The normalized pressure distribution
as selected by the optimizer has the form shown in Fig. 3.18. When these pressure points
60
Normalized Expansion Length
No
r
m
al
ize
d P
r
essu
r
e
Figure 3.18: Normalized pressure distribution selected from optimizer
are evaluated through
p = (pe ps)Pn +pe (3.21)
where pe is the exit pressure required to achieve the desired exit condition and ps is the
pressure found at the intersection of the Sauer constant Mach line and the centerline. The
resulting static pressure distribution along the centerline is shown in Fig. 3.19
61
Expansion Length
S
t
a
t
i
c 
P
r
e
ssur
e
,
 B
ar
Figure 3.19: Resulting static pressure distribution
3.2.1 Fitness Algorithm
The  tness algorithm analyzes the viability of the pressure distribution through ana-
lyzing the solution. For the producing a nozzle contour for a given exit Mach number the
 tness algorithm contains two parts. The  rst algorithm, shown in Fig. 3.20, analyzes the
initial MOC solution for viability and the second algorithm, shown in Fig.3.21, analyzes the
entire solution given that the initial solution is viable.
In addition to the above evaluation, the addition of solution smoothing terms were
added to the  tness evaluation. These functions allowed for better scoring solutions with
a more evenly distributed solution mesh and smoother transition along characteristic lines.
The resulting solutions are as far away from coalescing solutions as possible. The smoothing
 tness for the left running characteristics is
62
Form Kernel Mesh 
Find Intersection of 
Two Segments 
 AB x CD = xi,yi 
xC<xi>xD n=n+1 g=y
A-yC 
Yes 
yA<yB 
yC<yA  
& 
 yD<yB 
yC<yA 
 &  y
D<yB 
n=n+1 
g=xA-xC 
Yes 
?=???+ ? 
Figure 3.20: Evaluation of  tness through intersection of characteristics
d =
ncX
i
  
  xi+2 xi 2xi+1 xi
xi+2 xi
  
  (3.22)
The absolute value function, Eq. 3.22, allows for more evenly distributed characteristic to
have a lower  tness evaluation, where the x is the x location of the  nal characteristic point in
the kernel. Evaluation of the slope between characteristic points allows for a better transition
between points in the solution separating like characteristics further. The smoothing function
63
that evaluates the characteristic slopes is Eq. 3.23
if xC <xA then
mi =jxC xAj (3.23)
m =
X
mi (3.24)
The smoothing functions are added to the  tness value to allow their in uence.
Secondary  tness can only begin if there are no crossing characteristics from the kernel
solution. The program also skips the solution of the non-simple region if crossing character-
istics have been found; eliminating the processing time required to evaluate the non-simple
region.
Viable Kernel Found 
Evaluate Non-Simple Region 
NaN or Complex f=1.0e7 Yes 
Form Boundary 
Crossing   Characteristics 
Evaluation 
Nullify Pressure 
Distribution 
Check Error f=f+error 
Check Length f=f+length 
Yes 
Yes 
Figure 3.21: Evaluation of  tness for non-simple region and boundary
64
The  tness factors are weighted so that the optimizer will tend towards the most im-
portant factors  rst. These weights are presented in Table 3.2. The optimizer goals always
prefer solutions that have no crossing characteristics as its main objective. It was found that
a viable prescribed pressure distribution will allow for a solution that has no crossing charac-
teristics no matter the expansion length. Thus this length is found by meeting requirements
by the boundary, and continuity. The expansion length that satis es these requirements is
a unique solution for a selected pressure distribution however the expansion length is not a
unique solution for a desired exit Mach number. When nozzle length is of interest the use
of the optimizer also allows for the minimization of the expansion length.
The solution must also contain only real numbers. To eliminate in numerical issues that
could result in a solution with any variables being read as NaN (not a number) or complex
will result in a  tness value that is far greater than the desired result. This situation can
occur in the random selection of optimizer parameters.
Parameter Weight
Crossing in Kernel 3
Crossing in Non-Simple Region 3
Crossing in Boundary Region 1
Area Ratio Error 10
Expansion Length 0.01
Any NaN in solution 1E7
Any Complex in solution 1E7
Table 3.2: Fitness weights
3.3 Uniform Flow Nozzles
The following nozzle contours are nozzles designed for uniform  ow at a speci ed exit
Mach number. Table3.3 are the directly chosen expansion lengths and pressure points for
the given nozzles. For the following contours the chamber conditions and throat radius for
each nozzle was held constant. Figure 3.26 compares the contours.
65
Normalized Pressures for Viable Solutions
Mach 2.0 3.0 4.0 5.0
Length 1.1050 2.7463 5.8953 10.9632
P1 1.000 1.000 1.000 1.000
P2 0.9572 0.7943 0.6428 0.4576
P3 0.9034 0.6642 0.4145 0.2187
P4 0.8325 0.5655 0.2809 0.1047
P5 0.6898 0.4610 0.1836 0.0622
P6 0.5668 0.3494 0.1184 0.0372
P7 0.5136 0.2879 0.0856 0.0241
P8 0.4345 0.2273 0.0576 0.0163
P9 0.3774 0.1605 0.0429 0.0106
P10 0.3233 0.1305 0.0311 0.0076
P11 0.2677 0.1028 0.0211 0.0058
P12 0.1940 0.0737 0.0154 0.0039
P13 0.1664 0.0547 0.133 0.0031
P14 0.1350 0.0395 0.0107 0.0022
P15 0.1152 0.0291 0.0069 0.0017
P16 0.0958 0.0262 0.0056 0.0011
P17 0.0832 0.0230 0.039 0.0008
P18 0.0762 0.0201 0.0034 0.0007
P19 0.0561 0.0142 0.0028 0.0006
P20 0.000 0.000 0.000 0.000
*Speci c heat ratio,  = 1:23
**Pressures Values:P = (Pe Pt)Pn +Pe.
Table 3.3: Directly prescribed normalized pressure points
66
Figure 3.22: Mach 2 Nozzle Contour solved with 20 characteristics,  = 1:23
Figure 3.23: Mach 3 Nozzle Contour solved with 20 characteristics,  = 1:23
67
Figure 3.24: Mach 4 Nozzle Contour solved with 20 characteristics,  = 1:23
Figure 3.25: Mach 5 Nozzle Contour solved with 20 characteristics,  = 1:23
68
M=5
M=4
M=3
M=2
Figure 3.26: Comparison of the above contours
69
3.4 Nozzle Truncation and Performance Evaluation
Nozzles designed for uniform exit  ow can become excessively long. In order to design
a nozzle that produces acceptable performance without these lengths a perfect nozzle can be
truncated. The exit velocities of a truncated nozzle however will not be uniform or fully axial.
The thrust of the truncated nozzle cannot be simply calculated from Eq.(1.3). Instead the
mass  ow must be integrated along a characteristic to account for the  ow being non axial.
The integrated mass  ow is constant at any cross section of the nozzle, however the mass
 ow is not uniform across this cross section. Depending on the location of the truncation
point the fully axial  ow may only occur at the centerline.
The mass  ow along a characteristic can be found through integrating the mass  ow
along a characteristic line.
_m =
Z
 V
 
cos  sin tan( + )
 
rdr (3.25)
Similarly the thrust along a characteristic line can be found through
Ft =
Z  
P + V 2
 
cos  sin tan( + )
  
2 rdr  r2Pa (3.26)
The above thrust equation only accounts for the the thrust that is produced by a nozzle
truncated where a left running characteristic intersects the nozzle boundary. This does not
allow for generic truncation of the nozzle. If the desired truncation length is at a point that
does not lie on this intersection the calculated thrust would be not include the portion of
the boundary that is between characteristics. A scheme can easily be made to calculate the
the thrust associated with the portion of the boundary that does not end at a characteristic
intersection.
By de nition, the trust is the sum of all pressure forces along the boundary of thrust
chamber. Using this the remaining portion of the thrust can be found by integrating the
pressure along the remaining boundary.
70
Ft =
Z
PdA (3.27)
In terms of the nozzle radius, and linearly interpolating the pressure along additional
boundary section the thrust of the segment becomes
Ft =
Z re
ri
2 r(Pi kr)dr (3.28)
where k is the slope of the nozzle segment
This integral can easily be determined numerically and the thrust on this segment of
the nozzle can be calculated directly through
F = Pi (r2e r2i) 2 k3 (r3e r3i) (3.29)
The total thrust for the truncated nozzle is found by calculating the characteristic thrust
along the  rst characteristic upstream of the truncation point. The additional thrust from
the remaining boundary is then added to this to form the total thrust. With this method
the nozzle may be truncated at any distance and the trust can be calculated.
With the thrust calculated the speci c impulse can be found in the usual way as pre-
sented in Eq. (1.2).
3.4.1 Optimization of Truncated Nozzles
Optimization of a truncated perfect nozzles can be found in multiple ways. One such
way would to allow for an external optimization routine to select an exit Mach number for
a uniform exit  ow nozzle, a optimized pressure distribution would then be found through
an internal optimization routine. The external optimization would grade the  tness of the
solution by maximizing the thrust, or thrust coe cient of the this contour at a selected
truncation point.
71
This method for optimizing for maximum thrust at a given truncation length would
produce the optimum solution, however the approach would have enormous computational
expense due to to the shear number of calls to the objective function. The wall clock time
to  nd the solution would also be enormous, because the design space would be very large.
A more intelligent approach to this solution must be made to reduce both computational
expense and wall clock time to a solution. This can be done by mapping the thrust coe cient
at a given truncation length using pre-speci ed perfect contours.
M=2
M=3
M=4
M=5
Truncation Point
Figure 3.27: Example of contours that can be used for mapping thrust coe cient at a given
truncation length
This provides insight to the design space that contains the optimum case for a given
length. The exit Mach number for the perfect contour that results in this optimum solution
can then be found. Using the process above the MOC solution and contour can then be
found.
72
After the initial process of  nding the data to make such performance maps this approach
provides a very inexpensive way of  nding the optimum solution for a given truncation length.
Figure 3.28: Coe cient of thrust versus truncation length for shock free contours
73
Chapter 4
AXSMOC Nozzle Design with Chemical Kinetics
In order to design a nozzle that takes into consideration the e ect of the working  uid
being a nonequilibrium chemically reacting mixture the axis-symmetric method of character-
istics approach described in Chapter 3 can be coupled with the  nite reaction rate analysis
of the chemical kinetics. As was explained in Sec 2.4, analysis of  nite rate chemistry is not a
trivial undertaking; however, the nature of the method of characteristics allows for the MOC
solution to be coupled with such an analysis. This chapter will show how the the chemical
kinetics can be coupled with the AXSMOC and provides a method for implementing the
the coupled approach to existing isentropic solutions. These solutions are chemical kinetics
adjusted axis-symmetric method of characteristics solutions.
4.1 Combining Chemical Kinetics with AXSMOC
The axis-symmetric method of characteristics solutions presented thus far have been
completed considering isentropic  ow. When  nding a axis-symmetric method of character-
istics solution while considering the working  uid a nonequilibrium chemically reacting gas
mixture the computational expense is increased greatly due to the presence of  nite rate
chemistry.
The governing equations for a steady multidimensional adiabatic, inviscid  ow that
performs no work and has negligible body forces can be applied with two new considerations.
 The enthalpy of the system is the absolute enthalpy which is comprised of the sensible
enthalpy and the energy of formation
74
 The law of conservation of mass must be applied to each chemical species involved in
the nonequilbrium reactions
The governing equations are then
r 
 
 ~V
 
= 0 (4.1)
 d
~V
dt +rp = 0 (4.2)
d
dt
 
h+ V
2
2
 
= 0 (4.3)
r 
 
 i~V
 
=  i i = 1;:::;n (4.4)
Where  is the source function as presented in Eqn. 2.43
h =
nX
i=1
Ci (Ci(hi +h i) (4.5)
Noting that i denotes the individual chemical species involved, also note that hi is the
sensible enthalpy and h i is the energy of formation of the species.
These governing equations can be arranged to form a set of characteristic equations and
compatibility equations.
The relationship Ci =  i= can be substituted into 4.4 to form
Cir ( ~V) + ~V  rCi (4.6)
75
Substituting4.1 into the above equation will form the continuity equation in terms of global
continuity as well as speciese continuity.
 ~V  rCi =  i =  dCidt (4.7)
This gives way to
 uCix + vCiy =  i (4.8)
The momentum equation can be rewritten as
 ~V  d
~V
dt +
~V  rp = 0 (4.9)
Rewriting the above momentum equation so that the pressure term is a time derivative
yeilds
 ddtV
2
2 +
dp
dt = 0 (4.10)
Substituting the energy equation into the form of the momentum equation above forms
 dhdt  dpdt = 0 (4.11)
Due to the involvement of change in chemical composition the enthalpy derivative term
of the above equation must rewritten to include the species speci c enthalpy change per
time. With this substitution the above becomes can be rearranged to become
dp
dt =  
nX
i=1
(hi +h i) dCidt + 
nX
i=1
Cidhidt (4.12)
Note that the enthalpy of formation is constant for each species resulting in dh i=dt = 0
The sensible enthalpy hi is de ned as
76
hi =
Z T
To
cpidT (4.13)
The concentrations derivative can be rewritten in terms of the source function,  
 i =  dCidt (4.14)
resulting in the equation
dp
dt =
nX
i=1
(hi +h i) i + 
X
i=1
nCicpidTdt (4.15)
Substituting the species speci c equation of state.
p =  T
nX
i=1
CiRi (4.16)
Yields
dP
dt =
nX
i=1
(hi +h i) i + 
X
i=1
nCicpidTdt (4.17)
Noting that for a mixture of gases the speci c gas constant for the mixture is found through
R =
nX
i=1
CiRi (4.18)
The constant pressure speci c heat is found through
cp =
nX
i=1
Cicpi (4.19)
The resulting speci c heat ratio for the mixture is then found in the usual way
 = cpc
p R
(4.20)
77
The constant pressure speci c heat can be rewritten in terms of the speci c heat ration and
the mixtures speci c gas constant as
cp =  R  1 (4.21)
Using the mixture quantities and the de nition of sound speed Eqn. [4.17] can be rewrit-
ten as
dp
dt a
2d 
dt =
nX
i=1
[ RT (  1)(hi h i)] i (4.22)
For convenience let
  
nX
i=1
[ RT (  1)(hi h i)] i (4.23)
The governing equations for steady axis-symmetric nonequilibrium chemically reacting
gases become
 ux + vr +u x +v r +  vr = 0 (4.24)
 uux + vur +px = 0 (4.25)
 uvx + vvr +pr = 0 (4.26)
upx +vpr a2u x a2v r =  (4.27)
 uCix + vCir =  i (4.28)
With the governing equations containing the source function for  nite rate chemical
kinetics the characteristic and compatibility equations can be derived. The approach used
here for deriving the characteristic and compatibility equations is the approach outlined by
Ferri[22] and Zucrow[4].
78
Multiplying the governing equations above by an unknown factor  1:::4 to the global
governing equations and summing them yields
 1( ux + vr +u x +v r +  vr ) +:::
 2( uux + vur +px) +::: (4.29)
 3( uvx + vvr +pr) +:::
 4(upx +vpr a2u x a2v r  ) = 0
Factoring the derivative terms of u, v, p, and  the above equation can be rewritten in
the form
(  1 + u 2)
 
ux +  v 2  
1 + u 2
ur
 
+ ( u 3)
 
vx +   1 + v 3 u 
3
vr
 
+( 2 +u 4)
 
px +  3 +v 4 
2 +u 4
pr
 
(4.30)
+(u 1 a2u 4)
 
 x + v 1 a
2v 4
u 1 a2u 4 r
 
+  1 vr   4 = 0
Assuming that u, v, p, and  are continuous and that the slope of the characteristic
curves is dr=dx then the derivative terms can be written as
du
dx = ux +
dr
dxur (4.31)
du
dx = vx +
dr
dxvr (4.32)
du
dx = px +
dr
dxpr (4.33)
du
dx =  x +
dr
dx r (4.34)
From Eqn. (4.31) the slopes of the characteristic can be rewritten into 4 di erent and
equal expressions
79
dr
dx =
v 2
 1 +u 2 =
 1 +v 3
u 3 =
 3 +v 4
 2 +u 4 =
v 1 a2v 4
u 1 a2u 4 (4.35)
Yielding the compatibility equation in unknown form
 ( 1 +u 2)du+ u 3dv + ( 2 +u 4)dp+u( 1 a2 4)d + 1 vr dx  4 dx = 0 (4.36)
This compatibility equation is valid along each of the characteristic curves. However the
characteristic equation dr=dx has yet to be determined. By the elimination of the unknown
parameters in the system of equations in Eqn. (4.35) the characteristic equation can be found.
Solving each expression in Eqn. (4.35) for  1 to  4 the coe cient matrix for the system can
be formed
  
  
  
  
  
  
 
dr=dx udr=dx v 0 0
 1 0 udr=dx v 0
0 dr=dx  1 udr=dx v
udr=dx v 0 0  a2(udr=dx v)
  
  
  
  
  
  
 
(4.37)
To form a solution to the unknown parameters that is nonzero, the determinant of the
above coe cient matrix must be equal zero. This results in
(udrdx v)2
"
(udrdx v)2 a2
 
1 +
 dr
dx
 2!#
= 0 (4.38)
De ning   dr=dx and substituting into the above and letting S = u  v the above can
be rewritten as a fourth-order algebraic equation
S2 S2 a2(1 + 2) = 0 (4.39)
80
The fourth-order equation can is expressed above as the product of two second order
factors. The roots of this fourth-order equation be come the characteristic equation. The
 rst root is repeated and is found to be the di erential equation to a streamline.
 1;2 = drdx = v=u (4.40)
The other root two roots become the equation of for the Mach lines
 3;4 = (drdx) = tan(   ) (4.41)
It is important to note that the characteristic equations are identical for non-equilibrium
chemically reacting  ow as for isentropic  ow. However the compatibility equations are not.
Also the species continuity equation Eqn. (4.28) can be placed directly into characteristic
form. Allowing the concentrations to be found along a stream line. This is done in the
same numerical approach as with a one-dimensional non-equilibrium chemically reacting gas
calculation. Rearranging Eqn. (4.28).
 u(Cxi + vuCri) =  i (i = 1:::n) (4.42)
Rewriting the above equation to the following di erential form
dCi
dx =
 i
 u (4.43)
which can be applied along a stream line.
To form the compatibility equations the unknown parameters must be eliminated from
Eqn. (4.36). From the solution to the characteristic equations, the quadratic factor is zero
when following Mach lines. Rearranging and substituting this into the system of equations
in Eqn. (4.35). Reducing each equation to its simplest form results in
81
 2 = tan(   ) 3 (4.44)
 1 = (utan(   ) v) 3 (4.45)
 4 =  3
 1 + tan(   )2
(utan(   ) v)
 
(4.46)
 4 =  3 (utan(   ) v)a2 (4.47)
Substituting equations 4.44, 4.45, 4.47 into the unknown compatibility equation 4.36.
yields the compatibility equation in velocity component form
( v)du ( u)dv+
 
tan(   ) u(u tan(   ) v)a2
 
dp (u tan(   ) v)
  v
r  
 
a2
 
dx = 0
(4.48)
The compatibility equation can be written in the  nal form as
cot dVV  d  sin sin cos(   )dxr +   V 2a
f
dx
cos(   ) = 0 (4.49)
As can be seen from the above equation. The new compatibility equation is the same
as the compatibility equation for the axis-symmetric method of characteristics for isentropic
 ow with the addition of the source term which includes all the changes in chemistry that
a ect the solution. As shown earlier the characteristic equations remained the same.
Discretizing the compatibility equations yields the  nite di erence equations including
the non-equilibrium chemically reacting term
Solving for vc through the elimination of  C as was done before. The only real change
is that $ and M are adjusted to include the  nite rate chemistry term.
82
Method of Characteristic Equations
Isentropic Flow Non-equilibrium Chemically Reacting Flow
Characteristic Equations:
dr=dx = tan(   ) dr=dx = tan(   )
Compatibility Equations:
cot dVV  d  sin sin cos(   ) dxr = 0 cot dVV  d  sin sin cos(   ) dxr +   V2a
f
dx
cos(   ) = 0
Table 4.1: Method of Characteristics Equations
vc = 1cot A
va +
cot B
vB
fcot A(1 + $(xc xa))(1 + M(xc xB)) + B  Ag (4.50)
where
$ = tan A sin A sin Ay
A cos( A + A)
   
Av2A cos( A + A)
(4.51)
and
M = tan B sin B sin By
B cos( B  B)
   
Bv2B cos( B  B)
(4.52)
In turn the  ow angle of the characteristic may be obtained by solving the compatibility
equations for  c yeilding:
 c =  A + cot A
 V
c VA
VA
 
 $(xc xA)) (4.53)
The method may be performed in the same way as with the isentropic  ow case however
the  nite rate analysis must be performed throughout the solution simultaneously to update
the chemical species source function as the solution progresses.
To analyze the chemically reacting  ow using the intersection of a stream line with the
left and the left running characteristic instead of only following the left running and right
83
Streamline
Left Running
Characteristic
 
Solution Point
Right Running 
Characteristic
B
A
C
D
Figure 4.1: Solution can be developed through streamlines[4]
running characteristics through the solution, using the centerline pressure distribution and
the Sauer solution constant Mach number line as initial conditions.
The  nite di erence equations used to analyze the solution with chemistry now add the
streamline term for updating the species source function
yC yA
xC xA = tan( A + A) (4.54)
yC yB
xC xB = tan( B  B) (4.55)
yC yD
xC xD =
tan( D + C)
2 (4.56)
yD yB
xD xB =
yB yA
xB xA (4.57)
84
The location of point D is found after the location of C, and the properties at D
are found through linear interpolation. The location of point D provides the information to
calculate a length di erential between D and C the  nite chemistry scheme is then performed
along this di erential length allowing for the species source function to be updated at C.
Finding the streamline DC is very similar to the method used in determining the boundary.
However DC is the streamline that is associated with the intersection of the characteristic
lines, the same streamline is not being followed throughout the solution.
For the  rst iteration for  nding the location of point D a prediction of the  ow angle
by averaging the  ow angles of A and B can be made. Point D can be readily found by
iteration through
xD = yC yA 
1
2 (tan C + tan D)xC + sxD
s 12(tan C + tan D) (4.58)
and
rD = rC + (x xC)12(tan C + tan D) (4.59)
where
s = rA rBx
A xB
(4.60)
The  ow properties and mass fractions are now linearly interpolated between  AB. The
chemical kinetics can now be analyzed along the stream line with the length of  DC being
the di erential length of the kinetics problem. Because the line segment represents the
stream line  DC and the  ow is supersonic chemical kinetics can be analyzed using the
one-dimensional analysis as presented in Section 2.4.
4.2 Nozzle Contour Design Process with Chemistry Included
The approach in design of an axis-symmetric nozzle coupled with  nite rate chemistry
combines the the three-zone method described with in Chapter 3 with a two-zone down
85
stream approach to adjust for the chemical kinetics. The addition of the  nite rate chemistry
allows for the speci c heat ratio to be an independent variable instead of a constant.
4.2.1 Nozzle Contour Design Process
The following describes the process for axis-symmetric method of characteristics coupled
with  nite rate chemistry nozzle design. Providing the chamber conditions, propellant type,
mixture ratio, converging curvature, and isentropic exit uniform  ow exit Mach number.
1. Select Converging Section Geometry
2. Perform One-Dimensional Finite Rate Analysis
3. Form Isentropic Nozzle Solution
4. Perform MOC Coupled with Chemistry
5. Form Boundary
4.2.2 Converging Section Analysis
The converging section of the nozzle is analyzed with one-dimensional non equilibrium
chemically reacting as described in the above section. To initialize the solution the compo-
sition of the combustion exhaust gases must be known. To  nd composition, an equilibrium
composition calculation must be performed, using any number of methods to minimize Gibbs
free energy[34]. A modi ed version of the Chemical Equilibrium with Applications (CEA)
program developed by Gordon and McBride was used for this calculation[16][17].
The equilibrium calculation can be used to initiate the procedure using the assumption
that the gases are reacting at or very near equilibrium in the combustion chamber. This
analysis begins directly down stream of the combustion chamber at the start of the nozzle?s
converging section. To form the required pressure distribution, the converging geometry
must be found.
86
Converging Length
Ra
di
us
 ?
 ?
?  = 26.25?
r
sub
= 2.0
r
i     
= 3.0
r
sub
ri
l
Figure 4.2: Converging section geometry with dimensions normalized by throat radius
The design of the converging section has been simpli ed for this development. The
initial turning radius and the subsonic radius of curvature, along with the turning angles of
both curvatures, are held constant for any case. These curvatures were modeled after the
converging section of the Pratt and Whitney RL-10 [13] and have also been used in TDK
performance calculations[7].
To form the required geometry to produce choked  ow at the throat the length l is
allowed to change so that the correct expansion ratio is found. The expansion ratio is
determined by assuming that the speed of the gases at the beginning of the converging section
have a non-zero velocity. For the present development the converging section is assumed to
have Mach 0.15  ow at the beginning of the section. This value can be derived from the
design of the fuel injection system and combustion chamber design. For this development
the expansion ratio for the converging section is only dependent on propellant choice through
the speci c heat ratio. From the equilibrium calculation the speci c heat ratio and assumed
Mach number are used in the isentropic relation to determine the initial pressure ratio, and
87
expansion ratio. Because, the length l is not aligned axially, it is more convenient to evaluate
total converging section length.
x = ((rsub +ri) sin( ) +
pA
i 1 (rsub +ri) (1 cos( ))
tan (4.61)
As the required expansion ratio increases the required length l must increase to expand
the initial area. This directly results in a change of axial length of the converging section.
With the length of the converging section known the area of the converging section and its
derivative can be found easily through the following set of equations[7]
When x<x1 +ri +sin( )
A =
p
Ai ri
 
1 
p1 (x x
1)2
ri2
!2
(4.62)
dA
dx =
 2(x x1)
(r2i  (x x21))  (1=2))
pA (4.63)
When x+risin <xandx< rsubsin 
A =
p
Ai ri (1 cos ) (x x1 ri sin tan )2 (4.64)
dA
dx = 2
pAtan( ) (4.65)
When  rsubsin( )) <x(i) andx(i) < 0
A =
 
1 +rsub
 
1 
s
1 x
2
r2sub
!!2
(4.66)
dA
dx =
2xp
r2sub x2
pA (4.67)
88
As described by Nickerson, et. al.[7][35], the pressure distribution associated with this
geometry and chamber conditions can be found. The pressure ratio and its derivative can
be determined by using the average equilibrium expansion coe cient and the area of the
converging section for every position in x.
The average expansion coe cient can be found through iteration upon the following
formula
Cn+1e = Cne +
2
Cne +1
Cne
Cne 1  P e
Pc
2
Cne +1
Cne
Cne 1 1
Cne 1
 
1
Cne 1ln
 
2
Cne +1
 
+ CneCn
e +1
 (4.68)
Where P e is the pressure if  ow was in chemical equilibrium to the throat. A good
starting point for the iteration is 1.2 [7]
With the average expansion coe cient found the pressure ratio at each location of x
can now be found. The pressure ratio for the given location is iterated upon through the
expression, where P is the ratio P=Pc
Pj+1 = Pj + 2
(
Ce 1
Ce
 
1 P
Ce 1
Ce
j
  1
P
 1
Ce
j  
2
CeP
 1
j
) 1
 
8>
>><
>>>:
2
66
4
Ce 1
2
2
Ce+1
Ce+1
Ce 1
P
2
Ce
j
 
1 P
Ce 1
Ce
j
 
3
77
5
 12
A 1
9
>>>
=
>>>; (4.69)
The geometry illustrated in Fig. 4.2, was developed for a LOX-Hydrogen engine oper-
ating with a chamber pressure of 27.186 Bar. This is similar to that of the RL-10. The
pressure distribution prescribed for this engines converging section is depicted in Fig. 4.3.
With the equilibrium condition, geometry, and pressure distribution determined the
analysis of the  ow including  nite rate reactions may take place as described in Sec. 2.4.1.
Figure 4.6 shows the general procedure in which the chemistry analysis occurs.
89
Figure 4.3: Converging section pressure distribution for LOX H2 with Pc = 27:186 Bar
Find Equilibrium Composition & Combustion Temperature 
Solution Complete? 
Call CAP for Properties at 
Current Temp 
Yes 
Return  
Integrate to Determine Species 
Production Rates 
Update State 
Input:  Fuel/Oxidizer 
Mixture Ratio Chamber Pressure 
Obtain Reaction Rates for 
Equilibrium Species 
Thermodynamic Properties 
No 
Call CEA  
Default to 
C Cl F H N O 
Reaction  Unavailable 
Figure 4.4: Procedure of  nite rate chemistry analysis
90
The  nite rate chemistry analysis begins with a library search for the correct set of re-
actions that can be performed with the species in the equilibrium composition. The reaction
rate data available is listed but not limited to that presented in Appendix D. The  nite rate
chemistry application was developed as a generic tool provided that the correct rate data is
provided in the library. If the application cannot  nd a speci c reaction mechanism card it
will default to the C, Cl, F, H, N, O reaction card. This card contains the elementary reac-
tions of species that typically could be contained in a liquid propellant. Computationally, a
card holding the reaction mechanism with only the available species is preferred.
The propellants in Table 4.2 are some typical liquid propellants that have been used in
rockets. Reaction rate data is provided in Appendix D common propellants.
Propellant Chemistry
Elemental Breakdown Fuel Oxidizer
1 C , F, H, N C2H8N2 F2
2 C, F, H, N, O N2H4 C2H8N2 OF2
3 C, H, O CH4, C3H8, RP-1 O2
4 H, O H2 O2
5 F, H, H2 F2
6 F, H, N N2H4 F2
7 F, H, N, O OF2 N2H4
8 H, N, O H2 N2O4
*Chemical kinetic rates have been included into the program.
Table 4.2: Elemental breakdown of propellant fuel oxidizer pairs
Using the appropriate rates for the reaction mechanism the one-dimensional  ow analysis
of non-equilibrium chemically reacting gases is performed through numerical integration.
The scale of the rates of the chemical reactions compared to the scale of the rates of change of
the spacial and thermodynamic properties is very large. The extreme scale di erence makes
the governing equations a system of sti di erential equations. The numerical integration
technique used is required to be an implicit scheme so that the scheme remains stable and the
solution can be found[7][4][36]. The technique used to numerically integrate the governing
equations is fully outlined in the by Nickerson, Coats, Bartz[35].
91
4.2.3 Isentropic Nozzle Formation
Once the solution has been completed for the converging section, the end result provides
the  ow properties at the geometric throat. At this point the composition of the  uid and
its thermodynamic properties are known. The speci c heat ratio for this point will be used
for the  rst estimate of the speci c heat ratio for the Sauer solution to the constant Mach
number line. The line is formed and one dimensional- ow analysis is performed between the
geometric throat and the intersection of the Sauer constant Mach line with the centerline.
When the solution converges so that the speci c heat ratio is met for the displaced distance
from the geometric throat, this composition and thermodynamic properties are assumed to
be constant along the constant Mach line.
With the speci c heat ratio, location and Mach number known along the Sauer constant
Mach line known an axis-symmetric method of characteristics solution can performed for
the desired Mach number. If the centerline pressure distribution for this desired Mach
number and  is known the method can be carried out quickly, if not the centerline pressure
distribution can be determined and the appropriate solution formed as was done in Chapter
3.
4.2.4 Perform Chemically Coupled MOC
The isentropic nozzle solution will be the basis for the solution coupled with  nite rate
chemistry. To begin the solution it is necessary to use the the  rst set of points downstream
of the Sauer constant Mach line to begin the solution. The location of these points was
determined through the pressure distribution required to achieve correct expansion to the
desired Mach number, using these points helps initiate the solution.
The composition and new speci c heat ratio is now found at each initiation point through
the same means as the would be for an internal point. The solution can now be carried down
stream as described is Section4.1, with the exception of  nding the points along the center-
line. The points on the centerline can be found through symmetry. Here the magnitude of
92
L
C
Figure 4.5: Initial Point Selection
the slopes of the left running characteristic of A is the same as the right running character-
istic of B, but with opposite signs. This is also true for the  ow angles, and the r location of
the point. With this exception included the analysis can be completed for the kernel with
 nite rate chemistry included.
C
L
r
x
A
B
C
 ?
B
 ?
A
Figure 4.6: Points determined along the centerline through symmetry
93
Once the the chemically reacting kernel has been found the chemically reacting boundary
zone can be found. The process is similar to that of the isentropic case. However, the
addition of  nite rate complicates the procedure. The di culty in  nding the solution to
the boundary region is caused by the need to march down a streamline to update the  nite
rate chemistry in uence. Before the region was simple to  nd, assuming that the exit  ow
was uniform and axial, but using those assumptions the solution was marched upstream for
the isentropic case. To initiate the chemically reacting MOC scheme an assumption must be
made about the position of the next right running characteristic, the  ow angle at the  rst
point, and the composition of the working  uid. To initiate the chemically reacting boundary
r
x
Isentr opic Boundary
Chemically Reacting Kernel
i
j
B
A
C
Up
da
t
e
 
F
or
 Ch
em
i
st
r
y
Figure 4.7: Update the boundary region characteristic for chemistry
region, an initial right running characteristic is found using the isentropic boundary region
solution. Using a point outside the isentropic boundary an initial estimate of the  ow angle
is found. The value for  is assumed to be the value of  for the corresponding left running
characteristic at the edge of the kernel. Initializing the right running characteristic in this
way allows the method of characteristics with  nite rate chemistry to be used. The chemically
reacting MOC scheme is used to march downstream updating the location and the chemical
source function for each point along the right running characteristic. The next right running
characteristic is formed in the same way how ever the initial point should be shifted to the
94
next left running characteristic, and the assumed values are taken from the previous right
running characteristic. The process for forming the right running characteristics for the
chemically reacting boundary region is then repeated until ample characteristics have been
formed to produce a solution. This solution being overlayed on the isentropic boundary
produces a solution mesh as in Fig. 4.8 will show a jump in the non physical boundary mesh.
This has no a ect on the solution.
The boundary is formed in the same manner as that of the isentropic  ow nozzle. The
formation of the boundary is described in detail in Section 3.1.9.
Figure 4.8: Chemically reacting solution mesh for Mach 4 nozzle
4.3 E ects of Finite Rate Chemistry on Nozzle Contour
The addition of  nite rate chemistry to the method of characteristics for nozzle design
has visible e ects on the solution. The e ects are more noticeable as the nozzle length
increase.
Figure 4.9 shows the di erence between the isentropic kernel solution and the kernel
solution with the inclusion of  nite rate chemistry. Notice that the characteristic lines and
points have shifted due to the in uence of the chemically reacting system.
The speci c heat ratio is no longer constant through out the nozzle. The increase in
speci c heat ratio tends to be fairly linear along the centerline of the nozzle as depicted in
Fig 4.10. For this example, LOX Hydrogen was the propellant choice.
95
Figure 4.9: Comparison of isentropic kernel to kernel including  nite rate chemistry
Figure 4.10: Speci c heat ratio across expansion length
96
Referring to Eq 2.1, as  increases the area ratio to produce a given exit condition
decreases. This can be seen in the contour generated including  nite rate chemistry. Fig-
ure 4.11 shows the e ect of  nite rate chemistry on the contour of nozzle that was designed
for Mach 4  ow at the exit. The turning rate of the contour decreases along the length of the
nozzle, as the speci c heat ratio increases. This produces a contour that follows the same
contour as the isentropic case initially but the exit area is decreased.
Figure 4.11: Comparison of isentropic Mach 4 nozzle to nozzle including  nite rate chemistry
This result shows that the isentropic assumption does not have a drastic e ect on the
design of truncated nozzles. The contours are only slightly changed from the throat to about
10 time the radius of the throat in length. From this point on the isentropic assumption
begins to fail and the inclusion of  nite rate chemistry allows for a more accurate contour
to produce a shock free nozzle.
97
Chapter 5
Conclusions and Suggestions for Future Work
A method of designing axis-symmetric nozzles with included  nite-rate chemical kinetics
has been developed. The method allows for the physical boundary of the nozzle to be
determined based on a desired uniform exit Mach number for the isentropic case. The
solution is then coupled with and adjusted for  nite rate chemistry.
The contours designed in this way can have excessive and undesirable lengths, for use
on rocket engines. A method for analysis of truncated nozzle performance was explored.
The contour that produces maximum thrust for a given truncation length can be found.
A method of simplifying the optimization procedure by reduction of the design space was
developed through the use of performance maps. The performance maps can be extended
to include a larger number of contours to reduce the design space for the optimum nozzle
even further. The truncated nozzle analysis method is not limited to that of isentropic
 ow nozzles. The process can be extended to contours designed with  nite rate chemistry
included.
Improvement in contour geometry can be made through the incorporation of boundary
layer e ects could be added to the methodology developed. This would allow a more accurate
physical contour geometry. A simple method for adding a boundary layer to the current
development would be to displace the nozzle wall by a distance equal to the estimated
boundary layer thickness.
The current development?s analysis of the  nite rate chemistry in the subsonic portion
of the nozzle is limited to a one-dimensional analysis. This could be increased to a higher
degree of freedom. In doing so the analysis would consider striated  ows. The current
development uses a one dimensional analysis becomes the distances traveled by streamlines
98
in the converging section is approximately equal. With the traversed distance being relatively
uniform across the converging section it is a safe assumption that the composition of the
product gases are uniform at the sonic line.
Inaccuracies in the  nite rate analysis are largely due to provided reaction rate data. The
rate data may only be as good as the experimental process to obtain such rates. Techniques
for the measurement of reaction rates may become more accurate, and the reaction rate and
third-body in uence data should be updated accordingly.
The approach to the method of characteristics coupled with  nite rate chemistry pre-
sented could be extended to analysis of the exhaust plume. This method would allow for
the plume composition to be found. This can be useful for missile identi cation, or optimiz-
ing combustion and expansion of the gases to produce desired composition. Such a process
would allow reduction of unwanted exhaust chemicals.
The contours designed with the addition of  nite rate chemistry provide only a small de-
viation from the isentropic contour, when the truncation length is kept small. For propulsive
nozzles, the truncation length that would typically be selected will be within the range that
would see very little di erences when designed with or without chemistry. The additional
e ort may not be worth the slight changes in contour that are formed by the additional
analysis. However, if the composition of the exhaust gases needs to be known, as with plume
analysis, it is imperative to use the analysis with chemical kinetics included.
99
Bibliography
[1] Sutton, G. P. and Biblarz, O., Rocket Propulsion Elements, John Wiley & Sons, Inc.,
seventh ed., 2001.
[2]  ?Ostlund B. Muhammad-Klingmann, J., \Supersonic Flow Separation with Application
to Rocket Engine Nozzles," Applied Mechanics Reviews, ASME, Vol. 58, May 2005.
[3] Guentert, E. C. and Neumann, H. E., \Design of Axisymmetric Exhaust Nozzles By
Method of Characteristics Incorporating A Variable Isentropic Exponent," NASA TR
R-33, 1959.
[4] Zucrow, M. J. and Ho man, J. D., Gas Dynamics, Vol II., John Wiley and Sons, 1977.
[5] Sauer, R., \General Characteristics of the Flow Through Nozzles at Near Critical
Speeds," National Advisory Committee for Aeronautics, June 1947.
[6] Dunn, B., \Fuel Table," Tech. rep., 1997 [retrieved 10 March 2012].
[7] Nickerson, G. R., Dang, L. D., and Coats, D. E., \Engineering and Programming Manual
: Two-Dimensional Kinetic Reference Computer Program," NASA CR-178628, April
1985.
[8] Anderson, J. D., Introduction to Flight, McGraw Hill.
[9] Technologies, G. L., \Brief History of Rockets," http://www.grc.nasa.gov/WWW/K-12/
TRC/Rockets/history_of_rockets.html, Feb. 2010 [retrieved 24 January 2012].
[10] NASA, \Dr. Robert H. Goddard, American Rocketry Pioneer," http://www.nasa.gov/
centers/goddard/about/history/dr_goddard.html, June 2009 [retrieved 24 January
2012].
[11] Edwards, T., \Liquid Fuels and Propellants for Aerospace Propulsion: 1903-2003,"
Journal of Propulsion and Power, Vol. 19, No. 6, 2003, pp. 1089{1105.
[12] Huzel, D. K. and Huang, D. H., Modern Engieering for Design of Liquid-Propellant
Rocket Engines, AIAA.
[13] Pratt and Whitney, \RL10 Rocket Engine A National Historic Mechanical Engineering
Landmark," American Society of Mechanical Engineers.
[14] Pratt and Whitney, \Space Shuttle Main Engine," Tech. rep., [retrieved 10 March 2012].
100
[15] Williams, F. A., Combustion Theory, The Benjamin/Cummings Publishing Company,
Inc., 1985.
[16] Gordon, S. and McBride, B. J., \Computer Program for Calculation of Complex Chem-
ical Equilibrium Compositions, Rocket Performance, Incident and Re ected Shockes
and Chapman-Jouguet Detonations," Nasa sp-273, 1971.
[17] Gordon, S. and McBride, B. J., \Computer Program for Calculation of Complex Chem-
ical Equilibrium Compositions and Applications," NASA RP-1311, 1996.
[18] Anderson, J., Modern Compressible Flow: with Historical Perspective, McGraw-Hill.
[19] Rao, G. V. R., \Exhaust Nozzle Contour for Optimum Thrust," Jet Propulsion, Vol. 28,
No. 6, June 1958, pp. 377{382.
[20] Abbot, M. B., An Introduction to the Method of Characteristics, Thames and London.
[21] Hart eld, R. J. and Ahuja, V., \The Closing Boundary Condition for the Axis-
Symmetric Method of Characteristics," Joint Propulsion Conference, AIAA, San Diego,
CA, 2011.
[22] Ferri, A., \The Method of Characteristics," General Theory of High Speed Aerodynam-
ics, edited by W. Sears, The Macmillan Co., 1951.
[23] Rao, G., \Contoured Rocket Nozzles," Proc 9th, 1958.
[24] Haddad, A. and Mosst, J. B., \Aerodynamic Design for Supersonic Nozzles of Arbitrary
Cross Section," Journal of propultion, 1988.
[25] Wang, B. N. and Ho man, \Calculation of Annular Nozzle Trisonic Flow elds by the
Method of Characteristics," Tech. rep.
[26] Allman, J. G. and Ho man, J. D., \Design of Maximum Thrust Nozzle Contours by
Direct Optimization Methods," AIAA/SAE Joint Propulsion Conference, Vol. 14, Aug.
1978.
[27] Zucrow, M. J. and Ho man, J. D., Gas Dynamics, Vol I., John Wiley and Sons, 1976.
[28] Albarado, K. M., Hart eld, R. J., Hurston, B. W., and Jenkins, R. M., \Solid Rocket
Motor Performance Matching Using Pattern Search/Particle Swarm Optimization,"
Journal of Propulsion and Power, Vol. 47, 2011.
[29] Quarteroni, A., Sacco, R., and Saleri, F., Numerical Mathematics, Springer, 2nd ed.,
2007.
[30] Kushida, R., \Revision of CPIA 246, Section 6.2, Reaction Rate Data," Jpl 383cr-76-
211, 1976.
[31] Cohen, N. and Westberg, K. R., \Chemical Kinetic Data Sheets for High-Temperature
Chemical Reactions," Journal of Physics.
101
[32] Jr., K. W., Advanced Thermodynamics for Engineers, McGraw Hill.
[33] Zehe, M. J., Gordon, S., and McBride, B., \CAP: A Computer Code for Generating
Tabular Thermodynamic Functions from NASA Lewis Coe cients," NASA/TP- 2001-
210959, Oct. 2001.
[34] Anderson, J. D., Hypersonic and High-Temperature Gas Dynamics, AIAA, 2006.
[35] Nickerson, G. R., Coats, D. E., and Bartz, J. L., \Two-Dimensional Kinetic Reference
Computer Program," NASA CR-152999, Dec. 1973.
[36] Radhakrishnan, K., \Combustion Kinetics and Sensitivity Analysis Computations,"
Numerical Approaches to Combustion Modeling.
102
Appendices
103
Appendix A
Sauer?s Solution to Transonic Region
Transonic  ow occurs in a converging-diverging nozzle near the throat. Unlike in quasi-
one-dimensional  ow, the sonic line (line where the Mach number of the  uid is unity at all
points) is not a vertical line spanning the width of the throat. In fact the sonic line is curved.
The region close to this sonic line is the transonic  ow zone and contains both subsonic and
supersonic portions of the  ow. This causes di culties in the mathematical representation
of this region. To describe this zone, the coordinate system used should lie so that the x-axis
is the centerline of the nozzle and the origin is placed where the critical velocity  rst occurs
on the axis, see Fig A.1
Considering ~u;~v as the x and r components of the velocity and a as the sonic velocity
at the throat, the potential equation is written as
@U
@x
 
1 U2   1 + 1V 2
 
+ @V@r
 
1   1 + 1U2 V 2
 
 4 + 1@U@rUV + 
 
1   1 + 1  U2 +V 2 
 V
r = 0 (A.1)
where
U = ~ua , and V = ~va (A.2)
Knowing that the transonic  ow is only in the vicinity of the throat de ne
U = 1 +u; and V = v (A.3)
The quantities of u and v are to be considered small. Substituting Eqn (A.3) into Eqn (A.1)
results in the following approximation to be formed. The intermediate steps are found in
Sauer[5].
( + 1)u@u@x @v@r  vr = 0 (A.4)
The following relations must be de ned in order to continue. These relations take into
account the symmetry around the x-axis
104
r
x
e
h
Fl ow  Direction
Figure A.1: Transonic Flow Zone[5]
? = f0(x) +r2f2(x) +r4f4(x) +:::
u = @?@x = f00(x) +r2f04(x):::
v = @?@r = 2rf2(x) + 4r3f4(x) +::: (A.5)
The primes in the Eq. A.5 represent the  rst derivatives with respect to x of the appro-
priate function. From the expressions de ned the following can be found through Eqn (A.4).
( + 1)(f00 +r2f02 +:::)(f000 +y2f002 +:::) = 2f2 + 12r2f4 +:::+ 2 (f2 + 2r2f4 +:::) (A.6)
Rearranging Eqn (A.6) in order of powers of r and equating coe cients of the individual
powers to zero the following system of equations is obtained.
2(1 + )f2 = ( + 1)f00f000
2(6 + 2  )f4 = ( + 1)(f00f002 +f000f02) (A.7)
The function
f00 = uo(x)
characterizes the velocity distribution along the nozzle axis. This  ow velocity distribution
can be set given
f00 = uo(x) =  x (A.8)
105
Given uo(y) and understanding that the nozzle wall can be assumed to be a curve that is
revolved about the x-axis. This provides that the  ow is completely established in the nozzle.
With this being said Eqn (A.8) can be substituted into the system of equations in Eqn (A.7)
and solved for the functions f2 and f4 yielding the following equations
f2 =  + 14  2x (A.9)
f4 = ( + 1)
2
64  
3 (A.10)
Because the nozzles being considered are three dimensional = 1, substituting Eqn (A.9)
and (A.10) into Eqn (A.5) results in
u =  x+  + 14  2r2 +::: (A.11)
v =  + 12  2xr + ( + 1)
2
16  
3r3 +::: (A.12)
where higher order terms have been neglected. This result are the components of the
velocity as a function of x and r in the nozzle. From (A.11) and (A.12) the sonic line can
be obtained when the following requirement is met
(1 +u)2 +v2 = 1 (A.13)
When this condition met the following curve can be projected through the transonic
region
xK =  + 14  r2k (A.14)
This curve crosses through the axis of symmetry normally at the point x = 0;y = 0. At
this location its curvature is 1
 K =
 + 1
2  (A.15)
Considering where the  ow is tangent to the boundary, the  ow here is thus parallel
to the axis of the nozzle. At this point v = 0, this being so allows for the location of the
minimum area point to be found through
 = xs =  + 18  y2s (A.16)
From the location of the minimum area point and the curve of the sonic line the distance
behind the minimum area point where the  ow becomes sonic can be de ned through
 = (xk xs) = 18( + 1) y2s =  (A.17)
106
With the minimum area point and the sonic point at the wall of the nozzle known the
radius of curvature of the throat can be found. Following the assumption that the curvature
follows a circular arc the radius of curvature can be found to be
 s = 2( + 1) 2y
s
(A.18)
The preceding equations have been presented as if the velocity components of the  ow
 eld are known and the radius of curvature needs to be found. Typically the opposite
occurs. In use of Sauer?s solution in developing an automated nozzle design program neither
is known and no assumptions could be made about the  ow  eld beyond reasonable doubt.
This being said a radius of curvature of the subsonic portion of the throat is allowable,
especially with the use of an optimizer. Thus the  ow  eld was determined by an assumed
radius of curvature. Solving a few of the equations so that one of the dependent variables is
the radius of curvature is completed below.
 = sqrt 2( + 1) 
sys
(A.19)
 =  = ys8
r
2( + 1)ys 
s
(A.20)
The use of (A.11) and (A.12) along with the solution to (A.19) the Mach number within
the transonic region can be found.
107
Appendix B
Chemical Equilibrium with Applications
The development of a computer program for the calculation of complex chemical equilib-
rium began in the 1940?s through the NASA Lewis Research Center[16]. Sanford Gordon and
Bonnie McBride developed a series of software applications that allow for such calculation
as well as application of chemical equilibrium  ow such as rocket performance calculations.
The software is readily available at www.grc.nasa.gov/WWW/CEAWeb/. The fortran source
code can be found and modi ed as needed to  t the application at hand.
The approach used by Gordon and MacBride was a free-energy minimization technique,
formed in a generic sense as to beable to calculate the equilibrium compositions of any
speci ed mixture of propellants.
108
Appendix C
Thermodynamic Properties Database
As with applications such as CEA, the development at hand required the availability
of thermodynamic data of all species involved in the calculation performed. NASA Glen
Research Center maintains a database of over 2000 species[33]. The database has been
compiled into the Coe cients and Properties (CAP) program.
The database uses a 9-constant empirical representation of the thermodynamic data
of each species. The thermodynamic coe cients are generated by a least-squares  t to
measured thermodynamic functions of condensed and gas phase species[33]. The database
is contiually updated to re ect new species, improved measurements for current species, and
newer physical constants.
The properties available in the database include
 C p Molar Constant Pressure Speci c Heat
 G (T) Molar Gibbs energy at temperature T
 H (T) Molar Enthalpy
 K Equilibrium Constant
 M Molecular Weight
 S (T) Entropy at temperture T
The program also has the ability of changing the reference state. As well as compare
some of the above properties to the properties at the reference state.
For the development in hand the CAP program was used internally. The reference state
used was 298:15K. The source code can easily be found through www.grc.nasa.gov/WWW/
CEAWeb/. This is the same thermodynamic database used in the version of CEA employed
in the current development.
109
Appendix D
Reaction Mechanisms
The following tables provide the forward reaction rate coe cients and mechanisms that
can be used in this development.
O-H Reactions
Elementary Reaction 3rd Body A N B
1 2H +M  H2 +M M1 6.4E17 1.0 0.0
2 H +OH +M  H2O +M M2 8.4E21 2.0 0.0
3 O +H +M  OH +M M3 1.9E13 0.0 -1.79
4 O +O +M  O2 M4 3.62E18 1.0 0.0
5 O2 +H O +OH 0 2.2E14 0.0 16.8
6 H2 +O H +OH 0 1.8E10 -1.0 8.9
7 H2 +OH H2O +H 0 2.2E13 0.0 5.15
8 OH +OH H2O +O 0 6.3E12 0.0 1.09
* Forward Rate: K = AT NEXP( 1000B=RuT)
Table D.1: Simple reaction mechanism for Hydrogen/Oxygen combustion[7]
Third-Body Ratios of H,O system
Species M1 M2 M3 M4
H 25.0 12.5 12.5 12.5
H2 4.0 5.0 5.0 5.0
H2O 10.0 17.0 5.0 5.0
O 25.0 12.5 12.5 12.5
OH 25.0 12.5 12.5 12.5
O2 1.5 6.0 5.0.0 11.0
Table D.2: Third-body e ciency ratios for Hydrogen/Oxygen
110
Carbon-Hydrogen-Oxygen Nitrogen Reactions
Elementary Reaction 3rd Body A N B
1 2H +M  H2 +M M1 6.4E17 1.0 0.0
2 H +OH +M  H2O +M M2 8.4E21 2.0 0.0
3 O +H +M  OH +M M3 1.9E13 0.0 -1.79
4 O +O +M  O2 M4 3.62E18 1.0 0.0
5 N +O +M  NO +M M5 6.4E16 0.5 0.0
6 N +N  N2 +M M6 3.0E14 0.0 -0.99
7 CO +O +M  CO2 +M M7 1.0E14 0.0 0.0
8 O2 +H O +OH 0 2.2E14 0.0 16.8
9 H2 +O H +OH 0 1.8E10 -1.0 8.9
10 H2 +OH H2O +H 0 2.2E13 0.0 5.15
11 OH +OH H2O +O 0 6.3E12 0.0 1.09
12 CO +OH CO2 +H 0 1.5E7 -1.3 -0.765
13 N2 +O NO +N 0 7.6E13 0.0 75.5
14 O2 +N  NO +N 0 6.4E9 -1.0 6.25
15 CO +O CO2 0 2.5E6 0.0 3.18
16 CO2 +O CO +O2 0 1.7E13 0.0 52.7
* Forward Rate: K = AT NEXP( 1000B=RuT)
Table D.3: Simple reaction mechanism for CHON combustion[7]
Third-Body E eciency of CHON
Species M1 M2 M3 M4 M5 M6 M7
Ar 1.0 1.0 1.0 1.0 0.8 1.0 1.0
CO 1.5 3.0 4.0 4.0 1.0 1.0 1.0
CO2 6.4 4.0 5.0 8.0 3.0 2.0 5.0
H 25.0 12.5 12.5 12.5 10.0 10.0 1.0
H2 4.0 5.0 5.0 5.0 2.0 2.0 1.0
H2O 10.0 17.0 5.0 5.0 7.0 3.0 1.0
N 1.0 1.0 1.0 10.0 10.0 10.0 1.0
NO 1.5 3.0 4.0 4.0 1.0 1.0 1.0
N2 1.5 3.0 4.0 4.0 1.0 1.0 2.0
O 25.0 12.5 12.5 12.5 10.0 10.0 1.0
OH 25.0 12.5 12.5 12.5 10.0 10.0 1.0
O2 1.5 6.0 5.0.0 11.0 1.0 1.0 25.0
*Argon only present in combustion with air
Table D.4: Third-body e ciency ratios for CHON[7]
111
C,Cl,F,H,N,O Reactions
Elementary Reaction 3rd Body A N B
1 2H +M  H2 +M M1 6.4E17 1.0 0.0
2 H +OH +M  H2O +M M2 8.4E21 2.0 0.0
3 O +H +M  OH +M M3 1.9E13 0.0 -1.79
4 O +O +M  O2 M4 3.62E18 1.0 0.0
5 N +O +M  NO +M M5 6.4E16 0.5 0.0
6 N +N  N2 +M M6 3.0E14 0.0 -0.99
7 CO +O +M  CO2 +M M7 1.0E14 0.0 0.0
8 F +F +M  F2 +M M8 3.25E8 -1.0 -6.38
9 H +F +M  HF +M M9 7.5E12 0.0 -35.13
10 CL+CL+M  CL2 +M M10 2.23E14 0.0 -1.81
11 H +CL+M  HCL+M M11 2.6E13 0.0 -19.9
12 CL+F +M  CLF +M M12 3.0E16 0.5 0.0
13 O2 +H O +OH 0 2.2E14 0.0 16.8
14 H2 +O H +OH 0 1.8E10 -1.0 8.9
15 H2 +OH H2O +H 0 2.2E13 0.0 5.15
16 OH +OH H2O +O 0 6.3E12 0.0 1.09
17 CO +OH CO2 +H 0 1.5E7 -1.3 -0.765
18 N2 +O NO +N 0 7.6E13 0.0 75.5
19 O2 +N  NO +N 0 6.4E9 -1.0 6.25
20 CO +O CO2 0 2.5E6 0.0 3.18
21 CO2 +O CO +O2 0 1.7E13 0.0 52.7
22 F +H2  HF +F 0 9.2E13 0.0 1.08
23 H +F2  HF +F 0 8.8E13 0.0 2.40
24 Cl2 +H HCl +Cl 0 8.6E13 0.0 1.2
25 Cl +H2  HCl +H 0 1.45E13 0.0 4.4
26 HCl +F  HF +Cl 0 1.9E12 -0.68 0.6
27 Cl2 +F  Cl +ClF 0 6.2E12 -0.68 0.5
28 Cl +F2  F +ClF 0 7.6E12 -0.68 0.3
29 ClF +H HF +Cl 0 1.8E12 -0.58 3.2
30 ClF +F  HCl +F 0 5.6E12 -0.68 1.9
31 ClF +H2  HCl +HF 0 1.8e10 -0.5 46.34
32 F2 +HCl HF +ClF 0 1.8E10 -0.5 39.43
33 ClF +HCl HF +Cl2 0 1.8E10 -0.5 46.03
34 F2 +Cl2  ClF +ClF 0 1.8E10 -0.5 26.76
35 OH +F  HF +OH 0 2.9E12 -0.68 0.2
36 H2O +F  HF +OH 0 1.4E10 -0.68 .6
37 Cl +OH HCl +O 0 5.9E10 -0.5 39.43
38 HCL+OH H2O +Cl 0 2.25E12 0.0 1.02
Table D.5: Reaction mechanism for C,Cl,F,H,N,O system combustion[7]
112
Third-Body E eciency of C,Cl,F,H,N,O
Species M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 M12
Ar 1.0 1.0 1.0 1.0 0.8 1.0 1.0 0.0 0.0 0.0 1.0 1.0
CO 1.5 3.0 4.0 4.0 1.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0
CO2 6.4 4.0 5.0 8.0 3.0 2.0 5.0 0.0 0.0 0.0 1.0 1.0
H 25.0 12.5 12.5 12.5 10.0 10.0 1.0 1.0 20.0 0.0 1.0 1.0
H2 4.0 5.0 5.0 5.0 2.0 2.0 1.0 1.0 8.0 0.0 1.0 1.0
H2O 10.0 17.0 5.0 5.0 7.0 3.0 1.0 0.0 0.0 0.0 1.0 1.0
N 1.0 1.0 1.0 10.0 10.0 10.0 1.0 1.0 1.0 0.0 1.0 1.0
NO 1.5 3.0 4.0 4.0 1.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0
N2 1.5 3.0 4.0 4.0 1.0 1.0 2.0 1.0 0.0 0.0 1.0 1.0
O 25.0 12.5 12.5 12.5 10.0 10.0 1.0 0.0 0.0 0.0 1.0 1.0
OH 25.0 12.5 12.5 12.5 10.0 10.0 1.0 0.0 0.0 0.0 1.0 1.0
O2 1.5 6.0 5.0.0 11.0 1.0 1.0 25.0 0.0 0.0 0.0 1.0 1.0
HF 2.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 1.0 1.0
F 4.0 0.0 0.0 0.0 0.0 0.0 0.0 2.4 4.0 0.0 1.0 1.0
HCl 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 1.0 1.0
Cl2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 5.0 1.0 1.0
CL 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0
*Argon only present in combustion with air
Table D.6: Third-body e ciency ratios for C,Cl,F,H,N,O system[7]
113
Appendix E
One-dimensional Finite Rate Chemistry Data
The following plots are mole fractions of the chemical species in one-dimensional  -
nite rate chemistry analysis of the converging sections for LOX-Hydrogen and LOX-RP1
propellant mixtures.
E.1 LOX-Hydrogen
Analysis was made with LOX-Hydrogen propellant. Chamber pressure was set to 27.186
Bar and the mixture ratio was 6.0.
Figure E.1: Mass fraction of LOX Hydrogen system in converging section
114
Figure E.2: Mass fraction of H in converging section
Figure E.3: Mass fraction of H2 in converging section
115
Figure E.4: Mass fraction of H2O in converging section
Figure E.5: Mole fraction of O in converging section
116
Figure E.6: Mole fraction of O2 in converging section
Figure E.7: Mole fraction of OH
117