Implementing CEA calculations using Cantera
Contents
Implementing CEA calculations using Cantera¶
# this line makes figures interactive in Jupyter notebooks
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
import cantera as ct
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
# these lines are only for helping improve the display
import matplotlib_inline.backend_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('pdf', 'png')
plt.rcParams['figure.dpi']= 300
plt.rcParams['savefig.dpi'] = 300
plt.rcParams['mathtext.fontset'] = 'cm'
CEA (Chemical Equilibrium with Applications) is a classic NASA software tool developed for analyzing combustion and rocket propulsion problems. It was written in Fortran, but is available to run via a web interface.
Given rocket propellants, CEA can not only determine the combustion chamber equilibrium composition and temperature, but also calculate important rocket performance parameters.
Although CEA is extremely useful, it cannot (easily) be used within Python. Plus, we might want to Cantera is a modern software library for solving problems in chemical kinetics, thermodynamics, and transport, that offers a Python interface. Cantera natively supports phase and chemical equilibrium solvers. In particular, it can simulate finite-rate chemical reactions.
This article examines how we can use Cantera and Python to perform the calculations of CEA.
Fixed temperature and pressure¶
Given a fixed temperature and pressure, 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.
For example, say we have gaseous hydrazine (N2H4) as a propellant, with a chamber temperature of 5000 K and pressure of 50 psia. For this system, determine the equilibrium composition.
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: N2, H2, H, N, and HN.
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\).
We can perform the same equilibrium calculation in Cantera, but we need to construct an object that contains the appropriate chemical species. Cantera actually comes with a NASA database of gaseous species thermodynamic models,
in the nasa_gas.cti file.
# extract all species in the NASA database
full_species = {S.name: S for S in ct.Species.listFromFile('nasa_gas.cti')}
# extract only the relevant species
species = [full_species[S] for S in (
'N2H4', 'N2', 'H2', 'H', 'N', 'NH'
)]
gas = ct.Solution(thermo='IdealGas', species=species)
temperature = Q_(5000, 'K')
pressure = Q_(50, 'psi')
gas.TPX = to_si(temperature), to_si(pressure), 'N2H4:1.0'
gas.equilibrate('TP')
gas()
temperature 5000 K
pressure 3.4474e+05 Pa
density 0.055346 kg/m^3
mean mol. weight 6.6743 kg/kmol
phase of matter gas
1 kg 1 kmol
--------------- ---------------
enthalpy 4.2088e+07 2.8091e+08 J
internal energy 3.586e+07 2.3934e+08 J
entropy 29182 1.9477e+05 J/K
Gibbs function -1.0382e+08 -6.9294e+08 J
heat capacity c_p 3779.4 25225 J/K
heat capacity c_v 2533.6 16910 J/K
mass frac. Y mole frac. X chem. pot. / RT
--------------- --------------- ---------------
N2 0.8567 0.20411 -30.731
H2 0.013688 0.045315 -24.65
H 0.1121 0.74225 -12.325
N 0.017042 0.0081204 -15.365
NH 0.00046864 0.00020831 -27.691
[ +1 minor] 2.3328e-15 4.8586e-16
Comparing the results from CEA and Cantera, we see very good agreement between (most) thermodynamic properties and the species mole fractions.
But, the heat capacity \(C_p\) appears very different, and if we calculate the specific heat ratio,
we will see it also differs quite substantially:
gamma_ct = gas.cp_mole / gas.cv_mole
print(f'Cantera specific heat ratio: {gamma_ct: .4f}')
Cantera specific heat ratio: 1.4917
🤯
Well, 1.492 is quite different from 1.255, and this would lead to substantially different rocket performance parameters that depend on \(\gamma\).
So, what’s going on?
Well, the key lies in examining the actual definition of specific heat, following Gordon and McBride [GM94]:
In this derivative, while enthalpy and temperature change, and pressure is held constant, what happens to the species composition? We could assume the composition is “frozen” and remains fixed, or that the composition adjusts to a new equilibrium instantaneously. The “equilibrium” specific heat then has two components, a frozen contribution and reaction contribution:
where \(N_s\) is the number of species and \(N_g\) is the number of gas-phase species (so that \(N_g + 1\) refers to the first condensed-phase species, if present).
But, Cantera defines quantities like specific heat (and other thermodynamic quantities based on derivatives) at fixed composition, meaning Cantera’s specific heat is just the frozen contribution \(C_{p,f}\).
We can obtain the full equilibrium-based value of specific heat, but it requires determining additional thermodynamic derivatives. Following Gordon and McBride [GM94] again, we can obtain this system of linear equations:
where \({N_e}\) is the number of elements.
As a first pass, let’s assume that no condensed species are present. (This is fine for conditions in the combustion chamber, but for some systems the rapid expansion in the nozzle may drop below the dew point for some species.)
Then, the unknowns in that system of equations are \( \left( \frac{\partial \pi_i}{\partial \log T}\right)_P\), \( \left( \frac{\partial \log n}{\partial \log T}\right)_P\), \( \left( \frac{\partial \pi_i}{\partial \log P}\right)_T\), \( \left( \frac{\partial \log n}{\partial \log P}\right)_T\), with a total of \(2 \times N_e + 2\) unknowns.
For the current system of N2H4, \(N_e = 2\) and thus there are six unknowns.
Since this is a linear system of equations, we can solve it using linear algebra, via NumPy’s
linalg.solve function. Let’s set up a function to solve this system:
def get_thermo_derivatives(gas):
'''Gets thermo derivatives based on shifting equilibrium.
'''
# unknowns for system with no condensed species:
# dpi_i_dlogT_P (# elements)
# dlogn_dlogT_P
# dpi_i_dlogP_T (# elements)
# dlogn_dlogP_T
# total unknowns: 2*n_elements + 2
num_var = 2 * gas.n_elements + 2
coeff_matrix = np.zeros((num_var, num_var))
right_hand_side = np.zeros(num_var)
tot_moles = 1.0 / gas.mean_molecular_weight
moles = gas.X * tot_moles
condensed = False
# indices
idx_dpi_dlogT_P = 0
idx_dlogn_dlogT_P = idx_dpi_dlogT_P + gas.n_elements
idx_dpi_dlogP_T = idx_dlogn_dlogT_P + 1
idx_dlogn_dlogP_T = idx_dpi_dlogP_T + gas.n_elements
# construct matrix of elemental stoichiometric coefficients
stoich_coeffs = np.zeros((gas.n_elements, gas.n_species))
for i, elem in enumerate(gas.element_names):
for j, sp in enumerate(gas.species_names):
stoich_coeffs[i,j] = gas.n_atoms(sp, elem)
# equations for derivatives with respect to temperature
# first n_elements equations
for k in range(gas.n_elements):
for i in range(gas.n_elements):
coeff_matrix[k,i] = np.sum(stoich_coeffs[k,:] * stoich_coeffs[i,:] * moles)
coeff_matrix[k, gas.n_elements] = np.sum(stoich_coeffs[k,:] * moles)
right_hand_side[k] = -np.sum(stoich_coeffs[k,:] * moles * gas.standard_enthalpies_RT)
# skip equation relevant to condensed species
for i in range(gas.n_elements):
coeff_matrix[gas.n_elements, i] = np.sum(stoich_coeffs[i, :] * moles)
right_hand_side[gas.n_elements] = -np.sum(moles * gas.standard_enthalpies_RT)
# equations for derivatives with respect to pressure
for k in range(gas.n_elements):
for i in range(gas.n_elements):
coeff_matrix[gas.n_elements+1+k,gas.n_elements+1+i] = np.sum(stoich_coeffs[k,:] * stoich_coeffs[i,:] * moles)
coeff_matrix[gas.n_elements+1+k, 2*gas.n_elements+1] = np.sum(stoich_coeffs[k,:] * moles)
right_hand_side[gas.n_elements+1+k] = np.sum(stoich_coeffs[k,:] * moles)
for i in range(gas.n_elements):
coeff_matrix[2*gas.n_elements+1, gas.n_elements+1+i] = np.sum(stoich_coeffs[i, :] * moles)
right_hand_side[2*gas.n_elements+1] = np.sum(moles)
derivs = np.linalg.solve(coeff_matrix, right_hand_side)
dpi_dlogT_P = derivs[idx_dpi_dlogT_P : idx_dpi_dlogT_P + gas.n_elements]
dlogn_dlogT_P = derivs[idx_dlogn_dlogT_P]
dpi_dlogP_T = derivs[idx_dpi_dlogP_T]
dlogn_dlogP_T = derivs[idx_dlogn_dlogP_T]
# dpi_dlogP_T is not used
return dpi_dlogT_P, dlogn_dlogT_P, dlogn_dlogP_T
Using these derivatives, we can then calculate the specific heat, other relevant derivatives, and the ratio of specific heats:
The ratio of specific heats shows up via the speed of sound:
where the ratio of specific heats is
and
The constant volume specific heat is
def get_thermo_properties(gas, dpi_dlogT_P, dlogn_dlogT_P, dlogn_dlogP_T):
'''Calculates specific heats, volume derivatives, and specific heat ratio.
Based on shifting equilibrium for mixtures.
'''
tot_moles = 1.0 / gas.mean_molecular_weight
moles = gas.X * tot_moles
# construct matrix of elemental stoichiometric coefficients
stoich_coeffs = np.zeros((gas.n_elements, gas.n_species))
for i, elem in enumerate(gas.element_names):
for j, sp in enumerate(gas.species_names):
stoich_coeffs[i,j] = gas.n_atoms(sp, elem)
spec_heat_p = ct.gas_constant * (
np.sum([dpi_dlogT_P[i] *
np.sum(stoich_coeffs[i,:] * moles * gas.standard_enthalpies_RT)
for i in range(gas.n_elements)
]) +
np.sum(moles * gas.standard_enthalpies_RT) * dlogn_dlogT_P +
np.sum(moles * gas.standard_cp_R) +
np.sum(moles * gas.standard_enthalpies_RT**2)
)
dlogV_dlogT_P = 1 + dlogn_dlogT_P
dlogV_dlogP_T = -1 + dlogn_dlogP_T
spec_heat_v = (
spec_heat_p + gas.P * gas.v / gas.T * dlogV_dlogT_P**2 / dlogV_dlogP_T
)
gamma = spec_heat_p / spec_heat_v
gamma_s = -gamma/dlogV_dlogP_T
return dlogV_dlogT_P, dlogV_dlogP_T, spec_heat_p, gamma_s
derivs = get_thermo_derivatives(gas)
dlogV_dlogT_P, dlogV_dlogP_T, cp, gamma_s = get_thermo_properties(
gas, derivs[0], derivs[1], derivs[2]
)
print(f'Cp = {cp: .2f} J/(K kg)')
print(f'(d log V/d log P)_T = {dlogV_dlogP_T: .4f}')
print(f'(d log V/d log T)_P = {dlogV_dlogT_P: .4f}')
print(f'gamma_s = {gamma_s: .4f}')
speed_sound = np.sqrt(ct.gas_constant * gas.T * gamma_s / gas.mean_molecular_weight)
print(f'Speed of sound = {speed_sound: .1f} m/s')
Cp = 11104.47 J/(K kg)
(d log V/d log P)_T = -1.0400
(d log V/d log T)_P = 1.4722
gamma_s = 1.2549
Speed of sound = 2795.8 m/s
🎉 Success! These calculations agree very closely with those from CEA.
Adiabatic combustion¶
CEA also supports calculating the chamber temperature (along with composition) for adiabatic combustion, both with gaseous and liquid propellants.
Cantera’s equilibrium solver that we used above handles constant enthlapy and pressure equilibrium (HP) just fine with gaseous reactants, but how to
CEA has a database of reactants with assigned enthalpies, as described by Gordon and McBride [GM94]:
noncryogenic reactants are represented via enthalpy of formation (i.e., heat of formation) at the standard reference temperature of 298.15 K
cryogenic liquid reactants are represented via enthalpies given at their boiling points, which represent the standard enthalpy of formation minus the sensible heat (between 298.15 K and the boiling point), the heat of vaporization at the boiling point, and also the difference in enthalpy due to real gas effects at the boiling point.
For example, CEA’s thermodynamic database [MZG02] represents liquid dinitrogen tetroxide (N2O4), which is an oxidizer used with hydrazine, with
N2O4(L) Dinitrogen tetroxide. McBride,1996 pp85,93.
0 g 6/96 N 2.00O 4.00 0.00 0.00 0.00 1 92.0110000 -17549.000
298.150 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000
while cryogenic liquid hydrogen is given with
H2(L) Hydrogen. McBride,1996 pp84,92.
0 g 6/96 H 2.00 0.00 0.00 0.00 0.00 1 2.0158800 -9012.000
20.270 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000
The full format of species thermodynamic entries is given in the CEA User Manual [GM96], but for these reactants the key information includes
species name, given in the first line
elemental composition, given in a fixed column format in the second line
phase, given as an integer in the third-to-last entry of the second line (zero for gases, nonzero for condensed phases)
molecular weight, in the second-to-last entry of the second line
enthalpy at the boiling point, in J/mol, at the end of the second line
boiling point temperature, in K, at the beginning of the third line
In general, both CEA and Cantera represent the thermodynamic properties of gaseous and condensed species via more-sophisticated polynomial fits across multiple ranges of temperatures, but these problems only require initial enthalpy of reactants.
Let’s consider the Space Shuttle main engine (SSME), which used cryogenic liquid hydrogen and liquid oxygen at an oxidizer to fuel ratio of 6.0 and a chamber pressure of around 3000 psia.
This is a constant enthalpy and pressure problem (hp) in CEA:
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 Wed 27-Jan-2021 13:09:27
# Problem Type: "Assigned Enthalpy and Pressure"
prob case=_______________3446 hp
# Pressure (1 value):
p,psia= 3000
# 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 Mole Fractions.
# Heat will be expressed as siunits
output siunits
# Input prepared by this script:prepareInputFile.cgi
### IMPORTANT: The following line is the end of your CEA input file!
end
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 206.84
T, K 3598.76
RHO, KG/CU M 9.4113 0
H, KJ/KG -986.31
U, KJ/KG -3184.12
G, KJ/KG -62768.7
S, KJ/(KG)(K) 17.1677
M, (1/n) 13.614
(dLV/dLP)t -1.01897
(dLV/dLT)p 1.3291
Cp, KJ/(KG)(K) 7.3140
GAMMAs 1.1475
SON VEL,M/SEC 1588.1
MOLE FRACTIONS
*H 0.02543
HO2 0.00003
*H2 0.24740
H2O 0.68635
H2O2 0.00002
*O 0.00202
*OH 0.03659
*O2 0.00215
* THERMODYNAMIC PROPERTIES FITTED TO 20000.K
NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL OXIDANTS
The key results include the chamber pressure \(T_c\) of 3598.8 K, the specific heat ratio \(\gamma_s\) of 1.148, and the mean molecular weight of 13.614 kg/kmol.
To perform this calculation using Cantera, we need the reactant information:
H2(L) Hydrogen. McBride,1996 pp84,92.
0 g 6/96 H 2.00 0.00 0.00 0.00 0.00 1 2.0158800 -9012.000
20.270 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000
O2(L) Oxygen. McBride,1996 pp85,93.
0 g 6/96 O 2.00 0.00 0.00 0.00 0.00 1 31.9988000 -12979.000
90.170 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000
Cantera has fairly sophisticated ways of representing the thermodynamics of condensed phases, but in this case we actually do not need that—we just need a way of easily representing the elemental composition and enthalpy of the reactants, which is the data needed for constraining the equilibrium solver.
So, we can actually use the ideal gas thermodynamic model (ideal-gas) for the phase.
For each species, we can use the constant heat capacity (constant-cp) thermodynamic model,
with the reference temperature set to boiling point (for the cryogenic liquid propellants in this case; for non-cryogenic reactants, this would be 298.15 K), the reference enthalpy set to the assigned value, and the reference specific heat and entropy set to zero.
I’ve constructed a representative Cantera YAML input file, that describes separate phases for liquid hydrogen and liquid oxygen.
⚠️ Warning ⚠️ these phases are only valid at the specific cryogenic temperature specified, and should only be used for this specific purpose (as reactants).
h2o2_filename = 'h2o2_react.yaml'
print('Contents of ' + h2o2_filename + ':\n')
with open(h2o2_filename) as f:
file_contents = f.read()
print(file_contents)
Contents of h2o2_react.yaml:
phases:
- name: liquid_hydrogen
thermo: ideal-gas
elements: [H]
species: [H2(L)]
- name: liquid_oxygen
thermo: ideal-gas
elements: [O]
species: [O2(L)]
species:
- name: H2(L)
composition: {H: 2}
thermo:
model: constant-cp
T0: 20.270
h0: -9012.0 J/mol
s0: 0.0
cp0: 0.0
- name: O2(L)
composition: {O: 2}
thermo:
model: constant-cp
T0: 90.170
h0: -12979.0 J/mol
s0: 0.0
cp0: 0.0
To set up the system in Cantera, we create separate Solution objects for the liquid hydrogen and oxygen phases, and also a Solution containing the gas-phase products (actually, this could also include condensed species as well!).
Then, create a Mixture that contains all three objects, and specify the initial moles of hydrogen and oxygen based on the oxidizer-to-fuel ratio:
o_f_ratio = 6.0
temperature_h2 = Q_(20.270, 'K')
temperature_o2 = Q_(90.170, 'K')
pressure_chamber = Q_(3000, 'psi')
h2 = ct.Solution(h2o2_filename, 'liquid_hydrogen')
h2.TP = to_si(temperature_h2), to_si(pressure_chamber)
o2 = ct.Solution(h2o2_filename, 'liquid_oxygen')
o2.TP = to_si(temperature_o2), to_si(pressure_chamber)
molar_ratio = o_f_ratio / (o2.mean_molecular_weight / h2.mean_molecular_weight)
moles_ox = molar_ratio / (1 + molar_ratio)
moles_f = 1 - moles_ox
gas2 = ct.Solution('nasa_h2o2.yaml', 'gas')
# create a mixture of the liquid phases with the gas-phase model,
# with the number of moles for fuel and oxidizer based on
# the O/F ratio
mix = ct.Mixture([(h2, moles_f), (o2, moles_ox), (gas2, 0)])
# Solve for the equilibrium state, at constant enthalpy and pressure
mix.equilibrate('HP')
gas2()
derivs = get_thermo_derivatives(gas2)
dlogV_dlogT_P, dlogV_dlogP_T, cp, gamma_s = get_thermo_properties(
gas2, derivs[0], derivs[1], derivs[2]
)
print(f'Cp = {cp: .2f} J/(K kg)')
print(f'(d log V/d log P)_T = {dlogV_dlogP_T: .4f}')
print(f'(d log V/d log T)_P = {dlogV_dlogT_P: .4f}')
print(f'gamma_s = {gamma_s: .4f}')
speed_sound = np.sqrt(ct.gas_constant * gas2.T * gamma_s / gas2.mean_molecular_weight)
print(f'Speed of sound = {speed_sound: .1f} m/s')
gas:
temperature 3597.5 K
pressure 2.0684e+07 Pa
density 9.4137 kg/m^3
mean mol. weight 13.613 kg/kmol
phase of matter gas
1 kg 1 kmol
--------------- ---------------
enthalpy -9.8628e+05 -1.3426e+07 J
internal energy -3.1835e+06 -4.3338e+07 J
entropy 17175 2.3381e+05 J/K
Gibbs function -6.2775e+07 -8.5457e+08 J
heat capacity c_p 3795.4 51668 J/K
heat capacity c_v 3184.7 43354 J/K
mass frac. Y mole frac. X chem. pot. / RT
--------------- --------------- ---------------
H 0.0018901 0.025526 -8.7917
HO2 8.3393e-05 3.4395e-05 -40.623
H2 0.036632 0.24736 -17.583
H2O 0.908 0.68615 -33.499
H2O2 4.23e-05 1.693e-05 -49.415
O 0.0023965 0.0020392 -15.916
OH 0.045862 0.036711 -24.707
O2 0.0050891 0.0021651 -31.831
O3 2.257e-08 6.4014e-09 -47.747
Cp = 7330.18 J/(K kg)
(d log V/d log P)_T = -1.0190
(d log V/d log T)_P = 1.3305
gamma_s = 1.1474
Speed of sound = 1587.8 m/s
🎉 Success! We get an equilibrium temperature of 3597.5 K, which is just 0.036% off the value calculated by CEA. Similarly, the ratios of specific heats match within 0.009%, and the speed of sounds within 0.019%.
Rocket calculations¶
CEA also calculates performance quantities specific to rockets, such as the effective velocity (C-star, \(c^*\)), thrust coefficient (\(C_F\)), and specific impulse (\(I_{\text{sp}}\)).
For the above example, but choosing the rocket problem and specifying a nozzle area ratio
of 68.8, 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
*******************************************************************************
# Problem Type: "Rocket" (Infinite Area Combustor)
prob case=_______________3446 ro equilibrium
# Pressure (1 value):
p,psia= 3000
# 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
output short
output siunits
end
THEORETICAL ROCKET PERFORMANCE ASSUMING EQUILIBRIUM
COMPOSITION DURING EXPANSION FROM INFINITE AREA COMBUSTOR
Pin = 3000.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.7403 961.12
P, BAR 206.84 118.85 0.21521
T, K 3598.76 3381.67 1233.84
RHO, KG/CU M 9.4113 0 5.8080 0 2.9602-2
H, KJ/KG -986.31 -2161.66 -10544.8
U, KJ/KG -3184.12 -4208.00 -11271.8
G, KJ/KG -62768.7 -60217.2 -31727.0
S, KJ/(KG)(K) 17.1677 17.1677 17.1677
M, (1/n) 13.614 13.740 14.111
(dLV/dLP)t -1.01897 -1.01412 -1.00000
(dLV/dLT)p 1.3291 1.2605 1.0000
Cp, KJ/(KG)(K) 7.3140 6.6953 2.9097
GAMMAs 1.1475 1.1487 1.2539
SON VEL,M/SEC 1588.1 1533.2 954.8
MACH NUMBER 0.000 1.000 4.579
PERFORMANCE PARAMETERS
Ae/At 1.0000 68.800
CSTAR, M/SEC 2322.8 2322.8
CF 0.6601 1.8823
Ivac, M/SEC 2867.9 4538.6
Isp, M/SEC 1533.2 4372.3
MOLE FRACTIONS
*H 0.02543 0.02034 0.00000
HO2 0.00003 0.00002 0.00000
*H2 0.24740 0.24494 0.24402
H2O 0.68635 0.70506 0.75598
H2O2 0.00002 0.00001 0.00000
*O 0.00202 0.00123 0.00000
*OH 0.03659 0.02704 0.00000
*O2 0.00215 0.00137 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.8 m/s (based on combustion chamber conditions)
throat pressure of 118.85 bar and temperature of 3381.67 K
exit pressure of 0.21521 bar, temperature of 1233.84 K
at the nozzle exit, thrust coefficient = 1.8823, Isp = 4372.3 m/s
area_ratio = 68.8
pressure_throat_cea = Q_(118.85, 'bar').to('Pa')
temperature_throat_cea = 3381.67
pressure_exit_cea = Q_(0.21521, 'bar').to('Pa')
temperature_exit_cea = 1233.84
c_star_cea = 2322.8
thrust_coeff_cea = 1.8823
specific_impulse_cea = 4372.3
specific_impulse_vac_cea = 4538.6
C-star¶
We can calculate \(c^*\) directly using the combustion chamber state already obtained with Cantera:
def calculate_c_star(gamma, temperature, molecular_weight):
return (
np.sqrt(ct.gas_constant * temperature / (molecular_weight * gamma)) *
np.power(2 / (gamma + 1), -(gamma + 1) / (2*(gamma - 1)))
)
entropy_chamber = gas2.s
enthalpy_chamber = gas2.enthalpy_mass
mole_fractions_chamber = gas2.X
gamma_chamber = gamma_s
c_star = calculate_c_star(gamma_chamber, gas2.T, gas2.mean_molecular_weight)
print(f'c-star: {c_star: .1f} m/s')
print('Error in c-star: '
f'{100*np.abs(c_star - c_star_cea)/c_star_cea: .3e} %'
)
c-star: 2323.0 m/s
Error in c-star: 7.130e-03 %
Throat conditions¶
The nozzle flow from the combustion chamber to the throat is isentropic, and at the throat the flow velocity matches the sonic velocity. We need to iterate to determine the pressure and other properties.
From 1D isentropic flow assumptions, the equation
applies exactly, but only if \(\gamma_s\) remains constant from the chamber to the throat. This works for the frozen-flow assumption, but not for shifting equilibrium, where the gas composition will adjust with changing pressure and temperature.
We can use this equation to get a first estimate of throat pressure, \(p_{t,1}\), then equilibrate the gas mixture at \(s_c\) (chamber entropy) and \(p_{t,1}\).
The throat state is correct and converged when the velocity is sonic (i.e., equals the speed of sound). CEA checks for convergence using
where
using the properties at the current iteration. If the solution is not converged, we get an improved estimate for pressure:
where \(k\) is the iteration.
gas_throat = ct.Solution('nasa_h2o2.yaml', 'gas')
pressure_throat = pressure_chamber / np.power(
(gamma_chamber + 1) / 2., gamma_chamber / (gamma_chamber - 1)
)
# based on CEA defaults
max_iter_throat = 5
tolerance_throat = 0.4e-4
print('Throat iterations:')
mach = 1.0
num_iter = 0
residual = 1
while residual > tolerance_throat:
num_iter += 1
if num_iter == max_iter_throat:
break
print(f'Error: more than {max_iter_throat} iterations required for throat calculation')
pressure_throat = pressure_throat * (1 + gamma_s * mach**2) / (1 + gamma_s)
gas_throat.SPX = entropy_chamber, to_si(pressure_throat), mole_fractions_chamber
gas_throat.equilibrate('SP')
derivs = get_thermo_derivatives(gas_throat)
dlogV_dlogT_P, dlogV_dlogP_T, cp, gamma_s = get_thermo_properties(
gas_throat, derivs[0], derivs[1], derivs[2]
)
velocity = np.sqrt(2 * (enthalpy_chamber - gas_throat.enthalpy_mass))
speed_sound = np.sqrt(
ct.gas_constant * gas_throat.T * gamma_s / gas_throat.mean_molecular_weight
)
mach = velocity / speed_sound
residual = np.abs(1.0 - 1/mach**2)
print(f'{num_iter} {residual: .3e}')
temperature_throat = gas_throat.T
pressure_throat = Q_(gas_throat.P, 'Pa')
gamma_s_throat = gamma_s
print('Error in throat temperature: '
f'{100*np.abs(temperature_throat - temperature_throat_cea)/temperature_throat_cea: .3e} %'
)
print('Error in throat pressure: '
f'{100*np.abs(pressure_throat - pressure_throat_cea)/pressure_throat_cea: .3e~P} %'
)
Throat iterations:
1 9.420e-04
2 1.590e-06
Error in throat temperature: 2.640e-02 %
Error in throat pressure: 5.430e-03 %
Exit conditions¶
The conditions at the nozzle exit (or any location, really) can be determined with a given exit-to-throat area ratio \(A_e / A_t\), by an iterative approach.
First, calculate the area per unit mass flow rate at the throat:
where \(n_t = 1 / \text{MW}_t\) is the number of moles. Then, for a supersonic nozzle with an area ratio greater than two (\(A_e / A_t \geq 2\)), we can obtain an initial estimate for pressure ratio using an empirical formula:
where \(\gamma_s\) is evaluated using the throat state.
An improved estimate of pressure ratio for the next iteration can be found using:
where the derivative is
and the \(k\)th estimate of area ratio comes from
# this is constant
A_mdot_thr = gas_throat.T / (gas_throat.P * velocity * gas_throat.mean_molecular_weight)
gas_exit = ct.Solution('nasa_h2o2.yaml', 'gas')
gas_exit.SPX = gas_throat.s, gas_throat.P, gas_throat.X
# initial estimate for pressure ratio
pinf_pe = np.exp(gamma_s_throat + 1.4 * np.log(area_ratio))
p_exit = to_si(pressure_chamber) / pinf_pe
gas_exit.SP = entropy_chamber, p_exit
gas_exit.equilibrate('SP')
Ae_At = gas_exit.T / (gas_exit.P * velocity * gas_exit.mean_molecular_weight) / A_mdot_thr
print('Iter T_exit Ae/At P_exit P_inf/P')
num_iter = 0
print(f'{num_iter} {gas_exit.T:.3f} K {Ae_At: .2f} {gas_exit.P/1e5:.3f} bar {pinf_pe:.3f}')
max_iter_exit = 10
tolerance_exit = 4e-5
residual = 1
while np.abs(residual) > tolerance_exit:
num_iter += 1
if num_iter == max_iter_throat:
break
print(f'Error: more than {max_iter_exit} iterations required for exit calculation')
derivs = get_thermo_derivatives(gas_exit)
dlogV_dlogT_P, dlogV_dlogP_T, cp, gamma_s = get_thermo_properties(
gas_exit, derivs[0], derivs[1], derivs[2]
)
velocity = np.sqrt(2 * (enthalpy_chamber - gas_exit.enthalpy_mass))
speed_sound = np.sqrt(ct.gas_constant * gas_exit.T * gamma_s / gas_exit.mean_molecular_weight)
Ae_At = gas_exit.T / (gas_exit.P * velocity * gas_exit.mean_molecular_weight) / A_mdot_thr
dlogp_dlogA = gamma_s * velocity**2 / (velocity**2 - speed_sound**2)
residual = dlogp_dlogA * (np.log(area_ratio) - np.log(Ae_At))
log_pinf_pe = np.log(pinf_pe) + residual
pinf_pe = np.exp(log_pinf_pe)
p_exit = to_si(pressure_chamber) / pinf_pe
gas_exit.SP = entropy_chamber, p_exit
gas_exit.equilibrate('SP')
print(f'{num_iter} {gas_exit.T:.3f} K {Ae_At: .2f} {gas_exit.P/1e5:.3f} bar {pinf_pe:.3f}')
Iter T_exit Ae/At P_exit P_inf/P
0 1183.787 K 230.93 0.175 bar 1178.875
1 1234.196 K 80.36 0.215 bar 960.727
2 1234.177 K 68.80 0.215 bar 960.800
3 1234.177 K 68.80 0.215 bar 960.800
print(f'Exit temperature: {gas_exit.T: .2f} K')
print(f'Exit pressure: {Q_(gas_exit.P, "Pa").to("bar"): .5f~P}')
print()
print('Error in exit temperature: '
f'{100*np.abs(gas_exit.T - temperature_exit_cea)/temperature_exit_cea: .3e} %'
)
print('Error in exit pressure: '
f'{100*np.abs(Q_(gas_exit.P, "Pa") - pressure_exit_cea)/pressure_exit_cea: .3e~P} %'
)
Exit temperature: 1234.18 K
Exit pressure: 0.21528 bar
Error in exit temperature: 2.731e-02 %
Error in exit pressure: 3.329e-02 %
Those results look good! Now we can calculate thrust coefficient and specific impulse, using
CEA prints specific impulse with units of velocity, without the reference gravity term, so we will compute both versions for comparison.
derivs = get_thermo_derivatives(gas_exit)
dlogV_dlogT_P, dlogV_dlogP_T, cp, gamma_s = get_thermo_properties(
gas_exit, derivs[0], derivs[1], derivs[2]
)
velocity = np.sqrt(2 * (enthalpy_chamber - gas_exit.enthalpy_mass))
thrust_coeff = velocity / c_star
print(f'Thrust coefficient: {thrust_coeff: .4f}')
g0 = 9.80665
Isp = velocity
Ivac = Isp + gas_exit.T * ct.gas_constant / (velocity * gas_exit.mean_molecular_weight)
print(f'I_sp = {Isp: .1f} m/s')
print(f'I_vac = {Ivac: .1f} m/s')
print()
print('Error in Isp: '
f'{100*np.abs(Isp - specific_impulse_cea)/specific_impulse_cea: .3e} %'
)
print('Error in Ivac: '
f'{100*np.abs(Ivac - specific_impulse_vac_cea)/specific_impulse_vac_cea: .3e} %'
)
Thrust coefficient: 1.8822
I_sp = 4372.2 m/s
I_vac = 4538.5 m/s
Error in Isp: 2.903e-03 %
Error in Ivac: 2.513e-03 %
print('Actual specific impulse:')
print(f'I_sp = {Isp / g0: .1f} s')
print(f'I_vac = {Ivac / g0: .1f} s')
Actual specific impulse:
I_sp = 445.8 s
I_vac = 462.8 s
🎉 🚀 🛰
These comparisons look good!
The approach outlined here shows how we can perform calculations equivalent to CEA using Python and Cantera, at least for typical problems. CEA does support many additional options, including assuming frozen flow for the entire nozzle or after particular locations, calculations for a finite-area combustor (rather than an infinite-area combustor, as the work here assumes), and more.