Example 13 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 13 from RP-1311 [1] using the Python API.

This is a rocket problem assuming an infinite-area combustor (IAC), using N₂H₄(L), Be(a), and H₂O₂(L) as reactants. This problem demonstrates how to call the insert parameter for a rocket problem, adding BeO(L) to the initial guess for the equilibrium evaluation. The resulting equilibrium mixtures also include non-negligible amounts of condensed species.

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 and initial temperatures. Currently, the Python API requires species names to be in bytes format, so we use the b"" syntax to create byte strings. The initial reactant temperatures, reac_T, will be used later to compute the chamber enthalpy. The amounts of each are specified through the fuel_weights, oxidant_weights, and pct_fuel variables specified here. Setting the fuel_weights array equal to [0.8, 0.2, 0.0] means that N₂H₄(L) constitutes 80% of the fuel, and Be(a) makes up the other 20% of the fuel mixture by weight. Similarly, setting the oxidant_weights array equal to [0.0, 0.0, 1.0] means that H₂O₂(L) constitutes 100% of the oxidant. These values will be used in conjunction later with pct_weight to compute the overall weight fraction array of the reactant mixture.

reac_names = [b"N2H4(L)", b"Be(a)", b"H2O2(L)"]
fuel_weights = np.array([0.8, 0.2, 0.0])
oxidant_weights = np.array([0.0, 0.0, 1.0])
reac_T = np.array([298.15, 298.15, 298.15])  # Reactant temperatures (K)
pct_fuel = 67.0

Next we set a variable for our trace parameter, defining the threshold below which species concentrations are assumed to be 0.

trace = 1e-10

Next, set some states for the rocket analysis. We will pass these values into the RocketSolver later.

pc = cea.units.psi_to_bar(3000)  # Chamber pressure (bar)
pi_p = [3.0, 10.0, 30.0, 300.0]   # Pressure ratio

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. We will also create the list of insert species that we will pass into the RocketSolver next.

reac = cea.Mixture(reac_names)
prod = cea.Mixture(reac_names, products_from_reactants=True)
insert = [b"BeO(L)"]

Now initialize the RocketSolver and RocketSolution class instances. We pass in the optional arguments insert and trace to the RocketSolver class initialization.

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

Compute the array of weight fractions for the reactant mixture:

of_ratio = (100.0 - pct_fuel) / pct_fuel
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, reac_T)/cea.R

Now we can solve the solve() function:

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

Finally, querry the solution variables and print them out:

num_pts = solution.num_pts
T = solution.T
P = cea.units.bar_to_atm(solution.P)
rho = 1e-3*solution.density
enthalpy = cea.units.joule_to_cal(solution.enthalpy)
energy = cea.units.joule_to_cal(solution.energy)
gibbs = cea.units.joule_to_cal(solution.gibbs_energy)
entropy = cea.units.joule_to_cal(solution.entropy)
n = solution.n
M_1n = solution.M
MW = solution.MW
cp_eq = cea.units.joule_to_cal(solution.cp_eq)
Mach = solution.Mach
gamma_s = solution.gamma_s
v_sonic = solution.sonic_velocity
ae_at = solution.ae_at
c_star = cea.units.m_per_s_to_ft_per_s(solution.c_star)
Cf = solution.coefficient_of_thrust
Isp = solution.Isp/9.80665
Isp_vac = solution.Isp_vacuum/9.80665

print("P, atm         ", 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, g/cm^3", end=" ")
for i in range(num_pts):
    if i < num_pts-1:
        print("{0:10.3g}".format(rho[i]), end=" ")
    else:
        print("{0:10.3g}".format(rho[i]))

print("H, cal/g       ", 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, cal/g       ", 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, cal/g       ", 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, cal/g-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, cal/g-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.4f}".format(gamma_s[i]), end=" ")
    else:
        print("{0:10.4f}".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.4f}".format(ae_at[i]), end=" ")
    else:
        print("{0:10.4f}".format(ae_at[i]))

print("C*, ft/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.4f}".format(Cf[i]), end=" ")
    else:
        print("{0:10.4f}".format(Cf[i]))

print("Ivac, lb-s/lb  ", 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, lb-s/lb   ", 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("MOLE FRACTIONS")
print("")
trace_species = []
for prod in solution.mole_fractions:
    if np.any(solution.mole_fractions[prod] > trace):
        print("{0:15s}".format(prod), end=" ")
        for j in range(len(solution.mole_fractions[prod])):
            if j < len(solution.mole_fractions[prod])-1:
                print("{0:10.5g}".format(solution.mole_fractions[prod][j]), end=" ")
            else:
                print("{0:10.5g}".format(solution.mole_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, atm             204.138    125.281     68.046     20.414      6.805      0.680
T, K              3002.540   2851.000   2851.000   2449.698   2068.730   1398.106
Density, g/cm^3     0.0138    0.00891    0.00483    0.00169   0.000669   9.91e-05
H, cal/g           -234.01    -403.87    -611.80    -995.93   -1291.19   -1761.87
U, cal/g           -592.87    -744.29    -952.71   -1287.95   -1537.45   -1928.24
G, cal/g          -10072.6    -9745.9    -9953.9    -9023.0    -8069.9    -6343.1
S, cal/g-K           3.277      3.277      3.277      3.277      3.277      3.277
M, (1/n)            16.627     16.643     16.619     16.670     16.694     16.700
MW                  13.431     13.420     13.405     13.378     13.377     13.378
Cp, cal/g-K          0.974      0.000      0.000      0.813      0.747      0.665
Gamma_s             1.1515     0.9979     0.9974     1.1791     1.1916     1.2180
Son. vel., m/s     1314.87    1192.19    1192.73    1200.25    1108.03     920.76
Mach                 0.000      1.000      1.491      2.104      2.684      3.883

PERFORMANCE PARAMETERS

Ae/At               0.0000     1.0000     1.2363     2.4856     5.3385    30.0004
C*, ft/s           6386.75    6386.75    6386.75    6386.75    6386.75    6386.75
Cf                  0.0000     0.6124     0.9133     1.2971     1.5279     1.8368
Ivac, lb-s/lb        0.000    243.396    263.111    306.823    338.619    384.464
Isp, lb-s/lb         0.000    121.571    181.306    257.481    303.294    364.613

MOLE FRACTIONS

Be              8.5544e-06 4.1033e-06 7.5534e-06 2.5505e-07 1.9172e-09 1.4393e-16
Be(OH)2          0.0072678  0.0057344  0.0057215  0.0014065 0.00020281 4.5262e-07
Be2O            1.1221e-06 4.0985e-07 7.5446e-07 5.5565e-09 5.3694e-12  5.606e-22
Be2O2           1.8492e-07 9.0835e-08 1.6729e-07 3.1776e-09 9.8507e-12 3.6727e-20
Be3O3           3.2166e-06 1.9697e-06 3.6277e-06  7.568e-08 2.5733e-10 1.1253e-18
Be4O4           9.2887e-07 7.4489e-07 1.3719e-06  3.683e-08 1.6849e-10 1.5897e-18
BeH             1.3243e-06 4.8345e-07 6.5487e-07 1.1009e-08 4.1984e-11 6.6186e-19
BeH2            7.0102e-06 2.8949e-06 2.8856e-06 8.8386e-08 9.5162e-10 6.8426e-16
BeN             4.3183e-08 1.3137e-08 1.7812e-08 1.6559e-10 2.9079e-13 4.1644e-22
BeO              4.359e-07 1.8631e-07 3.4313e-07 8.1179e-09 3.7299e-11  6.566e-19
BeOH            0.00022402 0.00012482 0.00016916 1.1082e-05  2.335e-07 7.1452e-13
H                 0.007146  0.0055921  0.0075786  0.0028413 0.00062797 3.6751e-06
H2                 0.51468    0.51518    0.51377    0.51519    0.51626    0.51664
H2O               0.053371   0.054925   0.054802   0.059185   0.060446   0.060664
H2O2            4.3791e-09 2.1323e-09 2.1285e-09  2.004e-10 7.8656e-12 3.4764e-16
HNO             8.6849e-08 4.1263e-08 4.1208e-08 3.7031e-09 1.5209e-10 8.6972e-15
HNO2            1.9242e-10 8.0962e-11 8.0894e-11 4.7618e-12 1.0221e-13 7.4123e-19
HO2             2.0568e-09 9.3034e-10 1.2621e-09 7.7342e-11 1.3662e-12 2.9585e-18
N               4.1872e-07 1.9115e-07  2.593e-07  1.691e-08 3.7932e-10 1.8929e-15
N2                 0.22449    0.22435    0.22415    0.22374    0.22373    0.22376
N2H2            7.2231e-09 2.7375e-09  1.481e-09 9.5737e-11 4.5049e-12 1.3131e-15
N2O             7.9158e-09 3.8882e-09 3.8866e-09  3.857e-10 1.7835e-11 1.4552e-15
N3H              1.411e-10 4.5333e-11 2.4548e-11 9.1327e-13 2.0287e-14 5.2666e-19
NH              2.3329e-06 1.0891e-06 1.0871e-06 9.1901e-08 3.6387e-09 1.7282e-13
NH2             1.8652e-05 9.8459e-06 7.2321e-06  1.105e-06 1.2065e-07 2.2644e-10
NH2OH           2.6057e-09 1.0468e-09 5.6608e-10   4.32e-11 2.3068e-12 9.2313e-16
NH3              0.0003067 0.00021088 0.00011398 4.9816e-05 2.7364e-05 1.3079e-05
NO              1.6315e-05 1.0305e-05 1.3985e-05 2.5602e-06  2.032e-07 4.6237e-11
O               2.0862e-06 1.1814e-06 2.1769e-06 2.3209e-07 7.2162e-09  5.494e-14
O2              8.5536e-08 5.0299e-08 9.2727e-08  1.098e-08 3.6115e-10 3.0766e-15
OH               0.0002671 0.00019071 0.00025858 7.0365e-05 9.3363e-06 1.0156e-08
BeO(L)             0.19218    0.17345   0.031752          0          0          0
BeO(a)                   0          0          0          0    0.19869    0.19892
BeO(b)                   0   0.020211    0.16164     0.1975          0          0

TRACE SPECIES:
Be2             HNO3            N2H4            N2O3            N2O4            N2O5            N3              NH2NO2
NO2             NO3             O3              Be(L)           Be(OH)2(b)      Be(a)           Be(b)           Be3N2(L)
Be3N2(a)        Be3N2(b)        H2O(L)          H2O(cr)

[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)