ABS Hybrid Rocket Equilibrium Calculator

Chemical-equilibrium combustion properties for an ABS (acrylonitrile-butadiene-styrene) / oxidizer hybrid rocket motor, computed with Cantera. Feeds a lookup table of combustion properties (T_ad, gamma, MW, c*, ...) vs. (O/F, chamber pressure) to a separate transient chamber model, replacing that model's placeholder combustion assumption with real equilibrium chemistry.

Current Bambu Black basis

The checked-in defaults now match the planned physical test rig:

The inactive calorimetry block remains available as an alternative input path. The Hybrid model still flags its heat capacity, gasification heat, combustion efficiency, Prandtl number, and product viscosity as provisional.

Quick start

pip install cantera numpy scipy matplotlib pandas pyyaml shapely
python run.py --config config.yaml --of-min 1.1 --of-max 12 --n 90 --pressure-kpa 500 --pressure-reference gauge

Useful flags: - --no-soot - also runs a gas-only comparison sweep and plots T_ad with vs. without condensed carbon, so you can see how much soot formation actually matters. - --mechanism gri30 - also runs the sweep with full GRI-Mech 3.0 species and plots it against the curated nasa_gas.yaml mechanism (see below for why these should be close but not identical). - --area-ratio <value> - adds a frozen-flow Isp estimate at that nozzle expansion ratio to the results table.

Outputs land in output/: T_ad_vs_OF.png (primary deliverable), species_vs_OF.png, c_star_vs_OF.png, sweep_nominal.csv, combustion_table.npz, and (if requested) soot_comparison.png / mechanism_comparison.png.

Run the test suite (anchor, shape, soot, conservation, ABS-consumption, calorimetry sign-convention checks) with:

pip install pytest
python -m pytest tests/ -v

Cantera environment used to build this

Verified on Cantera 3.2.0, installed via pip install --user cantera (a plain pip install cantera failed here with a permissions error writing to a shared C:\Python312\share directory; --user sidesteps that). cantera.get_data_directories() returned the package's bundled .../site-packages/cantera/data directory, which includes (among others) nasa_gas.yaml, nasa_condensed.yaml, graphite.yaml, and gri30.yaml. Exactly what this tool needs. If you're on a different install, rerun run.py and check its startup banner: it always prints the Cantera version and data directories first, before doing anything else, per your explicit request.

Why not GRI-Mech 3.0? (mechanism/thermo-database choice)

The meeting notes suggested GRI-Mech 3.0. I pushed back on that, and here's the reasoning:

Equilibrium needs a thermodynamic database, not a kinetic mechanism. Chemical equilibrium is a Gibbs-free-energy minimization over a fixed set of species at fixed elemental abundances. Reaction rates (which is what a kinetic mechanism like GRI-Mech actually encodes: Arrhenius parameters, third-body efficiencies, pressure falloff) never enter the calculation at all. All that matters is (a) which species are allowed to exist in the products, and (b) each one's thermodynamic (NASA-polynomial) data.

GRI-Mech 3.0's species list is a byproduct of what's needed for methane/air ignition and flame-speed kinetics: 53 species, tuned and validated for that specific chemistry. It's missing or poorly covers much of the fuel-rich, soot-precursor, high-C/H chemistry that matters here: no solid carbon phase, thin coverage of larger hydrocarbon radicals (C3, C4H2, C2N2, etc.), and essentially no nitrile/HCN chemistry despite ABS containing nitrogen (from the acrylonitrile component).

Instead, this project builds the gas phase directly from Cantera's bundled nasa_gas.yaml thermodynamic database (748 species, NASA-7/9 polynomials, no kinetics at all) using a curated species list (config.yaml → mechanism.species) chosen to cover the C/H/N/O2 equilibrium products relevant to ABS combustion. CO, CO2, H2O, H2, the usual radical pool, plus the larger-hydrocarbon and nitrogen species (HCN, CN, C2N2, NH3, C2H2, C4H2, ...) that matter fuel-rich. Condensed carbon comes from a second phase (graphite.yaml) coupled in via a multiphase ct.Mixture; see "Soot / condensed carbon" below.

This project still supports --mechanism gri30 as an explicit side-by-side comparison (see mechanism_comparison.png) so you can show your supervisor exactly how much the species-list choice matters, quantitatively, rather than asserting it. On the validation anchor case (stoichiometric CH4/air @ 1 atm), gri30 gives T_ad = 2224.6 K and the curated nasa_gas.yaml mechanism gives T_ad = 2225.1 K. Both essentially exact against the textbook ~2225 K value, which is expected since methane/air combustion doesn't touch any of the species where the two mechanisms actually differ (soot precursors, nitriles, etc.).

Two curated-list species needed a name fix-up (config.yaml → mechanism.species_aliases), verified against the installed Cantera 3.2.0 data: - C2H2 → C2H2,acetylene (nasa_gas.yaml also separately lists the vinylidene isomer) - CH2O → HCHO,formaldehy

If you upgrade Cantera and species_setup.py raises a "species not found" error, re-verify these names against your installed nasa_gas.yaml.

Cantera for beginners

A few Cantera concepts this codebase leans on, if you haven't used it before:

• Solution: couples a thermodynamic model with a list of species. Ours has no kinetics (no reactions); see above, equilibrium doesn't need any.
• gas.equilibrate('HP'): finds the composition (and temperature) that minimizes Gibbs free energy at fixed enthalpy (H) and pressure (P). This is the adiabatic flame temperature calculation: H is conserved because the process is adiabatic (no heat added or removed), and P is conserved because the chamber is modeled as constant-pressure. It solves for composition and T simultaneously.
• ct.Mixture: when solid carbon can form, a single ideal-gas Solution can't represent it; there's no "mole fraction of a non-gas species" in a gas phase. ct.Mixture([(gas, n_gas), (carbon, n_carbon)]) couples the gas phase and the condensed graphite phase into one multiphase system and minimizes total Gibbs energy across both, with total elemental abundances (set via the gas phase's initial composition) held fixed. mix.equilibrate('HP', ...) is the multiphase equivalent. Important: after mix.equilibrate(...), the original gas/carbon objects you passed into the Mixture constructor are left synced to the equilibrium state; you read results (T, X, density, ...) directly off them, not off the Mixture object (which only tracks the extensive state: how many kmol of each phase).
• ct.Species / separate reactant and product phases: the custom ABS species is built from a small YAML mapping via ct.Species.list_from_yaml(...) and added to the reactant Solution. The product Solution uses the same base database without ABS because the placeholder ABS record supplies reactant enthalpy, not product thermodynamic data.

The custom ABS species: units and self-consistency

ABS isn't in any thermodynamic database, so it's defined as a small custom species (species_setup.py) using Cantera's constant-cp thermo model: no meaningful preheat happens in this application (reactants enter near 298.15 K), so a full temperature-dependent NASA polynomial for ABS would be unnecessary complexity.

species:
- name: ABS
  composition: {C: 5.381, H: 6.090, N: 0.240}
  thermo:
    model: constant-cp
    T0: 298.15
    h0: 66250 J/mol
    s0: 0 J/mol/K
    cp0: 0 J/mol/K

Composition and h0 must be self-consistent. The absolute stoichiometric coefficients (5.381, 6.090, 0.240) are tied to the measured Bambu formula unit. They cancel out in mole fractions and in the O/F-normalized energy balance. But h0 MUST be expressed per mole of exactly this formula unit. Change the composition without recomputing h0 for the same formula unit, and every downstream O/F and energy number is silently wrong; there is no way for Cantera to catch this for you. calorimetry.py and species_setup.py always derive/consume these together so they can't drift apart.

The units footgun: Cantera's internal molar basis is the kilomole, not the mole, even though the species YAML above is authored with quantity: mol for readability. gas.partial_molar_enthalpies and friends come back in J/kmol: a factor of 1000 away from the YAML's h0. src/utils.py provides j_per_mol_to_j_per_kmol / j_per_kmol_to_j_per_mol helpers so this conversion never happens as a bare * 1000 scattered through the code, and tests/test_units.py round-trips the conversion and verifies against a real Cantera Solution object (writes h0 = -6e4 J/mol, reads back partial_molar_enthalpies = -6e7 J/kmol) so a future refactor can't silently drop the factor of 1000.

s0 = 0 and cp0 = 0 are modelling conveniences, not measured values. Only the reactant enthalpy at the configured inlet temperature is used in the combustion energy balance. ABS is excluded from the allowed product phase, following NASA CEA's separate reactant and product-list formulation. A physical cp0 (about 1400 J/kg·K for ABS) is also available via abs_fuel.cp0_model: physical in config.yaml for studies with a reactant temperature away from the 298.15 K reference point.

Soot / condensed carbon

ABS burns sooty, and a gas-only phase can't represent solid carbon. species_setup.py loads Cantera's bundled graphite.yaml as a second Solution; verified (don't assume from memory) against the installed Cantera 3.2.0: phase name graphite, single species C(gr), NASA-7 thermo with a constant-volume equation of state (2.16 g/cm³). The two phases are coupled via ct.Mixture for a true multiphase equilibrium (see "Cantera for beginners" above).

Run with --no-soot to also compute the gas-only comparison and see soot_comparison.png. On the current Bambu Black/N2O sweep, allowing soot changes T_ad by very little near or lean of stoichiometric but matters substantially in fuel-rich states.

Solver robustness

The reactant phase first supplies the target enthalpy and total C, H, N, and O amounts. An element-equivalent seed initializes the separate product phase. The solver then root-finds temperature with scipy.optimize.brentq over [200, 5000] K, solving H_products(T) - H_reactants = 0. Each trial temperature performs a constant-TP multiphase equilibrium, trying Cantera's VCS solver and then its Gibbs solver with a larger iteration allowance. This preserves reactant enthalpy and elements without admitting pseudo-ABS as a product.

Every point is wrapped in try/except inside equilibrium_point(). A state with no complete-product enthalpy root in the supported product-thermo range returns converged=False and NaN fields rather than inventing incomplete-combustion data.

Bomb calorimetry → enthalpy of formation (calorimetry.py)

Converts a measured (or placeholder) bomb-calorimeter heat of combustion into Hf_min / Hf_nominal / Hf_max for ABS:

1. Bomb calorimeters measure heat of combustion at constant volume (ΔU_c, typically reported as a positive HHV in MJ/kg, with liquid product water).
2. Convert to constant pressure: ΔH_c = ΔU_c + Δn_gas·R·T, where for C_x H_y N_z + (x + y/4) O2 → x CO2 + (y/2) H2O(l) + (z/2) N2, Δn_gas = (x + z/2) − (x + y/4) = z/2 − y/4.
3. Hess's law: Hf(ABS) = [x·Hf(CO2,g) + (y/2)·Hf(H2O,l)] − ΔH_c (N2 taken as zero).
4. Reference Hf(CO2), Hf(H2O) are pulled directly from Cantera's own thermo data (nasa_gas.yaml / nasa_condensed.yaml), never hardcoded, so everything stays self-consistent with whatever Cantera version is actually installed.
5. Uncertainty on ΔU_c propagates straight through to hf_min/hf_max.

A sign-convention subtlety worth calling out explicitly, since the meeting notes glossed over this exact step: a bomb calorimeter reports a positive magnitude (HHV: "this much heat came out"), while enthalpies of formation are signed (negative for exothermic formation). It's very easy to combine a positive HHV number with the Hess's-law sum and get the sign wrong. The derivation is written out in full in calorimetry.py's module docstring, and tests/test_calorimetry.py cross-validates it against methane's well-known combustion/formation data (recovers Hf(CH4) ≈ −74.4 kJ/mol against the literature −74.6 kJ/mol) specifically to catch a regression here: get the sign wrong and Hf(ABS) comes out with the wrong sign and roughly double the expected magnitude.

Nitrogen (Washburn) correction: nitrogen in the sample can form nitric acid in the bomb rather than N2, releasing extra heat not attributable to the fuel's true heat of combustion. apply_nitrogen_acid_correction=True deliberately raises NotImplementedError rather than guessing a correction factor. Confirm with your calorimetry lab whether this has already been applied to their reported ΔU_c.

Model boundaries

The default lookup table spans O/F 1.1 through 40 and 100 through 1200 kPa absolute. The Hybrid solver rejects a final pressure or O/F state outside those bounds. Expand the chemistry table explicitly before modelling a higher-pressure or wider-O/F case.

Validation

• Anchor test (tests/test_anchor.py): stoichiometric methane/air reproduces T_ad ≈ 2225 K (gas-only: 2224.7 K; multiphase path with soot phase present but inactive: 2224.7 K, carbon_mass_fraction < 1e-6). Confirms the core solver path before ABS or soot are introduced at all, per the "build and validate one layer at a time" approach used throughout development.
• Shape check (tests/test_shape.py): T_ad peaks within 0.5-1.5x the stoichiometric O/F and the fuel-rich shoulder plateaus (slope shrinks approaching deep-rich) rather than spiking.
• Soot check (tests/test_soot.py): condensed carbon ~0 fuel-lean, present fuel-rich.
• Element and enthalpy conservation (tests/test_conservation.py): total C/H/N/O amounts and extensive enthalpy agree at fuel-rich, stoichiometric, and fuel-lean points. Current relative errors are near machine precision.
• Reactant-only ABS treatment (tests/test_abs_consumption.py): supported states contain no ABS product species, while a placeholder state without a complete-product enthalpy root is explicitly rejected.
• Units round-trip (tests/test_units.py): the mol/kmol factor-of-1000 conversion, cross-checked against a live Cantera Solution object.
• Calorimetry sign convention (tests/test_calorimetry.py): methane cross-check, described above.

One real bug this validation process caught during development, worth documenting: an earlier version of equilibrium_point() reconstructed a fresh ct.Mixture([(gas, 1.0), (carbon, 0.0)]) purely to "read back" the converged carbon amount after solving -- but the moles argument to the Mixture constructor resets the extensive state, so it silently read back zero carbon at every point regardless of the true solution. The element-conservation test (which checks carbon moles explicitly) and a manual fuel-rich sanity sweep both caught this; the fix threads the actual solving Mixture object through instead of rebuilding it (see equilibrium._solve_multiphase_hp's docstring for detail).

Cross-validation against NASA CEA / RPA

Structure for comparing against an independent tool: run the same (fuel formula, Hf, oxidizer, O/F, P) point through NASA CEA or RPA and compare T_ad, gamma, and c* against this tool's equilibrium_point() output for the same inputs. Paste your reference tool's numbers here once you have them. Deliberately left blank rather than inventing placeholder reference values:

Downstream lookup table (table.py)

table = table.build_table(of_grid, p_grid, hf_kJ_per_mol, oxidizer='N2O:1')
props = table.query(of=5.3, p=600e3)   # T_ad, gamma, MW, density, carbon_mass_fraction, c_star
table.save('output/combustion_table.npz')
loaded = table.CombustionTable.load('output/combustion_table.npz')

Precomputes equilibrium_point() over a (O/F, P) grid once (config: table.*), then interpolates (scipy.interpolate.RegularGridInterpolator) so the transient chamber model can query combustion properties every timestep without calling Cantera's Gibbs solver directly. Far too slow for that use case. Out-of-range queries extrapolate (so the transient model never crashes mid-timestep) but always emit a RuntimeWarning first. Never silent extrapolation.

Repo structure

abs_equilibrium/
  README.md
  config.yaml            # ABS composition, HF calorimetry inputs, oxidizer, pressure, sweep range, species list
  data/
    abs_species.yaml      # generated: exact ABS thermo used for the last run (Hf_nominal basis)
  src/
    utils.py              # config loading, mol/kmol unit-conversion helpers
    species_setup.py      # curated gas phase + custom ABS species + graphite phase
    calorimetry.py         # dU_bomb -> Hf, with uncertainty propagation
    equilibrium.py         # equilibrium_point() core + solver fallbacks
    sweep.py                # mixture-ratio sweep, HF sensitivity
    rocket.py               # c*, gamma, frozen Isp helpers
    table.py                 # lookup table + interpolator for the downstream transient model
    plotting.py
  tests/
  run.py                  # CLI entry point
  output/

CLI reference

python run.py --config config.yaml --of-min 1.1 --of-max 12 --n 90 \
              --pressure-kpa 500 --pressure-reference gauge \
              [--mechanism gri30] [--no-soot] [--area-ratio 8.0]

All flags override the corresponding config.yaml value for that run only.