Example 1 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 1 from RP-1311 [1] using the Python API. This example is a TP equilibrium problem, with H₂ and Air as reactants.

First import the required libraries:

import numpy as np
import cea

Use cea.units for unit conversions as needed.

# Unit conversions are handled with cea.units helpers.

Declare the product and reactnat species names. Currently, the Python API requires species names to be in bytes format, so we use the b"" syntax to create byte strings. Also note that prod_names is optional in general, but in this case, we explicitly define the list of species that we want to be in the final mixture (note for experienced CEA-users: this is akin to the only parameter in the legacy interface).

reac_names = [b"H2", b"Air"]
prod_names = [b"Ar",   b"C",   b"CO",  b"CO2", b"H",
              b"H2",   b"H2O", b"HNO", b"HO2", b"HNO2",
              b"HNO3", b"N",   b"NH",  b"NO",  b"N2",
              b"N2O3", b"O",   b"O2",  b"OH",  b"O3"]

Define the thermodynamic states at which we want to solve the equilibrium problem, in SI units.

pressures = cea.units.atm_to_bar(np.array([1.0, 0.1, 0.01]))
temperatures = np.array([3000.0, 2000.0])

Define the amounts of each reactant; in this case, a chemical equivalence ratio chem_eq_ratios is prescribed (\(r_{eq}\) in the RP-1311 [2]). The arrays fuel_moles and oxidant_moles correspond to the reac_names list, and sets the mole fraction of each that is part of the fuel and oxidant mixtures, respecttively. In this case, because we are using one fuel and one oxidizer, these are simply 1.0 to indicate which reactant is the fuel and which is the oxidizer.

fuel_moles = np.array([1.0, 0.0])
oxidant_moles = np.array([0.0, 1.0])
chem_eq_ratios = np.array([1.0, 1.5])

Now having defined all of the relevant inputs to the problem, we can begin creating the required CEA objects, starting with the Mixture.

reac = cea.Mixture(reac_names)
prod = cea.Mixture(prod_names)

Next we instantiate the EqSolver and EqSolution objects.

solver = cea.EqSolver(prod, reactants=reac)
solution = cea.EqSolution(solver)

Next, we convert the chem_eq_ratios to oxidant-to-fuel ratios. This also requires first converting the fuel_moles and oxidant_moles to fuel_weights and oxidant_weights, respectively.

fuel_weights = reac.moles_to_weights(fuel_moles)
oxidant_weights = reac.moles_to_weights(oxidant_moles)
of_ratios = len(chem_eq_ratios)*[0.0]
for i, eqrat in enumerate(chem_eq_ratios):
    of_ratios[i] = reac.chem_eq_ratio_to_of_ratio(oxidant_weights,
                                                fuel_weights,
                                                eqrat)

We will now initialize an array to store each of the solution variables for printing the output later.

n = len(chem_eq_ratios)*len(pressures)*len(temperatures)
of_ratio_out = np.zeros(n)
T_out = np.zeros(n)
P_out = np.zeros(n)
rho = np.zeros(n)
volume = np.zeros(n)
enthalpy = np.zeros(n)
energy = np.zeros(n)
gibbs = np.zeros(n)
entropy = np.zeros(n)
molecular_weight_M = np.zeros(n)
molecular_weight_MW = np.zeros(n)
gamma_s = np.zeros(n)
cp_eq = np.zeros(n)
cp_fr = np.zeros(n)
cv_eq = np.zeros(n)
cv_fr = np.zeros(n)
mole_fractions = {}
i = 0

Finally, we can loop through the defined pressures, temperatures, and oxidant-to-fuel ratios to solve the equilibrium problem at each state. We will also retrieve the solution variables and store them in the arrays we just initialized, and convert some units before storing. The key points to note here are:

1. The solve() requires a list of reactant weights, which we compute using the of_ratio_to_weights method of the cea.Mixture class.

2. The syntax of the solve() method is solver.solve(solution, cea.TP, temperature, pressure, weights), where weights is the list of reactant weights computed from the oxidant-to-fuel ratio.

for of_ratio in of_ratios:
    for p in pressures:
        for t in temperatures:
            weights = reac.of_ratio_to_weights(oxidant_weights, fuel_weights, of_ratio)
            solver.solve(solution, cea.TP, t, p, weights)

            # Store the output
            of_ratio_out[i] = of_ratio
            T_out[i] = t
            P_out[i] = cea.units.bar_to_atm(p)
            if solution.converged:
                rho[i] = solution.density*1.e-3
                volume[i] = solution.volume*1.e3
                enthalpy[i] = cea.units.joule_to_cal(solution.enthalpy)
                energy[i] = cea.units.joule_to_cal(solution.energy)
                gibbs[i] = cea.units.joule_to_cal(solution.gibbs_energy)
                entropy[i] = cea.units.joule_to_cal(solution.entropy)
                molecular_weight_M[i] = solution.M
                molecular_weight_MW[i] = solution.MW
                gamma_s[i] = solution.gamma_s
                cp_eq[i] = cea.units.joule_to_cal(solution.cp_eq)
                cp_fr[i] = cea.units.joule_to_cal(solution.cp_fr)
                cv_eq[i] = cea.units.joule_to_cal(solution.cv_eq)
                cv_fr[i] = cea.units.joule_to_cal(solution.cv_fr)

            if i == 0:
                for prod in solution.mole_fractions:
                    mole_fractions[prod] = np.array([solution.mole_fractions[prod]])
            else:
                for prod in mole_fractions:
                    mole_fractions[prod] = np.append(mole_fractions[prod], solution.mole_fractions[prod])

            i += 1

Finally, print everything out in a formatted manner consistent with the legacy CEA output format.

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

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

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

print("Density, g/cc   ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3e}".format(rho[i]), end=" ")
    else:
        print("{0:10.3e}".format(rho[i]))

print("Volume, cc/g    ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.3e}".format(volume[i]), end=" ")
    else:
        print("{0:10.3e}".format(volume[i]))

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

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

print("G, cal/g        ", end="")
for i in range(n):
    if i < n-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(n):
    if i < n-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(n):
    if i < n-1:
        print("{0:10.3f}".format(molecular_weight_M[i]), end=" ")
    else:
        print("{0:10.3f}".format(molecular_weight_M[i]))

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

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

print("Cp_eq, cal/g-K  ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4f}".format(cp_eq[i]), end=" ")
    else:
        print("{0:10.4f}".format(cp_eq[i]))

print("Cp_fr, cal/g-K  ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4f}".format(cp_fr[i]), end=" ")
    else:
        print("{0:10.4f}".format(cp_fr[i]))

print("Cv_eq, cal/g-K  ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4f}".format(cv_eq[i]), end=" ")
    else:
        print("{0:10.4f}".format(cv_eq[i]))

print("Cv_fr, cal/g-K  ", end="")
for i in range(n):
    if i < n-1:
        print("{0:10.4f}".format(cv_fr[i]), end=" ")
    else:
        print("{0:10.4f}".format(cv_fr[i]))

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

print()
print("TRACE SPECIES:")
max_cols = 10
nrows = (len(trace_species) + max_cols - 1) // max_cols
for i in range(nrows):
    print(" ".join("{0:10s}".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:

o/f                 34.296     34.296     34.296     34.296     34.296     34.296     22.853     22.853     22.853     22.853     22.853     22.853
P, atm               1.000      1.000      0.100      0.100      0.010      0.010      1.000      1.000      0.100      0.100      0.010      0.010
T, K              3000.000   2000.000   3000.000   2000.000   3000.000   2000.000   3000.000   2000.000   3000.000   2000.000   3000.000   2000.000
Density, g/cc    9.178e-05  1.499e-04  8.085e-06  1.496e-05  6.614e-07  1.487e-06  8.122e-05  1.297e-04  7.116e-06  1.296e-05  5.672e-07  1.293e-06
Volume, cc/g     1.090e+04  6.671e+03  1.237e+05  6.687e+04  1.512e+06  6.723e+05  1.231e+04  7.707e+03  1.405e+05  7.714e+04  1.763e+06  7.735e+05
H, cal/g           663.714   -203.569   1369.546   -191.833   2647.099   -164.346    718.654   -120.698   1550.122   -116.204   3208.235   -101.779
U, cal/g           399.856   -365.131   1070.016   -353.764   2280.923   -327.156    420.475   -307.350   1209.797   -303.012   2781.243   -289.088
G, cal/g           -7974.5    -5290.4    -8616.7    -5662.8    -9381.0    -6036.5    -8818.8    -5830.6    -9545.4    -6260.6   -10424.5    -6691.2
S, cal/g-K           2.879      2.543      3.329      2.735      4.009      2.936      3.179      2.855      3.698      3.072      4.544      3.295
M, (1/n)            22.594     24.600     19.903     24.544     16.281     24.412     19.994     21.293     17.518     21.275     13.962     21.219
MW                  22.594     24.600     19.903     24.544     16.281     24.412     19.994     21.293     17.518     21.275     13.962     21.219
Gamma_s             1.1312     1.2258     1.1206     1.2026     1.1318     1.1668     1.1337     1.2529     1.1196     1.2389     1.1296     1.2051
Cp_eq, cal/g-K      1.6816     0.4552     3.4398     0.5215     3.7168     0.6858     1.8253     0.4671     4.1929     0.5000     4.9230     0.6077
Cp_fr, cal/g-K      0.4250     0.4008     0.4282     0.4008     0.4368     0.4010     0.4812     0.4520     0.4853     0.4521     0.4968     0.4522
Cv_eq, cal/g-K      1.4371     0.3711     2.8458     0.4330     3.0554     0.5856     1.5585     0.3727     3.4492     0.4034     4.0091     0.5032
Cv_fr, cal/g-K      0.3370     0.3200     0.3284     0.3199     0.3148     0.3196     0.3818     0.3587     0.3719     0.3587     0.3544     0.3586

MOLE FRACTIONS

Ar               0.0070982  0.0077283  0.0062528  0.0077107  0.0051148  0.0076691  0.0061933  0.0065959  0.0054263  0.0065904  0.0043249  0.0065728
CO              0.00017073 1.0365e-05 0.00018417 2.0965e-05 0.00016776  4.081e-05 0.00017519 0.00015634 0.00016627  0.0001561  0.0001428 0.00015539
CO2             7.1077e-05 0.00025289 2.8816e-05 0.00024169 6.4686e-06 0.00022042 3.5772e-05 6.8336e-05 1.8562e-05 6.8384e-05 4.5151e-06 6.8501e-05
H                 0.040752 8.9523e-05    0.14289 0.00040886    0.31894  0.0018585    0.06025 0.00062121    0.18236  0.0019603    0.39312  0.0061596
H2                0.067277  0.0030607   0.082718  0.0063841   0.041209   0.013191    0.14706    0.14737    0.13472    0.14676   0.062606     0.1449
H2O                 0.2073    0.34207   0.095791    0.33714   0.011761    0.32637    0.22224    0.29508    0.11131     0.2945   0.014651    0.29261
HO2             1.0257e-05 1.0068e-07 5.0804e-06 1.0266e-07  6.887e-07 1.0243e-07 3.6481e-06 2.2423e-10 3.3005e-06 7.1074e-10 5.7072e-07 2.2616e-09
N               1.0577e-05 7.1874e-10 3.1348e-05   2.27e-09 8.9735e-05 7.1582e-09 9.9041e-06 6.6412e-10 2.9248e-05 2.0992e-09 8.2567e-05 6.6293e-09
NO                0.012303 0.00048345   0.013705  0.0007215  0.0096685  0.0010659  0.0056508  8.003e-06  0.0091234 2.5353e-05  0.0072944 8.0572e-05
N2                 0.58568    0.64414    0.51448    0.64255    0.42158    0.63891    0.51356    0.54995    0.44786    0.54948    0.35692    0.54799
O                 0.015397 2.1253e-05   0.057868 0.00010042    0.14261  0.0004705  0.0075519 3.8075e-07   0.041288  3.816e-06    0.11694 3.8402e-05
O2                0.018761  0.0010234   0.026501  0.0022849   0.016095  0.0050156  0.0045133 3.2847e-07   0.013491 3.2993e-06   0.010821 3.3412e-05
OH                0.045174  0.0011285   0.059534  0.0024353   0.032747  0.0051864   0.032758 0.00014029   0.054207 0.00044369   0.033096   0.001403

TRACE SPECIES:
C          HNO        HNO2       HNO3       NH         N2O3       O3

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

[2]

Gordon, S., McBride, B.J., "Computer program for calculation of complex chemical equilibrium compositions and applications. Part 1: Analysis", NASA RP-1311, 1994. [NTRS](https://ntrs.nasa.gov/citations/19950013764)