hybrid_engine_sim

A transient single-point hybrid rocket combustion-chamber simulation: couples the oxidizer tank, injector, solid-fuel grain regression, and nozzle into one ODE system integrated from ignition to oxidizer depletion, sweeping nozzle throat/exit geometry to find the largest total impulse.

This is a ground-up rewrite of githubcode/hybrid-rocket-system/src/: the "already-written transient single-point combustion chamber model" the original ABS equilibrium capstone brief referenced as the downstream consumer of abs_equilibrium's lookup table. That connection is now real: this project queries a precomputed abs_equilibrium combustion table instead of calling Cantera's gas.equilibrate('HP') live inside the transient loop.

The original githubcode/hybrid-rocket-system/ repo is untouched. This lives in its own sibling folder and doesn't modify or depend on anything in that repo.

Quick start

pip install cantera numpy scipy matplotlib pandas pyyaml CoolProp shapely

1. Build an abs_equilibrium combustion table first (if you haven't already): cd ../abs_equilibrium && python run.py. See that project's README.
2. python run.py. Runs a single burn at the configured throat/exit radii.
3. python run.py --sweep-nozzle. Also sweeps nozzle geometry (range set in config.yaml → nozzle.sweep) and reports the best (throat, exit) pair.

Outputs land in output/: regression_rate.png, thrust_curve.png, chamber_history.png, burn_history.csv, run_summary.json, and grain_geometry.png. The CLI keeps the geometry diagnostic, while the final web results omit it because the Geometry Designer already shows that plot.

Run tests with pytest tests/ -v (some tests skip automatically if the abs_equilibrium table hasn't been built yet).

What's kept, and what changed

Dropped entirely

• Excel_reader.py: unused Excel-column reader + a throwaway tank-fill-time calc (flagged by you as "a complete waste of time"; confirmed unused anywhere else).
• Test_Field.py: scratchpad duplicating Non_Equilibrium_Tank_Model_ZK.py's heat-transfer functions, plus unused line-intersection geometry helpers.
• "Heat Transfer FEA Model.py": disconnected wall-heat-transfer scratch work, never wired to anything.
• Regression_sample.py: early draft superseded by the ZK tank model.
• Tank_model_solver.py: confirmed broken: calls the ZK tank-physics functions with mismatched argument signatures (wrong count/order vs. how they're actually defined), and its own plotting section indexes an 11-element state array at positions 11 and 12 (sol['y'][11], [12]). Would raise IndexError if execution ever reached that point.
• Non_Equilibrium_Tank_Model_ZK.py, Tank_solver_LITE.py: a more physically detailed two-phase (separate liquid/vapor) tank model. Real effort, more realistic than the lumped model, but per your call: left behind, not ported. Tank_solver_LITE.py calls the ZK physics with matching signatures (unlike Tank_model_solver.py) but sets rtol=atol=1e99 on its solver, effectively disabling numerical error control. If you want this ported later, that'd be the starting point: fix the tolerances, keep the physics.

The abs_equilibrium connection

combustion.py loads a precomputed CombustionTable directly from abs_equilibrium/src/table.py (the .npz file built by abs_equilibrium/run.py). Inside the regression-rate root-find (still scipy.optimize.brentq, solving for the fuel mass flow rate that makes the Eilers correlation self-consistent), the flame temperature/gamma/molecular-weight needed each iteration comes from table.query(of, p) instead of a fresh Cantera Gibbs-energy minimization. The exact swap the original capstone brief asked for.

Implementation note on module loading: both projects independently have a sweep.py and a plotting.py. An early version of combustion.py used sys.path.insert(0, ...) to make abs_equilibrium's src/ importable, which made Python's import system resolve import sweep to abs_equilibrium's sweep module instead of this project's own. A very confusing AttributeError deep inside an unrelated function. Fixed by sys.path.append (not insert) instead, so this project's own same-named modules (already on sys.path at position 0 via run.py) always win the lookup first.

The chamber is now quasi-steady: a real fidelity trade-off

This is the one design decision in this rewrite that changes the simulation's physics, not just its code quality, so it's worth being explicit about.

The original tracked the chamber's gas mass and internal energy as ODE states (plus one state per chemical species) and inverted them back to temperature/pressure via Cantera's full equation of state at every step. That only works with a live gas object that knows the species thermodynamics; which is exactly what querying a precomputed table instead of a live equilibrium solver gives up.

Instead, engine_sim.py treats the chamber as quasi-steady: at every timestep, chamber pressure is whatever value makes the throat's choked mass-flow capacity equal the current total propellant mass flow rate (mdot_ox + mdot_fuel); a root-find over pressure, with the regression-rate root-find nested inside it (see engine_sim.chamber_state_at). This means:

• The ODE state vector shrank from [radius, mc, Uc, mo, Uo, *mY_species] (5 + N species) to just [radius, mo, Uo]. There's no in-chamber composition to accumulate anymore.
• Chamber-filling / ignition-transient pressure dynamics are gone. The original had some of this (imperfectly) through its chamber energy balance; this rewrite assumes the chamber's residence time is short enough that it's always close to the equilibrium state for the current instantaneous O/F, not modeling a lag.
• This is a standard, well-understood simplification in hybrid-rocket internal ballistics education/modeling when full in-chamber composition/energy tracking isn't the focus. But it is a real trade-off, not a pure cleanup, and worth revisiting if ignition-transient chamber pressure buildup matters for your application.
• The upside: every ODE evaluation is now two cheap nested 1-D root-finds against table lookups, not a live Cantera Gibbs-energy minimization; which is the entire point of the exercise.

Test-rig pressure and geometry inputs

All exposed pressure values use kPa. target_peak_chamber_pressure_kPa is interpreted using pressure_reference: gauge or absolute; gauge input uses ambient_pressure_kPa_absolute for the internal absolute state. The default is the planned 500 kPa gauge rig condition. The fixed nozzle throat radius remains 8 mm, and injector sizing runs complete transient burns until peak pressure is within 1 percent of the selected target.

The default grain is 100 mm long with a 31.82 mm outer diameter and a 4.5 mm circular bore radius. A selected Geometry Designer .npz stores both dimensions and the remaining-fuel quench margin as authoritative metadata. A dimension mismatch is rejected before simulation, and older geometry files must be upgraded in the Designer. The quench margin is shown in millimetres and defaults to 1 mm. Geometry diagnostic plots show regression in mm, areas in mm², and remaining fuel volume in cm³ while the solver remains in SI units.

The Bambu Black density default is 1050 kg/m³ from the Bambu ABS Technical Data Sheet.

Provisional regression inputs

The Eilers regression-rate correlation needs the Prandtl number and dynamic viscosity of the combustion products; neither of which abs_equilibrium's table carries (its curated NASA-thermo species set has no transport data at all; see abs_equilibrium/README.md, "Why not GRI-Mech 3.0?"). Only a full kinetic mechanism like GRI-Mech 3.0 ships transport polynomials, which is exactly what that project deliberately avoids for the equilibrium calculation itself. Rather than reworking abs_equilibrium's species database to add transport data (a much bigger change to an already-finished, tested project), this project uses:

• propellant.prandtl_products_placeholder (config.yaml, default 0.7): a standard approximation for hot combustion gases used in several published hybrid-rocket regression-rate studies, not computed from the actual gas state.
• propellant.viscosity_products_placeholder_pa_s (config.yaml, default 9e-5 Pa s): a retained numerical placeholder for the dynamic viscosity of the combustion-products mixture used in the boundary-layer Reynolds number.

Both are clearly marked in config.yaml. Replace with computed/measured values if you need better regression-rate accuracy.

Resolved inherited finocyl artifact

The former ripple near 6 mm regression radius came from a discontinuous legacy formula switch followed by a degree-120 Chebyshev fit. The built-in finocyl now uses the same uniform normal-offset geometry solver as custom ports. Its port area and fuel volume are monotonic, its burn surface is continuous through the old switch, and full burnout means no fuel cross-section remains.

Verification pass: what was checked, and what it found

The geometry code was put through a dedicated verification pass. Every check was made against something independent of the code under test: a closed-form formula, a physical conservation law, a second implementation, or published data. Because re-running the same function and getting the same answer proves nothing. The checks now live in tests/test_grain_shapes.py and tests/test_custom_geometry_integration.py.

Two real bugs were found and fixed.

1. Burning surface area included the motor casing (grain_shapes.py). Surface area was taken as the perimeter of the port after clipping it to the chamber. Once the port grew out to the casing, that perimeter ran along the chamber wall. So inert steel was being counted as burning fuel. It inflated burning surface by 13-42% depending on shape and broke the dV/dr = -S conservation identity by the same amount. Now measured as the burning front itself restricted to inside the chamber, which cannot include the casing by construction. Conservation error dropped from 13-42% to under 0.03%.

Impact on results you may already have: with the default 4.63 L tank this bug was dormant; the burn is tank-limited, so the oxidiser runs out before the port ever reaches the casing (margin was 0.6-5.8 mm depending on shape). Numbers produced with the default config were unaffected. On a grain-limited motor (e.g. 4x tank volume) the same bug inflated total impulse by 1.6-4.3%.

1. Straight-bore port cross-section grew linearly instead of quadratically (grain_geometry.py, inherited from the original project). Port area was built with np.linspace between the start and end areas, ramping area linearly with regression depth. But a circle's area goes as r². It was exact at both ends and overstated by up to 18.8% mid-burn. Port area sets the oxidiser mass flux, and the Eilers correlation scales as G^0.8, so this fed through to roughly a 15% error in regression rate for any straight-bore run. The burning-surface curve was fine (2πRL genuinely is linear in R), which is likely why it went unnoticed.

Also fixed: burnout radius was quantised to a scan grid (0.2-0.9% overshoot, now 0.003% via bisection); geometry curves raised ValueError instead of clamping when the ODE solver stepped a hair past burnout; v_fun meant port volume in one module and remaining fuel volume in the other while sharing a plot label. Now both mean remaining fuel; nothing in the simulation reads v_fun, so this affects only the geometry plot. Uploaded .npz geometries were accepted with NaNs, negative areas or non-monotonic data and silently poisoned the whole burn (now validated and rejected with a clear message); and the webapp's preview plots forced a GUI matplotlib backend that would fail on a headless server.

Checked and found correct (no change needed): nozzle critical pressure ratio and choked mass flow against textbook closed forms (exact); the area-Mach solver against an independent implementation in abs_equilibrium (agree to 2e-9) and the textbook Me = 2.94 at Ae/At = 4, γ = 1.4; c* against P_c·A_t/ṁ (exact); N2O/ABS chemistry against published hybrid figures (optimum O/F 6.8, peak T_ad 3093 K, c 1450 m/s, γ 1.24-1.37); lookup-table interpolation against live Cantera at off-grid pressures (≤0.17%, so n_p: 4 is adequate); Eilers regression-rate magnitude (3.1 mm/s) and its mdot^0.8 / L^-0.2 exponent scaling (exact); CoolProp N2O saturation pressure against published data (≤0.54%); and the Dyer injector blend, whose plain 50/50 average is exactly* the published κ-weighted form because κ ≡ 1 for a saturated tank.

The chamber mass-balance solve no longer substitutes a fixed pressure when it cannot bracket a solution. It expands only within the combustion table's absolute-pressure bounds and fails clearly if no physical root exists there. Final O/F states outside the table are also rejected instead of clamped.

Sanity-checked magnitudes (default config)

A production run with the 100 mm grain, 31.82 mm outer diameter, 4.5 mm circular bore radius, Bambu Black density, 8 mm throat radius, and 500 kPa gauge target sized the injector to 2.476039 mm². It reached 500.000 kPa gauge and stopped at the 1 mm quench margin after 8.288 s. Chamber pressure fell to 401.502 kPa gauge as the tank cooled, so thrust fell from 137.753 N to 110.177 N. The run produced 978.496 N·s total impulse and 4.565612 mm/s peak regression rate. Reynolds number crossed below the Eilers cited range after about 0.160 s, leaving 1.93 percent of the active burn inside the cited range. Red regression-plot segments are therefore model extrapolations, not numerical solver failures.

Repo structure

hybrid_engine_sim/
  README.md
  config.yaml
  src/
    grain_geometry.py     # built-in finocyl / straight-bore + custom-geometry loader
    grain_shapes.py        # generic polygon-offset engine: star, wagon wheel, custom outlines
    mass_flux.py
    oxidizer_tank.py
    combustion.py       # <- the abs_equilibrium integration point
    nozzle.py
    engine_sim.py        # ODE system + quasi-steady chamber closure
    sweep.py               # single-burn runner + nozzle geometry sweep
    plotting.py
  tests/
    test_combustion.py
    test_engine_sim.py
    test_grain_shapes.py                  # geometry engine vs closed forms + conservation law
    test_custom_geometry_integration.py    # custom geometry actually drives the burn
  run.py                  # CLI entry point
  output/

CLI reference

python run.py [--config config.yaml] [--sweep-nozzle] [--t-final 20.0]

Website integration is deliberately not built in this round. The CLI + config.yaml pattern here is the same one abs_equilibrium's webapp already wraps, so adding a page for this later follows the exact same recipe once you're ready.