Example 12 from RP-1311

Note

The python script for this example is available in the source/bind/python/cea/samples directory of the CEA repository.

Here we describe how to run example 12 from RP-1311 [1] using the Python API.

This is a rocket problem assuming an infinite-area combustor (IAC), with the composition frozen at the throat, using CH₆N₂(L) and N₂O₄(L) as reactants.

First import the required libraries:

import numpy as np
import cea

Use cea.units for unit conversions and cea.R inline when normalizing enthalpy.

Declare the reactants and set their amounts. Currently, the Python API requires species names to be in bytes format, so we use the b"" syntax to create byte strings. The amounts of each are specified through the fuel_weights, oxidant_weights, and of_ratio variables specified here. Setting the fuel_weights array equal to [1.o, 0.0] means that CH₆N₂(L) constitutes 100% of the fuel, and similarly, setting the oxidant_weights array equal to [0.0, 1.0] means that N₂O₄(L) constitutes 100% of the oxidant. These values will be used in conjunction later with of_ratio to compute the overall weight fraction array of the reactant mixture.

reac_names = [b"CH6N2(L)", b"N2O4(L)"]
fuel_weights = np.array([1.0, 0.0])
oxidant_weights = np.array([0.0, 1.0])
of_ratio = 2.5

Next, set some states for the rocket analysis. We will pass these values into the RocketSolver later. Note that here we define a variable n_frz that we will use later to declare which rocket station should begin the frozen composition.

pc = cea.units.psi_to_bar(1000)  # Chamber pressure (bar)
pi_p = 68.0457   # Pressure ratio
supar = [5.0, 10.0, 25.0, 50.0, 75.0, 100.0, 150.0, 200.0]  # Supersonic area ratio
n_frz = 2

Now we will instantiate the reactant and product Mixture objects. To create the product Mixture, we pass the list of reactant names along with the flag products_from_reactants=True, which will return the full set of possible product species.

reac = cea.Mixture(reac_names)
prod = cea.Mixture([b"CO", b"CO2", b"H", b"HNO", b"HNO2", b"HO2",
                    b"H2", b"H2O", b"H2O2", b"N", b"NO", b"NO2",
                    b"N2", b"N2O", b"O", b"OH", b"O2", b"HCO", b"NH",
                    b"CH4", b"NH2", b"NH3", b"H2O(L)", b"C(gr)"])

Now initialize the RocketSolver and RocketSolution class instances.

solver = cea.RocketSolver(prod, reactants=reac)
solution = cea.RocketSolution(solver)

Compute the array of weight fractions for the reactant mixture:

weights = reac.of_ratio_to_weights(oxidant_weights, fuel_weights, of_ratio)

Compute the chamber enthalpy value based on the reactants weights and temperatures. We will pass this in later when we call solve(). Note that this value is normalized by R here.

hc = reac.calc_property(cea.ENTHALPY, weights, 298.15)/cea.R

Now we can solve the solve() function:

solver.solve(solution, weights, pc, pi_p, supar=supar, iac=True, hc=hc, n_frz=n_frz)

Finally, querry the solution variables and print them out:

num_pts = solution.num_pts
T = solution.T
P = solution.P
rho = solution.density
enthalpy = solution.enthalpy
energy = solution.energy
gibbs = solution.gibbs_energy
entropy = solution.entropy
n = solution.n
M_1n = solution.M
MW = solution.MW
cp_eq = solution.cp_eq
Mach = solution.Mach
gamma_s = solution.gamma_s
v_sonic = solution.sonic_velocity
ae_at = solution.ae_at
c_star = solution.c_star
Cf = solution.coefficient_of_thrust
Isp = solution.Isp
Isp_vac = solution.Isp_vacuum

print("P, bar         ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(P[i]), end=" ")
    else:
        print("{0:10.3f}".format(P[i]))

print("T, K           ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(T[i]), end=" ")
    else:
        print("{0:10.3f}".format(T[i]))

print("Density, kg/m^3", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(rho[i]), end=" ")
    else:
        print("{0:10.3f}".format(rho[i]))

print("H, kJ/kg       ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.2f}".format(enthalpy[i]), end=" ")
    else:
        print("{0:10.2f}".format(enthalpy[i]))

print("U, kJ/kg       ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.2f}".format(energy[i]), end=" ")
    else:
        print("{0:10.2f}".format(energy[i]))

print("G, kJ/kg       ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.1f}".format(gibbs[i]), end=" ")
    else:
        print("{0:10.1f}".format(gibbs[i]))

print("S, kJ/kg-K     ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(entropy[i]), end=" ")
    else:
        print("{0:10.3f}".format(entropy[i]))

print("M, (1/n)       ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(M_1n[i]), end=" ")
    else:
        print("{0:10.3f}".format(M_1n[i]))

print("MW             ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(MW[i]), end=" ")
    else:
        print("{0:10.3f}".format(MW[i]))

print("Cp, kJ/kg-K    ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(cp_eq[i]), end=" ")
    else:
        print("{0:10.3f}".format(cp_eq[i]))

print("Gamma_s        ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(gamma_s[i]), end=" ")
    else:
        print("{0:10.3f}".format(gamma_s[i]))

print("Son. vel., m/s ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.2f}".format(v_sonic[i]), end=" ")
    else:
        print("{0:10.2f}".format(v_sonic[i]))

print("Mach           ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(Mach[i]), end=" ")
    else:
        print("{0:10.3f}".format(Mach[i]))

print("")
print("PERFORMANCE PARAMETERS")
print("")

print("Ae/At          ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(ae_at[i]), end=" ")
    else:
        print("{0:10.3f}".format(ae_at[i]))

print("C*, m/s        ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.2f}".format(c_star[i]), end=" ")
    else:
        print("{0:10.2f}".format(c_star[i]))

print("Cf             ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(Cf[i]), end=" ")
    else:
        print("{0:10.3f}".format(Cf[i]))

print("Isp, vac., m/s ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(Isp_vac[i]), end=" ")
    else:
        print("{0:10.3f}".format(Isp_vac[i]))

print("Isp, m/s       ", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3f}".format(Isp[i]), end=" ")
    else:
        print("{0:10.3f}".format(Isp[i]))

print()
print("MASS FRACTIONS")
print("")
trace_species = []
for prod in solution.mass_fractions:
    if np.any(solution.mass_fractions[prod] > 5e-6):
        print("{0:15s}".format(prod), end=" ")
        for j in range(len(solution.mass_fractions[prod])):
            if j < len(solution.mass_fractions[prod])-1:
                print("{0:10.5g}".format(solution.mass_fractions[prod][j]), end=" ")
            else:
                print("{0:10.5g}".format(solution.mass_fractions[prod][j]))
    else:
        trace_species.append(prod)

print()
print("TRACE SPECIES:")
max_cols = 8
nrows = (len(trace_species) + max_cols - 1) // max_cols
for i in range(nrows):
    print(" ".join("{0:15s}".format(trace_species[j]) for j in range(i * max_cols, min((i + 1) * max_cols, len(trace_species)))))

This results in the following output to the terminal:

P, bar              68.948     39.814      1.013      2.068      0.790      0.228      0.090      0.052      0.035      0.020      0.014
T, K              3382.319   3203.673   1627.679   1868.526   1549.626   1202.295    983.296    870.192    796.165    700.162    637.715
Density, kg/m^3      5.846      3.606      0.181      0.321      0.148      0.055      0.026      0.017      0.013      0.008      0.006
H, kJ/kg            199.89    -426.81   -3389.58   -2960.26   -3525.96   -4113.12   -4463.73   -4637.96   -4749.30   -4890.46   -4980.28
U, kJ/kg           -979.48   -1530.94   -3950.56   -3604.24   -4060.04   -4527.49   -4802.62   -4937.86   -5023.70   -5131.77   -5200.07
G, kJ/kg          -36833.8   -35504.5   -21211.4   -23419.1   -20493.1   -17277.3   -15230.0   -14165.9   -13466.7   -12556.7   -11962.8
S, kJ/kg-K          10.949     10.949     10.949     10.949     10.949     10.949     10.949     10.949     10.949     10.949     10.949
M, (1/n)            23.845     24.125     24.125     24.125     24.125     24.125     24.125     24.125     24.125     24.125     24.125
MW                  23.845     24.125     24.125     24.125     24.125     24.125     24.125     24.125     24.125     24.125     24.125
Cp, kJ/kg-K          5.130      4.982      1.757      1.807      1.738      1.638      1.562      1.519      1.490      1.451      1.426
Gamma_s              1.138      1.135      1.244      1.236      1.247      1.266      1.283      1.294      1.301      1.312      1.319
Son. vel., m/s     1158.38    1119.55     835.41     892.06     816.21     724.39     659.43     622.86     597.49     562.57     538.37
Mach                 0.000      1.000      3.207      2.818      3.344      4.054      4.631      4.994      5.266      5.672      5.979

PERFORMANCE PARAMETERS

Ae/At                0.000      1.000      8.342      5.000     10.000     25.000     50.000     75.000    100.000    150.000    200.000
C*, m/s            1707.92    1707.92    1707.92    1707.92    1707.92    1707.92    1707.92    1707.92    1707.92    1707.92    1707.92
Cf                   0.000      0.656      1.569      1.472      1.598      1.720      1.788      1.821      1.842      1.868      1.885
Isp, vac., m/s       0.000   2105.782   2888.724   2770.175   2925.428   3078.093   3165.016   3206.992   3233.386   3266.348   3287.032
Isp, m/s             0.000   1119.549   2679.355   2514.019   2729.780   2937.008   3054.053   3110.577   3146.170   3190.720   3218.748

MASS FRACTIONS

CO                0.076919   0.067456   0.067456   0.067456   0.067456   0.067456   0.067456   0.067456   0.067456   0.067456   0.067456
CO2                0.15207    0.16694    0.16694    0.16694    0.16694    0.16694    0.16694    0.16694    0.16694    0.16694    0.16694
H               0.00043944 0.00033645 0.00033645 0.00033645 0.00033645 0.00033645 0.00033645 0.00033645 0.00033645 0.00033645 0.00033645
HNO             1.8511e-05 9.9871e-06 9.9871e-06 9.9871e-06 9.9871e-06 9.9871e-06 9.9871e-06 9.9871e-06 9.9871e-06 9.9871e-06 9.9871e-06
HNO2             8.301e-06 4.7208e-06 4.7208e-06 4.7208e-06 4.7208e-06 4.7208e-06 4.7208e-06 4.7208e-06 4.7208e-06 4.7208e-06 4.7208e-06
HO2             0.00014312 9.0662e-05 9.0662e-05 9.0662e-05 9.0662e-05 9.0662e-05 9.0662e-05 9.0662e-05 9.0662e-05 9.0662e-05 9.0662e-05
H2               0.0031271  0.0026704  0.0026704  0.0026704  0.0026704  0.0026704  0.0026704  0.0026704  0.0026704  0.0026704  0.0026704
H2O                0.28474    0.29296    0.29296    0.29296    0.29296    0.29296    0.29296    0.29296    0.29296    0.29296    0.29296
H2O2            2.2385e-05 1.3449e-05 1.3449e-05 1.3449e-05 1.3449e-05 1.3449e-05 1.3449e-05 1.3449e-05 1.3449e-05 1.3449e-05 1.3449e-05
N               5.0012e-06 2.5164e-06 2.5164e-06 2.5164e-06 2.5164e-06 2.5164e-06 2.5164e-06 2.5164e-06 2.5164e-06 2.5164e-06 2.5164e-06
NO                0.021824   0.017486   0.017486   0.017486   0.017486   0.017486   0.017486   0.017486   0.017486   0.017486   0.017486
NO2             4.8048e-05 3.1177e-05 3.1177e-05 3.1177e-05 3.1177e-05 3.1177e-05 3.1177e-05 3.1177e-05 3.1177e-05 3.1177e-05 3.1177e-05
N2                 0.38097    0.38301    0.38301    0.38301    0.38301    0.38301    0.38301    0.38301    0.38301    0.38301    0.38301
N2O             9.9328e-06  6.023e-06  6.023e-06  6.023e-06  6.023e-06  6.023e-06  6.023e-06  6.023e-06  6.023e-06  6.023e-06  6.023e-06
O                0.0050373  0.0037866  0.0037866  0.0037866  0.0037866  0.0037866  0.0037866  0.0037866  0.0037866  0.0037866  0.0037866
OH                0.034904    0.02888    0.02888    0.02888    0.02888    0.02888    0.02888    0.02888    0.02888    0.02888    0.02888
O2                0.039711   0.036317   0.036317   0.036317   0.036317   0.036317   0.036317   0.036317   0.036317   0.036317   0.036317

TRACE SPECIES:
HCO             NH              CH4             NH2             NH3             H2O(L)          C(gr)

[1]

McBride, B.J., Gordon, S., "Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications II. Users Manual and Program Description: Users Manual and Program Description - 2", NASA RP-1311, 1996. [NTRS](https://ntrs.nasa.gov/citations/19960044559)