Add RocketCEA integration for propellant thermochemistry

README.md

@@ -1,8 +1,61 @@
-# Welcome to openrocketengine
+# OpenRocketEngine
+
 ![Python package](https://github.com/cmflannery/openrocketengine/workflows/Python%20package/badge.svg)
-[openrocketengine](https://github.com/cmflannery/openrocketengine) is an open source project designed to help with the design and development of liquid rocket engines.
 
-<!-- References -->
-[1]: http://soliton.ae.gatech.edu/people/jseitzma/classes/ae6450/bell_nozzle.pdf "GATech: Bell Nozzles"
-[2]: https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19710019929.pdf "Design of Liquid Propellant Rocket Engines"
-[3]: https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19720026079.pdf "Liquid Propellant Rocket Combustion Instability, NASA SP-194"
+Tools for liquid rocket engine design and analysis.
+
+## Installation
+
+Requires a Fortran compiler for RocketCEA (NASA CEA thermochemistry).
+
+```bash
+# macOS
+brew install gcc
+
+# Linux (Debian/Ubuntu)
+sudo apt-get install gfortran
+
+# Then install
+pip install openrocketengine
+```
+
+## Quick Start
+
+```python
+from openrocketengine import EngineInputs, design_engine, plot_engine_dashboard
+from openrocketengine.units import kilonewtons, megapascals
+
+# Design from propellant selection (thermochemistry auto-calculated)
+inputs = EngineInputs.from_propellants(
+    oxidizer="LOX",
+    fuel="RP1",
+    thrust=kilonewtons(100),
+    chamber_pressure=megapascals(7),
+    mixture_ratio=2.7,
+)
+
+# Compute performance and geometry
+performance, geometry = design_engine(inputs)
+
+print(f"Isp (sea level): {performance.isp}")
+print(f"Isp (vacuum): {performance.isp_vac}")
+print(f"Throat diameter: {geometry.throat_diameter}")
+
+# Visualize
+plot_engine_dashboard(inputs, performance, geometry)
+```
+
+## Features
+
+- **Type-safe**: Runtime type checking with beartype
+- **Units handling**: Built-in `Quantity` class prevents unit errors
+- **Fast**: Numba-accelerated isentropic flow calculations
+- **Visualization**: Engine cross-sections, performance curves, dashboards
+- **NASA CEA**: Accurate thermochemistry via RocketCEA
+- **Nozzle contours**: Rao bell and conical nozzle generation with CSV export
+
+## References
+
+1. [GATech: Bell Nozzles](http://soliton.ae.gatech.edu/people/jseitzma/classes/ae6450/bell_nozzle.pdf)
+2. [Design of Liquid Propellant Rocket Engines](https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19710019929.pdf)
+3. [Liquid Propellant Rocket Combustion Instability, NASA SP-194](https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19720026079.pdf)

openrocketengine/__init__.py

@@ -54,6 +54,15 @@
     plot_performance_vs_altitude,
 )
 
+# Propellants and thermochemistry
+from openrocketengine.propellants import (
+    CombustionProperties,
+    get_combustion_properties,
+    get_optimal_mixture_ratio,
+    is_cea_available,
+    list_database_propellants,
+)
+
 __all__ = [
     # Version
     "__version__",
@@ -80,4 +89,10 @@
     "plot_nozzle_contour",
     "plot_performance_vs_altitude",
     "plot_engine_dashboard",
+    # Propellants
+    "CombustionProperties",
+    "get_combustion_properties",
+    "get_optimal_mixture_ratio",
+    "is_cea_available",
+    "list_database_propellants",
 ]

openrocketengine/engine.py

@@ -33,9 +33,11 @@
 )
 from openrocketengine.units import (
     Quantity,
+    kelvin,
     kg_per_second,
     meters,
     meters_per_second,
+    pascals,
     seconds,
     square_meters,
 )
@@ -131,6 +133,108 @@ def effective_ambient_pressure(self) -> Quantity:
             return self.ambient_pressure
         return self.exit_pressure
 
+    @classmethod
+    def from_propellants(
+        cls,
+        oxidizer: str,
+        fuel: str,
+        thrust: Quantity,
+        chamber_pressure: Quantity,
+        mixture_ratio: float | None = None,
+        exit_pressure: Quantity | None = None,
+        lstar: Quantity | None = None,
+        ambient_pressure: Quantity | None = None,
+        contraction_ratio: float = 4.0,
+        contraction_angle: float = 45.0,
+        bell_fraction: float = 0.8,
+        name: str | None = None,
+    ) -> "EngineInputs":
+        """Create EngineInputs from propellant names, automatically computing thermochemistry.
+
+        This factory method uses RocketCEA (NASA CEA) to determine the combustion
+        properties (chamber temperature, molecular weight, and gamma) from the
+        specified propellant combination.
+
+        Args:
+            oxidizer: Oxidizer name (e.g., "LOX", "N2O4", "N2O", "H2O2")
+            fuel: Fuel name (e.g., "RP1", "LH2", "CH4", "Ethanol", "MMH")
+            thrust: Sea-level thrust
+            chamber_pressure: Chamber pressure
+            mixture_ratio: O/F mass ratio. If None, uses optimal ratio for max Isp.
+            exit_pressure: Nozzle exit pressure. Defaults to 1 atm (101325 Pa).
+            lstar: Characteristic length. Defaults to 1.0 m (typical for biprop).
+            ambient_pressure: Ambient pressure for performance calc. Defaults to exit_pressure.
+            contraction_ratio: Chamber/throat area ratio. Default 4.0.
+            contraction_angle: Convergent section half-angle [deg]. Default 45.
+            bell_fraction: Bell length as fraction of 15° cone. Default 0.8.
+            name: Optional engine name.
+
+        Returns:
+            EngineInputs with thermochemistry computed from propellant combination.
+
+        Example:
+            >>> inputs = EngineInputs.from_propellants(
+            ...     oxidizer="LOX",
+            ...     fuel="RP1",
+            ...     thrust=kilonewtons(100),
+            ...     chamber_pressure=megapascals(7),
+            ...     mixture_ratio=2.7,
+            ... )
+            >>> print(f"Tc = {inputs.chamber_temp}")
+        """
+        from openrocketengine.propellants import (
+            get_combustion_properties,
+            get_optimal_mixture_ratio,
+        )
+
+        # Default exit pressure to 1 atm
+        if exit_pressure is None:
+            exit_pressure = pascals(101325)
+
+        # Default L* to 1.0 m
+        if lstar is None:
+            lstar = meters(1.0)
+
+        # Get chamber pressure in Pa for CEA
+        pc_pa = chamber_pressure.to("Pa").value
+
+        # Find optimal mixture ratio if not specified
+        if mixture_ratio is None:
+            mixture_ratio, _ = get_optimal_mixture_ratio(
+                oxidizer=oxidizer,
+                fuel=fuel,
+                chamber_pressure_pa=pc_pa,
+                metric="isp",
+            )
+
+        # Get combustion properties from CEA
+        props = get_combustion_properties(
+            oxidizer=oxidizer,
+            fuel=fuel,
+            mixture_ratio=mixture_ratio,
+            chamber_pressure_pa=pc_pa,
+        )
+
+        # Generate name if not provided
+        if name is None:
+            name = f"{oxidizer}/{fuel} Engine"
+
+        return cls(
+            thrust=thrust,
+            chamber_pressure=chamber_pressure,
+            chamber_temp=kelvin(props.chamber_temp_k),
+            exit_pressure=exit_pressure,
+            molecular_weight=props.molecular_weight,
+            gamma=props.gamma,
+            lstar=lstar,
+            mixture_ratio=mixture_ratio,
+            ambient_pressure=ambient_pressure,
+            contraction_ratio=contraction_ratio,
+            contraction_angle=contraction_angle,
+            bell_fraction=bell_fraction,
+            name=name,
+        )
+
 
 # =============================================================================
 # Output Data Structures

openrocketengine/examples/propellant_design.py

@@ -0,0 +1,173 @@
+#!/usr/bin/env python
+"""Propellant-based engine design example for OpenRocketEngine.
+
+This example demonstrates the simplified workflow where you specify
+propellants and the library automatically determines combustion properties.
+
+No need to manually look up Tc, gamma, or molecular weight!
+"""
+
+from openrocketengine import design_engine, plot_engine_dashboard
+from openrocketengine.engine import EngineInputs
+from openrocketengine.nozzle import full_chamber_contour, generate_nozzle_from_geometry
+from openrocketengine.units import kilonewtons, megapascals
+
+
+def main() -> None:
+    """Run the propellant-based design example."""
+    print("=" * 70)
+    print("OpenRocketEngine - Propellant-Based Design")
+    print("=" * 70)
+    print()
+    print("Using NASA CEA (via RocketCEA) for thermochemistry calculations")
+    print()
+
+    # =========================================================================
+    # Design a LOX/RP-1 Engine (like Merlin)
+    # =========================================================================
+
+    print("-" * 70)
+    print("Design 1: LOX/RP-1 Engine (Kerolox)")
+    print("-" * 70)
+    print()
+
+    # Just specify propellants, thrust, and pressure - that's it!
+    lox_rp1 = EngineInputs.from_propellants(
+        oxidizer="LOX",
+        fuel="RP1",
+        thrust=kilonewtons(100),  # 100 kN
+        chamber_pressure=megapascals(7),  # 7 MPa (~1000 psi)
+        mixture_ratio=2.7,  # Typical for LOX/RP-1
+        name="Kerolox-100",
+    )
+
+    print(f"Engine: {lox_rp1.name}")
+    print("  Propellants: LOX / RP-1")
+    print(f"  Mixture Ratio: {lox_rp1.mixture_ratio}")
+    print(f"  Chamber Temp: {lox_rp1.chamber_temp.to('K').value:.0f} K (auto-calculated!)")
+    print(f"  Gamma: {lox_rp1.gamma:.3f} (auto-calculated!)")
+    print(f"  Molecular Weight: {lox_rp1.molecular_weight:.1f} kg/kmol (auto-calculated!)")
+    print()
+
+    perf1, geom1 = design_engine(lox_rp1)
+    print("Performance:")
+    print(f"  Isp (SL): {perf1.isp.value:.1f} s")
+    print(f"  Isp (Vac): {perf1.isp_vac.value:.1f} s")
+    print(f"  Thrust Coeff: {perf1.thrust_coeff:.3f}")
+    print(f"  Mass Flow: {perf1.mdot.value:.2f} kg/s")
+    print()
+    print("Geometry:")
+    print(f"  Throat Diameter: {geom1.throat_diameter.to('m').value * 100:.1f} cm")
+    print(f"  Exit Diameter: {geom1.exit_diameter.to('m').value * 100:.1f} cm")
+    print(f"  Expansion Ratio: {geom1.expansion_ratio:.1f}")
+    print()
+
+    # =========================================================================
+    # Design a LOX/Methane Engine (like Raptor)
+    # =========================================================================
+
+    print("-" * 70)
+    print("Design 2: LOX/Methane Engine (Methalox)")
+    print("-" * 70)
+    print()
+
+    lox_ch4 = EngineInputs.from_propellants(
+        oxidizer="LOX",
+        fuel="CH4",
+        thrust=kilonewtons(200),
+        chamber_pressure=megapascals(10),  # Higher pressure
+        mixture_ratio=3.2,
+        name="Methalox-200",
+    )
+
+    print(f"Engine: {lox_ch4.name}")
+    print(f"  Chamber Temp: {lox_ch4.chamber_temp.to('K').value:.0f} K")
+    print(f"  Gamma: {lox_ch4.gamma:.3f}")
+    print()
+
+    perf2, geom2 = design_engine(lox_ch4)
+    print("Performance:")
+    print(f"  Isp (SL): {perf2.isp.value:.1f} s")
+    print(f"  Isp (Vac): {perf2.isp_vac.value:.1f} s")
+    print()
+
+    # =========================================================================
+    # Design a LOX/LH2 Engine (like RS-25/SSME)
+    # =========================================================================
+
+    print("-" * 70)
+    print("Design 3: LOX/LH2 Engine (Hydrolox)")
+    print("-" * 70)
+    print()
+
+    lox_lh2 = EngineInputs.from_propellants(
+        oxidizer="LOX",
+        fuel="LH2",
+        thrust=kilonewtons(50),  # Smaller for demo
+        chamber_pressure=megapascals(15),  # High pressure like SSME
+        mixture_ratio=6.0,  # Typical for LOX/LH2
+        name="Hydrolox-50",
+    )
+
+    print(f"Engine: {lox_lh2.name}")
+    print(f"  Chamber Temp: {lox_lh2.chamber_temp.to('K').value:.0f} K")
+    print(f"  Molecular Weight: {lox_lh2.molecular_weight:.1f} kg/kmol (low = high Isp!)")
+    print()
+
+    perf3, geom3 = design_engine(lox_lh2)
+    print("Performance:")
+    print(f"  Isp (SL): {perf3.isp.value:.1f} s")
+    print(f"  Isp (Vac): {perf3.isp_vac.value:.1f} s  <- Highest!")
+    print()
+
+    # =========================================================================
+    # Comparison Summary
+    # =========================================================================
+
+    print("=" * 70)
+    print("COMPARISON SUMMARY")
+    print("=" * 70)
+    print()
+    print(f"{'Engine':<20} {'Isp(SL)':<10} {'Isp(Vac)':<10} {'MW':<8} {'Tc (K)':<10}")
+    print("-" * 70)
+
+    for name, inputs, perf in [
+        ("LOX/RP-1", lox_rp1, perf1),
+        ("LOX/CH4", lox_ch4, perf2),
+        ("LOX/LH2", lox_lh2, perf3),
+    ]:
+        print(
+            f"{name:<20} "
+            f"{perf.isp.value:<10.1f} "
+            f"{perf.isp_vac.value:<10.1f} "
+            f"{inputs.molecular_weight:<8.1f} "
+            f"{inputs.chamber_temp.to('K').value:<10.0f}"
+        )
+
+    print()
+    print("Note: Lower molecular weight (MW) = higher Isp")
+    print("      LH2 engines have best Isp but require large tanks (low density)")
+    print()
+
+    # =========================================================================
+    # Generate Dashboard for LOX/RP-1 Engine
+    # =========================================================================
+
+    print("Generating visualization for LOX/RP-1 engine...")
+
+    nozzle = generate_nozzle_from_geometry(geom1)
+    contour = full_chamber_contour(lox_rp1, geom1, nozzle)
+
+    fig = plot_engine_dashboard(lox_rp1, perf1, geom1, contour)
+    fig.savefig("kerolox_engine_dashboard.png", dpi=150, bbox_inches="tight")
+    print("  Saved: kerolox_engine_dashboard.png")
+
+    print()
+    print("=" * 70)
+    print("Done!")
+    print("=" * 70)
+
+
+if __name__ == "__main__":
+    main()
+