NASA
Reference
Publication
1311
October 1994

Computer Program for Calculation ofComplex Chemical Equilibrium
Compositions and Applications
I. Analysis

Sanford Gordon and Bonnie J. McBridee.N,.;:':":'";"/
~•.~

"

y

.'~'

National Aeronautics andSpace Administration

NASA
Reference
Publication
1311
1994

Computer Program for Calculation ofComplex Chemical Equilibrium
Compositions and Applications
I. Analysis

Sanford Gordon
Sanford Gordon and Associates
Cleveland, Ohio
Bonnie J. McBride
Lewis Research Center
Cleveland, Ohio

National Aeronautics andSpace Administration
Office of Management
Scienti1ic and Technical
Information ProgramPreface
This report presents the latest in a number of versions of chemical
equilibrium and applications programs developed at the NASA Lewis ResearchCenter over more than 40 years. These programs have changed over the years toinclude additional features and improved calculation techniques and to takeadvantage of constantly improving computer capabilities. The minimization-offree-energy approach to chemical equilibrium calculations has been used in all
versions of the program since 1967.
The two principal purposes of this report are presented in two parts. The
first purpose, which is accomplished here in part I, is to present in detail anumber of topics of general interest in complex equilibrium calculations. These
topics include mathematical analyses and techniques for obtaining chemical
equilibrium; formulas for obtaining thermodynamic and transport mixture
properties and thermodynamic derivatives; criteria for inclusion of condensed
phases; calculations at a triple point; inclusion of ionized species; and various
applications, such as constant-pressure or constant-volume combustion, rocket
performance based on either a finite- or infinite-chamber-area model, shock wave
calculations, and Chapman-Jouguet detonations.
The second purpose of this report, to facilitate the use of the computercode, is accomplished in part II, entitled "Users Manual and Program
Description." Vacious aspects of the computer code are discussed, and a numberof examples are given to illustrate its versatility.

iii

Contents
Page
1

1. Introduction
2. Equations Describing Chemical Equilibrium·
2.1 Units
2.2 Equation of State
2.3 Minimization of Gibbs Energy
2.3.1 Gibbs Iteration Equations
2.3.2 Reduced Gibbs Iteration Equations
2.4 Minimization of Helmholtz Energy
2.4.1 Helmholtz Iteration Equations
2.4.2 Reduced Helmholtz Iteration Equations
2.5 Thermodynamic Derivatives From Matrix Solutions
2.5.1 Derivatives With Respect to Temperature
2.5.2 Derivatives With Respect to Pressure
2.6 Other Thermodynamic Derivatives

3
3
3
4
6
7
8
9
9
10
10
11
12

3. Procedure for Obtaining Equilibrium Compositions
3.1 Initial Estimates
3.2 Magnitude of Species Used During Iteration
3.3 Convergence
3.4 Tests for Condensed Phases
3.5 Phase Transitions and Special Derivatives
3.6 Singularities
3.7 Iteration Procedure and Tests for Ions

13
13
13
14
15
16
17
17

4. Thermodynamic Data
4.1 Assigned Enthalpies
4.2 Least-Squares Coefficients

19
19
19

5. Thermal Transport Property Data
5.1 Data for Individual Species
5.2 Mixture Property Data
5.2.1 Viscosity and Frozen Thermal Conductivity
5.2.2 Reaction Thermal Conductivity
5.2.3 Specific Heat for Gases Only
5.2.4 Prandtl Number

21
21
21
21
22
23
23

6. Theoretical Rocket Performance
6.1 Assumptions

25
25
v

6.2 Parameters26
6.2.1 Conservation Equations
26
6.2.2 Velocity of Flow
26
6.2.3 Force26
6.2.4 Specific Impulse
26
27
6.2.5 Mach Number
6.2.6 Characteristic Velocity
27
6.2.7 Area per Unit Mass Flow Rate
27
27
6.2.8 Coefficient of Thrust
6.2.9 Area Ratio
27
6.3 Procedure for Obtaining Equilibrium Rocket Performance27
for lAC Model
6.3.1 Combustion Conditions28
6.3.2 Exit Conditions28
6.3.3 Throat Conditions28
6.3.4 Discontinuities at Throat
28
6.3.5 Assigned Subsonic or Supersonic Area Ratios
29
6.3.6 Empirical Formulas for Initial Estimates of PinlPe .•••. 29
6.3.7 Analytic Expression for Improved Estimatesof PinlP e • . . . . . . . . . . . . . . • • . . • • . . . . . . • . . . . . • 29
6.4 Procedure for Obtaining Equilibrium Rocket Performancefor FAC Model
30
6.4.1 Initial Estimates for Pinf and AjAt • • • • • . • . . • • • • • . . 30
6.4.2 Improved Estimates for Pinf and AjAt • • • . . . . . . . . . . . 31
6.5 Procedure for Obtaining Frozen Rocket Performance31
6.5.1 Exit Conditions
31
6.5.2 Throat Conditions
32
6.5.3 Thermodynamic Derivatives for Frozen Composition32

7. Incident and Reflected Shocks
7.1 Incident Shocks
7.1.1 Iteration Equations
7.1.2 Corrections and Convergence7.1.3 Initial Estimates of TiT} and PiP}
7.2 Reflected Shocks
7.2.1 Iteration Equations
7.2.2 Corrections and Convergence
7.2.3 Initial Estimates of Ts/Tz and PsIPz

33
33
34
35
35
36
36
37
37

8. Chapman-Jouguet Detonations
8.1 Iteration Equations
8.2 Corrections and Convergence
8.3 Initial Estimates of TzlT} and pzlP}

39
39
40
40

9. Input Calculations

41

Appendix-Symbols

45

References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

vi

Chapter 1

Introduction
Knowing the chemical equilibrium compositions ofa chemical system permits one to calculate theoretical
thermodynamic properties for the system. These properties can be applied to a wide variety of problems inchemistry and chemical engineering. Some applications
are the design and analysis of equipment such ascompressors, turbines, nozzles, engines, shock tubes, heatexchangers, and chemical processing equipment.
For more than 40 years the NASA Lewis ResearchCenter has been involved in developing methods andcomputer programs for calculating complex chemical
equilibrium compositions and thermodynamic propertiesof the equilibrium mixtures and in applying theseproperties to a number of problems (Gordon et aI., 1959,
1962, 1963, 1970, 1976, 1984, 1988; Huff et aI., 1951;
Svehla and McBride, 1973; and Zeleznik and Gordon,
1960, 1962a,b, 1968). Earlier versions of the chemical
equilibrium computer program (Zeleznik and Gordon,
1962a; and Gordon and McBride, 1976) have had wide
acceptance. Since the last publication this program has
been under continuous revision and updating, including
several substantial additions. One addition is an option
for obtaining the transport properties of complex
mixtures (Gordon et al., 1984) by methods simpler thanthose of Svehla and McBride (1973). A second addition
is an option to calculate rocket performance for a rocket
motor with a finite-area combustor (Gordon and
McBride, 1988). The present report documents these andother additions and revisions to the program since 1976.
The revised program is called CEA (Chemical
Equilibrium and Applications).
The program can now do the following kinds ofproblems:(1) Obtaining chemical equilibrium compositions
for assigned thermodynamic states. These states may be
specified by assigning two thermodynamic state functions as follows:(a) Temperature and pressure, tp(b) Enthalpy and pressure, hp(c) Entropy and pressure, sp(d) Temperature and volume (or density), tv(e) Internal energy and volume (or density), uv(f) Entropy and volume (or density), sv(2) Calculating theoretical rocket performance fora finite- or infinite-area combustion chamber(3) Calculating Chapman-Jouguet detonations(4) Calculating shock tube parameters for both incident and reflected shocks
Some problems handled by the program use just one
combination of assigned states-namely, the tp, hp, sp,tv, uv, and sv problems. For example, the tp problem,
which consists of a schedule of one or more assigned
temperatures and pressures, might be used to construct
Mollier diagrams. The hp problem gives constantpressure combustion properties and the uv problem gives
constant-volume combustion properties. Other problems
make use of more than one combination of assigned
thermodynamic states. For example, the shock anddetonation problems use hp and tp; the rocket problem
uses hp or tp and also sp.
This report consists of two parts. Part I, containing
the analysis, includes(1) The equations describing chemical equilibrium
and the applications previously mentioned (i.e., rocket
performance, shocks, and Chapman-Jouguet detonations)(2) The reduction of these equations to forms
suitable for mathematical solution by means of iterative
procedures(3) Equations for obtaining thermodynamic andtransport properties of mixtures
Part II, a program description and users manual,discusses the modular form of the program and briefly
1describes each subroutine. It also discusses the preparation of input, various permitted options, output tables,and error messages. A number of examples are also
given to illustrate the versatility of the program.
In addition to the work in chemical equilibrium
calculations and applications over the past 40 years,
progress in computer programs, data generation, and data
fitting has also been made at NASA Lewis for the thermodynamic and thermal transport properties of individual species required for the equilibrium calculations.
Some examples of this effort are Burcat et al. (1985),
McBride et al. (1961, 1963, 1967, 1992, 1993a,b),

2

Svehla (1962), and Zeleznik and Gordon (1961). In
addition to data calculated by us, other thermodynamic
data included in our files are taken from sources such as
Chase et al. (1985), Garvin et al. (1987), Gurvich et al.
(1989), and Marsh et al. (1988). Files distributed withthe computer program are described in part II.
Various versions of the equilibrium program ormodifications of the program have been incorporated
into a number of other computer codes. An example
is Radhakrishnan and Bittker (1994) for kinetics calculations.

Chapter 2

Equations Describing Chemical Equilibrium
Chemical equilibrium is usually described by eitherof two equivalent formulations-equilibrium constants or
minimization of free energy. Reports by Zeleznik andGordon (1960, 1968) compare the two formulations.
Zeleznik and Gordon (1960) shows that, if a generalized
method of solution is used, the two formulations reduceto the same number of iteration equations. However,with the minimization-of-free-energy method eachspecies can be treated independently without specifying
a set of reactions a priori, as is required with equilibrium
constants. Therefore, the minimization-of-free-energy
formulation is used in the CEA program.
The condition for equilibrium can be stated in termsof any of several thermodynamic functions, such as theminimization of Gibbs or Helmholtz energy or the
maximization of entropy. If one wishes to use temperature and pressure to characterize a thermodynamic state,
Gibbs energy is most easily minimized inasmuch astemperature and pressure are its natural variables.
Similarly, Helmholtz energy is most easily minimized if
the thermodynamic state is characterized by temperatureand volume (or density).
Zeleznik and Gordon (1960) presents equations
based on minimization of Gibbs energy. Some of theseequations are repeated and expanded herein forconvenience. In addition, a set of equations based onminimization of Helmholtz energy is also presented.
However, because only ideal gases and pure condensed
phases are being considered, the general notation of
Zeleznik and Gordon (1960) is not used.

2.1 Units
The International System of Units (SI) used in thisreport is

Physicalquantity

Unit

SymbolLength

meter

m

Mass

kilogramkg

Time

seconds

Temperature

kelvinK

Force

newtonN

Pressure

newton per
square meter

Work, energy

joule

N/m 2

J

The numerical values of a number of fundamental
contants are taken from Cohen and Taylor (1987). For
example, the value of the gas constant R taken from thisreference is 8314.51 J/(kg-mole)(K). In those sections
dealing with the computer program, other units are usedin addition to or instead of SI units.

2.2 Equation of StateIn this report we assume that all gases are ideal andthat interactions among phases may be neglected. The
equation of state for the mixture is
PV:::: nRT(2.1a)

P = nRTp(2.1b)

or3

where P is pressure (in newtons per square meter), V
specific volume (in cubic meters per kilogram), n moles
per unit mass of mixture (in kilograms-mole perkilogram), T temperature (in kelvin), and p density (inkilograms per cubic meter). Symbols used in this report
are defined in the appendix. For a reacting chemicalsystem the number of moles n is generally not constant.
Equation (2.1) is assumed to be correct even whensmall amounts of condensed species (up to several
percent by weight) are present. In this event thecondensed species are assumed to occupy a negligible
volume relative to the gaseous species. An example
given in part II of this report illustrates the validity ofthis assumption. In the variables V, n, and p the volume
and mole number refer to gases only, but the mass is forthe entire mixture including condensed species. The
word "mixture" is used in this report to refer to mixturesof reaction products as distinguished from mixtures ofreactants, which are referred to as "total reactants."
On the basis of this definition, n can be written as(2.2)

where nj is the number of kilogram-moles of species j
per kilogram of mixture and the index NG refers to thenumber of gases in the mixture. The molecular weight ofthe mixture M is defined as
1M=-(2.4a)

MW

Molecular weight is given the symbol MW in equation (2.4a) to differentiate it from M. The two different
definitions of molecular weight, M and MW, give different results only in mixtures of products containing condensed as well as gaseous species. Only M is given in
the output, but MW may be obtained from M by meansof

MW=M(l -1: Xj)(2.4b)

j=NG+lwhere x· is the mole fraction of species j relative to all
species In the mixture. Some additional discussion of thedifferences in these molecular weights is given in part IIof this report.

2.3 Minimization of Gibbs Energy
For a mixture of NS species the Gibbs energy per
kilogram of mixture g is given by(2.3a)(2.5)

nor equivalently as

where the chemical potential per kilogram-mole ofspecies j is defined to be
NS

Enflj

M = ~j=-=-l _ _(2.6)(2.3b)

NG

Enj
j=l

where M· is the molecular weight of species j and theindex N~ refers to the number of species in the mixture.
In the CEA computer program, among the NS species,
gases are indexed from 1 to NO and condensed species
from NG + 1 to NS.
More conventionally, molecular weight is definedas4

The condition for chemical equilibrium is the minimization of free energy. This minimization is usually subject
to certain constraints, such as the following mass-balance
constraints:NS

E aijn

j

j=l

-

b~

= 0

(i

1,... ,Q)

(2.7a)

orb., - b~,=:0(i =: 1,... ,~)(2.7b)

where the stoichiometric coefficients aij are the numberof kilogram-atoms of element i per kilogram-mole ofspecies j, the index ~ is the number of chemical elements (if ions are considered, the number of chemicalelements plus one), bi is the assigned number ofkilogram-atoms of element i per kilogram of total
reactants (see eq. (9.5», andNS(i = 1,...,f)b i = ,Eailnj
j .. l(2.7c)

is the number of kilogram-atoms of element i per kilogram of mixture.
Defining a term G to bee
G = g + ,E Ai(b; - b~)

(2.8)

i=l

where Ai are Lagrangian multipliers, the condition forequilibrium becomes

aG =

,t

~ (fL + APii)an +jjt(hi - hnaAi = 0

state. For a gas the standard state is the hypothetical
ideal gas at the standard-state pressure. For a pure solid
or liquid the standard state is the substance in thecondensed phase at the standard-state pressure.
Historically, the defined standard-state pressure has been
1 atmosphere (101 325 Pa). Most early tabulations ofthermodynamic data were based on this pressure.
However, in 1982 the International Union of Pure andand Applied Chemistry (Cox, 1982) recommended thatthe standard-state pressure should be defined as 1 bar(l05 Pa). Most recent compilations have used 1 bar as
the standard pressure (e.g., Chase et al. (1985), Garvin
et al. (1987), Gurvich et al. (1989), Marsh et al. (1988),
and McBride et al. (l993a,b». The unit of pressure inequation (2.11) should be consistent with the unit ofpressure in the thermodynamic data being used.
The term JlJ and other thermodynamic terms thatappear later in the text, such as C;')' HJ, SJ, C:'j' andare all functions of temperature. However, including
T as part of the symbol notation, such as Hj(]) or HO(]),
is done only when needed for clarity. For example,
sensible heat for species j between a temperature T anda temperature of 298.15 K can be written as HJ(]) HJ(298.15).
Equations (2.7a) and (2.10) permit the determination of equilibrium compositions for thermodynamic
states specified by an assigned temperature To andpressure Po. That is, in addition to equations (2.7a) and(2.10), we have the pair of trivial equations

VJ,

(2.9)
(2. 12a)Treating the variations on) and OAi as independent gives(2. 12b)~~j +L Aflij ~ 0(j ~ 1,...,NS)(2.10)

i=l

and also gives the mass-balance equation (2.7b).
From the assumptions in section 2.2 the chemical
potential can be written asn.

Il·J =r;0lJoj+ RT In-l.. + RT In P

n

U = NG + 1,...,NS)

However, the thermodynamic state can be specified byassigning any two state functions. For example, thethermodynamic state corresponding to a constantpressure combustion is specified, instead of byequations (2.12), by(j = 1,...,NG)h -- h0(2.13a)(2.11)

P = Po(2.13b)where ~l for gases (j = 1 to NG) and for condensed
phases (j> NG) is the chemical potential in the standard

where h is the specific enthalpy of the mixture and ho aconstant equal to the specific enthalpy of the reactants(see eq. (9.7». The expression for h is5NS

h = Enjlj O(2.14)

j=1

H/

is the standard-state molar enthalpy for
where
species j at temperature T.
For assigned entropy and pressure (such as for an
isentropic compression or expansion to a specified
pressure), the thermodynamic state is specified by(2.15a)pansion of the appropriate equations with all terms
truncated that contain derivatives higher than the- first.
The correction variables used are Aln nj (j = 1,...,NG),
Aln n, 1t i = -'AdRT, andI1nj (j = NO + 1,...,NS),

MnT. As Zeleznik and Gordon (1968) points out, it is
no restriction to start each iteration with the estimate forthe Lagrangian multipliers equal to zero inasmuch as
they appear linearly in equation (2.10). After making
dimensionless those equations containing thermodynamic
functions, the Newton-Raphson equations obtained fromequations (2.10), (2.7), (2.13a), and (2.15a) are(2.15b)where S is the specific entropy of the mixture and So theassigned specific entropy, or the specific entropy of thetotal reactant (see eq. (9.22». The expression for s is
NS

S

=

:E nJSj-~(j = 1,... ,NO)

(2.18)(j = NO + 1,...,NS)(2.19)

RT(2.16)

j=l

where

s. =
J

l

S~ - R In nj - R In P
J

SjO

n

(j = NO + 1,... ,NS)(j = 1,...,NG)NG

NS

L a/gn AIn n + L akj An = b; - b(2.17)

j

j

j

j=l

(k = 1,... ,Q)

and S ~ is the standard-state molar entropy for species j.
Equation (2.17) is similar to equation (2.11), and thesame discussion concerning standard-state pressure thatapplied to equation (2.11) also applies to equation (2.17).
NG

2.3.1 Gibbs Iteration Equations

L nj I1ln nj6(2.20)

NG-

n 11In n = n -j=1

The equations required to obtain composition are
not all linear in the composition variables and therefore
an iteration procedure is generally required. In the iteration procedure to be described it will be convenient totreat n as an independent variable. A descent NewtonRaphson method is used to solve for corrections to initial
estimates of compositions nj , Lagrangian multipliers ~,moles of gaseous species n, and (when required)temperature T. This method involves a Taylor series ex-ij=NG+lEn(2.21)jj=l+

F;J)l1ln T= _ho~-_h
(E
li.1
n

R

RT

(2.22)

NS
L (NO
L a.;rl'"~)fJ JAJ 1t + L_J

H? An + ~NG
L n."?) Amn

RT

RT

Q

i=1

So - SRj+lii=NG+l_J_-JJjRT'=1(2.23)+ n

where Cp~j is the standard-state specific heat at constantpressure for species j at temperature T.

(2.27)

2.3.2 Reduced Gibbs Iteration Equations
For problems with assigned thermodynamic states
tp, hp, or sp, various combinations of equations (2.18) to(2.23) could be used to obtain corrections to estimates.
However, for chemical systems containing many species,
it would be necessary to solve a large number ofsimultaneous equations. This large number of equations
can be reduced quite simply to a much smaller number
by algebraic substitution. The expression for Aln nj
obtained from equation (2.18) is substituted into equations (2.20) to (2.23). When equation (2.19) written withsigns reversed is included, the resulting reduced
equations are

t(E
i=1

j=1

~

ai}:tlf j )1t i +R

Sj Anj +j=NG+lR(Ej=1nfj)Aln n
R

NO
NO
S
~
~n .. 1-1.
- LJ n. + L-...LLJ.

j=1

J

j=1

R 2T

(2.28)

Equations (2.24) to (2.28) are given in table 2.1 ina form that permits a direct comparison with other setsof simultaneous equations presented in later sections(tables 2.2 to 2.4). (Note that tables 2.1 to 2.4 appear at
the end of this document.) In a previous report (Gordon
and McBride, 1976), some special script symbols wereused for the sake of compactness in preparing tables 2.1
to 2.4, as for example,

3-C == h(2.29)

RT(2.24)

H~

Q

L
i=1

aij1t

i

+ _J Aln T == ~RT(j == NO + 1,...,NS)

These script symbols are no longer used in this report.
The correction equations required for several typesof constant-pressure problems are summarized asfollows, where i = 1,... ,~ and j = NG + 1,...,NS:RT(2.25)

NO

= n -

Ln

j

j=1

Type ofproblemEquationsrequiredCorrection
variables

Assigned temperatureand pressure, tp(2.24), (2.25),(2.26)

1tj , &!j' ~ln n

Assigned enthalpy
and pressure, hp(2.24), (2.25),(2.26), (2.27)

1tj , ~nj' ~ln n,

Assigned entropy andpressure, sp(2.24), (2.25),(2.26), (2.28)

1tj • ~nj' ~ln n,~In T(2.26)~ln T

7

After obtaining the correction variables shown
above, the corrections for gaseous species Aln nj(j = 1•...,NG) are then obtained from equation (2.18).
Section 3.3 discusses controlling the size of corrections
before they are applied to obtain improved estimates.

2.4 Minimization of Helmholtz Energy
The equations presented in this section are similar
to those in section 2.3. Whatever differences appear are
due to the different forms of the chemical potential
IJ.j (j = 1,...,NG). In section 2.3 pressure was one of theassigned thermodynamic states, and consequently Gibbs
energy was minimized. In this section volume (ordensity) is one of the assigned thermodynamic states,
and consequently Helmholtz energy is minimized.
The two energies (Gibbs and Helmholtz) have thefollowing thermodynamic relationship:

f = g - PV

(2.30)

where f is the Helmholtz energy per kilogram of mixture. After substituting Gibbs energy g as given byequation (2.5), equation (2.30) becomes

Treating fmj and OAi as independent again gives, as insection 2.3, equations (2.7) and (2.10). Now, however,instead of equation (2.11),
n~'Jr

IJ.. =
J(j = 1,...,NO)).L? + RTln-V

J

1
IJ.;

(j = NO + 1,...,NS)

where R' = Rxl0-5•
Equations (2.7) and (2.10), with ~j given byequation (2.35), permit the determination of equilibriumcompositions for thermodynamic states specified by anassigned temperature To and volume Vo; that is, inaddition to equations (2.7) and (2.10), we have the pairof trivial equations(2.36a)T= To(2.36b)NS

f = L IJ.jnj - PV

(2.31)

j.:lThe chemical potential ~j can be expressed as a thermodynamic derivative in several ways (Kirkwood andOppenheim, 1961). One way is given by equation (2.6).
Another expression is(2.32)

Analogous to equation (2.13) for a constantpressure combustion process, we can set down thefollowing conditions for constant-volume combustion:

,

,

= Uo

(2.37a)v = Vo(2.37b)Uwhere u' is the specific internal energy of the mixtureand U a constant equal to the specific internal energy of
the reactants. The expression for u' is

o

If
f

F =f +

L 1; (bi - b:)NS(2.33)

i=1

the condition for equilibrium based on the minimizationof Helmholtz energy subject to mass-balance constraintsis8(2.35)

U'= L njUj

O(2.38)

1'.:1

where Uj O is the standard-state molar internal energy forspecies j.

Analogous to equation (2.15), for assigned entropy
and volume (such as for an isentropic compression or
expansion to a specified volume), the thermodynamic
state is specified byS
L n (s.)-.1 - 1 Alnn. + L -'-An.
R
R
NO

j=1

NS

j

J

j=NG+lJ

NS+(2.39a)n'CO'J Aln T = s - s
L~
R
G R
_0_ _

(2.44)'",1(2.39b)

V=V:o

Iteration equations are derived in section 2.4.1 that
permit solution of the composition variables for constantvolume problems.

where C:'j is the stand~d-state specific heat at constantvolume for species j at temperature T.

2.4.2 Reduced Helmholtz Iteration Equations

2.4.1 Helmholtz Iteration Equations
Correction equations are obtained in a manner
similar to that described in section 2.3 .1. In this case,however, the expression for flj is equation (2.35) rather
than equation (2.11). Because n does not appear
expliCitly as a variable in equation (2.35), Mn n nolonger appears as a correction variable. The NewtonRaphson equations obtained from equations (2.10), (2.7),(2.37a), and (2.39a) areUOj~.
a ..1t.- -AInT= __J
i=1 J} IRTRT

Equations (2.40) to (2.44) may be reduced to amuch smaller set of working correction equations byeliminating Aln nj , obtained from equation (2.40), fromequations (2.42) to (2.44). When equation (2.41) writtenwith the signs reversed is included, the resulting reduced
set of equations is

Q

Alnn.JE(j =' 1,...,NG)(k = 1,... ,~)(2.45)(2.40)
U~Q-Eaij1t; i=1_J

RT

AInT = -~RT(j =' NG + 1,... ,8 )(j = NG+ 1,...,NS)(2.41)
NONS~ akl"j Aln nj +J=1

E akjAnj = b; - b(k = 1,... ,~)kj=NG+l(2.42)
NO

n.U~

NS

E .-L...!- Aln n. + E
j=1

RT

J

j=NO+l

U~
_I

RT

G

'=1

v,j

R

Q

i",1

_I

J

j=l

RT

j",NG+l

RT

J

J.E n CO J Aln T = u~ - ujNSU~ An.
L (NG
L a'Jl.u~J
IrJ
1t; + L

An.

NS

+

(2.46)

RTl(2.43)(2.47)

9

t (E
i",l

j=l

attj Sj - R) Tt i +

R

+

E

Sj b.njRj"'NG+l[E

njC:' j +R

j"'l:::: So - S

R

E
j=l

+

R)] L\In T

njUjO (Sj -

E

2

R T

niSj - R)~j

i"'l

R 2T

(2.48)

Equations (2.45) to (2.48) are given in table 2.2 in aform that permits direct comparison with the iteration
equations in table 2.1 and the derivative equations intables 2.3 and 2.4.
The correction equations required for several typesof constant-volume problems are summarized as follows,
where i = l,... ,~ and j = NG + 1,... ,NS:
Correction
variablesType ofproblemEquationsrequired

Assigned temperatureand volume, tv(2.45), (2.46)'Itj ' Anj

Assigned internalenergy and volume,uv(2.45), (2.46),(2.47)'Itj , An • Aln T

Assigned entropy andvolume, sv(2.45), (2.46),(2.48)'Itj ' Any Ain T

j

After obtaining the correction variables thecorrections for gaseous species L\In nj (j=l,...,NG) arethen obtained from equation (2.40). (See section 3.3 fordiscussion on controlling size of corrections before
applying them to obtain improved estimates.)

2.5 Thermodynamic Derivatives From
Matrix Solutions
All thermodynamic first derivatives can be expressed in terms of any three independent first derivatives. The Bridgman tables, as tabulated, for example, inLewis and Randall (1961), express first derivatives interms of (aVIan 1" (av/ap)T' and (ah/a1)p == cpo We use
the logarithmic form of the volume derivatives because
it gives an indication of the extent of chemical reaction
occurring among the reaction species. These derivatives
may have more than one value depending on what isassumed occurs to composition in a thermodynamic
process from one condition to another. If, for example,
10

composition is assumed to reach its equilibrium value
instantaneously, the derivatives are referred to as"equilibrium" derivatives. If, on the other hand, reaction
times are assumed to be infinitely slow, composition
remains fixed (frozen) and the derivatives are referred toas "frozen." Special subscripts are used to differentiatethese different conditions only for c, thermal
conductivity, and Prandtl number. The equilibrium valueof cp may be expressed as the the sum of a "frozen"
contribution and a "reaction" contribution as follows:(2.49a)whereNS

cp ,! =

L njC~j(2.49b)i=l

andcp~ = ~nJ~(::;L J~f+1 I;"(:~L (2.49c)+

The expressions in equations (2.49b) and (2.49c) wereobtained by differentiating equation (2.14). (Subscripts
different from those used in equation (2.49a) are used intransport property calculations. See section 5.2.3).
From equation (2.1)

V)

aIn( aInT p

1 +

am V)
( aInP T-1 +(alnn)
aInT(2.50)
p(aInn)(2.51)

alnP T

2.5.1 Derivatives With Respect to Temperature
The derivatives of nj and n with respect to temperature are needed to evaluate equations (2.49c) and (2.50).
These may be obtained by differentiating equations(2.10), (2.7), and (2.3), which gives the following:

Ea..(5) _(alnn)aInnj ) _( aIn T p i"'lfJ

aIn T p-t aJ5) ::,aInT RT
Hj O

i"'l

p

aIn T pH~J

RT(j = 1,...,NG)(2.52)(j :::: NG + 1,...,NS)(2.53)EEj ) = 0algnj(_BIn_ni ) +akj(_a_ni"lBIn T P i'''NO+laIn T p(k = 1,...,0)Cp,e =HO (aLe ~NO
L aijnr - (a)~ + L --L----..!!L )

R

i=1

I

(2.54)+'=1E(2.55)As in the case of the iteration correction equations
previously discussed, the derivative equations can be
reduced to a much smaller number of simultaneous
equations by eliminating (BIn nj/aln 1)p, obtained fromequation (2.52), from equations (2.54) and (2.55). Whenequation (2.53) written with the sign reversed isincluded, the resulting reduced number of temperature
derivative equations is

RT

NSaIn T pj=NG+lRT aln T p

NG 11l0)
NS
CO
NO
(H~2
L
3::L (aIn n) + L n + L n~

j'=1 RTnj( aIn nj ) _ n( BIn n) = 0
j=1
aIn T pBIn T pj
_IlOJ

aIn T p

j=1

j

p, j

R

j }

j=1 R 2 T 2(2.59)

In equation (2.59) only the temperature derivatives
obtained directly from solution of equations (2.56) to(2.58) are required. Furthermore, all the coefficients inequation (2.59) are exactly the coefficients appearing in
the reduced-enthalpy equation (2.27). The second-lastterm in equation (2.59) is the frozen contribution to
specific heat (eq. (2.49b»; the remainder of the terms are
the reaction contribution (eq. (2.49c».

2.5.2 Derivatives With Respect to Pressure
The derivative (aIn n/aln P)T can be obtained in amanner similar to that described for obtaining derivatives
with respect to the temperature. Differentiating equations (2.10), (2.7), and (2.3) gives(2.56)(5)Ea ..
i=1 J} BIn T p

H~

__
1

RT

(j = NO + 1,...,NS)

a1t. )'
LLa.n.i=1 j=1
aIn T
D

NO

NO

(

J}

1

j

-1

(j = 1,...,NG)(2.57)e . (an.
-Ea.
-aInP' ) T =0;=1

(j = NG+ 1,...,NS)(2.60)(2.61)J},Ilo-E3::Lp_ Ea .. (~) _(aInn))( alnnaInP T i.. l J} aInP T aInP T(2.58)

j=1 RT

Equations (2.56) to (2.58) are given in table 2.3 ina form that simplifies the comparison of derivative equations with the iteration equations in tables 2.1 and 2.2.
The values of the derivatives obtained from equations (2.56) to (2.58) could be used in equation (2.52) toobtain derivatives for the gaseous species (BIn niBIn 1)p(j = 1,...,NG) and all of the temperature derivatives
could then be used to evaluate cp from equations (2.49).
However, there is an alternatIve and much simpler
procedure for obtaining cpo Substituting BIn nialn 1)pobtained from equation (2.52) into equations (2.49) anddividing by R yield

E
i"l

E

akjnj( aIn nj ) +akj( Bnj ) = 0
Bin P T j=NG+1
aIn P T(k = 1,...,0)(2.62)nj(alnnj ) - n( aInn) = 0
j=1
aInP T
alnP T(2.63)

E

Equations (2.60) to (2.63) can be reduced to asmaller set by eliminating (BIn n/BIn P)T' obtained fromequation (2.60), from equations (2.62) and (2.63). Whenequation (2.61) written with the sign reversed isincluded, the results are

11

L L ak/li/lj _a_NOQ

11:. )(I

aln P T;=1 j=1NO(k = 1,...,0)

=

j --

T

j=1

p~(aln V)2

aIn P T

j=NG+l

L ak.Jn (aInn)
L aiflj
aln P
NO

+

where

(a n.I )+ L alcj _ _NS

TaInT paln V)(2.70)( aInP T(2.64)

j=1

Using the symbols

La..(arc.)__ = 0Q(j = NO + 1,...,NS)

I

;=1

IJ

aInP

(2.65)
Y5 E

T

aInP)
(aln
p 5

(2.71)

anda'it. )
LLa..n. aIn 'P = En.
QNG;=1 j=1

(

I}

NO

J

T

(2.72)(2.66)

j=l}equation (2.69) may be written as
Equations (2.64) to (2.66) are given in table 2.4 forcomparison with tables 2.1, 2.2, and 2.3. The only
derivative obtained from solution of equations (2.64)
to (2.66) that is used is (aInn/aInP)T (see eq. (2.51).

2.6 Other Thermodynamic Derivatives
As stated previously, all thermodynamic first
derivatives can be expressed in terms of the threethermodynamic first derivatives discussed in the previous
sections-namely, cp' aIn waIn 1)p, and(aln V/alnP)T(see Bridgman tables in Lewis and Randall (1961».
Velocity of sound a, a frequently used parameter, isdefined by
2

a

(ap)

p(aInP)p(aInP)= ap S = P aln p s-P aIn V s(2.67)

From Bridgman tables

aInP)
(aln
V s

(2.68)V) + pv(aIn V)2
aInP T T aIn T p

c (aInP

aIn V)C ( aInP T
V

12(2.69)(2.73)(aIn V)
alnp T

Using the equation of state given in equation (2.1),
we obtain from equation (2.67) the familiar expression
for velocity of sounda = JnRTy s(2.74)

Note that the y s defined by equation (2.71) is requiredin equation (2.74) and not the specific heat ratio 'Y
defined in equation (2.72).
In section 3.5 an alternative expression is derived
for 'Ys for the special situation of triple phases, where theexpressions in equations (2.68), (2.69), and (2.73) are nolonger valid.
Gordon and Zeleznik (1962) gives numerous first
derivative relations that are of interest in rocket
perfonnance calculations. One of these derivatives,
which is used in section 6.3.4, is(alnP)aInT

This may be written as

alnP)
(aIn
V 5

Y

Y8 = -(2.75)5nR(aln V)

alnT

p

A numerical example involving this derivative is given
in Zeleznik and Gordon (1968).

Chapter 3

Procedure for Obtaining Equilibrium
Compositions
In principle, obtaining equilibrium compositions bymeans of the Newton-Raphson iteration procedure discussed in sections 2.3.1, 2.3.2, 2.4.1, and 2.4.2 should
offer no difficulties. However, a number of practicalitems require detailed attention in order to avoid numerical difficulties: initial estimates, tests for condensed
phases, phase transitions and triple points, convergence,
accidental singularities, special handling of ions, andconsideratio.n of trace species.

3.1 Initial Estimates
An extremely simple procedure is used in thisreport to assign estimates for composition. For the first
iteration of the first point in a schedule of points, weassign n = 0.1, which is equivalent to an estimate of 10
for molecular weight. Then the number of kilogrammoles of each gaseous species per kilogram of mixture
is set equal to O.lING, where NG is the number of gaseous species being considered. The number of moles ofeach condensed species is set equal to zero. For hp, sp,uv, and sv problems an arbitrary initial estimate of
T = 3800 K is used by the program unless a different
estimate is included in the input.
Admittedly, this simple procedure will often givepoor initial estimates. However, for a general chemistry
program, we find this technique preferable to the alternative of devising numerous special routines for obtaining
good estimates for numerous possible chemical systems.
Furthermore, the estimating technique is used only forthe first point in any schedule of points. For all points
after the first the results of a preceding point serve as
initial estimates.
Because no attempt is made to obtain good initial
estimates, the question arises whether convergence canbe "guaranteed."section 3.3.

Thisquestion is

discussedm

3.2 Magnitude of Species Used During
Iteration
Both the linear and logarithmic composition variables are used for gaseous species during the composition
iteration process. Only the linear variable is used forcondensed species. Corrections to compositions forgases are in the form of logarithmic variables .Aln nj' andtherefore the logarithmic values of gaseous compositions
In nj are continuously updated from iteration to iteration.
The linear values of the compositions nj are obtained bytaking the antilogarithm of In nj' However, to save computer time during iteration, nj are calculated only for
those species whose mole fractions are greater than acertain specified size. In a previous version of the
program (Gordon and McBride, 1976) this specified size
had only one value, namely njln = 10-8 (or In njln
= -18.420681). A program variable SIZE was defined as
SIZE = 18.420681. Thus, antilogarithms of In nj wereobtained only for gases meeting the following condition:
In nj/n ~ -SIZE (njln > 10-8 .) For gaseous species not
passing this test nj were set equal to zero. In addition,
the maximum number of iterations permitted was 35.
In the present CEA program several additionalvariables relating to SIZE have been added to handle
more demanding situations. Two of the present variables
relating to the mole fraction size for which antilogarithms are obtained are TSIZE for non-ionized species and
ESIZE for ionized species. TSIZE may be modified forany of the following reasons: inclusion of species in thecalculations with mole fractions smaller than 10-8 (bymeans of an input parameter TRACE); a change of13

components; a singular matrix; or the chemical system
under consideration containing a chemical element thatdiffers in magnitude from the largest of the otherelements by more than 105. The purpose of changing
TSIZE in the last case is to ensure that not all species
containing the trace element will be eliminated during
iteration. To aid in testing for trace elements, a parameter BRATIO is defined to be the ratio of the elementswith the lowest to highest kilogram-atoms per kilogramof mixture. The following, then, are the conditions underwhich several parameters relating to species size are set
to various values. These conditions are based on extensive practical experience with obtaining solutions formany types of problems:(1) TSIZE = SIZE until convergence or a singularmatrix occurs.
(2) TSIZE = XSIZE if TRACE ':f:. 0 after first convergence, or if a singular matrix or new components
occur.
(3) ITN == maximum number of iterations.

can cause large corrections. The first situation occurs in
the early stages of the calculation and is due to poorestimates. The second may occur at later stages of the calculation when the iteration process sometimes attempts tomake extremely large increases in moles of species thatare present in small amounts. An example of this second
situation is given in Zeleznik and Gordon (1962a). In
both of these cases a control factor 'A is used to restrict
the size of the corrections to In nj (j = 1,...,NG) and nj(j = NG + 1,...,NS) as well as to In n and In T obtainedby solving the equations in tables 2.1 and 2.2.
The numerical value of 'A is determined by empirical rules that experience has shown to be satisfactory.
For T and n, corrections are limited to a factor of eO A
= 1.4918. For gas-phase species two different correction
controls are calculated that depend on the magnitude of
the mole fractions. The logarithm of each mole fraction
is compared with a parameter called SIZE whose default
value is -In 10-8 = 18.420681. If In(nj/n) > -SIZE,
corrections to nj are limited to a factor of e2 = 7.3891.
For these limitations on corrections to T, n, and nj/n thevalue of a control factor A.l may be calculated as

Default:(1) SIZE = -In 10-8 = 18.420681(2) XSIZE = -In 10-11 = 25.328436(3) ESIZE = -In 10-14 = 32.236191(4) ITN =50(5) TRACE = 0

2

(3.1)

For those gaseous species for which In(n/n)~ -SIZE and .Aln nj ~ 0, a control factor 'A2 is defined as

Nondefault:(1) If TRACE t:- 0, ITN = 50 + NS/2
(2) If TRACE < 10-8 , XSIZE = -In TRACE and
ESIZE = -In (TRACExl0-3) =XSIZE + 6.9077553(3) If BRATIO < 10-5 , SIZE = In l000IBRATIO
and XSIZE = SIZE + In 1000 (6.9077553)(4) If singular matrix, XSIZE = TSIZE = 80.
The use of ESIZE to control the size of ionizedspecies permitted to be present during iteration is
discussed in the section 3.7.

3.3 Convergence
The problem of convergence is discussed in
Zeleznik and Gordon (1962a) and Gordon and McBride(1976). Zeleznik and Gordon (1962a) points out that theiteration equations sometimes give large corrections that,
if used directly, could lead to divergence. Two situations
14n.
-1n-.1 - 9.2103404

A = minl--n-1
2
aln n - aln n
j

(3.2)

This prevents a gaseous species with a small mole
fraction from increasing to a mole fraction greater than
10-4. The control factor 'A to be used in equations (3.4)
is defined in terms of Al and ~ as
(3.3)
A value for A is determined for each iteration.
Whenever current estimates of composition and/or temperature are far from their equilibrium values, A will beless than 1. Whenever they are close to their equilibrium
values, 'A will equal 1. New estimates for composition

and temperature are then obtained from the correction
equations

When TRACE(k)In+ A(i)(~ln n.)W
nj(i+1) = In n~i)
J
J

1t i

(j = 1,...,NG)*' 0, an additional test is used:(1+1)- 1t i(1+1)(i = 1,... ,~)~ 0.001(3.6d)

1tin)i+l)= n),J + A(i)(~nj)(I)(j = NG + 1,...,NS)

In n (i+l) = In n (i) + AW(AIn n)(i)where the superscript i represents the ith estimate.
The iteration procedure is continued until corrections to composition satisfy the following criteria:(j1,... ,NG)

(j = NO + 1,...,NS)!!J~ln nj ~ O.5xlO- 5
NS

En
j=1

(3.5)

j

For those chemical elements for which b~ >1.0xl0-6 , the convergence test for mass balance is(i = 1,... ,~)where the subscript "max" refers to the chemical element
i with the largest value of b~. When temperature is avariable, the convergence test for temperature is(3.6b)For a constant-entropy problem (sp, sv, or rocket), thefollowing convergence test on entropy is also required:s - slT
I~ O.5xlO- 4(3.6c)

where the superscript refers to the kth iteration. The
convergence tests in equations (3.5) and (3.6) ensure
accuracy to five places in composition when expressed
as mole fractions.
As pointed out in section 3.2 the maximum numberof iterations permitted by the CEA program is
50 + NS/2. For most of the hundreds of different kindsof problems that have been solved by the program,convergence has been obtained in fewer than thisnumber. For the first point, which starts with arbitrary
initial estimates, a typical number of iterations is 8 to 20.
For succeeding points, which use compositions of apreviously calculated point for initial estimates, a typicalnumber of iterations is 3 to 10. In the sample problems
given in part II, for example, the number of iterations forthe first point was 9 to 15. In some special cases a
problem may be singular, and solutions to the correction
matrix are then unobtainable. Techniques for handling
this situation are discussed in section 3.6.

3.4 Tests for Condensed Phases
For the first point in a schedule of points, unless
INSERT records are used (see part II of this report), theprogram considers only gaseous species during the iteration to convergence. For each point after the first, theprogram uses the results of a previous point for its initial
estimate. After every convergence the program automatically checks for the inclusion or elimination of condensed species.
The test is based on the minimization of Gibbs
energy. At equilibrium, equation (2.9) is satisfied (i.e.,= 0). The requirement for a condensed species j,which was not previously included as a possible species,
to now be included is that its inclusion will decrease
Gibbs energy; that is, from equation (2.9)

oG

J.1~)
-BG = ( _J
anj

RT c

-

E
Q

1tflij <

0

(3.7)

i=1

where the subscript c refers to a condensed species.
Equation (3.7) is identical to the vapor pressure test used
in Zeleznik and Gordon (l962a) when data for the gas
15phase as well for the condensed phase are available. The
vapor pressure test is

n,.p)
0) (0
(RT c RT + I nrn- < 0
~.

_J

_

~.

_J

(3.8)

g

where the subscript g refers to the gas phase of the same
species as the solid phase referred to by the subscript c.
Equations (3.7) and (3.8) are equivalent. The last term inequation (3.7) is identical for the gas and condensed
phases of the same species. From equation (2.10) thisterm for the gas phase equals ('ri'!RT)g' which, using thedefinition of 'rij in equation (i.l1), leads directly toequation (3.8). However, the advantage of equation (3.7)
over (3.8) is that equation (3.8) can be used only when
the gas phase of the species corresponding to the condensed phase to be tested is present, whereas equation(3.7) can always be used. Using equation (3.7) eliminates
the need for the extensive programming required in
Zeleznik and Gordon (l962a) to accommodate
AI 20 3(s,I), for which gas-phase data are not available.
At most, only one new condensed species is included after each convergence. In the event that several
condensed species pass the test required by equation (3.7), only that species giving the largest negativechange to Gibbs energy is included as a possible species
and convergence to a new equilibrium composition is
obtained. This process is repeated until all condensed
species required by equation (3.7) are included.
If, after convergence, the concentration of acondensed species is negative, the species is removed
from the list of currently considered species, and convergence to a new equilibrium composition is obtained.

3.5 Phase Transitions and Special
Derivatives
The calculation method is based on the assumptionthat condensed phases are pure. Therefore, the possibility
exists of encountering phase transition between solid andliquid (melting points) or between two stable solid
phases. Such transitions constitute triple points because
three phases of the same species coexist, one gaseousand two condensed. Such triple points are characterized
by a definite vapor pressure and temperature, independent of the relative proportions of each phase. This
characterization is shown by the fact that the iteration
16

equations of table 2.1 become singular for an assigned
temperature and pressure and the inclusion of two
condensed phases of the same species. At a triple point,
for a specified system pressure, the relative amounts ofthe phases can be determined only if either the enthalpy
or the entropy is assigned.
The program can obtain equilibrium compositions
containing either one or two condensed phases of a
species with or without the corresponding gas phase.
When temperature is assigned, no problems arise as towhich one of the two or more condensed phases of a
species is to be considered. However, when temperature
is a variable, such as in combustion or rocket problems,
several possibilities need to be considered. For example,
if a liquid phase is being considered by the program andthe temperature at convergence is below the melting
point, two possibilities exist. First, the solid phase might
be substituted for the liquid phase and a new convergence obtained, or second, both liquid and solid might
be considered simultaneously and a new convergenceobtained. Similar possibilites exist when convergence isobtained with a solid phase above the melting point.
These possibilities and methods of treating them are
discussed in detail in Zeleznik and Gordon (1962a); thediscussion is not repeated herein. In brief, the following
criteria are used by the program to determine whether toswitch one condensed phase of a species to another orwhether to consider both simultaneously:(1) Liquid present at temperature T < Tm:(a) If Tm - T > 50 K, switch solid for liquid.
(b) If Tm - T ~ 50 K, include solid and liquid.
(2) Solid present at temperature T> Tm:(a) If T - Tm > 50 K, switch liquid for solid.
(b) If T - Tm ~ 50 K, include solid and liquid.
Similar tests apply for two solid phases of a species.
The unusual situation of constant temperatureduring an isentropic compression or expansion process
can occur when two phases of the same species coexist.
The transition temperature remains constant while one
phase is being converted to the other. Under thiscircumstance derivatives with respect to temperature are
not defined. Thus, several of the derivatives previously
discussed, cp' c v' and (dIn Vlaln T)p' cannot be obtainedin this case. As a consequence equation (2.73) cannot be
used to obtain "'(s'In Gordon (1970) the following expression wasderived for the situation of constant temperature andentropy:

YS,T ==P)aln( aln p

1
S,T

am V)
(aln P T(3.9)

The velocity of sound for this situation is no longer
given by equation (2.74) but rather byas,T= JnRTY S.T(3.10)

These derivatives will be used in connection withdiscontinuities in a rocket throat (see section 6.3.4).

3.6 Singularities
The iteration method used in this report has
successfully handled numerous chemical systems undera wide variety of thermodynamic conditions. Nevertheless, special procedures are required to take care ofseveral situations involving singularities that would
otherwise cause the iteration method to fail. One suchsituation is a singularity in the coefficient matrix that iscaused when all species that contain one of the elements
are temporarily eliminated during the iteration process.
Another situation resulting in a singularity occurs when
two rows of the coefficient matrix are identical. This
happens when the ratio of the assigned elements in these
two rows is equal to the ratio of the stoichiometriccoefficients of these two elements in every gaseous
species being considered during the current iteration that
contains both elements. That is,(j = 1,...,NG)(3.11)One example for which equation (3.11) applies is
stoichiometric hydrogen and oxygen at low temperaturesand pressures, where the only species with nj > 10-8 is
H20(g). Another example is stoichiometric lithium andfluorine at low temperatures, where the only species are
LiF, Li 2F 2, and Li 3F3. When either of these twosituations occurs, the program will automatically reset all
species currently set to zero to nj = 10-6. This resetfeature can be tried twic y and will usually take care ofthese causes of singularity.
If the restart procedure just described is not successful and the situation represented by equation (3.11)
is the cause of the singularity, the program will attemptan additional procedure: (1) selecting one of the larger
species containing both elements to be one of the newcomponents; (2) reducing the number of components andthe size of the coefficient matrix by 1; and (3) switchingcomponents for the rest of the chemical system.
Another singularity is caused when several
condensed species are currently being considered andone of these species can be formed by some linear
combination of the others. An example of this is achemical system containing iron and oxygen where thecurrent iteration is simultaneously considering FeO(s),Fe203(s), and Fe30 4 (s). When this occurs, the programautomatically removes one of the species involved withthe singularity, prints an error message with thisinformation, and restarts.

3.7 Iteration Procedure and Tests
for Ions
The program is capable of calculating equilibrium
properties of plasmas (mixtures containing ionizedspecies) if the plasma is considered ideal. Ideal is here
meant to imply that no coulombic interactions are
considered. Plasma textbooks, such as Oriem (1964),
point out that effects of coulombic interactions do need
to be considered in plasmas (by means of the DebyeHuckel approximation, e.g.). However, special programming is needed to take care of these effects.
Therefore, because the program does not consider theseeffects, the results of calculations when ions are
considered will be valid only for those conditions where
the ionic density is so small that the coulombic effects
are unimportant.
For ions to be considered, the charge-balance
equationNG'Laejnj ==j=l0(3.12)

is required, where aej indicates the excess or deficiencyof electrons in the ion relative to the neutral species. For
example, in a mole of an ionized species, aej = -3 for
Ar+++ and +1 for 02' To prevent difficulties in matrix

17

solutions, the program automatically removes the chargebalance equation when, for each ionized species being
considered, In nj/n < -ESIZE (defined in section 3.2).
There are situations when all the previous
convergence tests have been passed but the ion balanceis still incorrect. A special iteration procedure wasdeveloped to obtain the correct ion balance for thesesituations. It consists of obtaining a value of the
Lagrangian multiplier for ions divided by RT, based on
the assumption that the magnitude of the ionized species
is small relative to the un-ionized species. The initial
estimate for 1t e is taken to be the value in storage for thecurrent point or from a previous point. The iteration
procedure consists of the following steps:(1) Corrections to 1te are obtained fromNG

- Lael"jj=lNG

L (aej?nj
j=l18(3.13)(2) The test for convergence is(3.14)(3) If this convergence test is not met, newestimates for the composition of ionized species are
obtained from
= (In n·Y + a . A 1t
(In n.)i+l
J
J
eJ
e

(3.15)

where the superscript i refers to the ith iteration. The
previous sequence of steps is repeated until equation (3.14) is satisfied. The CEA program allows for a
maximum of 80 ion-balance iterations, but generally
convergence is reached in about 2 to 10 iterations.

Chapter 4

Thermodynamic Data
Thermodynamic data are included with the current
program for reaction products and reactants. The numberof these products and reactants changes over time as newones are added to the data base. The data are selectedfrom a number of sources, but the principal current
sources are Chase et aL (1985), Cox et aL (1989),
Gurvich et al. (1989), and Marsh et al. (1988). McBrideet al. (1993a) documents the sources and the data for 50
reference elements plus electron gas and deuterium,
which are presently included in the thermodynamic data
set. These elements are discussed in section 4.1. The data
for the atomic gases as well as for a number of diatomic
and polyatomic gases were calculated at NASA Lewis byusing the PAC91 computer program described in
McBride and Gordon (1992).
The thermodynamic data provided with the CEA
program are in the form of least-squares coefficients (tobe described). These data, in formatted form, areprocessed by subroutine UTHERM and stored for further
use in unformatted form. Subroutine UTHERM and theformat for the coefficient data are both described inpart II of this report.

4.1 Assigned Enthalpies
For each species heats of formation (and, whenapplicable, heats of transition) were combined withsensible heats to give assigned enthalpies BO(T). By
definition
W(1)== W(298.15)+ [(W(1) - Jr(298.15)](4.1)

In general, HO(T) *" LifHO(T) for T *" 298.15 K.
A set of reference elements must be specified inorder for heats of formation to be unambiguously related
to specific reactions. Included among the species forwhich thermodynamic data are on a data file are thefollowing 50 reference elements plus deuterium andelectron gas: Ag, AI, Ar, B, Ba, Be, Bi, Br2 , C, Ca, Cd,
C1 2, Co, Cr, Cu, 02' e-, F 2, Fe, Ge, H2, He, Hg, 12, K,
Kr, Li, Mg, Mn, Mo, N2 , Na, Nb, Ne, Ni, 02' P, Pb, Rb,
S, Si, Sn, Sr, Ta, Th, Ti, U, V, W, Xe, Zn, and Zr. The
thermodynamic data for these elements are documented
in McBride et al. (l993a). For all reference elements
Li HO(298.15) = HO(298.15) = O.
f
Assigned enthalpies for a number of reactants arein a reactant data file. For noncryogenic reactants
assigned enthalpies (heats of formation) are given at298.15 K. For cryogenic liquids assigned enthalpies are
given at their boiling points. These enthalpies are usuallyobtained by subtracting the following quantities from theheat of formation of the gas phase at 298.15 K: thesensible heat between 298.15 K and the boiling point,
the difference in enthalpy between the ideal gas and thereal gas at the boiling point, and the heat of vaporization
at the boiling point.

4.2 Least-Squares Coefficients
For each reaction species the thermodynamic
functions specific heat, enthalpy, and entropy asfunctions of temperature are given in the form of leastsquares coefficients. The general form of these equations
is as follows:

We have arbitrarily assumed HO(298.15) =~'fHO(298.15).
Equation (4.1) then becomesW(1) = Af H'{298.15) + [W(1) - H'(298.15)]

(4.2)

(4.3)

19

Jr

fC;dT

RT

RT

So

Co

-

(4.4)

Jr

-

RT

2

T
T
-a1T- 2 + a2 T- 1 In T + a3 +a+a42

T3

Ii = f R~dT

(4.5)

Two sets of least-squares coefficients have beenused in the chemical equilibrium program to represent
thermodynamic data. The "old" form has been used forthe last 25 years (Gordon and McBride, 1976; andMcBride et aI., 1993b). It consists of a fourth-orderpolynomial for the C'j/R function plus integration
constants for enthalpy and entropy as follows:

Sr + aJr =a + a2T + a3T2 + a4T3 + a
RT
1
2
3
4
5

6

(4.7)

T

The "new" form for these functions will be used torepresent the thermodynamic data with the program
described in part II. This form consists of seven termsfor C'j/R and corresponding terms for enthalpy andentropy as well as the integration constants a8 and a9 as
follows:+ a T2 + a T 3 + a T4
5

20

6

7

(4.9)

T- + a8+ a -+ a 6 4
7 5

So

-

R

T- 2

= -a 1

2

-a T
2

-1

T

5 3

(4.10)+ a3 lnT+a4 T(4.11)For gases the temperature intervals for both the oldand new functional forms are fixed. These intervals are
200 to 1000 K. and 1000 to 6000 K for the old form
(i.e., the fourth-order polynomial form for C;). The newform uses these same intervals plus an additional intervalfrom 6000 to 20 000 K for some gases. For the condensed species each phase has its own set of coefficients.
If possible, the old form uses the same two temperatureintervals for condensed species as for the gases, but theintervals are usually limited by transition points. Furthermore, there are two intervals only if the 1000 K common
point is within the species temperature range. Otherwise,
there is just one. By contrast, the new functional form
has a flexible number of intervals in order to fit the
selected data more accurately.
Generally, the three functions C;fR, HOIRT, andSOIR are fit simultaneously. The fit is constrained tomatch the functions exactly at T = 298.15 K. Thus, theleast-squares coefficients reproduce heats of formation at
T = 298.15 K exactly. For temperature intervals that do
not contain T = 298.15 K, the fit is constrained to givethe same functional values at the common temperaturepoint for any contiguous intervals. When phase transitions occur, the fit is constrained so that the differencein Gibbs energy is zero between the phases at thetransition point.

Chapter 5

Thermal Transport Property Data
The CEA program provides an option for calculating mixture viscosities and thennal conductivities. The
fonnulas for these mixture properties are given in section 5.2. Viscosity and thennal conductivity data for
individual gaseous species, which are required for the
mixture calculations, are included with the current program. Thermal transport properties for condensed species
are not included inasmuch as there is no feasible method
for calculating thennal transport properties for a multiphase mixture. The thennal transport data were takenfrom Svehla (1995).

5.2 Mixture Property Data

5.1 Data for Individual Species

where Aeq' Afr , and Are are the equilibrium, frozen, andreaction thermal conductivities of the mixture,
respectively. The mixture viscosity, on the other hand,
has only one term.

The thennal transport property data provided with
the CEA program are in the fonn of least-squares coefficients. The data for each species were fitted to thefollowing fonn, which is also used in Gordon et al.
(1984):

'll}

In
B
C
=AlnT+-+-+D
In A
T
T2

(5.1)

Thermal conductivity and specific heat of a mixture
each consists of two parts, the so-called "frozen" and"reaction" contributions. This was discussed for specificheat in section 2.5. Equation (2.49a) shows specific heatas the sum of these two parts. Analogously, thermal
conductivity can be written as(5.2)

5.2.1 Viscosity and Frozen Thermal Conductivity
As pointed out in Gordon et al. (1984), mostapproximate mixture methods have the following fonn
for mixture viscosities and frozen thermal conductivities(also the form used in the CEA program):NM

11 mix =Xi'lli

L

NM!

i=1

where fl is viscosity and A is thennal conductivity. A
binary interaction parameter flij is also included for somepairs of species in the same fonn as equation (5.1).
The coefficients in equation (5.1) were generated togive viscosity in units of micropoise (flP) and thermal
conductivity in units of microwatts per centimeter-kelvin(flW/cm-K). The order and format of the transport data
coefficients are given in part II of this report.

Xi

(5.3)

+ I,XJ-<Pij
j=l

rt:i

and
NM

Aft =

x;A t

L

NM

i=1

Xi +

I ;Xj.ij

(5.4)

j=l
fi;i

21

where NM is the number of gaseous species for thermaltransport property calculations (NM ~ 50), x; is the mole
fraction of species i relative to the NM gaseous species
used for thermal transport property calculations, 11; is theviscosity of species i, Tlmix is the viscosity of the mixture, Ai is the thermal conductivity of species i, 4>ij is theviscosity interaction coefficient between species i and jin equation (5.3) (<pij -:;:. 4>ji), and 'IIij is the interaction
coefficient between species i and j in equation (5.4) ('I'ij-:;:. 'IIj).
Several forms for the interaction coefficients 4>ij and'IIij available in the literature are compared in Gordon
et ale (1984), and the following were used:/)112<1>.. -_ -1[ 1 + (1'1
- i )II2(M- i )1 4]2( 2M
'J
4
1'1
Mi
Mi + M
j

i

(5.5)

i

(5.6)

For some pairs of species an interaction parameter 11ij isavailable from the literature. As discussed in section 5.1,
these values have been least squared and are included in
the transport property data file. When this parameter isavailable, 4>ij is obtained from the following equation
rather than trom equation (5.5):

where NR is the total number of chemical reactions ofgaseous species and li.,Hi is the heat of reaction forreaction i at temperature T expressed as
NM!J.,H: (1) = L ex ufI;(il,2,...,NR) (5.9)

k;1where (J,ik are the stoichiometric coefficients written forthe chemical reactions involving species A k as follows:
l,2,...,NR)(i(5.10)

Equations (5.8) to (5.10) are also applicable to ionizing
gases. In this case the heat of reaction is replaced by theionization potential, and ions and electrons are species.
Note that the stoichiometric coefficients (J,ik discussed here are defined differently from the stoichiometric coefficients aij discussed in chapter 2. In thecoefficients (J,ik the subscripts refer to species k inreaction i, and the coefficient may be positive or
negative. In the coefficients aij the subscripts refer tochemical element i in species j, and the coefficient isalways positive.
The Ar,i required in equation (5.8) are found bysolving the following set of linear equations:
NR

1'1.

«1> .. = - '
IJ
1'1 M

9

2M.
J

i

(5.7)

LgijA rj(i

1,2,...,NR)(5.11)j;l+ Mj

where the gij are given by
The interaction parameter llij also appears in the expression for the reaction contribution to thermal conductivity.
The same values of 4>ij that are obtained from equations (5.5) or (5.7) are also used in equation (5.6).

5.2.2 Reaction Thermal Conductivity

and where

The reaction contribution to thermal conductivity is
obtained in the same manner as discussed by Butler andBrokaw (1957) and Brokaw (1960) and as used by
Svehla and McBride (1973). The following equation isused, which is applicable when local equilibrium exists
in a mixture of reacting or ionizing gases:
NR

6.,HioA =R"-A.
re~
RT T,',;1

22(5.12)(5.8)(5.13)The A:~ factor in the previous equation is a collisioncross-section ratio that is used in Svehla and McBride(1973) but is not used in this report. In its place we havesubstituted the value of 1.1, which, as may be seen in
Hirschfelder et al. (1954) (table I-N for the LennardJones potential and table VII-E for the modified

Buckingham 6-exp potential), is a typical value (to two
figures) for this parameter. The parameter l1kQ inequation (5.13) was briefly discussed in section 5.2.1 andis available for some pairs of species with the transport
properties of individual species. For most pairs ofspecies, however, it is calculated from equation (5.7)
where <PkQ is defined, with appropriate subscripts, by
equation (5.5). The thermal conductivity of the gas
mixture as given by equation (5.2) is then obtained byadding the results of equations (5.4) and (5.8).

An alternative method of calculating the reaction
contribution to specific heat is given in Svehla andBrokaw (1966) that differs from but is equivalent toequation (2.49c) or the reaction part of equation (2.59).
It is analogous to the method for obtaining the reaction
contribution to thermal conductivity (eq. (5.8», namely(5.16)

5.2.3 Specific Heat for Gases Only
Some functions, such as Prandtl number, involve a
number of other parameters (e.g., specific heat, viscosity,and thermal conductivity). In this event, for consistency,
the parameters involved should be based on similarassumptions. Inasmuch as viscosity and thermal conductivity are calculated for a maximum of 50 gases (with nocondensed species permitted), specific heat, when usedwith these properties, should also be based on the samegases. The following equation, which is similar toequation (2.49a) but with different subscripts, is used toexpress specific heat for transport property calculations:(5.14)
When no condensed species are present and thesame gaseous species are included in the transport
property calculations as in the thermodynamic property
calculations, equations (2.49a) and (5.14) will produce
the same numerical results. Otherwise, they may yield
different results.
The equation for cp,fr is similar to equation (2.49b)
except for being restricted to gases and having a
different limit for the number of gases involved; that is,

Note that the upper case Xi in equation (5.16) are not thesame variables as the lower case Xi' The Xi are found bysolving the following set of linear equations:
NR

A,Er/

j=1

RT

Ld~i==-i == 1,2,...,NR(5.17)

where the dij are given by(5.18)
Note that the terms in equation (5.18) are the same as
those in equation (5.12) except that the RTIPDkQ groupis not present.

5.2.4 Prandtl Number
Prandtl numbers have application in heat transfer
calculations. Frozen and equilibrium Prandtl numbers
may be calculated from previously calculated properties
by means of(5.19a)NM~x.Co.

L...J

I

p,l

;=1

Cp,fr

==

NM

LX!!;

(5.15)

and(5. 19b)i=1

23

Chapter 6

Theoretical Rocket Performance
Before the publication of Gordon and McBride(1988), the Chemical Equilibrium Calculations (CEC)
computer program described in Gordon and McBride(1976) could calculate theoretical rocket performance
only for an infinite-area combustion chamber (lAC).
Calculation of rocket performance for a finite-area combustor (FAC), presented in Gordon and McBride (1988),
was added as an option to the Chemical Equilibriumwith Transport Properties (CET) program in 1988.
Figure 6.1 presents schematic cross sections of FAC
and lAC rocket engines. Various points at which calculaInfinite-area chamber, inf
I
I
I

I
I
\

Exit, e

\
\
\

I
I
I

\
\

I

\

\ Combustor end, c
\

I

" \

Propellants

~

I
\, I

Exhaustproducts- - - - - 1I " " - - - - - - - - - - I
I
I

I
I
I

Injector, inj

Throat, t(a)Infinite-area chamber, inf-------------------------~(b)

Figure 6.1.-Schematics of rocket combustion chamber crosssections. (a) Finite-area combustion chamber. (b) Infinite-area
combustion chamber.

tions are made in the CEA program to obtain rocket
performance are indicated in these figures. Combustion
and throat parameters are always calculated first automatically. For the lAC model, only one combustion
point is calculated, namely, at infinite area (int). However, for the FAC model, two combustion points arecalculated, namely, at the combustion chamber inlet (orequivalently at the injector face, inj) and at the combustor end, c. In addition to these two combustion points for
the FAC, a combustion calculation for an infinite-area
combustor, indicated in figure 6.1(a) by the dashed line,
is also made. The results at this fictitious point are used
as an aid in an iteration procedure to obtain combustor
end conditions, as discussed in section 6.4. In addition,
the pressure at this point is used in calculating c* (seesection 6.2.6). Throat conditions are indicated by thesubscript t; other nozzle exit conditions, either subsonicor supersonic, are indicated by the subscript e. Nozzle
conditions are assigned as an option and may be in theform of assigned area ratios, pressure ratios, or both.

6.1 Assumptions
The calculation of theoretical rocket performance
involves a number of assumptions. For the same propellant and operating conditions theoretical performance canvary depending on which assumptions are used. For thisreport most of the assumptions are the same for both thelAC and FAC models. These assumptions are onedimensional form of the continuity, energy, andmomentum equations; zero velocity at the combustion
chamber inlet; complete combustion; adiabatic combustion; isentropic expansion in the nozzle; homogeneous
mixing; ideal-gas law; and zero temperature lags andzero velocity lags between condensed and gaseous
species. The chamber in the FAC model is assumed tohave a constant cross-sectional area. In this chamber

25

combustion is a nonisentropic, irreversible process.
During the burning process part of the energy released isused to raise the entropy, and the pressure drops.
Expansion in the nozzle is assumed to be isentropic.
Combustion conditions are obtained with theassumption of chemical equilibrium of the combustion
products. For the lAC model the CEA program provides
the option of calculating either equilibrium or frozen
theoretical rocket performance. Equilibrium performanceis based on the assumption of instantaneous chemical
equilibrium during expansion in the nozzle. Frozen performance is based on the assumption that composition
remains frozen at the combustion composition during
expansion. For the FAC model only equilibrium performance is permitted.
Assuming the same velocity (either zero or otherwise) at the combustion chamber inlet, identical thermodynamic results are obtained for the combustion inlet
condition for both the lAC and FAC models.

6.2 Parameters

Equation (6.2) applies only for constant-area, onedimensional flow.

6.2.2 Velocity of Flow
The combustion chamber inlet is indicated by thesubscripts inf for the lAC model and inj for the FAC
model. Then using these subscripts instead of 1 andusing e instead of 2 in equation (6.3) and assuming thevelocity at the combustion chamber inlet to be negligible
relative to the exit velocity result in equation (6.3)
becoming
for IAC model(6.5)

for PAC model

where h is in units of joules per kilogram and u is inunits of meters per second.

6.2.3 Force

6.2.1 Conservation Equations
Rocket performance, as well as other fluid dynamic
problems in the program, is based on the following
conservation equations, which are consistent with theassumptions previously discussed:

From the momentum principle of fluid mechanics
the external force on a body in a steadily flowing fluid
is due to the change of momentum of the fluid and to
the increase in pressure forces acting on the body. For
rocket applications this is expressed asmU

e
F=-+(P-P)~ea e(1) Continuity:

(6.6)

gc

(6.1)

(2) Momentum:(6.2)

The conversion factor gc has been introduced to allow
for various units. For some commonly used systems ofunits, such as the cgs system or the International System(Goldman and Bell, 1986), gc :;:; 1. However, in the
English Technical System, commonly used by engineers,

gc :;:; 32.1740 (lbmllbf)(ftls 2)(3) Energy:2~h2 + -

2

= hI

(6.3)

Equation (6.1) describes the condition of constant mass
flow rate, which will be given the symbol m; that is,

m = pAu

26(6.4)

6.2.4 Specific Impulse
Specific impulse is defined as force per unit mass
flow rate. From equation (6.6)
F
1=-

m

ue
gc

(Pe - Pa)A

+ - - - - -e

m

(6.7)

In rocket literature specific impulse is often expressed in
English Technical System units of pounds force perpound mass per second. However, for those systems ofunits previously mentioned for which gc = 1, I is bothdimensionally and numerically equal to velocity.
In this report when the exit pressure is equal to theambient pressure, specific impulse will be given thesymbol I sp . From equation (6.7)(6.8)
When the ambient pressure is assumed to be zero(vacuum conditions), specific impulse will be given thesymbol I vac . From equations (6.7) and (6.8)(6.9)

6.2.5 Mach Number
Mach number is defined as the ratio of velocity offlow to velocity of sound:vII = !!:.

(6.10)a

Velocity of flow is given by equation (6.5). Velocity ofsound is given by equation (2.74) (or eq. (3.10».
6.2.6 Characteristic Velocity
Characteristic velocity is given the symbol c* andis defined as

PmAtgcc*m(6.11)

6.2.7 Area per Unit Mass Flow Rate
From equation (6.4)

A

m

1

(6.12)The terms p and u e are obtained from equations (2.1b)and (6.5), respectively.

6.2.8 Coefficient of Thrust
The coefficient of thrust IS defined in terms ofpreviously defined parameterscF = ~
c*(6.13)

6.2.9 Area Ratio
The ratio of area at any exit nozzle station to areaat the throat is obtained from the values of equation (6.12) at these two points
Ae(A/m)eAt(A/m)t(6.14)

6.3 Procedure for Obtaining
Equilibrium Rocket Performancefor lAC Model
The procedure for obtaining equilibrium performance differs somewhat for the lAC and FAC models.
The principal difference is due to only one combustion
point being required for the lAC model (point inf in
fig. 6. 1(b» but two combustion points being required for
the FAC model, namely, at the inlet and exit of the finite
chamber (points inj and c in fig. 6.1(a». The procedure
for the lAC model is discussed first, inasmuch as it issomewhat simpler.
For the lAC model the procedure consists of first
determining combustion properties and then determining
exhaust properties at the throat and at other assigned
stations, if any, in the nozzle exit. Combustion and throat
conditions are always obtained first automatically by the
CEA program. To obtain other nozzle conditions (either
subsonic or supersonic), the desired points must be
specified as part of the input in the form of assigned area
ratios or pressure ratios.
For the FAC model the procedure involves first
determining combustion properties at the combustor
inlet. This station is also referred to as the "injector (inj)station." Conditions at the end of the combustor c and at
the throat are then both determined by means of an
iteration loop that also includes the fictitious point inf

27

(fig. (6.1 (a», as discussed in section 6.4. The two
combustion stations inj and c and the throat are always
calculated first automatically by the CEA program. As isthe case for the lAC model other nozzle exit points, if
desired, must be specified.

is equal to the velocity of sound. The second procedure
is used in this report. Throat pressure is determined by
iteration.
The initial estimate for the pressure ratio at thethroat is obtained from the approximate formula6.3.1 Combustion Conditions(6.15)For both the FAC model at station inj and the lAC
model at station inf, the same combustion temperatureand equilibrium compositions are obtained by the program for an assigned chamber pressure and the reactant
enthalpy (in the program HP = .TRUE.). From the
combustion compositions and temperature the combustion entropy and other combustion properties are determined. For the lAC model the combustion entropy sinf
is assumed to be constant during isentropic expansion in
the nozzle. However, combustion in a finite-area rocket
chamber (the FAC model) is a nonisentropic process.
The entropy increases during the combustion processfrom si!1i to Sc while pressure drops from P inj to Pc' For
the FAC model the entropy at the end of the combustion
chamber s c is held constant during isentropic expansion
in the nozzle.

6.3.2 Exit Conditions
Exit conditions include the throat conditions andassigned area ratios A/At or pressure ratios PinclPe or
Pin/Pe' Throat conditions are always detennined
automatically by the program. Other exit conditions, on
the other hand, are optional and, if included, will becalculated after throat calculations are completed.
For an assigned pressure ratio equilibrium compositions and exit temperature are detennined for thepressure P corresponding to the assigned pressure ratio
and for the combustion entropy (sinf for lAC and Sinj for
FAC). For throat and assigned area ratios iteration
procedures are used to detennine the correct pressure
ratios. These procedures are described in sections 6.3.3
to 6.3.7.
After equilibrium compositions and temperature are
obtained for an assigned pressure ratio or area ratio, all
the rocket parameters for that point can be determined.

6.3.3 Throat Conditions
Throat conditions can be detennined by locating thepressure or pressure ratio for which the area ratio is aminimum or, equivalently, for which the velocity of flow

28

Equation (6.15) is found in many references on rocket
propulsion, such as Sutton and Ross (1976), but is exact
only when 'Ys is constant from combustion point tothroat. Because the value of 'Ys is not yet known at thethroat, the value of 'Ys from the combustion point is usedby the CEA program in equation (6.15). It generally
gives an excellent initial estimate.
Equilibrium properties for sinf and for the value ofP t calculated from equation (6.15) are obtained as for(usingany exit point. From these properties(using eq. (2.74» are calculated andeq. (6.5» andthe following test for convergence is made:u;a;~ O.4xlO- 4(6.16)

This criterion is equivalent to ensuring that at the throat
the Mach number is within 1±O.2x10-=-4.
If the convergence test is not met, an improved
estimate of the throat pressure ratio is obtained from theiteration formula2

1 + Y s ul( )Pt,k+l == ( P 1 +Yst,k(6.17)

where the subscript k indicates the kth iteration. A
maximum of four iterations on throat pressure ispermitted. Usually, no more than two are required.

6.3.4 Discontinuities at Throat
Gordon (1970) gives a special procedure forobtaining throat conditions when the velocity of soundis discontinuous at the throat. This type of discontinuity
may occur when a transition point, such as a melting
point, is being calculated at the throat. The solution ofthis problem requires the following equation, whichpermits estimating the throat pressure at the melting
point, where the solid phase just begins to appear:

lnP = InP + (aIn P) (In T - In 1)
t
aIn T
m

(6.18)

S

where the derivative is given by equation (2.75).

and(~: ~ 2) (6.22)

6.3.5 Assigned Subsonic or Supersonic Area Ratios
An iteration process is used in the CEA program toobtain the pressure ratios corresponding to the assigned
area ratios. In general, initial estimates for pressure ratios
in the iteration scheme are obtained from empiricalequations derived from plots of previously calculated
data. However, if the current as well as the previous
assigned supersonic area ratios are greater than 2, ananalytic expression is used to obtain the initial estimate.
The same analytic expression is also used in all cases toobtain improved estimates for the pressure ratios.

6.3.6 Empirical Formulas for Initial Estimates of
Pin/PeInitial estimates of pressure ratios corresponding to
subsonic area ratios are obtained from the following
empirical formulas:

In equation (6.22) the value of Ys is that determined forthroat conditions.

6.3.7 Analytic Expression for Improved Estimates of
Pin/Pe
The equilibrium properties obtained for the initialand each subsequently improved estimate of PinlPe areused in equations (6.23) and (6.24) to obtain the next
improved estimates. From table I of Gordon and
Zeleznik (1962) the following derivative can be obtained:

Pinf
aln-

Pe

aIn~

At

s

This derivative is used in the following correction
formula to obtain an improved estimate for PinlPe:

A e + lO.587(ln A e )3 + 9.454 In A e
At
At
At

(~: ~ 1.09)(6.23)

A(6.19)Pinf) -_(I nPint)
- +
(InPe
Pe
k+l

k

and

Pinf
0.9 InPt

Pinf(6.24)

In-=----------AeeP
A e + lO.587(ln A e )3 +9.454ln
AtAt(1.0001 < ~: < 1.09)

At(6.20)

When an assigned supersonic area ratio requires aninitial estimate of pressure ratio to be obtained from anempirical formula, the following formulas are used:

where the subscript k refers to the kth estimate andwhere the area ratio with no iteration subscript is theassigned value. The iteration procedure is continued until(I nPint)
Pe(Pm)- Ink+lPe k:5;0.00004(6.25)

29

with a maximum of 10 iterations permitted. Generally,
convergence is reached within two to four iterations.

6.4 Procedure for Obtaining
Equilibrium Rocket Performancefor FAC ModelIn this report the finite-area combustor is assumed
to have a constant cross section. For this constant-area
combustion chamber the momentum equation (6.2) canbe written as(6.26)
Substituting rizlA = pu (from the continuity eq. (6.4» into
equation (6.26) gives= (p mu)(p mu)..
A
A+(6.27)+mJc

IPinj - (P + pu )cl ~ 2xlO-5
2(6.29)Pinj
If equation (6.29) is satisfied, the calculations forthe finite-area combustor are complete for points c and t.
Calculations are then continued if other values ofpressure ratio and/or area ratio have been specified in the
input. If equation (6.29) is not satisfied, an improved
estimate for Pinf is obtained as described in section 6.4.2,
and the procedure for the points inf, t, and c is repeateduntil equation (6.29) is satisfied.
An iteration procedure similar to that described for
option 1 is also used for option 2. However, thecontraction ratio is not known for option 2. Therefore,
the iteration procedure involves starting with anestimated value for A/At as well as for P inj and thenobtaining improved estimates for both P inj and A/At.
Not surprisingly, more iterations are required for
option 2 than for option 1, which requires improved
estimates for P inj only. As for option 1 iteration iscomplete when equation (6.29) is satisfied.

6.4.1 Initial Estimates for Pinf and A/At

When velocity at the injector face is negligible, equations (6.26) and (6.27) reduce to

A curve given in figure 3-14 of Sutton and Ross
(1976) relates PinlPinj with A/At for an assumed valueof y = 1.2. The following empirical equation was derived
by fitting three selected points read from the curve:(6.28)

A

1.0257 - 1.2318~

An iteration procedure is required to satisfy
equation (6.28). This procedure is somewhat different foreach of the two possible input options available for FAC.
In option 1 the contraction ratio A/At is assigned. In
option 2 the mass flow rate per unit combustor area mlA c
is assigned. The iteration procedure for option I issimpler and therefore will be described first. Four points
shown in figure 6.1 (inj, inf, c, and t) are involved in theiteration procedure. Thermodynamic parameters at the
injector face are obtained by a combustion calculation.
This point needs to be calculated only once. Startingwith an estimated value for Pinf' calculations are thenmade for points inf, t, and c (the assigned contraction
ratio) exactly as is done for the lAC model (for theinfinite-area combustion point, the throat, and an
assigned subsonic area ratio). A check is then made tosee if equation (6.28) is satisfied to within the tolerance30

At(6.30)

A

1 - 1.26505~

At

Equation (6.30) is used only to obtain an initial estimate
for Pinf'
For option 1 the assigned value of the contraction
ratio A/At is used in equation (6.30). For option 2 an
initial estimate of A/At is required. This initial estimate
is obtained from

(6.31)

Equation (6.31) was derived by starting with therelationship for characteristic velocity c* = PinrAtlrit(eq. (6.11) with gc = 1), multiplying both sides by A c'and using a representative but approximate value of
c* = 2350 mls. When equation (6.31) is substituted into
equation (6.30), a quadratic in Pinf is obtained. The
solution of this quadratic gives the initial estimate for
Pinf' This estimate for Pinf' when used in equation (6.31),
gives the initial estimate for A/At. If the input value ofrit/A c is so large that equation (6.31) calculates a valueless than 1, the program will stop the calculations andprint out the error message INPUT VALUE OF
MDOT/A = (value) IS TOO LARGE. GIVES CR
ESTIMATE LESS THAN 1.

tion rates) during expansion is simpler than thatassuming equilibrium composition. The reason is that .
equilibrium compositions need be determined only for
combustion conditions. After obtaining combustion conditions in the identical way described for equilibrium
rocket performance, the remainder of the procedure is asfollows.

6.5.1 Exit Conditions
Improved estimates of the exit temperaturecorresponding to some assigned Pin/Pe are obtained bymeans of the following iteration formulas:(6.34)

6.4.2 Improved Estimates for Pinf and A/At
where
For option 1 an improved estimate for Pinf
(Pinfnew) is obtained by assuming that the ratio of thisdesired Pinf,new value to the current value of Pinf is equal
to the ratio of the assigned value of P inj (Pinj,a) to thecurrent value of P inj (obtained by means of eq. (6.28)).
This assumption leads to(6.32)( dIn Te) k =

s

- s
inC

e,k

(6.35)

Cp,e.)cand where k refers to the kth estimate. The initial estimate of an exit temperature is the value of temperature
for the preceding point. The iteration procedure iscontinued until(6.36)

Equation (6.32) often gives such an excellent improved
estimate for Pinf that it need be used only once to obtain
convergence (eq. (6.29». For option 2 an improved
estimate for A/At is required in addition to the one for
Pinf and is obtained from

m

Ac
At

-

At

m
-

(6.33)

Ac

For some test cases involving hydrogen and oxygen as
propellant, approximately four iterations involving equation (6.33) were required for convergence (eq. (6.29».

6.5 Procedure for Obtaining Frozen
Rocket Performance
The procedure for obtaining rocket performance
assuming that composition is frozen (infinitely slow reac-

The maximum number of iterations permitted by the program is eight, although the convergence criterion ofequation (6.36) is generally reached in two to four
iterations.
Phases are also considered to be frozen. Therefore,
the program will calculate frozen rocket performance forassigned schedules of pressure andlor area ratios only
until an exit temperature is reached that is 50 K below
the transition temperature of any condensed species
present at the combustion point. If a calculated exit temperature is more than 50 K below the transition temperature, this point and all subsequent points in the schedule
are ignored by the program and do not appear in the
output. In addition, the following message is printed:
CALCULATIONS WERE STOPPED BECAUSE NEXT
POINT IS MORE THAN 50 DEG BELOW TEMPRANGE OF A CONDENSED SPECIES.
After an exit temperature has been determined, all
the rocket performance parameters (section 6.2) can be
determined.

31

6.5.2 Throat Conditions
Calculations for frozen throat conditions are similar
to those for equilibrium conditions. That is, equation (6.15) is used to get initial estimates for PinlPt ,equation (6.17) is used to get improved estimates for
PinlPt , and equation (6.16) is used as the convergencecriterion. With composition and phases frozen, there isno possibility of discontinuities at the throat, in contrast
to equilibrium compositions.

sition, has already been given in equation (2.49b). Someother derivatives based on frozen composition are as
follows: From equations (2.50) and (2.51), respectively.

(6.37)(6.38)

6.5.3 Thermodynamic Derivatives for Frozen
Composition
From equation (2.70)
The thermodynamic derivatives discussed in previous sections were based on the assumption that in any
thermodynamic process, from one condition to another,composition reaches its equilibrium values instantaneously. If, on the other hand, reaction times are assumed
to be infinitely slow, composition remains fixed (frozen).
In this event expressions for derivatives become simpler(see Zeleznik and Gordon (1968) for further discussion).
An expression for specific heat, based on frozen compo-

32(6.39)

From equations (2.73) and (6.38), it is clear that forfrozen composition
Y8 = Y(6.40)

Chapter 7

Incident and Reflected Shocks
The solution to the conservation equations that
describe conditions for incident and reflected shocks ismost conveniently obtained in terms of assigned
temperature before and after shock. Therefore, theoreticalvalues for shock parameters, such as velocities and gas
compositions, are generally found tabulated as functionsof assigned temperatures (or temperature ratios). However, in shock tube experiments shock velocities are
generally known, rather than temperature ratios, thus
requiring interpolation for purposes of comparison. It istherefore useful to have a calculating scheme that calculates shock properties in terms of assigned velocities, andthat is the scheme used in this report.
No consistent nomenclature has been found in theliterature for shock parameters. For example, in Gaydon
and Hude (1963) u is the velocity of the gases relative
to the velocity of the shock front w, and v is the actual
velocity of the gases in the shock tube. In Glass and Hall
(1959), however, the converse is given. We will use essentially the nomenclature of Gaydon and Hude (1963).
Figure 11.6 in Gaydon and Hude (1963) shows theusual incident and reflected conditions in bothlaboratory-fixed and shock-fixed coordinates. Velocities
can be expressed as
or

it =

w-v

(7.4)
and
(7.5)
In equation (7.5) it is assumed that the shocked gases atthe end of the shock tube are brought to rest (""v5 = 0).
The asterisk is used within equation (7.4) to differentiate it from 112 in equation (7.3). The gas properties
are the same for condition 2, but the relative velocities
are different. All quantities that appear in equations (7.2)
to (7.5) are positive, and henceforth the arrows will be
dropped.

It;

7.1 Incident Shocks
For a constant-area shock the conservation equations describing incident conditions are those given byequations (6.1) to (6.3) with Al = A 2. For iterative
purposes it is convenient to reduce these equations to aform similar to that given in Zeleznik and Gordon(1962b) as follows:

(7.1)
(7.6)

For incident shock waves equation (7.1) gives
(7.2)

(7.7)

and
(7.3)

In equation (7.2) it is assumed that the test gas is at rest(v 1 = 0). For reflected waves equation (7.1) gives

It is also convenient to use the symbols p* and h* forthe right sides of equations (7.6) and (7.7), respectively.
These equations then become33

(7.8)

P*(7.13)(7.9)

h* - h2 = 0

In equation (7.9), h2 isequation (2.14).

defined by means of(7.14)

7.1.1 Iteration Equations
Applying the Newton-Raphson method to equations (7.8) and (7.9) divided by R and using for the
independent variables the logarithm of temperature ratio
and pressure ratio across the shock give

(7.15)(7.10)

T[(BlnaIn V)T2+ Ma(h*-h2 )pa(h*-h2 )R~In~ +

R

P

P

Bln~

2

P;l -

1]
(7.16)

1',
aln~

I

T1

PI

h2 - h*(7.17)(7.11)

R

where

( Aln

~2) = In(~2)
Ik.

-In(i)(7. 12a)I k+1 1 kh2 - h*(7.18)R( Alni) In( ~2) -In( ~2)
=

1 k

I 1+1

(7.12b)1 k.

and the subscript k stands for the kth iteration.
The partial derivatives and right sides in equations (7.10) and (7.11) are given as follows:

34

where from equations (2.1 b) and (2.3a)

PI = PI M 1 T2
P2
P2 M 2 T1(7.19)

7.1.2 Corrections and Convergence
An arbitrary control factor A is applied to thecorrections obtained from solution of equations (7.10)
and (7.11) before using the corrections in equations (7.12) to obtain improved estimates. The purposeof the control factor is to limit the size of the corrections
in order to prevent divergence. Different maximum
corrections are permitted during the course of theiteration process. Depending on the iteration number thiscontrol factor permits a maximum correction of 1.5 to
0.05 times the previous estimates of PiPt and TiTt.
This is the same as permitting a maximum absolutecorrection of 0.4054651 1 to 0.04879016 on In PiP! andIn TiTt. As an illustration, the control factor A,
corresponding to the maximum correction of 0.40546511
is obtained by means of the following equations:A = 0.40546511

lL\ht ::1 < 0.00005(7.24a)T
Aln ~ < 0.00005
T1(7.24b)7.1.3 Initial Estimates ofTiTl and P21P1
Formulas for temperature ratio and pressure ratio
across the incident shock, assuming constant 1, can befound in texts such as Gaydon and Hude (1963). These
formulas, slightly rearranged, are
P2

2Yl~At~-Yl +1

PI

Yl + 1(7.25)(7.20)

3(7.26)

0.40546511(7.21)

lL\ht~1(7.22)A schedule of numerical values used in equations (7.20)
and (7.21) (in addition to the value of 0.40546511) andthe criteria for when they are used are discussed inpart II of this report.
Improved estimates are obtained by using equation (7.22) with equation (7.12) to giveIn(~2)
1 k+l= In(~2)
I k+

A(L\ht ~2)(7.23a)1 kht(~L = ht(~:L + A(L\ht ~J. (7.23b)

If'Y is constant over the temperature range Tt to T2 ,equations (7.25) and (7.26) give exact theoretical values
for PiPt and TiTt. When composition is assumed to befrozen across the shock, 12 does not vary greatly from 11'
and equations (7.25) and (7.26) generally give excellent
first estimates. Improved estimates are then obtained as
described in sections 7.1.1 and 7.1.2. However, if
equilibrium composition is assumed across the shock,equation (7.26) generally gives estimates of temperature
ratio that are too high. In Gordon and McBride (1976)
an empirical equation is used to lower the T2 estimatesobtained from equation (7.26). This procedure has now
been replaced by a new estimating scheme that gives
considerably better estimates, especially at high Mach
numbers. The initial estimate for P 2 is still obtained fromequation (7.25). However, the initial estimate for T2 isobtained by carrying out a combustion calculation using
the estimate for P2 for the assigned pressure and thevalue obtained from the following equation for theassigned enthalpy:2

U1

ho=h2 =h 1 +2The iteration process is continued until corrections
obtained from equations (7.10) and (7.11) meet thefollowing criteria:(7.27)The properties for condition 1 are for the unshocked
gases.

35

7.2 Reflected ShocksP= ~ _P'

The reflected shock conservation equations could be
written like equations (6.1) to (6.3), but it is more convenient to use the relations in equations (7.4) and (7.5),
which giveP2

(7.35)(7.28)

h5 -h'(7.29)

I 2 h
1)2
h5 + -W
2 R = 2 + -(v
2 2 + WrR

(7.30)For iterative purposes it is convenient to reduce equations (7.28) to (7.30) to a form similar to equations (7.6)
and (7.7) as follows:(7.36)
where( ~ln ~5l : : In(~5) - In(~5)(7.37)( ~ln i) : In( i) -In( i)(7.38)

2 k(7.31)

2 k+l

2 k

Vi[1 _p 2]

h5 = h2 + 2+a(p, - Ps )2vi
P_2~ _ M

For convenience, the right sides of equations (7.31)
and (7.32) will be given the symbols P' and h'. The
reflected shock parameters may now be solved by thesimultaneous solutions of(7.33)

Paln-sa(p, - ~:)
T

aln 2

Ps

~

(aln V) _ Ps (7.39)(PS)2 alnP
--1
P;

T,S

P2

Ps

2
M2v2

---

RTz

-

P2

(alnV)(:: -lr alnT

(7.40)
P5(7.34)

7.2.1 Iteration Equations
Applying the Newton-Raphson method to equations (7.33) and (7.34) divided by R and using forindependent variables the logarithm of temperature ratio
and pressure ratio across the shock give

36

RT2

P2

T2

h' - h s = 0

2 k(7.32)---->-_ _

P

2 k+l

The partial derivatives and right sides for equations (7.35) and (7.36) are as follows:--.2 - 1
P2P'-~ ::: 0
P2

2 k

Ps
2
M 2V 2

P2
1-----RT.2

Ps

P2

- 1

(7.41)

7.2.2 Corrections and Convergence{h' ~ h
5

Ps

)

-

2
V2

P

aln~

R ( :: _P2

r

1

alnP ~5T5 [( Bin V)+ M

a( hI - h5 )2-..:l..._R----"-- _V2salnT P,s -T2

1](7.42)

P2

7.2.3 Initial Estimates of T s/T2 and P S/P2
A value of TsiTz = 2 is generally a satisfactory
initial estimate. An estimate for PsIPz in terms of Ts/Tz
can be obtained by inverting equation (2.36) in Gaydon
and HurIe (1963). For Ts/Tz = 2 this gives

Ps( aln V)R ( -Ps- 1 )2 alnT P,5
P2

Ts
aln-

The same control factors, correction equations, andtests for convergence discussed in the previous section
for In PiPl and In TiTl (eqs. (7.20) to (7.24)) areapplicable for In PsIPz and In Ts/Tz.

(aln)

P2

TS cp ,5R

(7.43)_PS)2 _ (Ys+l)P- - -s _ 2-0
(-P
Y 1 P
2

hs - hi(7.44)

s -(7.45)

2

Only one solution of equation (7.45) is positive.

R

37

Chapter 8

Chapman-Jouguet Detonations
The method used for obtaining Chapman-Jouguet
detonation parameters is described in Zeleznik andGordon (l962b). There are three steps in the calculation
procedure. First, an initial estimate of the detonationpressure and temperature is made. Second, the estimateof these parameters is improved by using a recursion
formula. Third, the correct values are obtained by meansof a Newton-Raphson iteration procedure. The required
equations are derived in Zeleznik and Gordon (l962b)
and are summarized herein for convenience in slightly
modified form.
The same conservation equations (6.1) to (6.3) for
continuity, momentum, and energy that apply for shock
also apply here with the additional constraint thatu2 = a2· For iterative purposes, as was true for the shock
equations, it is convenient to reduce the three conservation equations to two:
- y
s;J.

(pz _ 1)PIh tl - hz -- 0

(8.4)

8.1 Iteration Equations
Applying the Newton-Raphson method to equations (8.3) and (8.4), dividing by R, and using forindependent variables the logarithm of temperature ratio
and pressure ratio across the detonation give

(8.5)

(8.1)

a(h hz)ll(8.2)

a(h" -hTz)p-----''---_R_....t.... AIn ~ + ----''---_R_....t.... AIn ~
P
PI
T
T1
dln-z
aIn-z
-

PI

hz - h t '
R

TI

(8.6)
The pressure ratio in equation (8.1) is the reciprocalof the pressure ratio in equation (7.6) because, as pointed
out in Zeleznik and Gordon (1962b), the NewtonRaphson iteration encounters fewer problems in the formof equation (8.1).
For convenience in writing the iteration equations
the symbols P' and h" are used to represent the rightsides of equations (8.1) and (8.2), respectively. Theseequations then become
P " _ PI

Pz

= 0

(8.3)

where AIn PiPl and AIn TiTl are defined by equation (7.12).
The partial derivatives appearing in equations (8.5)
and (8.6) can be evaluated if Ys is taken to be independent of temperature and pressure. This assumption isreasonable for moderate ranges of temperatures andpressures. To within the accuracy of this assumption thepartial derivatives and right sides of equations (8.5)
and (8.6) are39

8.3 Initial Estimates of T2/T1 andP2/P l(8.7)

Pz= YS,2~(aln V)aln T P,2

(8.8)

A good initial estimate for pressure ratio is not asimportant as a good estimate for temperature ratio. For
a number of chemical systems that were investigated, an
initial estimate of (PiPl)O :::: 15 has been found to besatisfactory. An initial estimate for temperature ratio isfound by calculating the flame temperature T 2
corresponding to the enthalpy

h

2

(8.9)

= h

+
1

i RT1(PPZ)

4MI(8.13)

1 0

Improved initial estimates for the assumed value of(PiPl)O and the estimated value of (TiTl)O corresponding to h2 in equation (8.13) can be obtained by using therecursion formulas(8.10)

(T
T
Z

)

1 hI

(T2)

3

R (P

Z

)

= T 1 0 -"4 M 1cp ;1. PI 0
2

+

(P

Ry S,2 'hI - 1 Z )2M1 cp2'k+l
PI k+1,(8.15)(8.11)
whereh2 - h"R_ YS,2T2 ~(P2)Z _1]
2M2

l PI

(8.12)

8.2 Corrections and Convergence
The discussion of convergence controls and corrections for shock calculations applies equally well for
Chapman-Jouguet detonations. Equations (7.20) to (7.22)
give the control factor A; equation (7.23) gives the new
estimates for P 21P 1 and TiTl; and equation (7.24) is thetest for convergence.

40(8.16)(8.17)The quantities M2 , 'Ys ,2' and cp ,2 in equations (8.14)
to (8.16) are the equilibrium values for (PiP 1)0 and(TiTl)O' Repeating the use of equations (8.14) to (8.17)
three times in the program generally provides excellent
initial estimates for the Newton-Raphson iteration.

Chapter 9

Input Calculations
A number of options are provided in the programfor specifying input, and these options are discussed indetail in part II. Part of the input concerns reactants. The
total reactant may be composed of a number of reactants,
each of which may be specified as an oxidant or a fuel.
If the total reactant contains more than one oxidant,
these oxidants may be combined into a total oxidant byspecifying the relative proportion of each oxidant.
Similarly, if the total reactant contains more than one
fuel, the fuels may be combined into a total fuel by
specifying the relative proportion of each fuel. The
overall properties of the total reactant (such as elemental
compositions hi, enthalpy ho' internal energy
molecular weight M o' density Po' positive and negativeoxidation states yt- and V-, specific heat co' and entropyso) can then be calculated by specifying the relative
amounts of total oxidant and total fuel. This method isparticularly convenient if calculations are to be obtained
for a number of oxidant-fuel ratios.
To obtain assigned properties for the total reactant,each reactant j may be specified as either a fuel or anoxidant even though a reactant may ordinarily be thoughtof as inert (e.g., N2). If reactants are not divided into
fuels and oxidants, they are treated like fuels with
olf = O. Letting the superscript k equal I for oxidant(provided that there is an oxidant) and 2 for fuel, thekilogram-atoms of the ith element per kilogram of total
oxidant or total fuel isw5k)n(k) jj- --NRE--AC--

M~k)
J

'(""

LJ

(j = 1,...,NREAC; k = 1,2)

W~k)
J

j=l

(9.2)

Mt)

wjt)

where
is the weight of the jth reactant andisthe molecular weight of the jth reactant. If, however,amounts of oxidants and fuels are specified in terms ofkilogram-moles Nt),un,(j = 1,...,NREAC; k = 1,2)

(9.3)

Ml)

can be calculated from the atomic weights ofThe
the chemical elements M i as
(j = 1,...,NREAC; k = 1,2)

The b~k) can be combined by means of thefollowing e'quation to give the kilogram-atoms per
kilogram of total reactant:b~2) + (!!-)b~l)~
f l(i = 1,... ,~; k = 1,2)

(9.1)

(9.4)

b~ = --~"----l0(i1,... ,~)

(9.5)

1+-

f

where NREAC is the total number of reactants andn}k)is the number of kilogram-moles of the jth reacJant per
kilogram of total oxidant (k = 1) or total fuel (k = 2). If
the amounts of oxidants and fuels are specified in termsof weights,
is obtained bynt)

where olf is the oxidant-to-fuel weight ratio.
Formulas analogous to equations (9.1) and (9.5) canbe used to obtain other properties of the total oxidant,
total fuel, and total reactant (such as enthalpies, internal
41

energies, molecular weights, and densities). The energies
given in the next four equations are in the fonn in which
they appear in the matrix arrays (tables 2.1 to 2.3).

(6) Molecular weight of total reactant, kg/kg-mole:(l) Specific enthalpy of total oxidant (k = 1) or
total fuel (k = 2) h(k) divided by RT, (kg-molelkg) (k):(9.11)NREAC (H~)(k)LJ+ 1] M (l)M (2)

(7r~)

+

M(I)

~k)

n}

(k = 1,2)

(9.6)

RT

j=l

7

(If MCl) = 0, M o = MC2); and if Af~) = 0, M o = MCl).)(7) Density of total oxidant or fuel:kwhere (HjOr ) is the molar enthalpy of the jth reactant,(Jlkg-mole)(k)](2) Specific enthalpy of total reactant ho divided by
RT, kg-molelkg:
h (2) + (!!.-) h (I)
RT
f RT

p(k)

=

1

_

(9.12)NREAC n~k) M~k)

L

J

j=l

}

pt)

(8) Density of total reactant:(k = 1,2)(9.7)

7 1] p(l)p~)

1 + !!.-

+

f

Po =

(3) Specific internal energy of total oxidant(k = 1) or total fuel (k = 2) (u')(k) divided by RT,(kg-molelkgi k):(U')(k)(k = 1,2)

RT

(9.8)

k

where (UjOr ) is the molar internal energy of the jth
reactant, (Jlkg-mole)Yk).
(4) Specific internal energy of total reactant
divided by RT, kg-molelkg:

un(U')C2) + (!!.-) (U,)(I)
RT
f RT

(9.9)

p(l) +(7)p~)(9.13)In the main program of CEA several alternative
expressions are given for specifying the relative amounts
of total oxidant to total fuel and for relating them to off.
Two expressions are given for equivalence ratio-achemical equivalence ratio r and a fuel-to-air (or fuel-tooxidant) equivalence ratio <1>. If the chemical equivalence
ratio is specified, it will be necessary to make use ofoxidation states. Letandbe positive andnegative oxidation states of the ith element in its
commonly occurring compounds. At least one of thesewill be zero. Thus, for example, the negative oxidation
state for chlorine is -1 and its positive oxidation state iszero. The oxidation states per kilogram of total oxidant
or total- fuel arevtVi

1 + !!..

f(k = 1,2)(9.14)(k = 1,2)(9.15)(5) Molecular weight of total oxidant (k = 1) or
total fuel (k = 2):
~

I __

M(k)

NREAC
"

LJ nj(k)j=l42

V-(k) = "(k = 1,2)(9.10)

LJ

v.- b~k)JJ

i=lThe positive and negative oxidation states for the total
reactant are thenV<(2) +(7)v<(I)(9.16)
1+~f

V-(2) +

(7oxidizer (e.g., H20 2), then for off-stoichiometricconditions r and <I> differ for the same reaction. For
example, for CH30H + 1/2 02' r = 2 and <P = 3.
One of the options in the SHOCK problem requires
reactant compositions relative to the total reactant. Thesemay be obtained as follows:)v-(I)o(9.17)

1+-f

The chemical equivalence ratio is now defined asn~k)(9.18)

for k=2

= _J_
1+~

f

The equivalence ratio <1>, commonly usedengineering practice, is defined as

III

1
a

(9.19)(j = 1,...,NREAC)(9.20)

where mj is kilogram-moles of jth reactant per kilogramof total reactant.
Specific heat and entropy for the total reactant inkilogram-moles per kilogram can be calculated by using
the mj' C;i' and Sj of the j reactants:

where i/a is the fuel-to-air weight ratio. As Gordon(1982) points out, the two equivalence ratios r and <I> arealways identical for the stoichiometric condition. In
addition, they are always identical when all the positive
valence atoms are in the fuel and all the negative valence
atoms are in the oxidizer (eg., CH4 + 02)' However, ifsome negative valence atoms are in the fuel (e.g.,
CH 30H) or if some positive valence atoms are in the

Co
R

So
R

NREAC

CO.

j=l

R

NREAC

mj(~l

L mj---.!!:LL
j=l(9.21)(9.22)

43

Appendix-Symbols
A,B,C,Dtransport property coefficients, eq. (5.1)

A

area, m2

A/At

ratio of nozzle exit area to throat area,eq. (6.14)

A *kQ

collision cross-section ratio, eq. (5.13)

Ak
a

C;J

J/(kg-mole)/K)

C~J

molar heat capacity at constant volume
for standard state for species j,
J/(kg -mole)/K)

chemical formula of species k,dimensionless, eq. (5.10)

Cospecific heat at constant pressure of total
reactant, Jlkg-K, eq. (9.21)

velocity of sound, mis, eqs. (2.67)
and (2.74)

cp

specific heat at constant pressure of
mixture, equilibrium or frozen, Jlkg-K

ai(i = 1...9) least-squares coefficients, eqs. (4.6)

cp,eto (4.11)
a··
I)

(k)

a··I)

molar heat capacity at constant pressure
for standard state for species or reactant j,

stoichiometric coefficients, kilogramatoms of element i per kilogram-mole ofspecies j, (kg-atom)!Ckg-mole)j' eq. (2.7a)
stoichiometric coefficients, kilogramatoms of element i per kilogram-mole ofreactant j (oxidant if k = 1, fuel if k = 2),(kg-atom)~k)/(kg-mole)Y), eq. (9.1)

Cp,eqequilibrium specific heat at constantpressure, Jlkg-K, eq. C2.49a)
equilibrium specific heat at constantpressure for transport properties,
Jlkg-K, eq. (5.14)

cp ,!

frozen specific heat at constant pressure,Jlkg-K, eq. (2.49b)

cp,fr

frozen specific heat at constant pressure
for transport properties, Jlkg-K, eq. (5.15)

as,Tvelocity of sound, mis, defined by,eq. (3.10)

cp,rreaction specific heat at constant pressure,Jlkg-K, eq. (2.49c)b.I

kilogram-atoms of element i per kilogramof mixture, (kg-atom)ilkg, eq. (2.7c)

cp,rereaction specific heat at constant pressure
for transport properties, Jlkg-K, eq. (5.16)b(k)assigned kilogram-atoms of element i per
kilogram of total oxidant (k = 1) or total
fuel (k = 2), (kg-atom)~k)lkg(k)' eq. (9.1)

I

bC?l

CF

assigned kilogram-atoms of element i per
kilogram of total reactant, (kg-atom)/kg,
eq. (9.5)

Cv

specific heat at constant volume of
mixture, equilibrium or frozen, Jlkg-K,
eq. (2.70) or (6.39)

C*

characteristic velocity, mis, eq. (6.11)

DkQ

binary diffusion coefficient, fi 2/s,
eq. (5.13)

d··I)

matrix coefficient in eq. (5.17)

coefficient of thrust, eq. (6.13)

45

Helmholtz energy of mixture withconstraints, defined by eq. (2.33), Jlkg; orforce, N, eq. (6.6)

F%F

percent of total fuel in total reactant by
weight or massfHelmholtz energy of mixture, Jlkg,
eqs. (2.30) and (2.31)

G

Gibbs energy of mixture with constraints,
Jlkg, defined by eq. (2.8)

g

Gibbs energy of mixture, Jlkg, eq. (2.5)

gcconversion factor, eqs. (6.6) to (6.8)

H?
(H?)(k)M·J

molecular weight of species j, eqs. (2.3b)and (2.4)

Mo

molecular weight of total reactant,kglkg-mole, eq. (9.11)

Alk)molecular weight of total oxidant (k = 1)
or total fuel (k = 2), (kglkg-molei k), eq.
(9.10)

A-I:k)

molecular weight of reactant j (oxidant if
k = 1; fuel if k = 2, (kglkg-mole»)k),eqs. (9.4)]

MW

molecular weight of mixture, kg/kg-mole,eqs. (2.4)

molar standard-state enthalpy forspecies j, J/(kg-mole)jJ(

Mach number, eq. (6.10)

molar enthalpy of reactant j (k = 1 foroxidant, k = 2 for fuel), J/(kg-mole»)k)'eq. (9.6)

mj

kilogram-moles of reactant j per kilogramof total reactant, (kg-mole)jlkg, eq. (9.20)

rhmass flow rate, kg/s, eq. (6.4)!1 HO(nfheat of formation at temperature T,
Jlkg-mole, eq. (4.2)

M!1,H°(nheat of reaction at temperature T,
J/(kg-mole), eq. (5.9)

nkilogram-moles of mixture per unit mass,
kg-molelkg, eq. (2.2)

hspecific enthalpy of mixture, Jlkg,'eq. (2.14)

n·kilogram-moles of species j per kilogramof mixture, (kg-mole)/kg, eq. (2.2)h(k)

specific enthalpy of total oxidant (k = 1)
or total fuel (k = 2), (J/kgi k), eq. (9.6)n(k)

kilogram-mole per kilogram of reactant j(oxidant if k = 1; fuel if k = 2),(kg-mole)jk)/kg(k), eq. (9.2) or (9.3)

J

J

J

k)

J

J

kilogram-moles of reactant j (oxidant if
k = 1; fuel if k = 2), (kg-mole»)k)specific enthalpy of total reactants, J/kg,
eq. (9.7)

offoxidant-to-fuel weight (or mass) ratioh'term defined by right side of eq. (7.32)

Ppressure, N/m 2h"term defined by right side of eq. (8.2)

Pin/Peh*term defined by right side of eq. (7.7)

ratio of combustion pressure to exit
pressure for lAC modelI

specific impulse, N/(kg/s) or mis,eq. (6.7)

PoP'

assigned or initial pressure, N/m 2
term defined by right side of eq. (7.31)

I sp

specific impulse with exit and ambientpressures equal, N/(kg/s) or mis, eq. (6.8)

P"

term defined by right side of eq. (8.1)

p*term defined by right side of eq. (7.6)

I vacvacuum specific impulse, N/(kg/s) or mis,eq. (6.9)

Pr

Prandtl number, eqs. (5.19)

M

molecular weight of mixture, kglkg-mole,eqs. (2.3)

Runiversal gas constant, 8314.51 J/(kg-mole)(K)M·latomic weight of chemical element i,(kglkg-atom)i' eq. (9.4)

requivalence ratio based on oxidation
states, eq. (9.18)

ho

46

s·J

molar entropy of species j,
J/(kg-mole)jCK), eq. (2.17)

Ws

incident shock front velocity, mis,eq. (7.2)

S?J

molar standard-state entropy for species j,
J/(kg -mole)jCK)

Xi

unknowns in eq. (5.16)

x·1

S

specific entropy of mixture, J/kg-K,
eq. (2.16)

Sospecific entropy of total reactants, J/kg-K,
eq. (9.22)

mole fraction of species i relative to all
species (gaseous and condensed) in eq.
(2.4a) and relative to NM gaseous species
used for thermal transport property
calculations in eqs. (5.3), (5.4), (5.12),(5.15), (5.16), and (5.18)

T

temperature, Kur;Jmolar standard-state internal energy forspecies j, J/kg-mole, eq. (2.38)

u ik

stoichiometric coefficients of species k inreaction i for transport property
calculations, eq. (5.10)(U?)Ck)molar internal energy of reactant j(oxidant if k = 1; fuel if k = 2),(Jlkg-mole»)k), eq. (9.8)

Ukparameter defined by equation (8.16)'Yspecific heat ratio, eq. (2.72)'Ys

isentropic exponent, eqs. (2.71) and (2.73)

11i

viscosity of species i, ~, eq. (5.3)

11ij

viscosity interaction parameter, JlP,
eq. (5.7)

11 mixA

viscosity of mixture, ~, eq. (5.3)

J

velocity (in shock problems, velocity
relative to incident or reflected shock
front), mis, eqs. (6.5) and (7.1)

uu'specific internal energy of mixture, Jlkg,
eq. (2.38)

u'Ck)

specific internal energy of total oxidant(k = 1) or total fuel (k = 2), (Jlkg)Ck),

eq. (9.8)

control factor, eqs. (3.3) and (7.22);
thermal conductivity, eq. (5.1)

Ai

Lagrangian multipliers, eqs. (2.8) and(2.33); thermal conductivity of species i,~W/cm-K, eq. (5.4)

A1,A 2,A3,A4

control factors, eqs. (3.1), (3.2), (7.20),
and (7.21)

Aeq

equilibrium thermal conductivity of
mixture, JlW/cm-K, eq. (5.2)

Afrfrozen thermal conductivity of mixture,
flW/cm-K, eq. (5.4)

Ar,ivariables in eqs. (5.8) and (5.11)

negative oxidation state of total oxidant(k = 1) or total fuel (k = 2), eq. (9.15)

Arereaction thermal conductivity of mixture,~W/cm-K, eq. (5.8)

actual velocity of gases in shock tube,
mis, eq. (7.1)

fljchemical potential of species j,
J/(kg-mole)j' eqs. (2.11) and (2.35)

Mk)

weight of jth reactant (oxidant if k = 1;
fuel if k = 2), kgY\ eq. (9.2)

Jl?J

standard-state chemical potential of
species j, J/(kg-mole)j' eq. (2.11)

Wshock front velocity, mis, eq. (7.1)

1t e

Lagrangian multiplier for ions divided bywR

reflected shock front velocity, mis,eq. (7.4),
Uospecific internal energy of total reactant,Jlkg, eq. (9.9)

Vspecific volume, m3/kg, eq. (2.1a)yt

positive oxidation state of total reactant,eq. (9.16)

yytCk)

negative oxidation state of total reactant,eq. (9.17)
positive oxidation state of total oxidant(k = 1) or total fuel (k = 2), eq. (9.14)

y-Ck)

v

J

RT

47

1t.
l

p

Lagrangian multiplier for element i
divided by RT, eqs. (2.18), (2.19), (2.40),
and (2.41)

Ppressures

entropydensity of mixture, kg/m3, eq. (2.1b)

stoichstoichiometricdensity of total reactant, kg/m3, eq. (9.13)

T

temperaturetthroat0

total reactant; zeroth iteration

1,2,5

stations

density of total oxidant (k = 1) or total
fuel (k = 2), (kg/m 3 k ), eq. (9.12)

i

density of reactant? (oxidant if k = 1; fuelif k = 2), (kg/m3
eq. (9.12)))k ,frozen thermal conductivity interaction
coefficient between species i and j,
eqs. (5.4) and (5.6)
<p

<Pij

equivalence ratio based on f/a ratios,eq. (9.19)

Superscripts:osymbol for standard state; an assigned or
initial conditionk

1 if oxidant, 2 if fuelIndices:number of chemical elements (if ions are
considered, number of chemical elements
plus one)viscosity interaction coefficient between
species i and j, eqs. (5.5) and (5.7)

Subscripts:

NG

number of gaseous species in thermodynamic data for current chemical system

NMnumber of gaseous species for transport
property calculations

gas

NR

combustor station assuming infinite
chamber area

number of reactions for transport property
calculations

NREAC

number of reactantsNS

number of species, gaseous and condensed, in thermodynamic data for current
chemical systemaambient; assignedcend of combustion chamber; condensedeexit; electrons

g

inf
inj

combustor inlet or injector station for
finite-area combustorkiteration kmmelting48

References
Brokaw, R.S., 1960, "Thermal Conductivity of Gas
Mixtures in Chemical Equilibrium," Journal ofChemical Physics, VoL 32, No.4, pp. 1005-1006.

Gaydon, A.G. and Hurle, LR., 1963, The Shock Tube inHigh-Temperature Chemical Physics, Reinhold
Publishing Corp.

Burcat, A., Zeleznik, F.J., and McBride, B.J., 1985,"Ideal Gas Thermodynamic Properties for the
Phenyl, Phenoxy and o-Biphenyl Radicals," NASA
TM-83800.

Glass, I.L and Hall, J.G., 1959, "Handbook of Supersonic Aerodynamics, Section 18. Shock Tubes,"
NAVORD Rep. 1488, Vol. 6, Bureau of Naval
Weapons.

Butler, J.N. and Brokaw, R.S., 1957, "Thermal Conductivity of Gas Mixtures in Chemical Equilibrium,"
Journal of Chemical Physics, Vol. 26, No.6,
pp. 1636-1643.

Goldman, D.T. and Bell, R.J., 1986, "The International
System of Units (SI)/' NBS Special Publication 330,
1986 Ed.

Chase, M.W., Jr., Davies, C.A., Downey, J.R., Jr.,
Frurip, D.J., McDonald, R.A., and Syverud, A.N.,
1985, "JANAF Thermochemical Tables," Third ed.,Pts. I-II, Journal of Physical and Chemical
Reference Data, Vol. 14, SuppL 1, American
Chemical Society, Washington, DC.
Cohen, E.R. and Taylor, B.N., 1987, "The 1986
CODATA Recommended Values of the Fundamental
Physical Constants," National Bureau of Standards,
Journal ofResearch, Vol. 92, Mar.-Apr., pp. 85-95.
Cox, J.D., ed., 1982, "Notation of States and Processes,
Significance of the Word . Standard' in Chemical
Thermodynamics, and Remarks on Commonly Tabulated Forms of Thermodynamic Functions," Pureand Applied Chemistry, Vol. 54, pp. 1239-1250.

Gordon, S., Zeleznik, F.J., and Huff, V.N., 1959, "A
General Method for Automatic Computation of
Equilibrium Compositions and Theoretical Rocket
Performance of Propellants," NASA TN D-132.
Gordon, S. and Zeleznik, F.J., 1962, "Thermodynamic
Extrapolation of Rocket Performance Parameters,"
American Rocket Society Journal, VoL 32, No.8,
pp. 1195-1202.
Gordon, S. and Zeleznik F.J., 1963, "A General IBM
704 or 7090 Computer Program for Computation ofChemical Equilibrium Compositions, Rocket Performance, and Chapman-Jouguet Detonations,
Suppl. 1: Assigned Area-Ratio Performance," NASATN D-1737.

Cox, J.D., Wagman, D.O., and Medvedev, V.A., 1989,
CODATA Key Values for Thermodynamics, Hemisphere Publishing Corp., New York.

Gordon, S., 1970, "Calculation of Theoretical Equilibrium Nozzle Throat Conditions When Velocity ofSound Is Discontinuous," American Institute of
Aeronautics and Astronautics Journal, Vol. 9, No.1,
pp. 179-182.

Garvin, D., Parker, V.B., and White, H.J., Jr., 1987,
CODATA Thermodynamic Tables, Hemisphere
Publishing Corp., New York.

Gordon, S. and McBride, B.J., 1976, Computer Programfor Calculating Complex Chemical Equilibria,
Rocket Peifonnance, Incident and Reflected Shocks,
49and Chapman-Jouguet Detonations, NASA SP-273,(Revision of 1971 publication).

Ideal Gas State," Journal of Chemical Physics.,
Vol. 35, No.6, pp. 2198-2206.

Gordon, S., 1982, "Thermodynamic and Transport
Combustion Properties of Hydrocarbons With Air,
I-Properties in SI Units," NASA TP-1906.

McBride, B.J., Heimel, S., Ehlers, J.G., and Gordon S.,
1963, Thermodynamic Properties to 6000 Kfor 210
Substances Involving the First 18 Elements, NASA
SP-3001.

Gordon, S., McBride, B.J., and Zeleznik, F.J., 1984,"Computer Program for Calculation of Complex
Chemical Equilibrium Compositions and Applications, Supplement I-Transport Properties," NASA
TM-86885.
Gordon, S. and McBride, BJ., 1988, "Finite Area
Combustor Theoretical Rocket Performance," NASA
TM-100785.

McBride, BJ. and Gordon S., 1967, "FORTRAN IV
Program for Calculation of Thermodynamic Data,"
NASA TN D-4097.
McBride, B.J. and Gordon, S., 1992, Computer Programfor Calculating and Fitting Thermodynamic Functions, NASA RP-1271.

Griem, H.R., 1964, Plasma Spectroscopy, McGraw-Hill
Book Co., Inc.

McBride, B.J., Gordon, S., and Reno, M.A., 1993a,"Thermodynamic Data for Fifty Reference Elements," NASA TP-3287.

Gurvich, L.V., Veyts, I.V., and Alcock, C.B., 1989,
Thermodynamic Properties ofIndividual Substances,Vol. 1, Pts. 1-2, Fourth ed., Hemisphere Publishing
Corp., New York.

McBride, B.J., Gordon, S., and Reno, M.A., 1993b,"Coefficients for Calculating Thermodynamic andTransport Properties of Individual Species," NASA
TM-4513.

Hirschfelder, J.D., Curtiss, C.P., and Bird, R.B., 1954,
Molecular Theory ofGases and Liquids, John Wiley& Sons, Inc.

Radhakrishnan, K. and Bittker, D.A., 1994, "LSENS, A
General Chemical Kinetics and Sensitivity Analysis
Code for Homogeneous Gas-Phase Reactions. II.
Code Description and Usage," NASA RP-1329.

Huff, V.N., Gordon, S., and Morrell, V.E., 1951, "General Method and Thermodynamic Tables for Computation of Equilibrium Composition and Temperatureof Chemical Reactions," NACA Report 1037.

Sutton, G.P. and Ross, D.M., 1976, Rocket Propulsion
Elements, Fourth ed., John Wiley & Sons, Inc.

Kirkwood, J.G. and Oppenheim, I., 1961, Chemical
Thermodynamics, McGraw-Hill Book Co., Inc.

Svehla, R.A., 1962, "Estimated Viscosities and Thermal
Conductivities of Gases at High Temperatures,"NASA TR R-132.

Lewis, G.N. and Randall, M., 1961, Thermodynamics,
Second ed., McGraw-Hill Book Co., Inc., pp.
667-668.
Marsh, K.N., Das, A., Frenkel, M. Gadalla, N.M., andWilhoit, R.C., 1988, TRC Thermodynamic Tables,
Non-Hydrocarbons, Vols. I-VIII; and Hydrocarbons,
Vols. I-XII, Thermodynamics Research Center,
Texas A&M University, College Station, TX.
McBride, B.J. and Gordon S., 1961, "Thermodynamic
Functions of Several Triatomic Molecules in the

50

Svehla, R.A. and Brokaw, R.S., 1966, "Heat Capacity
and the Lewis Number of a Reacting Gas,"
American Institute of Aeronautics and Astronautics
Journal, Vol. 4, No.1, pp. 182-184.
Svehla, R.A. and McBride, BJ., 1973, "FORTRAN IV
Computer Program for Calculation of Thermodynamic and Transport Properties of Complex
Chemical Systems," NASA TN D-7056.
Svehla, R.A., 1995, "Transport Coefficients for theNASA Lewis Chemical Equilibrium Program,"NASA TM-4647.

Zeleznik, F.J. and Gordon, S., 1960, "An Analytic
nvestigation of Three General Methods ofCalculating Chemical-Equilibrium Compositions,"
NASA TN D-473.
Zeleznik, F.J. and Gordon, S., 1961, "Simultaneous
Least-Squares Approximation of a Function and Its
First Integrals With Application to Thermodynamic
Data," NASA TN D-767.
Zeleznik, F.J. and Gordon, S., 1962a, "A General IBM
704 or 7090 Computer Program for Computation ofChemical Equilibrium Compositions, Rocket Performance, and Chapman-Jouguet Detonations,"·
NASA TN D-1454.
Zeleznik, P.J. and Gordon, S., 1962b, "Calculation of
Detonation Properties and Effect of Independent
Parameters on Gaseous Detonations," American
Rocket Society Journal, Vol. 32,No. 4,pp.60~15.
Zeleznik, FJ. and Gordon, S., 1968, "Calculation ofComplex Chemical Equilibria," Industrial andEngineering Chemistry, Vol. 60, No.6, pp. 27-57.

51

lJI.

tv

TABLE 2.1-ITERATION EQUATIONS TO DETERMINE EQUILffiRIUM COMPOSmONS FOR EITHER ASSIGNED TEMPERATURE AND PRESSURE,ENTHALPY AND PRESSURE, OR ENTROPY AND PRESSURE
Equation

Variable~.. .

alJa2jnJ...

~alJnJ

L a2ja2jnJ
i-I

.. .

NO

NO

L a.)alJ'li
i-I

E atN)J-l.. .

a 1,NG+la2,NG+l...

1tl(2.24)NONO

L alJaVnj
i-I

L

NO

L
J-l

(2.25)

NO

NO

E
~ni
J-l

NO

C(2.28)ES

E alJn)::.l.R
1-1

NO

NO

a1,NG+l

. ..

alNs

a2,NG+l...

a2NS

E a,p,fJi-Ial,NG+l...

al,NG+l0. ..

L atf'.tJ
i-I

E a~.tJ/-1...

Right side

. ..

0

a~s

0

0

...

H~O+1
RT

. ..

/-1

Ir.

)-1

S

E
tlzJ nl -!.
j-l
R

NO

H

C

...

Ea,fJ ~
i-I...

L a,f',I-I

NO

aColumn not used for assigned temperature and volume.
bRow used only for assigned enthalpy and pressure.
CRow used only for assigned entropy and pressure.

S

i

SNQ+l

R

. ..

RT
SNS

R

~

i*1

RT

NO

~

J-1

RT

E alJni

E a,f) itT
j-l(b~ - bJ + E atf'J ..!:.L

0

H~Q+l
RT

~NQ+l

~s

~NS

RT

RT

)-1

J-1

H

NO

C

RT

If'.

NO

NO

If'.

NO

E n.::.L
RT

1-1

J

NO
S
En::.l.
J-l 'j R

NS

C'.-

NO

Ir. Ir.

R

i-I'j RT RTh

E n -pJ+ L n ::.L - i
i-I

'j

~

NO

n-En +En-.-L)-1 'ji-I 'j RT

i-I

H';s

NO

(b~ - b2 ) + L ~nJ..!:.L

E ni - n L
nJ::.L
/-1
RT

0

~.

~nJ

L

NO

Eav"'JNO(b~ - b1) + L avn/:.L
J-l
RT

HO

NO

0

0

HO
NO
J
L alJnJ J-l
RT

E a2jnJ-l..
/-1
RT

NO

a~s

NO

.. .

L aljnJi-I
NONO

E a2jnJ::L
RT
NO

~T

~n

.1.nNSNO

a2NS

E aljnJ::L
i-I
RT

...

NO

E
alJni
i-I
NO

b(2.27).1.nNG+1

(a)

NO

alNS

(2.26)

j-l

1t~

NS

C'.-

NO

~S

i-I

R

i-I

RT R

E nJ..:..El. + E nJ..i::.l.

So

Ir. ~
En ::L..!:.L

_0
___

h

+ NO

RT __ KIi-I
NO'j RT RT
NO

S~...Q-!...+n-Ln+ Ln ::.l....!:1.
RR
J-l 'jj-l 'j R RT

TABLE 2.2-ITERATION EQUATIONS TO DETERMINE EQUILIBRIUM COMPOSmONS FOR EITHER ASSIGNED TEMPERATURE AND VOLUME,
INTERNAL ENERGY AND VOLUME, OR ENTROPY AND VOLUME
Equation

Variable
1t1

NO(2.45)

E a Va1JnJ
J-l(2.46)NO

L a1Ja2JnJ
i-I

NO

NO

L aQfLJf'J
J-lal,NG+la2,NG+lC(2.48)

NO

NO

l?,

NO

R

aColumn not used for assigned temperature and volume.

~ow used only for assigned internal energy and volume.
CRow used only for assigned entropy and volume.

VI

v.>

...

...

E
a2JnJ::.L
J-1
RT
j-l

...

...

(S ) L ~n) (S-'. - 1)

L aVn) J.. - 1
J-l
R

M NS

al,NG+l

. ..

alNS

L a1JnJ~
J-1
RT

aZ,NG+l

...

aZNS

L

L atpfflJ
J-1

aQ,NG+l

. ..

aQNS

aQ,NG+l

0

. ..

0

E
avaQflJJ-l
L ¥tflJ
J-l

NO

aZNS

l?,

. ..

NO

...
...

L
)-1

NO

. ..

0

aQNS

UO
aljnJ::LRT(S )L aqnj ::1 - 1
i-I
R

U~O+I
RT
SNO+l

R

. ..

. ..

u?

(b~ - b1)

RT
SNS

R

j.L.

J-1

RT

NO

j.L

J-l

RT

NO

1.1.

i-1

RT

(b~ - b1) +L ~nJ.!:1...

NOUO
L aPJnJ-.L
j-1
RT(b~ - b~ +L aqnJ..:::.L

U~O+1

j.LNGoI

RT

RT

j-1

U~s
RT

u~s

NO

+E a1JnJ...2UO
a2JnJ::L
RT

NO

0

Right side

MnT
(a)
NO

NO

L aQ
P0J
i-1

NO

MNG+l

1tQ

NO

...

N

NO

L alJni-.L
j-l
RT

...

NO

E a Va1JnJ

L a1JaVnJ
J-1

alNS

b(2.47)1t2

NS

I.I.NS

RT

CO

NO

uO uo

u'

,

RT

RT

NO

l?,j.L

J-l

J RT RT+L n..:L.!:1...

En)::.::J. +L ni::L::..L
)-1
R
j-1
RT RT

....Q.. - ~

NS
CO
NO
UO(S
)
Ln::.::J. + En::..L...l - 1
j-l J R)-1 J RT R

S
(S)
.J!-.!+En::1-1..2
NO

R

R

j-1

1.1.

J R

RT

Ut

+:-.

TABLE 2.3-EQUATIONS FOR EVALUATING DERIVATIVES WI1H RESPECT TO
LOGARITHM OF TEMPERATURE AT CONSTANT PRESSURE
Equation

Variable

aln T

...

~

a~G+1
aln T

E
aVDlJnJ
J-1

E DVa2JnJ., .

E
alJa"nJ
J-l

NO

NO

E
~alJnJ
J-1

E
~~nJ
J-1

...

NO

...

NO

(2.56)

NO

E a"alJnJ(2.57)(2.58)

ax 2

a'lt 1
aIn T
NO

aIn T

aIn n
aIn T

Right sideal,NG+l

"

.

a INS

E
alJnJ
J-l

-E
alJnJ::L
J-1
RT

a2,NG+l

. ..

aZNS

E
~nJ
J-1

E alJa"nJ
J-1

aG,NG+l

...

aINS

E
alJnJ
J-1

-E alJnJ::L
RT-H~G+1--NOJ-1

E DfJa2JnJ

.

~s

"

aIn T

NO

E
~a"nJ
J-1

NO

NO

NO

NO

'-1

'-1

al,NG+la2,NG+l...

aG,NG+l

0

. ..

0

0

aINS

aZNS

...

aINS

0

. ..

0

0

NO

NO

.. .

:E
alJnJ
i-I

0

...

0

0

E olJnJ
J-1

E
~nJ
J-l

NO

NO

HO

NO

HD

NO

HD

-E
~nJJJ-1
RT

j-1

RT

H~s
RT
NO

HD

-E
n ::J~ RTi-I

TABLE 2.4--EQUATIONS FOR EVALUATING DERIVATIVES WITH RESPECT TO
LOGARITHM OF PRESSURE AT CONSTANT TEMPERATURE
Equation

Variable

l1lt 1

Cht 2

aIn P

aIn P

...

~
aIn P

~(]+1
aIn P

E
aVaVnJj-l

E a Va2JnJ.. .

E QlJaPJnJ

al,NG+lNO

NOE~alJnJj-l

E~a2JnJ
j-l...

E~afJnJ
j-1

a2,NG+lNO

a4,NG+lNO(2.64)

NO(2.65)(2.66)

Vt
Vt

NO

J-1

~s

aIn n
aIn P

Right side

...

a iNS

E
aVnJ
J..l

E
allnJ
j-1

NO

NO

...

a 2NS

E
~nJ
j-l

E
~nJ
j-t

"

.

atNS

. ..

NO

NO

J-l

E
alJ~nJJ-l...

E alJa"nJ
J-1al,NG+la2,NG+l...

at,NG+l

0

"

.

"INS

ams

...

a~s

0

...

E
a'lnJ
j-1

NO

0

"

NO

NO

E
~nJ
J-l

NO

NO

E
aPJallnJ
j-l

E
"lJnJ
J-1

aIn P

NO

NO

NO

E
alJnJ
J-l

E
aPJnJ
j-l

0

0

0

. ..

0

0

0

.

0

0

NO

E
nJ
j-l

I

REPORT DOCUMENTATION PAGE

Form Approved,
OMB No. 0704-0188

Public reporting burden for this collection of Information is estimated to average 1 hour per response, including the time for reviewing Instructions, searching existing data sources,
gathering and maintaining the data needed, and completing and reviewing the collection of Information. Send comments regarding this burden estImate or any other aspect of thiscollection of Informatlon, Including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson
Davis Highway, Suite 1204, Arlington. VA 22202-4302, and to the Office of Management and BUdget, Paperwork Reduction Project (0704-0188). Washington, DC 20503.

2. REPORT DATE

1. AGENCY USE ONLY (Leave blank)

October 1994

1

3. REPORT TYPE AND DATES COVEREDReference Publication

1

5. FUNDING NUMBERS

4. TITLE AND SUBTITLE

Computer Program for Calculation of Complex Chemical Equilibrium
Compositions and Applications
I. Analysis

WU-505--62-52

6. AUTHOR(S)Sanford Gordon and Bonnie J. McBride

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)8. PERFORMING ORGANIZATIONREPORT NUMBER

National Aeronautics and Space Administration
Lewis Research Center
Cleveland, Ohio 44135-3191

E-8017

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)10. SPONSORING/MONITORINGAGENCY REPORT NUMBER

National Aeronautics and Space Administration
Washington, D.C. 20546-0001

NASA RP-1311

11. SUPPLEMENTARY NOTES

Sanford Gordon, Sanford Gordon and Associates, Cleveland, Ohio, and Bonnie J. McBride, NASA Lewis Research Center.
Responsible person, Bonnie J. McBride, organization code 2670, (216) 433-5870.
12a. DISTRIBUTION/AVAILABILITY STATEMENT

12b. DISTRIBUTION CODE

Unclassified - Unlimited
Subject Categories 20 and 25

13. ABSTRACT (Maximum 200 words)

This report presents the latest in a number of versions of chemical equilibrium and applications programs developed at
the NASA Lewis Research Center over more than 40 years. These programs have changed over the years to include
additional features and improved calculation techniques and to take advantage of constantly improving computer capabilities. The minimization-of-free-energy approach to chemical equilibrium calculations has been used in all versions of
the program since 1967. The two principal purposes of this report are presented in two parts. The first purpose, which isaccomplished here in part I, is to present in detail a number of topics of general interest in complex equilibrium calculations. These topics include mathematical analyses and techniques for obtaining chemical equilibrium; formulas for
obtaining thermodynamic and transport mixture properties and thermodynamic derivatives; criteria for inclusion ofcondensed phases; calculations at a triple point; inclusion of ionized species; and various applications, such as constantpressure or constant-volume combustion, rocket performance based on either a finite- or infinite-chamber-area model,
shock wave calculations, and Chapman-Jouguet detonations. The second purpose of this report, to facilitate the use of thecomputer code, is accomplished in part II, entitled "Users Manual and Program Description." Various aspects of thecomputer code are discussed, and a number of examples are given to illustrate its versatility.

15. NUMBER OF PAGES

14. SUBJECT TERMS

Chemical equilibrium; Combustion products; Combustion temperatures; Computer
program; Thermodynamic mixture properties; Thermal transport properties; Rocket
performance
17. SECURITY CLASSIFICATIONOF REPORT

18. SECURITY CLASSIFICATION
OF THiS PAGE

Unclassified

Unclassified

NSN 7540-01-280-5500

19. SECURITY CLASSIFICATION
OF ABSTRACT

61
16. PRICE CODE

A04
20. LIMITATION OF ABSTRACT

Unclassified
Standard Form 298 (Rev. 2-89)
Prescribed by ANSI Std. Z39-18
298-102