Example 4 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 4 from RP-1311 [1] using the Python API. This is a UV equilibrium problem, which is otherwise identical to Example 3. For that reason, the comments in this example will be limited to the minor differences.

import numpy as np
import cea

# Options
trace = 1e-15

# Species
reac_names = [b"Air", b"C7H8(L)", b"C8H18(L),n-octa"]
omit_names = [b"CCN", b"CNC", b"C2N2", b"C2O",
              b"C3H4,allene", b"C3H4,propyne", b"C3H4,cyclo-", b"C3",
              b"C3H5,allyl", b"C3H6,propylene", b"C3H6,cyclo-", b"C3H3,propargyl",
              b"C3H6O", b"C3H7,n-propyl", b"C3H7,i-propyl", b"Jet-A(g)",
              b"C3O2", b"C4", b"C4H2", b"C3H8O,2propanol",
              b"C4H4,1,3-cyclo-", b"C4H6,butadiene", b"C4H6,2-butyne", b"C3H8O,1propanol",
              b"C4H8,tr2-butene", b"C4H8,isobutene", b"C4H8,cyclo-", b"C4H6,cyclo-",
              b"(CH3COOH)2", b"C4H9,n-butyl", b"C4H9,i-butyl", b"C4H8,1-butene",
              b"C4H9,s-butyl", b"C4H9,t-butyl", b"C4H10,isobutane", b"C4H8,cis2-buten",
              b"C4H10,n-butane", b"C4N2", b"C5", b"C3H8",
              b"C5H6,1,3cyclo-", b"C5H8,cyclo-", b"C5H10,1-pentene", b"C10H21,n-decyl",
              b"C5H10,cyclo-", b"C5H11,pentyl", b"C5H11,t-pentyl", b"C12H10,biphenyl",
              b"C5H12,n-pentane", b"C5H12,i-pentane", b"CH3C(CH3)2CH3", b"C12H9,o-bipheny",
              b"C6H6", b"C6H5OH,phenol", b"C6H10,cyclo-", b"C6H2",
              b"C6H12,1-hexene", b"C6H12,cyclo-", b"C6H13,n-hexyl", b"C6H5,phenyl",
              b"C7H7,benzyl", b"C7H8", b"C7H8O,cresol-mx", b"C6H5O,phenoxy",
              b"C7H14,1-heptene", b"C7H15,n-heptyl", b"C7H16,n-heptane", b"C10H8,azulene",
              b"C8H8,styrene", b"C8H10,ethylbenz", b"C8H16,1-octene", b"C10H8,napthlene",
              b"C8H17,n-octyl", b"C8H18,isooctane", b"C8H18,n-octane", b"C9H19,n-nonyl",
              b"Jet-A(L)", b"C6H6(L)", b"H2O(s)", b"H2O(L)"]

Here, the fixed value of \(u/R\) and \({\rho}\) (density) are pulled from the output of Example 3.

# Thermo states
u_R = -45.1343
density = 14.428  # kg/m^3

# Mixture states
fuel_weights = np.array([0.0, 0.4, 0.6])
oxidant_weights = np.array([1.0, 0.0, 0.0])
T_reac = np.array([700.0, 298.15, 298.15])
of_ratios = np.array([17.0])

# Mixtures
reac = cea.Mixture(reac_names)
prod = cea.Mixture(reac_names, products_from_reactants=True, omit=omit_names)

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

# Initialize variable arrays
n = len(of_ratios)
of_ratio_out = np.zeros(n)
T_out = np.zeros(n)
P_out = np.zeros(n)
rho = 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 = {}
trace_species = []
i = 0

Note that the prescribed value of \({\rho}\) (density) is converted to specific volume before calling solve().

# Equilibrium solve
for of_ratio in of_ratios:
    weights = reac.of_ratio_to_weights(oxidant_weights, fuel_weights, of_ratio)
    solver.solve(solution, cea.UV, u_R, 1.0/density, weights)

    # Store the output
    of_ratio_out[i] = of_ratio
    if solution.converged:
        P_out[i] = solution.P
        T_out[i] = solution.T
        rho[i] = solution.density
        enthalpy[i] = solution.enthalpy
        energy[i] = solution.energy
        gibbs[i] = solution.gibbs_energy
        entropy[i] = solution.entropy
        molecular_weight_M[i] = solution.M
        molecular_weight_MW[i] = solution.MW
        gamma_s[i] = solution.gamma_s
        cp_eq[i] = solution.cp_eq
        cp_fr[i] = solution.cp_fr
        cv_eq[i] = solution.cv_eq
        cv_fr[i] = 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

print("P, bar         ", 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, kg/m^3", 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("H, kJ/kg       ", 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, kJ/kg       ", 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, kJ/kg       ", 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, kJ/kg-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, kJ/kg-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, kJ/kg-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, kJ/kg-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, kJ/kg-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] > trace):
        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 = 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)))))

Resulting in the output to the terminal:

P, bar              99.974
T, K              2418.538
Density, kg/m^3  1.443e+01
H, kJ/kg           317.645
U, kJ/kg          -375.270
G, kJ/kg          -19437.2
S, kJ/kg-K           8.168
M, (1/n)            29.021
MW                  29.021
Gamma_s             1.2257
Cp_eq, kJ/kg-K      1.6093
Cp_fr, kJ/kg-K      1.4176
Cv_eq, kJ/kg-K      1.3121
Cv_fr, kJ/kg-K      1.1311

MOLE FRACTIONS

Ar               0.0088617
CN                5.84e-14
CO               0.0016763
CO2                0.11537
COOH            5.1528e-08
H               2.7535e-05
H2              0.00025052
H2O                0.10277
H2O2            9.8756e-07
HCHO,formaldehy 1.7093e-11
HCN             1.0535e-11
HCO             7.3708e-10
HCOOH           6.4546e-09
HNC             1.0139e-12
HNCO            1.2303e-09
HNO             4.2176e-07
HNO2            1.8492e-06
HNO3            1.1294e-09
HO2             7.9882e-06
N               1.1482e-08
N2                 0.73547
N2H2            5.1653e-13
N2O             3.6436e-06
N2O3            2.4361e-11
N2O4            3.0948e-15
N3              1.2321e-12
N3H              3.848e-13
NCO             8.3006e-11
NH              2.9308e-09
NH2             1.7069e-09
NH2NO2          2.4338e-15
NH2OH           1.0202e-11
NH3             3.8872e-09
NO               0.0067747
NO2             2.3462e-05
NO3             1.9237e-10
O               0.00015498
O2                0.026252
O3               1.209e-08
OH               0.0023475

TRACE SPECIES:
(HCOOH)2        C               C10H8,naphthale C2              C2H             C2H2,acetylene  C2H2,vinylidene C2H3,vinyl
C2H4            C2H4O,ethylen-o C2H5            C2H5OH          C2H6            C3H3,1-propynl  C3H3,2-propynl  C3H6O,acetone
C3H6O,propanal  C3H6O,propylox  C4H2,butadiyne  C4H6,1butyne    C4H6,2butyne    C6H14,n-hexane  C7H16,2-methylh CH
CH2             CH2CO,ketene    CH2OH           CH3             CH3CHO,ethanal  CH3CN           CH3CO,acetyl    CH3COOH
CH3N2CH3        CH3O            CH3O2CH3        CH3OCH3         CH3OH           CH3OOH          CH4             CNCOCN
CNN             HCCN            HCCO            HO(CO)2OH       N2H4            N2O5            NCN             O(CH)2O
OCCN            OHCH2COOH       C(gr)           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)