Thermochemical Performance - Rocket Propulsion notes

This module covers the thermochemistry of rocket propulsion, including combustion fundamentals, propellant chemistry, and methods for calculating combustion chamber conditions and performance parameters.

# this line makes figures interactive in Jupyter notebooks
%matplotlib inline
from matplotlib import pyplot as plt

import numpy as np
from scipy.optimize import root_scalar
from scipy import constants

from pint import UnitRegistry
ureg = UnitRegistry()
Q_ = ureg.Quantity

# for convenience:
def to_si(quant):
    '''Converts a Pint Quantity to magnitude at base SI units.
    '''
    return quant.to_base_units().magnitude

# define some constants
Ru = Q_(constants.R, 'J/(K*mol)')
g0 = Q_(constants.g, 'm/s^2')

# from previous modules

def get_area_ratio(pressure_ratio, gamma):
    '''Calculates area ratio based on specific heat ratio and pressure ratio.
    pressure ratio: chamber / exit
    area ratio: exit / throat
    '''
    return (
        np.power(2 / (gamma + 1), 1/(gamma-1)) * 
        np.power(pressure_ratio, 1 / gamma) *
        np.sqrt((gamma - 1) / (gamma + 1) /
                (1 - np.power(pressure_ratio, (1 - gamma)/gamma))
                )
        )

def root_area_ratio(pressure_ratio, area_ratio, gamma):
    ''' Returns zero for a given area ratio, pressure ratio, and gamma.
    pressure ratio: chamber / exit
    area ratio: exit / throat
    '''
    return area_ratio - get_area_ratio(pressure_ratio, gamma)

def get_thrust_coeff(pressure_ratio, gamma):
    '''Calculates thrust coefficient for optimum expansion.
    pressure ratio: chamber / exit
    area ratio: exit / throat
    '''
    return np.sqrt(
        2 * np.power(gamma, 2) / (gamma - 1) * 
        np.power(2 / (gamma + 1), (gamma + 1)/(gamma - 1)) * 
        (1 - np.power(1.0 / pressure_ratio, (gamma - 1)/gamma))
        )

def get_cstar(Tc, MW, gamma):
    '''Calculates cstar using chamber properties.'''
    return (
        np.sqrt(Ru * Tc / (gamma * MW)) * 
        (2 / (gamma + 1))**(-0.5*(gamma + 1)/(gamma - 1))
        )

Overview¶

The key performance parameters for a rocket engine are:

$c^*$ — characteristic velocity (combustion efficiency)
$C_F$ — thrust coefficient (nozzle efficiency)
$I_{\text{sp}}$ — specific impulse
$T$ — thrust

These parameters depend on the following thermochemical properties:

Chamber temperature: $T_c$
Molecular weight of combustion products: $\overline{\textrm{MW}}$
Ratio of specific heats: $\gamma = c_p/c_v$
Heat of combustion: $q_R$

where $c^*$ particularly depends on these properties as the primary thermochemical/energetic performance parameter of a rocket:

c^* = \frac{p_c A_t}{\dot{m}} = \sqrt{\frac{\mathcal{R}_u T_c}{\textrm{MW} \, \gamma} \left( \frac{\gamma + 1}{2} \right)^{\frac{\gamma+1}{\gamma-1}} } \;.

(1)

Up to this point, we have been provided or assumed these values, but they are actually functions of:

The propellant combination (fuel and oxidizer types)
The oxidizer-to-fuel ratio (O/F or $r$ )

Figure 1:Schematic of liquid bipropellant rocket engine showing: propellant injection (LH₂ and LO₂), combustion chamber (where $T_c$ , $\overline{MW}$ , $\gamma$ are established), throat section (where $c^*$ is evaluated), and nozzle expansion (where $C_F$ is determined). Exit composition $\vec{X}_e$ may differ from chamber composition depending on whether flow is frozen or in shifting equilibrium.

Propulsion System Categories¶

The relationship between chamber temperature and propellant properties differs by propulsion type:

Nuclear Thermal/Electrothermal Propulsion¶

$T_c$ is a design parameter (set by heater power)
Temperature is independent of propellant chemistry
Goal: choose propellant to optimize equilibrium mole fractions $X_i$ , $\gamma$ , $\overline{\textrm{MW}}$

Chemical Propulsion¶

$T_c$ is determined by combustion chemistry
$T_c$ , $X_i \rightarrow \gamma$ , $\overline{\textrm{MW}}$ all depend on:
- Propellant combination
- O/F ratio

In either case, once $T_c$ , $\gamma$ , $\overline{\textrm{MW}}$ are determined, then $c^*$ , $C_F$ , and $I_{\text{sp}}$ can be found.

But, where do the chamber properties come from? Given a propellant or combination of propellants and heating (due to combustion, nuclear reactions, or electricity), the gas in the chamber and moving into the nozzle will form a mixture of chemical species at the state of chemical equilibrium. At this state, the forward and reverse rates of all chemical reactions are balanced, and the species remain at a fixed composition (as long as the temperature and/or pressure remain constant).

To continue, we need to define some terminology around combustion and chemical reactions.

The amounts of the chemical species at the equilibrium state are unknown, but can be determined based on the initial conditions using methods based on reaction equilibrium constants or the minimization of Gibbs free energy. If the temperature is not known/fixed, then it is also an unknown and must be found.

Fortunately, we can use software tools such as NASA’s CEA or Cantera to find the equilibrium state for us. We’ll focus on CEA here.

Combustion Terminology¶

Reactants and Products¶

Reactants = Fuel + Oxidizer

\text{Reactants} \longrightarrow \text{Products} + \text{Energy}

(2)

For a general hydrogen-oxygen reaction:

\text{H}_2 + \tfrac{1}{2}\text{O}_2 \longrightarrow n_1\text{H}_2\text{O} + n_2\text{H}_2 + n_3\text{O}_2 + n_4\text{H} + n_5\text{O} + n_6\text{OH} + n_7\text{H}_2\text{O}_2 + \cdots

(3)

Propellant Types¶

Fuels¶

Most rocket fuels are hydrocarbons (H-C compounds) or their derivatives:

Category	Formula Type	Examples
Hydrogen	H₂	LOX, GOX
Pure hydrocarbon	$\text{C}_x\text{H}_y$	Methane (CH₄), RP-1
Nitrogen-containing	$\text{C}_x\text{H}_y\text{O}_z\text{N}_i$	Hydrazine (N₂H₄), UDMH (H₂NN(CH₃)₂)
Alcohols	$\text{C}_x\text{H}_{y-1}\text{OH}$	Methanol (CH₃OH), Ethanol (C₂H₅OH)

Common liquid fuels:

H₂ — Liquid hydrogen (highest $I_{sp}$ , cryogenic)
RP-1 — Refined kerosene (storable, moderate $I_{sp}$ )
CH₄ — Liquid methane (cryogenic, gaining popularity)
C₂H₅OH — Ethanol
N₂H₄ — Hydrazine (storable, hypergolic with certain oxidizers)

Solid fuels:

Polybutadiene (PB) — Rubber-like binder
HTPB — Hydroxyl-terminated polybutadiene

Metal additives:

Al, Mg — Added to solid propellants for energy density

\text{Al} + \text{O}_2 \longrightarrow \text{Al}_2\text{O}_3 + \text{energy}

(4)

Oxidizers¶

Liquid oxidizers:

LOX (LO₂): Liquid oxygen (cryogenic, most common)
H₂O₂: Hydrogen peroxide
F₂: Fluorine (extremely toxic, rarely used)
N₂O: Nitrous oxide
N₂O₄ (NTO): Nitrogen tetroxide

Solid oxidizers:

AP: ammonium perchlorate (NH₄ClO₄)
AN: ammonium nitrate (NH₄NO₃)

Classification:

Cryogenic: LOX, LH₂, LCH₄
Storable: N₂O₄, N₂H₄, RP-1

“Storable” means that the propellant can be stored as a liquid under typical conditions, without requiring cooling like cryogenics.

Products¶

Fuels containing carbon and hydrogen (and oxygen, for alcohols) will form primarily carbon dioxide (CO₂) and water vapor (H₂O) as products in the exhaust:

\text{RP-1} + \text{O}_2 \longrightarrow \text{CO}_2 +\text{H}_2 \text{O}

(5)

If the propellants are hydrogen and oxygen, then the only product is water vapor:

\text{H}_2 + \text{O}_2 \longrightarrow \text{H}_2 \text{O}

(6)

These are referred to as the products of complete combustion. However, in reality, this is an ideal outcome and instead the equilibrium mixture will include a large number of other intermediate species (e.g., CO, H, OH).

Mole and Mass Fractions¶

Definitions¶

One mole contains Avogadro’s number of molecules/atoms:

1 \text{ mol} = 6.022 \times 10^{23} \text{ molecules/atoms}

(7)

Mole fraction of species $i$ :

X_i = \frac{n_i}{\sum n_i} = \frac{\text{moles of species } i}{\text{total moles}} \;,

(8)

where

\sum X_i = 1 \;.

(9)

Mass fraction of species $i$ :

Y_i = \frac{m_i}{\sum m_i} = \frac{\text{mass of species } i}{\text{total mass}} \;,

(10)

where

\sum Y_i = 1 \;.

(11)

Conversion Between Mole and Mass Fractions¶

Since molecular weight $\textrm{MW} = \text{mass}/\text{moles}$ :

Y_i = \frac{X_i \cdot \textrm{MW}_i}{\sum(X_i \cdot \textrm{MW}_i)}

(12)

Average molecular weight of the mixture:

\boxed{\overline{\textrm{MW}} = \sum X_i \cdot \textrm{MW}_i}

(13)

Stoichiometry¶

Stoichiometric mixture: The exact ratio of fuel to oxidizer needed for complete combustion, where all fuel and oxidizer are consumed and converted to fully oxidized products.

Example 1: Hydrogen-Oxygen¶

\text{H}_2 + \tfrac{1}{2}\text{O}_2 \longrightarrow \text{H}_2\text{O}

(14)

The stoichiometric molar ratio is 1 mol H₂ : 0.5 mol O₂.

Example 2: Methane-Oxygen¶

1.0 \, \text{CH}_4 + b \, \text{O}_2 \longrightarrow c \, \text{CO}_2 + d \, \text{H}_2\text{O}

(15)

To find the unknown stoichiometric coefficients $b$ , $c$ , and $d$ , we can apply conservation of elements:

Element	Balance	Result
C	$1 = c$	$c = 1$
H	$4 = 2d$	$d = 2$
O	$2b = 2c + d = 2 + 2 = 4$	$b = 2$

The stoichiometric reaction for methane and oxygen is then:

\text{CH}_4 + 2\text{O}_2 \longrightarrow \text{CO}_2 + 2\text{H}_2\text{O}

(16)

which can be achieved via an oxidizer-to-fuel molar flow rate of

\frac{\dot{n}_{\text{O}_2}}{\dot{n}_{\text{CH}_4}} = 2 \;.

(17)

Mole and Mass Fraction Example¶

Problem: For H₂ + ½O₂ → H₂O injected stoichiometrically, find the mole and mass fractions of the reactants.

Figure 2:Schematic of propellant injection showing H₂ and O₂ streams entering the combustion chamber at stoichiometric proportions.

Solution:

Stoichiometric reaction: $\text{H}_2 + \tfrac{1}{2}\text{O}_2 \rightarrow \text{H}_2\text{O}$

Mole fractions of reactants:

X_{\text{H}_2} = \frac{n_{\text{H}_2}}{n_{\text{H}_2} + n_{\text{O}_2}} = \frac{1}{1 + \tfrac{1}{2}} = 0.667

(18)

X_{\text{O}_2} = 0.333

(19)

Mass fractions:

Using $Y_i = \frac{X_i \cdot MW_i}{\sum(X_i \cdot MW_i)}$ :

Y_{\text{O}_2} = \frac{0.333 \times 32 \text{ kg/kmol}}{0.333(32) + 0.667(2)} = \frac{10.66}{10.66 + 1.33} = 0.89

(20)

Y_{\text{H}_2} = 0.11

(21)

Oxidizer-to-Fuel Ratio¶

The oxidizer-to-fuel ratio (O/F ratio or mixture ratio $r$ ) is defined as:

\boxed{r = \frac{\dot{m}_{ox}}{\dot{m}_F} = \frac{m_{ox}}{m_F}}

(22)

Example: H₂/O₂ Stoichiometric Ratio¶

For $\text{H}_2 + \tfrac{1}{2}\text{O}_2 \rightarrow \text{H}_2\text{O}$ :

Molar ratio:

\frac{\dot{n}_{ox}}{\dot{n}_F} = \frac{1/2}{1} = \frac{1}{2}

(23)

Mass ratio:

r = \frac{\dot{m}_{ox}}{\dot{m}_F} = \left(\frac{\dot{n}_{ox}}{\dot{n}_F}\right)\left(\frac{MW_{ox}}{MW_F}\right) = \left(\frac{1}{2}\right)\left(\frac{32 \text{ kg/kmol}}{2 \text{ kg/kmol}}\right) = 8.0

(24)

\boxed{r_{\text{stoich}} = (O/F)_{\text{stoich}} = 8.0 \quad \text{for H}_2/\text{O}_2}

(25)

Mass Flow Rate Relations¶

The total propellant mass flow rate is:

\dot{m} = \dot{m}_F + \dot{m}_{ox} = \dot{m}_F\left(1 + \frac{\dot{m}_{ox}}{\dot{m}_F}\right)

(26)

\boxed{\dot{m} = \dot{m}_F(1 + r)}

(27)

Solving for individual flow rates:

\boxed{\dot{m}_F = \frac{\dot{m}}{1 + r}} \qquad \boxed{\dot{m}_{ox} = \frac{\dot{m} \cdot r}{1 + r}}

(28)

Since thrust is a function of mass flow rate:

\text{Thrust} = f(\dot{m})

(29)

Combustion Chamber Modeling¶

Figure 3:Combustion chamber schematic showing fuel and oxidizer injection, residence time $t_{res}$ , chamber conditions ( $T_c$ , $P_c$ , $X_i$ ), and 1D isentropic expansion through the nozzle. The expansion can follow either frozen or shifting equilibrium assumptions.

Assumptions for Ideal Combustion Analysis¶

In the combustion chamber:

Adiabatic: $\dot{Q} \approx 0$ (no heat loss through walls)
Constant pressure: $P = P_c$
Chemical equilibrium: Fast reactions

t_{res} \gg t_{rxn}

(30)

where $t_{res}$ is the residence time in the chamber and $t_{rxn}$ is the characteristic reaction time.

In the nozzle:

1D flow
Isentropic expansion
Fully reacted gas

The nozzle expansion can be modeled as:

Frozen flow: Composition fixed at chamber values ( $X_i$ = constant)
Shifting equilibrium: Composition adjusts continuously to local $T$ , $P$

Chemical Equilibrium¶

Equilibrium Concept¶

At equilibrium, there is no net change in species concentrations — forward and reverse reaction rates are equal.

Figure 4:Graph showing approach to equilibrium: reactant concentration decreases while product concentration increases, both reaching steady-state values as $t \rightarrow \infty$ .

Example: H₂/O₂ Equilibrium¶

Initial state:

0.5 mol O₂
1.0 mol H₂
$T = 90$ K
$P = 20$ MPa

Final equilibrium state (at constant $P$ , $t \rightarrow \infty$ ):

$T = 3785$ K

Species	Mole Fraction
H₂O	0.733
H₂	0.11
H	0.025
O	0.012
O₂	0.031
OH	0.091

Thermodynamic Basis¶

For an adiabatic closed system ( $Q = 0$ , $W = 0$ ):

E_2 - E_1 = Q - W = 0

(31)

\Rightarrow U_2 = U_1

(32)

For an adiabatic open system (steady flow, no shaft work):

h_1 + \cancel{q} = h_2 + \cancel{W_s}

(33)

\Rightarrow h_1 = h_2

(34)

This is the basis for adiabatic flame temperature calculations.

Calculating Equilibrium Composition¶

Method 1: Analytical (Equilibrium Constants)¶

For the reaction:

\text{H}_2 + \tfrac{1}{2}\text{O}_2 \rightarrow n_{\text{H}_2\text{O}}\text{H}_2\text{O} + n_{\text{H}_2}\text{H}_2 + n_{\text{O}_2}\text{O}_2 + n_{\text{H}}\text{H} + n_{\text{O}}\text{O} + n_{\text{OH}}\text{OH} + \cdots

(35)

Governing equations:

First Law (energy conservation)
H atom conservation: Total H atoms in = Total H atoms out
O atom conservation: Total O atoms in = Total O atoms out

Equilibrium constants for dissociation reactions:

\tfrac{1}{2}\text{H}_2 + \tfrac{1}{2}\text{O}_2 \leftrightarrow \text{OH}

(36)

\text{H}_2 \leftrightarrow 2\text{H}

(37)

\text{O}_2 \leftrightarrow 2\text{O}

(38)

For the OH formation reaction:

\frac{n_{\text{OH}}^2}{n_{\text{H}_2} \cdot n_{\text{O}_2}} = K_p(T_c)

(39)

where $K_p$ is the equilibrium constant (function of temperature only).

Method 2: Software Tools¶

Modern approach uses software that minimizes Gibbs free energy:

NASA CEA (Chemical Equilibrium with Applications) — Free, widely used
Cantera — Open-source, Python interface

These tools solve the equilibrium problem numerically given:

Reactant species and amounts
Pressure and either temperature or enthalpy constraint

Thermochemical Calculation Types¶

Fixed temperature and pressure¶

Let’s first consider the problem where the pressure and temperature of the combustion/heating chamber are fixed and known, and determine the equilibrium composition of chemical species. This problem is relevant to an isothermal process, or where temperature is a design variable, such as in nuclear thermal or electrothermal rockets.

Example: Electrothermal Arcjet¶

Figure 5:Arcjet thruster schematic showing N₂H₄ propellant heated by electric arc to temperature $T_c$ , then expanded through nozzle. Electrical energy input raises gas temperature independent of chemical heat release.

For example, say we have an arcjet operating on gaseous hydrazine (N₂H₄) as a propellant, with a chamber temperature of 5000 K and pressure of 50 psia. The nozzle area ratio is 100, and the arcjet will operate in vacuum. For this system, determine the equilibrium composition, the average molecular weight, ratio of specific heats $\gamma$ , and then use these to get $c^*$ , $C_F$ , and $I_{\textrm{sp}}$ .

In CEA, this is a TP problem, or fixed temperature and pressure problem. We should expect that, at such high temperatures, the equilibrium state will have mostly one- and two-atom molecules, based on the elements present: N₂, H₂, H, N, and HN.

Decomposition reaction:

\text{N}_2\text{H}_4 \rightarrow a\,\text{N}_2 + b\,\text{H}_2 + c\,\text{H} + d\,\text{N} + e\,\text{NH} + \cdots \;.

(40)

The CEA plaintext input file looks like:

prob tp
 
p,psia= 50  t,k= 5000

reac
name N2H4 mol 1.0

output siunits
end

and the output is (with the repeated input removed):

*******************************************************************************

         NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA2, FEBRUARY 5, 2004
                   BY  BONNIE MCBRIDE AND SANFORD GORDON
      REFS: NASA RP-1311, PART I, 1994 AND NASA RP-1311, PART II, 1996

 *******************************************************************************

               THERMODYNAMIC EQUILIBRIUM PROPERTIES AT ASSIGNED

                           TEMPERATURE AND PRESSURE

             REACTANT                    WT FRACTION      ENERGY      TEMP
                                          (SEE NOTE)     KJ/KG-MOL      K  
 NAME        N2H4                         1.0000000         0.000      0.000

 O/F=    0.00000  %FUEL=  0.000000  R,EQ.RATIO= 0.000000  PHI,EQ.RATIO= 0.000000

 THERMODYNAMIC PROPERTIES

 P, BAR            3.4474
 T, K             5000.00
 RHO, KG/CU M    5.5368-2
 H, KJ/KG         42058.0
 U, KJ/KG         35831.8
 G, KJ/KG       -103744.4
 S, KJ/(KG)(K)    29.1605

 M, (1/n)           6.677
 (dLV/dLP)t      -1.04028
 (dLV/dLT)p        1.4750
 Cp, KJ/(KG)(K)   11.1350
 GAMMAs            1.2548
 SON VEL,M/SEC     2795.1

 MOLE FRACTIONS

 *H               0.74177
 *H2              0.04573
 *N               0.00806
 *NH              0.00021
 *N2              0.20422

So, CEA not only provides the equilibrium composition in terms of mole fraction ( $X_i$ ), but also the mean molecular weight of the mixture $MW$ ; thermodynamic properties and derivatives density $\rho$ , enthalpy $h$ , entropy $s$ , $\left(\partial \log V / \partial \log P\right)_T$ , $\left(\partial \log V / \partial \log T\right)_P$ , specific heat $C_p = \partial h / \partial T)_P$ , the ratio of specific heats ( $\gamma$ ), and the sonic velocity (i.e., speed of sound) $a$ . Some of these quantities are not particularly interesting to us right now, but we can use these quantities to find our quantities of interest.

area_ratio = 100
Tc = Q_(5000, 'K')
Pc = Q_(50, 'psi')

# output from CEA
MW = Q_(6.677, 'kg/kmol')
gamma = 1.2548

First, we need to find the exit pressure of the nozzle based on the area ratio:

# initial guesses for Pc / Pe
root = root_scalar(root_area_ratio, x0=1000, x1=2000, args=(area_ratio, gamma))
Pc_Pe = root.root
Pe = Pc / Pc_Pe
print(f"Exit pressure = {Pe.to('psi'): .2e~P}")

cstar = get_cstar(Tc, MW, gamma)
print(f"Cstar = {cstar.to('m/s'): .1f~P}")

CF0 = get_thrust_coeff(Pc_Pe, gamma)
# ambient pressure is zero in vacuum
CF = CF0 + (1/Pc_Pe) * area_ratio
print(f"C_F = {CF: .3f}")

Exit pressure =  2.54e-02 psi
Cstar =  3786.6 m/s
C_F =  1.884

Isp = CF * cstar / g0
print(f"Isp = {Isp.to('s'): .1f~P}")

Isp =  727.4 s

Problem Type 2: Constant H-P Problem (Adiabatic Combustion)¶

For chemical rockets, the temperature in the combustion chamber is unknown, and is a function of the propellant combination, oxidizer/fuel ratio, and chamber pressure.

Figure 6:Schematic for constant H-P problem: Propellants at known conditions enter chamber, combustion occurs at constant pressure, exit conditions ( $T_c$ , $X_i$ ) are determined by energy balance.

Example: Space Shuttle Main Engine (SSME)¶

Consider the Space Shuttle Main Engine (SSME, also known as the RS-25), which uses liquid oxygen and liquid hydrogen as propellants, at an O/F ratio of 6.0 and with a chamber pressure of 2870 psia. The engine nozzle has an area ratio of 68.8.

Find the chamber temperature, mole fractions of equilibrium products, specific heat ratio, average molecular weight, $c^*$ , $C_F$ , and specific impulse at sea level.

In CEA, this is an HP type problem (assigned enthalpy and pressure), where the enthalpy of the initial mixture is known based on the propellants and combustion is assumed to be adiabatic.

The CEA plaintext input file (with some comments removed) looks like:

prob hp
 
# Pressure (1 value):
p,psia= 2870
  
# Oxidizer/Fuel Wt. ratio (1 value):
o/f= 6.0
 
reac
fuel H2(L)             wt%=100.0000
oxid O2(L)             wt%=100.0000

output short
output siunits

end

and the output is (with the repeated input removed):

*******************************************************************************

         NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA2, FEBRUARY 5, 2004
                   BY  BONNIE MCBRIDE AND SANFORD GORDON
      REFS: NASA RP-1311, PART I, 1994 AND NASA RP-1311, PART II, 1996

 *******************************************************************************

         THERMODYNAMIC EQUILIBRIUM COMBUSTION PROPERTIES AT ASSIGNED

                                   PRESSURES

 CASE = _______________

             REACTANT                    WT FRACTION      ENERGY      TEMP
                                          (SEE NOTE)     KJ/KG-MOL      K  
 FUEL        H2(L)                        1.0000000     -9012.000     20.270
 OXIDANT     O2(L)                        1.0000000    -12979.000     90.170

 O/F=    6.00000  %FUEL= 14.285714  R,EQ.RATIO= 1.322780  PHI,EQ.RATIO= 1.322780

 THERMODYNAMIC PROPERTIES

 P, BAR            197.88
 T, K             3594.37
 RHO, KG/CU M    9.0105 0
 H, KJ/KG         -986.31
 U, KJ/KG        -3182.40
 G, KJ/KG        -62790.6
 S, KJ/(KG)(K)    17.1948

 M, (1/n)          13.608
 (dLV/dLP)t      -1.01921
 (dLV/dLT)p        1.3335
 Cp, KJ/(KG)(K)    7.3661
 GAMMAs            1.1472
 SON VEL,M/SEC     1587.2

 MOLE FRACTIONS

 *H               0.02575
 HO2              0.00003
 *H2              0.24744
 H2O              0.68555
 H2O2             0.00002
 *O               0.00207
 *OH              0.03694
 *O2              0.00220

  * THERMODYNAMIC PROPERTIES FITTED TO 20000.K

 NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL OXIDANTS

From the output, we see that the chamber temperature is 3594.37 K, the average molecular weight is 13.608 g/mol, and the specific heat ratio $\gamma$ is 1.1472. We also see that the equilibrium mixture contains mostly H $_2$ , H $_2$ O, OH, H, and with some trace amounts of O, O $_2$ , and H $_2$ O $_2$ . We can use this information to find the quantities of interest:

area_ratio = 68.8
Pc = Q_(2870, 'psi')
# sea level
Pa = Q_(1, 'atm')

# output from CEA
Tc = Q_(3594.37, 'K')
MW = Q_(13.608, 'kg/kmol')
gamma = 1.1472

# initial guesses for Pc / Pe
root = root_scalar(root_area_ratio, x0=500, x1=1000, args=(area_ratio, gamma))
Pc_Pe = root.root
Pe = Pc / Pc_Pe
print(f"Exit pressure = {Pe.to('psi'): .2f~P}")

cstar = get_cstar(Tc, MW, gamma)
print(f"Cstar = {cstar.to('m/s'): .1f~P}")

print('at sea level:')
CF0 = get_thrust_coeff(Pc_Pe, gamma)
CF = CF0 + (1/Pc_Pe + Pa/Pc) * area_ratio
print(f"C_F = {CF.to_base_units(): .3f~P}")

Isp = CF * cstar / g0
print(f"Isp = {Isp.to('s'): .1f~P}")

Exit pressure =  3.81 psi
Cstar =  2322.5 m/s
at sea level:
C_F =  2.350
Isp =  556.5 s

Problem Type 3: Rocket Problem¶

CEA can also calculate performance quantities specific to rockets, such as the effective velocity (C-star, $c^*$ ), thrust coefficient ( $C_F$ ), and specific impulse ( $I_{\text{sp}}$ ).

Figure 7:Complete rocket analysis showing stations: chamber (c), throat (t), and exit (e). CEA calculates properties at each station and performance parameters. Note: Rocket problem type assumes optimal expansion ( $P_e = P_a$ ) for vacuum $I_{sp}$ ; actual $C_F$ requires specifying $P_a$ .

For the above example, but choosing the rocket problem and specifying a (supersonic) nozzle area ratio of 68.8, and leaving other options at their defaults, CEA provides this output:

 
 *******************************************************************************

         NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA2, FEBRUARY 5, 2004
                   BY  BONNIE MCBRIDE AND SANFORD GORDON
      REFS: NASA RP-1311, PART I, 1994 AND NASA RP-1311, PART II, 1996

 *******************************************************************************



  
 ### CEA analysis performed on Tue 03-Mar-2026 13:16:02
  
 # Problem Type: "Rocket" (Infinite Area Combustor)
  
 prob case = _______________1929 ro equilibrium
  
 # Pressure (1 value):
 p,psia= 2870
 # Supersonic Area Ratio (1 value):
 supar= 68.8
  
 # Oxidizer/Fuel Wt. ratio (1 value):
 o/f = 6.0
  
 # You selected the following fuels and oxidizers:
 reac
 fuel H2(L)             wt%=100.0000
 oxid O2(L)             wt%=100.0000
  
 # You selected these options for output:
 # short version of output
 output short
 # Proportions of any products will be expressed as Mass Fractions.
 output massf
 # Heat will be expressed as siunits
 output siunits
  
 # Input prepared by this script:/var/www/sites/cearun/cgi-bin/CEARUN/prepareInpu
 tFile.cgi
  
 ### IMPORTANT:  The following line is the end of your CEA input file!
 end





              THEORETICAL ROCKET PERFORMANCE ASSUMING EQUILIBRIUM

           COMPOSITION DURING EXPANSION FROM INFINITE AREA COMBUSTOR

 Pin =  2870.0 PSIA
 CASE = _______________

             REACTANT                    WT FRACTION      ENERGY      TEMP
                                          (SEE NOTE)     KJ/KG-MOL      K  
 FUEL        H2(L)                        1.0000000     -9012.000     20.270
 OXIDANT     O2(L)                        1.0000000    -12979.000     90.170

 O/F=    6.00000  %FUEL= 14.285714  R,EQ.RATIO= 1.322780  PHI,EQ.RATIO= 1.322780

                 CHAMBER   THROAT     EXIT
 Pinf/P            1.0000   1.7401   960.27
 P, BAR            197.88   113.71  0.20607
 T, K             3594.37  3378.30  1234.47
 RHO, KG/CU M    9.0105 0 5.5605 0 2.8331-2
 H, KJ/KG         -986.31 -2160.57 -10543.0
 U, KJ/KG        -3182.40 -4205.62 -11270.4
 G, KJ/KG        -62790.6 -60249.7 -31769.4
 S, KJ/(KG)(K)    17.1948  17.1948  17.1948

 M, (1/n)          13.608   13.735   14.111
 (dLV/dLP)t      -1.01921 -1.01432 -1.00000
 (dLV/dLT)p        1.3335   1.2644   1.0000
 Cp, KJ/(KG)(K)    7.3661   6.7421   2.9102
 GAMMAs            1.1472   1.1484   1.2539
 SON VEL,M/SEC     1587.2   1532.5    955.0
 MACH NUMBER        0.000    1.000    4.578

 PERFORMANCE PARAMETERS

 Ae/At                      1.0000   68.800
 CSTAR, M/SEC               2322.1   2322.1
 CF                         0.6599   1.8827
 Ivac, M/SEC                2867.0   4538.3
 Isp, M/SEC                 1532.5   4371.9


 MASS FRACTIONS

 *H               0.00191  0.00151  0.00000
 HO2              0.00008  0.00004  0.00000
 *H2              0.03665  0.03595  0.03486
 H2O              0.90756  0.92389  0.96514
 H2O2             0.00004  0.00002  0.00000
 *O               0.00243  0.00146  0.00000
 *OH              0.04616  0.03386  0.00000
 *O2              0.00516  0.00326  0.00000

  * THERMODYNAMIC PROPERTIES FITTED TO 20000.K

 NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL OXIDANTS

The key properties include:

C-star of 2322.1 m/s (based on combustion chamber conditions)
throat pressure of 113.71 bar and temperature of 3378.30 K
exit pressure of 0.20607 bar, temperature of 1234.47 K
at the nozzle exit, thrust coefficient = 1.8827, Isp = 4371.9 m/s

The reported specific impulse $I_{\textrm{sp}}$ assumes that the ambient pressure is the same as the exit pressure (i.e., optimally expanded nozzle), since CEA does not know the actual ambient pressure. (It also reports the specific impulse in vacuum, $I_{\textrm{vac}}$ .) The value is given in m/s and needs to be divided by the gravity constant to get our usual units of seconds:

Isp_cea = Q_(4371.9, 'm/s')
Isp = Isp_cea / g0
print(f"Isp = {Isp.to('s'): .1f~P}")

Isp =  445.8 s

The next thing to note is that CEA defaulted to assuming shifting equilibrium mixture composition from the combustion chamber and throughout the nozzle; in other words, it is assuming that at every location the equilibrium composition adjusts immediately based on the local temperature and pressure in the nozzle. This is based on the assumption that chemical reactions occur infinitely fast, relative to the timescale of the flow velocity.

An alternative assumption is frozen flow, which is the other extreme: the composition remains the same throughout the nozzle, based on the composition in the combustion chamber. This assumption essentially states that the timescales due to high flow velocity in the nozzle are shorter than the chemical timescales of the reactions.

Neither of these assumptions are truly correct; in the real system, changes in composition will be driven by finite-rate chemical reactions, but it is extremely challening to model these time-dependent effects.

CEA can use the frozen flow approximation by selecting the Frozen rocket problem option. This produces:

           THEORETICAL ROCKET PERFORMANCE ASSUMING FROZEN COMPOSITION

 Pin =  2870.0 PSIA
 CASE = _______________

             REACTANT                    WT FRACTION      ENERGY      TEMP
                                          (SEE NOTE)     KJ/KG-MOL      K  
 FUEL        H2(L)                        1.0000000     -9012.000     20.270
 OXIDANT     O2(L)                        1.0000000    -12979.000     90.170

 O/F=    6.00000  %FUEL= 14.285714  R,EQ.RATIO= 1.322780  PHI,EQ.RATIO= 1.322780

                 CHAMBER   THROAT     EXIT
 Pinf/P            1.0000   1.7688  1110.26
 P, BAR            197.88   111.87  0.17823
 T, K             3594.37  3276.93   980.79
 RHO, KG/CU M    9.0105 0 5.5875 0 2.9742-2
 H, KJ/KG         -986.31 -2182.76 -9907.46
 U, KJ/KG        -3182.40 -4184.90 -10506.7
 G, KJ/KG        -62790.6 -58528.7 -26771.8
 S, KJ/(KG)(K)    17.1948  17.1948  17.1948

 M, (1/n)          13.608   13.608   13.608
 Cp, KJ/(KG)(K)    3.7951   3.7416   2.7469
 GAMMAs            1.1919   1.1952   1.2861
 SON VEL,M/SEC     1617.9   1546.9    877.9
 MACH NUMBER        0.000    1.000    4.812

 PERFORMANCE PARAMETERS

 Ae/At                      1.0000   68.800
 CSTAR, M/SEC               2289.4   2289.4
 CF                         0.6757   1.8450
 Ivac, M/SEC                2841.2   4365.9
 Isp, M/SEC                 1546.9   4224.0

 MASS FRACTIONS

 *H              0.00191   HO2             0.00008   *H2             0.03665
 H2O             0.90756   H2O2            0.00004   *O              0.00243
 *OH             0.04616   *O2             0.00516

  * THERMODYNAMIC PROPERTIES FITTED TO 20000.K

 NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL OXIDANTS

With the key properties:

C-star of 2289.4 m/s (based on combustion chamber conditions)
exit pressure of 0.17823 bar
at the nozzle exit, thrust coefficient = 1.8450, Isp = 4224.0 m/s

Isp_cea = Q_(4224.0, 'm/s')
Isp_eq = Isp_cea / g0
print(f"Isp (equilibrium) = {Isp_eq.to('s'): .1f~P}")

print(f"Difference: {100*np.abs(Isp_eq - Isp)/Isp: .2f~P}%")

Isp (equilibrium) =  430.7 s
Difference:  3.38%

The shifting equilibrium assumption typically overpredicts the real performance by 1–4%, while the frozen composition assumption underpredicts performance by 1–4%.

Varying oxidizer-to-fuel ratio¶

The thermochemical performance of a rocket, represented by $c^*$ , is a strong function of the oxidizer-to-fuel ratio, or $r =$ O/F, and thus the overall performance ( $I_{\textrm{sp}}$ ) also depends on the O/F ratio. For liquid rockets, this ratio can be controlled directly since the propellants are injected separately. CEA can help us identify the optimal O/F ratio.

For example, for the SSME considered above using liquid hydrogen and liquid oxygen, we can have CEA sweep over a range of O/F values, from 2 to 8 in increments of 0.25. This gives an input file of

 prob case=_______________8942 ro equilibrium
  
 # Pressure (1 value):
 p,psia= 2870
 # Supersonic Area Ratio (1 value):
 supar= 68.8
  
 # Oxidizer/Fuel Wt. ratio (25 values):
 o/f= 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, 5.25, 5.5, 5
 .75, 6, 6.25, 6.5, 6.75, 7, 7.25, 7.5, 7.75, 8
  
 reac
 fuel H2(L)             wt%=100.0000
 oxid O2(L)             wt%=100.0000
  
 output short
 # Proportions of any products will be expressed as Mole Fractions.
 # Heat will be expressed as siunits
 output siunits
  
 end

The output of CEA will be quite long, since there are 25 separate cases; the output is in the file cea_OF_ratio_out.txt. We can use Python to extract the desired information:

# Click the + to see the full contents of the output file
filename = 'cea_OF_ratio_out.txt'
print('Contents of ' + filename + ':')
with open(filename) as f:
    file_contents = f.read()
    print(file_contents)

mixture_ratios = np.arange(2, 8.1, 0.25)

with open('cea_OF_ratio_out.txt', 'r') as f:
    lines = f.readlines()

cstars = []
specific_impulses = []
specific_impulses_vac = []

for line in lines:
    # ignore blank lines
    if not line.strip():
        continue
    
    words = line.split()
    if words[0] == 'CSTAR,':
        cstars.append(float(words[3]))
    elif words[0] == 'Isp,':
        specific_impulses.append(float(words[3]))
    elif words[0] == 'Ivac,':
        specific_impulses_vac.append(float(words[3]))
cstars = Q_(np.array(cstars), 'm/s')
specific_impulses = Q_(np.array(specific_impulses), 'm/s') / g0
specific_impulses_vac = Q_(np.array(specific_impulses_vac), 'm/s') / g0

# just checking that things line up
assert len(mixture_ratios) == len(cstars)
assert len(mixture_ratios) == len(specific_impulses)
assert len(mixture_ratios) == len(specific_impulses_vac)

fig, axes = plt.subplots(3, 1)

axes[0].plot(mixture_ratios, to_si(cstars))
axes[0].set_ylabel('$c^*$ (m/s)')
axes[0].grid(True)

axes[1].plot(mixture_ratios, to_si(specific_impulses))
axes[1].set_ylabel('$I_{\\mathrm{sp}}$ (s)')
axes[1].grid(True)

axes[2].plot(mixture_ratios, to_si(specific_impulses_vac))
axes[2].set_xlabel('O/F ratio')
axes[2].set_ylabel('$I_{\\mathrm{vac}}$ (s)')
axes[2].grid(True)

plt.tight_layout()
plt.show()

Loading...

This allows us to determine the mixture ratio for optimal performance, although in practice something slightly off-optimal might be used for other reasons (e.g., combustion stability, reducing peak temperature).

Frozen vs. shifting equilibrium¶

We have already seen how assuming the mixture composition either remains frozen from the combustion chamber or follows a shifting equilibrium model leads to different calculations of performance. Frozen equilibrium underpredicts, while shifting equilibrium overpredicts. We can use these together to estimate the true actual performance as somewhere in between, such as 40% of the difference between the two results.

We already have results for shifting equilibrum, and can similarly extract the specific impulse for frozen equilibrium calculations saved in CEA_OF_ratio_frozen_out.txt:

# Click the + to see the full contents of the output file
filename = 'CEA_OF_ratio_frozen_out.txt'
print('Contents of ' + filename + ':')
with open(filename) as f:
    file_contents = f.read()
    print(file_contents)

with open('CEA_OF_ratio_frozen_out.txt', 'r') as f:
    lines = f.readlines()

specific_impulses_frozen = []

for line in lines:
    # ignore blank lines
    if not line.strip():
        continue
    
    words = line.split()
    if words[0] == 'Isp,':
        specific_impulses_frozen.append(float(words[3]))
specific_impulses_frozen = Q_(np.array(specific_impulses_frozen), 'm/s') / g0

# just checking that things line up
assert len(mixture_ratios) == len(specific_impulses_frozen)

specific_impulses_act = specific_impulses_frozen + 0.4*(specific_impulses - specific_impulses_frozen)

plt.plot(mixture_ratios, to_si(specific_impulses), label='Shifting equilibrium')
plt.plot(mixture_ratios, to_si(specific_impulses_frozen), label='Frozen')
plt.plot(mixture_ratios, to_si(specific_impulses_act), '--', label='Actual')
plt.ylabel('$I_{\\mathrm{sp}}$ (s)')
plt.xlabel('Mixture ratio')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

Loading...