<jupyter_start><jupyter_text>Implementing CEA calculations using Cantera<jupyter_code># 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
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'png')
plt.rcParams['figure.dpi']= 150
plt.rcParams['savefig.dpi'] = 150<jupyter_output><empty_output><jupyter_text>[CEA](https://www1.grc.nasa.gov/research-and-engineering/ceaweb/) (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](https://cearun.grc.nasa.gov/index.html).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](https://cantera.org/) is a modern software library for solving problems in chemical kinetics,thermodynamics, and transport, that offers a Python interface. Cantera natively supports phase and chemicalequilibrium 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 pressureGiven 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, suchas in nuclear thermal or electrothermal rockets.For example, say we have gaseous hydrazine (N2H4) as a propellant, with a chambertemperature of 5000 K and pressure of 50 psia. For this system, determine the equilibrium composition.In [CEA](https://cearun.grc.nasa.gov/), 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-atommolecules, based on the elements present: N2, H2, H, N, and HN.The CEA plaintext input file looks like:```textprob tp p,psia= 50  t,k= 5000reacname N2H4 mol 1.0output siunitsend```and the output is (with the repeated input removed):```text*******************************************************************************         NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA2, FEBRUARY 5, 2004                   BY  <NAME> AND <NAME>      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.<jupyter_code># 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()<jupyter_output>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
       [...]<jupyter_text>Comparing the results from CEA and Cantera, we see very good agreement between (most) thermodynamicproperties and the species mole fractions.But, the heat capacity $C_p$ appears **very** different, and if we calculate the specific heat ratio,$$\gamma = \frac{C_p}{C_v} \;,$$we will see it also differs quite substantially:<jupyter_code>gamma_ct = gas.cp_mole / gas.cv_mole
print(f'Cantera specific heat ratio: {gamma_ct: .4f}')<jupyter_output>Cantera specific heat ratio:  1.4917<jupyter_text>🤯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 {cite}`cea_analysis`:$$C_p = \left( \frac{\partial h}{\partial T} \right)_P \;.$$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:$$\begin{align}C_{p,e} &= C_{p,f} + C_{p,r} \\&= \sum_{j=1}^{N_s} n_j C_{p,j}^{\circ} + \sum_{j=1}^{N_g} n_j \frac{H_j^{\circ}}{T} \left( \frac{\partial \log n_j}{\partial \log T}\right)_P + \sum_{j=N_g+1}^{N_s} \frac{H_j^{\circ}}{T} \left( \frac{\partial n_j}{\partial \log T}\right)_P \;,\end{align}$$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 {cite}`cea_analysis` again, we can obtain this system of linear equations:$$\begin{align}\sum_{i=1}^{N_e} \sum_{j=1}^{N_g} a_{kj} a_{ij} n_j \left( \frac{\partial \pi_i}{\partial \log T}\right)_P + \sum_{j=N_g+1}^{N_s} a_{ij} \left( \frac{\partial n_j}{\partial \log T}\right)_P + \sum_{j=1}^{N_g} a_{kj} n_j \left( \frac{\partial \log n}{\partial \log T} \right)_P &= -\sum_{j=1}^{N_g} \frac{a_{kj} n_j H_j^{\circ}}{RT} \;, \quad k=1, \ldots, {N_e} \\\sum_{i=1}^{N_e} a_{ij} \left( \frac{\partial \pi_i}{\partial \log T}\right)_P &= - \frac{H_j^{\circ}}{RT} \;, \quad j = N_g + 1, \ldots, N_s \\\sum_{i=1}^{N_e} \sum_{j=1}^{N_g} a_{ij} n_j \left( \frac{\partial \pi_i}{\partial \log T}\right)_P &= -\sum_{j=1}^{N_g} \frac{n_j H_j^{\circ}}{RT} \\\sum_{i=1}^{N_e} \sum_{j=1}^{N_g} a_{kj} a_{ij} n_j \left( \frac{\partial \pi_i}{\partial \log P}\right)_T + \sum_{j=N_g + 1}^{N_s} a_{kj} \left( \frac{\partial n_j}{\partial \log P}\right)_T + \sum_{j=1}^{N_g} a_{ij} n_j \left( \frac{\partial \log n}{\partial \log P}\right)_T &= \sum_{j=1}^{N_g} a_{kj} n_j \;, \quad k=1, \ldots, {N_e} \\\sum_{i=1}^{N_e} a_{ij} \left( \frac{\partial \pi_i}{\partial \log P}\right)_T &= 0 \;, \quad j = N_g + 1, \ldots, N_s \\\sum_{i=1}^{N_e} \sum_{j=1}^{N_g} a_{ij} n_j \left( \frac{\partial \pi_i}{\partial \log P}\right)_T &= \sum_{j=1}^{N_g} n_j \;,\end{align}$$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:<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Using these derivatives, we can then calculate the specific heat, other relevant derivatives, and the ratio of specific heats:$$\begin{align}\frac{C_{p,e}}{R} &= \sum_{i=1}^{N_e} \left( \sum_{j=1}^{N_g} \frac{a_{ij} n_j H_j^{\circ}}{RT} \right) \left( \frac{\partial \pi_i}{\partial \log T}\right)_P + \sum_{j=N_g+1}^{N_s} \frac{H_j^{\circ}}{RT} \left( \frac{\partial n_j}{\partial \log T}\right)_P \\&+ \left( \sum_{j=1}^{N_g} \frac{n_j H_j^{\circ}}{RT} \right) \left( \frac{\partial \log n}{\partial \log T}\right)_P + \sum_{j=1}^{N_s} \frac{n_j C_{p,j}^{\circ}}{R} + \sum_{j=1}^{N_g} \frac{n_j (H_j^{\circ})^2}{R^2 T^2} \\\left( \frac{\partial \log V}{\partial \log T}\right)_P &= 1 + \left( \frac{\partial \log n}{\partial \log T}\right)_P \\\left( \frac{\partial \log V}{\partial \log P}\right)_T &= -1 + \left( \frac{\partial \log n}{\partial \log P}\right)_T \;.\end{align}$$The ratio of specific heats shows up via the speed of sound:$$\begin{align}a^2 &= \left( \frac{\partial P}{\partial \rho}\right)_s = -\frac{P}{\rho} \left( \frac{\partial \log P}{\partial \log V} \right)_s \\&= n R T \gamma_s\end{align} \;,$$where the ratio of specific heats is$$\gamma_s = \left( \frac{\partial \log P}{\partial \log \rho} \right)_s = - \frac{\gamma}{ \left( \frac{\partial \log V}{\partial \log P}\right)_T}$$and $$\gamma \equiv \frac{C_p}{C_v} \;.$$The constant volume specific heat is$$C_v \equiv \left( \frac{\partial u}{\partial T}\right)_V = C_p + \frac{ \frac{PV}{T} \left( \frac{\partial \log V}{\partial \log T}\right)_P^2}{ \left( \frac{\partial \log V}{\partial \log P}\right)_T} \;.$$<jupyter_code>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')<jupyter_output>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<jupyter_text>🎉 Success! These calculations agree very closely with those from CEA.  Adiabatic combustionCEA 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 {cite}`cea_analysis`:- 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 {cite}`NASA_thermo` represents liquid dinitrogen tetroxide (N2O4), which is an oxidizer used with hydrazine, with```textN2O4(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```textH2(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{cite}`cea_manual`, 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 lineIn general, both CEA and Cantera represent the thermodynamic properties of gaseous andcondensed 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  <NAME> AND <NAME>      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$ of1.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.000O2(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](https://cantera.org/documentation/dev/sphinx/html/yaml/phases.htmlsec-yaml-fixed-stoichiometry), 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 neededfor constraining the equilibrium solver.So, we can actually use the [ideal gas thermodynamic model](https://cantera.org/documentation/dev/doxygen/html/d7/dfa/classCantera_1_1IdealGasPhase.htmldetails) (`ideal-gas`) for the phase.For each species, we can use the [constant heat capacity](https://cantera.org/science/science-species.htmlconstant-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](https://cantera.org/tutorials/yaml/defining-phases.html), 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).<jupyter_code>h2o2_filename = 'h2o2_react.yaml'
print('Contents of ' + h2o2_filename + ':\n')
with open(h2o2_filename) as f:
    file_contents = f.read()
    print(file_contents)<jupyter_output>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<jupyter_text>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:<jupyter_code>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')<jupyter_output>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[...]<jupyter_text>🎉 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 calculationsCEA 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 ratioof 68.8, CEA provides this output:```*******************************************************************************         NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA2, FEBRUARY 5, 2004                   BY  <NAME> AND <NAME>      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<jupyter_code>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<jupyter_output><empty_output><jupyter_text>C-starWe can calculate $c^*$ directly using the combustion chamber state already obtained with Cantera:<jupyter_code>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} %'
      )<jupyter_output>c-star:  2323.0 m/s
Error in c-star:  7.130e-03 %<jupyter_text>Throat conditionsThe 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$$\frac{p_c}{p_t} = \left( \frac{\gamma_s + 1}{2} \right)^{\frac{\gamma_s}{\gamma_s - 1}}$$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 compositionwill adjust with changing pressure and temperature.We can use this equation to get a first estimate of throat pressure, $p_{t,1}$, then equilibratethe 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$$\left| \frac{u_t^2 - a_t^2}{u_t^2} \right| = \left| 1 - \frac{1}{M_t^2} \right| \leq 0.4 \cdot 10^{-4} \;,$$where $$\begin{align}M_t &= \frac{u_t}{a_t} \\u_t &= \sqrt{2 \left( h_c - h_t \right) } \\a_t &= \sqrt{ \gamma_s R T_t } \;,\end{align}$$using the properties at the current iteration.If the solution is not converged, we get an improved estimate for pressure:$$p_{t, k+1} = \left( p \frac{1 + \gamma_s M^2}{1 + \gamma_s} \right)_{t, k} \;,$$where $k$ is the iteration.<jupyter_code>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} %'
      )<jupyter_output>Throat iterations:
1   9.420e-04
2   1.590e-06
Error in throat temperature:  2.640e-02 %
Error in throat pressure: 5.430e-03 %<jupyter_text>Exit conditionsThe 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:$$\left( \frac{A}{\dot{m}} \right)_t = \frac{1}{\rho_t u_t} = \frac{T_t n_t \mathcal{R}}{p_t u_t} \;,$$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:$$\log \frac{p_c}{p_e} = \gamma_s + 1.4 \log \frac{A_e}{A_t} \;,$$where $\gamma_s$ is evaluated using the throat state.An improved estimate of pressure ratio for the next iteration can be found using:$$\left( \log \frac{p_c}{p_e} \right)_{k+1} = \left( \log \frac{p_c}{p_e} \right)_k + \left[ \left( \frac{\partial \log \frac{p_c}{p_e} }{\partial \log \frac{A_e}{A_t} } \right)_s \right]_k \times \left[ \log \frac{A_e}{A_t} - \left( \log \frac{A_e}{A_t} \right)_k \right] \;,$$where the derivative is$$\left( \frac{\partial \log \frac{p_c}{p_e} }{\partial \log \frac{A_e}{A_t} } \right)_s = \left( \frac{\gamma_s u^2}{u^2 - a^2} \right)_e$$and the $k$th estimate of area ratio comes from$$\left( \frac{A_e}{A_t} \right)_k = \left( \frac{T_e n_e \mathcal{R}}{p_e u_e} \right)_k \frac{1}{ \left(A/\dot{m}\right)_t } \;.$$<jupyter_code># 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}')
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} %'
      )<jupyter_output>Exit temperature:  1234.18 K
Exit pressure: 0.21528 bar

Error in exit temperature:  2.731e-02 %
Error in exit pressure: 3.329e-02 %<jupyter_text>Those results look good! Now we can calculate thrust coefficient and specific impulse, using$$\begin{align}C_F &= \frac{v_e}{c^*} \\I_{\text{sp}} &= \frac{v_e}{g_0} \\I_{\text{vac}} &= I_{\text{sp}} + \frac{p_e}{A_e}{\dot{m}} = I_{\text{sp}} + \frac{T_e \mathcal{R}}{v_e \overline{M}} \;.\end{align}$$CEA prints specific impulse with units of velocity, without the reference gravity term,so we will compute both versions for comparison.<jupyter_code>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} %'
      )
print('Actual specific impulse:')
print(f'I_sp = {Isp / g0: .1f} s')
print(f'I_vac = {Ivac / g0: .1f} s')<jupyter_output>Actual specific impulse:
I_sp =  445.8 s
I_vac =  462.8 s
{'fraction_non_alphanumeric': 0.12993243412955569, 'fraction_numerical': 0.05656443875116168, 'mean_word_length': 3.003821399839099, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 48, 'https://': 8, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 5, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '57750', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 15616, 'n_tokens_neox': 13843, 'n_words': 4028, 'path': 'book/thermochemistry/cea_cantera.ipynb', 'repo_name': 'kyleniemeyer/rocket-propulsion', 'size': '39861'}