Different output from running Python script and cea binary directly

[minh-dng]

I compared the same problem running through the binary and the Python distribution in bin/cea. You could see the numerical differences between the two program (despite the rounding of Python output), e.g. (C*).

Just want to double check that I haven't made any mistake in my Python program.

The input file is

reac
 fuel Acrylonitrile  C 3   H 3   N 1     wt%=39.92
 h,kj/mol=140.     t(k)=298.15
 fuel Butadiene  C 4   H 6     wt%=47.37
 h,kj/mol=90.50    t(k)=298.15
 fuel Styrene  C 8   H 8     wt%=12.71
 h,kj/mol=103.4     t(k)=298.15
 oxid N2O wt%=100 t(k)=298.15
prob rocket eql o/f=5.098591, p,bar=1.01325, supar=4.000000 ,pi/p=1
end

The input code to cea is

import cea
import numpy as np

acrylonitrile_args = {
    "name": "Acrylonitrile",
    "formula": {"C": 3, "H": 3, "N": 1},
    "enthalpy": 140,
    "temperature": 298.15,
    "molecular_weight": 53.0626,
}
"""https://webbook.nist.gov/cgi/cbook.cgi?ID=C107131&Mask=2"""

butadiene_args = {
    "name": "Butadiene",
    "formula": {"C": 4, "H": 6},
    "enthalpy": 90.50,
    "temperature": 298.15,
    "molecular_weight": 54.0904,
}
"""https://webbook.nist.gov/cgi/cbook.cgi?ID=C106990&Mask=7"""

styrene_args = {
    "name": "Styrene",
    "formula": {"C": 8, "H": 8},
    "enthalpy": 103.4,
    "temperature": 298.15,
    "molecular_weight": 104.1491,
}
"""https://webbook.nist.gov/cgi/cbook.cgi?ID=C100425&Mask=F"""

reac_names = [
    cea.Reactant(**acrylonitrile_args),
    cea.Reactant(**butadiene_args),
    cea.Reactant(**styrene_args),
    "N2O",
]
T_reactant = np.array([298.15] * 4)  # Reactant temperatures (K)
fuel_weights = np.array([0.3992, 0.4737, 0.1271, 0.0])
oxidant_weights = np.array([0.0] * 3 + [1.0])
of_ratio = 5.098591

pi_p = 1.0  # Pressure ratio
supar = 4.0  # Supersonic area ratio
of = 5.098591  # O/F ratio
pc = 1.01325  # Chamber pressure (bar)

reac = cea.Mixture(reac_names)
prod = cea.Mixture(reac_names, products_from_reactants=True)

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

weights = reac.of_ratio_to_weights(oxidant_weights, fuel_weights, of_ratio)
hc = reac.calc_property(cea.ENTHALPY, weights, T_reactant) / cea.R

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

# see below for the printing stuff

The output of cea binary

THEORETICAL ROCKET PERFORMANCE ASSUMING EQUILIBRIUM
COMPOSITION DURING EXPANSION FROM INFINITE AREA COMBUSTOR
 
Pc =   1.01 bar
 
pi/p =     1.0000
supar =     4.0000
 
 
            REACTANT     MOLES    ENERGY      TEMP
                                  KJ/MOL         K
     Acrylonitrile      39.920  140000.000    298.15
     Butadiene          47.370   90500.000    298.15
     Styrene            12.710  103400.000    298.15
     N2O               100.000   81600.000    298.15
 
O/F =    5.09859    % Fuel =  16.39723    r, Eq. Ratio =  0.00000    phi, Eq. Ratio =  0.00000
 

                        CHAMBER       THROAT         EXIT         EXIT
 Pinf/P                  1.0000       1.7327        1.000       21.448
 P, bar                  1.0132       0.5848       1.0132       0.0472
 T, K                   3054.12      2899.49      3054.12      2138.78
 Density, kg/m^3      0.9993E-1    0.6161E-1    0.9993E-1    0.7041E-2
 H, kJ/kg              1873.352     1334.037     1873.354     -707.805
 U, kJ/kg               859.344      384.915      859.345    -1378.776
 G, kJ/kg             -32116.28    -30934.65    -32116.29    -24510.47
 S, kJ/kg-K              11.129       11.129       11.129       11.129
 
 M, (1/n)               25.0427      25.4001      25.0427      26.5032
 (dln(V)/dln(P))t       -1.0299      -1.0241      -1.0299      -1.0016
 (dln(V)/dln(T))p        1.6003       1.5139       1.6003       1.0431
 Cp, kJ/(kg-K)           5.6400       5.2020       5.6400       1.9123
 Gamma_s                 1.1375       1.1365       1.1375       1.2150
 Son. Vel., m/s         1073.99      1038.57      1073.99       902.89
 Mach                    0.0000       1.0000          NaN       2.5164
 
 PERFORMANCE PARAMETERS
 
 Ae/At                                1.0000          NaN       4.0000
 C*, m/s                             1583.45      1583.45      1583.45
 Cf                                   0.6559          NaN       1.4349
 Ivac, m/s                           1952.44          NaN      2567.38
 Isp, m/s                            1038.57          NaN      2272.07
 
MOLE FRACTIONS
 
 CO                    0.226498     0.222595     0.226498     0.206825
 CO2                   0.050113     0.057965     0.050113     0.085919
 H                     0.042405     0.034232     0.042405     0.004395
 H2                    0.052549     0.051875     0.052549     0.058441
 H2O                   0.087324     0.097099     0.087324     0.123134
 N                     1.359E-5     6.549E-6     1.359E-5     1.947E-8
 N2                    0.487266     0.495331     0.487266     0.519707
 NO                    0.007722     0.005618     0.007722     1.445E-4
 O                     0.011919     0.007993     0.011919     6.210E-5
 O2                    0.007806     0.005955     0.007806     5.550E-5
 OH                    0.026376     0.021326     0.026376     0.001316

whereas the output of cea Python package

P, bar               1.013      0.584      1.013      0.046
T, K              2994.246   2836.530   2994.246   2000.870
Density, kg/m^3      0.103      0.064      0.103      0.007
H, kJ/kg           1550.02    1026.02    1550.02    -953.04
U, kJ/kg            566.87     106.96     566.87   -1579.59
G, kJ/kg          -31453.1   -30238.8   -31453.1   -23007.0
S, kJ/kg-K          11.022     11.022     11.022     11.022
M, (1/n)            25.322     25.661     25.322     26.552
MW                  25.322     25.661     25.322     26.552
Cp_eq, kJ/kg-K       5.156      4.604      5.156      1.684
Cp_fr, kJ/kg-K       1.566      1.560      1.566      1.506
Cv_eq, kJ/kg-K       4.415      3.961      4.415      1.360
Cv_eq, kJ/kg-K       1.237      1.236      1.237      1.192
Gamma_s              1.139      1.140      1.139      1.237
Son. vel., m/s     1058.25    1023.72    1058.25     880.53
Mach                 0.000      1.000        nan      2.541

PERFORMANCE PARAMETERS

Ae/At                0.000      1.000        nan      4.000
C*, m/s            1557.81    1557.81    1557.81    1557.81
Cf                  0.0000     0.6572        nan     1.4363
Isp, vac., m/s       0.000   1921.485        nan   2517.464
Isp, m/s             0.000   1023.721        nan   2237.435

MOLE FRACTIONS

CO                 0.22434    0.22018    0.22434    0.20474
CO2               0.055364   0.063269   0.055364   0.088546
H                  0.03523   0.027642    0.03523  0.0019052
H2                0.052091   0.051538   0.052091   0.061562
H2O               0.094997    0.10437   0.094997    0.12204
N               9.3531e-06 4.2282e-06 9.3531e-06 3.0869e-09
N2                 0.49336    0.50104    0.49336    0.52072
NO               0.0065091  0.0044611  0.0065091   3.67e-05
O                0.0087627  0.0054229  0.0087627 8.5008e-06
O2               0.0063217  0.0043898  0.0063217 7.2603e-06
OH                0.023015    0.01769   0.023015 0.00042668

---

# printing stuff
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
M_1n = solution.M
MW = solution.MW
cp_eq = solution.cp_eq
cp_fr = solution.cp_fr
cv_eq = solution.cv_eq
cv_fr = solution.cv_fr
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_eq, 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("Cp_fr, kJ/kg-K ", end=" ")
for i in range(num_pts):
    if i < num_pts - 1:
        print("{0:10.3f}".format(cp_fr[i]), end=" ")
    else:
        print("{0:10.3f}".format(cp_fr[i]))

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

print("Cv_eq, kJ/kg-K ", end=" ")
for i in range(num_pts):
    if i < num_pts - 1:
        print("{0:10.3f}".format(cv_fr[i]), end=" ")
    else:
        print("{0:10.3f}".format(cv_fr[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.4f}".format(Cf[i]), end=" ")
    else:
        print("{0:10.4f}".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("MOLE FRACTIONS")
print("")
trace_species = []
for prod in solution.mole_fractions:
    if np.any(solution.mole_fractions[prod] > 5e-6):
        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)

[markleader]

Thanks for reporting this and for comparing the binary and Python results so carefully.

I took a closer look, and I do not think this is a solver divergence. The difference comes from the enthalpy units used in the Python API. Right now, the Python Reactant API interprets enthalpy using J/kg, while the .inp file used kJ/mol.
Changing the values in your Python input to match these default units should make the Python and binary results consistent today.

I am also working on improving this in the Python API so the units are more explicit. The plan is to add an enthalpy_units argument going forward. Initially this will be optional, but recommended, so existing code keeps working. After that transition, I expect to require users to pass enthalpy units explicitly rather than relying on defaults, which should avoid this kind of confusion entirely.

Thanks again for flagging it. This was a helpful catch.