diff --git a/README.md b/README.md index 281193c..2494bd9 100644 --- a/README.md +++ b/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. - -[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) diff --git a/openrocketengine/__init__.py b/openrocketengine/__init__.py index cb1a4b5..fcae48e 100644 --- a/openrocketengine/__init__.py +++ b/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", ] diff --git a/openrocketengine/engine.py b/openrocketengine/engine.py index 197438d..ec8a870 100644 --- a/openrocketengine/engine.py +++ b/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 diff --git a/openrocketengine/examples/propellant_design.py b/openrocketengine/examples/propellant_design.py new file mode 100644 index 0000000..d134f6f --- /dev/null +++ b/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() + diff --git a/openrocketengine/nozzle.py b/openrocketengine/nozzle.py index 26d250a..099fa16 100644 --- a/openrocketengine/nozzle.py +++ b/openrocketengine/nozzle.py @@ -375,10 +375,12 @@ def full_chamber_contour( y_cone = np.array([]) # Generate convergent circular arc (transition to throat) + # Arc center is at (0, Rt + R1), tangent to throat at bottom + # Arc goes from tangent point with cone (angle = theta_c) to throat (angle = 0) n_arc = num_convergent_points - len(x_cone) - theta_range = np.linspace(math.pi - theta_c, math.pi, n_arc) - x_arc = R1 * np.cos(theta_range) # Goes from negative to 0 - y_arc = Rt + R1 + R1 * np.sin(theta_range) # Connects to throat + theta_range = np.linspace(theta_c, 0, n_arc) + x_arc = -R1 * np.sin(theta_range) # Negative (upstream of throat) + y_arc = Rt + R1 * (1 - np.cos(theta_range)) # From y_tan down to Rt # Shift nozzle contour (it starts at x=0 at throat) x_nozzle = nozzle_contour.x diff --git a/openrocketengine/propellants.py b/openrocketengine/propellants.py new file mode 100644 index 0000000..7483a8f --- /dev/null +++ b/openrocketengine/propellants.py @@ -0,0 +1,475 @@ +"""Propellant thermochemistry module for OpenRocketEngine. + +This module provides combustion thermochemistry calculations using NASA CEA +via RocketCEA. It computes chamber temperature, molecular weight, gamma, +and other properties needed for rocket engine performance analysis. + +Example: + >>> from openrocketengine.propellants import get_combustion_properties + >>> props = get_combustion_properties( + ... oxidizer="LOX", + ... fuel="RP1", + ... mixture_ratio=2.7, + ... chamber_pressure_pa=7e6, + ... ) + >>> print(f"Tc = {props.chamber_temp_k:.0f} K") +""" + +from dataclasses import dataclass +from typing import Literal + +from beartype import beartype +from rocketcea.cea_obj import CEA_Obj + +# ============================================================================= +# Data Structures +# ============================================================================= + + +@beartype +@dataclass(frozen=True, slots=True) +class CombustionProperties: + """Thermochemical properties from combustion analysis. + + These properties are needed to compute rocket engine performance + using isentropic flow equations. + + Attributes: + chamber_temp_k: Adiabatic flame temperature in chamber [K] + molecular_weight: Mean molecular weight of combustion products [kg/kmol] + gamma: Ratio of specific heats (Cp/Cv) [-] + specific_heat_cp: Specific heat at constant pressure [J/(kg·K)] + characteristic_velocity: Theoretical c* [m/s] + oxidizer: Oxidizer name + fuel: Fuel name + mixture_ratio: Oxidizer-to-fuel mass ratio [-] + chamber_pressure_pa: Chamber pressure [Pa] + source: Data source ("rocketcea" or "database") + """ + + chamber_temp_k: float | int + molecular_weight: float | int + gamma: float | int + specific_heat_cp: float | int + characteristic_velocity: float | int + oxidizer: str + fuel: str + mixture_ratio: float | int + chamber_pressure_pa: float | int + source: str + + +# ============================================================================= +# Propellant Name Mapping +# ============================================================================= + +# Map common names to RocketCEA names +OXIDIZER_NAMES: dict[str, str] = { + "LOX": "LOX", + "LO2": "LOX", + "O2": "LOX", + "OXYGEN": "LOX", + "N2O4": "N2O4", + "NTO": "N2O4", + "N2O": "N2O", + "NITROUS": "N2O", + "NITROUSOXIDE": "N2O", + "H2O2": "H2O2", + "HTP": "H2O2", + "PEROXIDE": "H2O2", + "MON25": "MON25", + "MON3": "MON3", + "IRFNA": "IRFNA", + "RFNA": "IRFNA", + "CLF5": "CLF5", + "F2": "F2", + "FLUORINE": "F2", +} + +FUEL_NAMES: dict[str, str] = { + "LH2": "LH2", + "H2": "LH2", + "HYDROGEN": "LH2", + "RP1": "RP1", + "RP-1": "RP1", + "KEROSENE": "RP1", + "JET-A": "Jet-A", + "JETA": "Jet-A", + "CH4": "CH4", + "METHANE": "CH4", + "LCH4": "CH4", + "C2H5OH": "Ethanol", + "ETHANOL": "Ethanol", + "C3H8O": "IPA", + "IPA": "IPA", + "ISOPROPANOL": "IPA", + "MMH": "MMH", + "UDMH": "UDMH", + "N2H4": "N2H4", + "HYDRAZINE": "N2H4", + "A50": "A-50", + "A-50": "A-50", + "AEROZINE50": "A-50", +} + + +def _normalize_propellant_name(name: str, is_oxidizer: bool) -> str: + """Normalize propellant name to RocketCEA format.""" + normalized = name.upper().replace(" ", "").replace("-", "") + lookup = OXIDIZER_NAMES if is_oxidizer else FUEL_NAMES + + if normalized in lookup: + return lookup[normalized] + + # Try original name (RocketCEA might accept it) + return name + + +# ============================================================================= +# Fallback Database (When CEA Not Available) +# ============================================================================= + +# Tabulated data for common propellant combinations at typical conditions +# Format: (oxidizer, fuel): {MR: (Tc_K, MW, gamma, cstar_m/s)} +# Data at approximately 1000 psia (6.9 MPa) chamber pressure +# Sources: Sutton & Biblarz, various NASA reports + +_PROPELLANT_DATABASE: dict[tuple[str, str], dict[float, tuple[float, float, float, float]]] = { + ("LOX", "LH2"): { + 4.0: (3015, 12.0, 1.20, 2290), + 5.0: (3250, 13.5, 1.18, 2360), + 6.0: (3400, 14.8, 1.16, 2390), + 7.0: (3470, 16.0, 1.15, 2380), + 8.0: (3450, 17.0, 1.14, 2340), + }, + ("LOX", "RP1"): { + 2.0: (3450, 21.5, 1.21, 1750), + 2.3: (3550, 22.5, 1.19, 1780), + 2.5: (3600, 23.0, 1.18, 1790), + 2.7: (3620, 23.3, 1.17, 1800), + 3.0: (3580, 24.0, 1.16, 1780), + }, + ("LOX", "CH4"): { + 2.5: (3400, 19.5, 1.19, 1820), + 3.0: (3530, 20.5, 1.17, 1850), + 3.2: (3560, 21.0, 1.16, 1860), + 3.5: (3570, 21.5, 1.15, 1850), + 4.0: (3520, 22.5, 1.14, 1820), + }, + ("LOX", "Ethanol"): { + 1.0: (2800, 20.0, 1.24, 1650), + 1.3: (3100, 21.0, 1.22, 1720), + 1.5: (3250, 21.5, 1.20, 1750), + 1.8: (3350, 22.0, 1.19, 1760), + 2.0: (3380, 22.5, 1.18, 1750), + }, + ("N2O4", "MMH"): { + 1.5: (3000, 21.0, 1.24, 1680), + 1.8: (3150, 21.5, 1.22, 1720), + 2.0: (3220, 22.0, 1.21, 1730), + 2.2: (3260, 22.5, 1.20, 1730), + 2.5: (3250, 23.0, 1.19, 1710), + }, + ("N2O4", "UDMH"): { + 1.8: (3050, 21.5, 1.23, 1690), + 2.0: (3150, 22.0, 1.22, 1710), + 2.2: (3200, 22.5, 1.21, 1720), + 2.5: (3220, 23.0, 1.20, 1710), + 2.8: (3180, 23.5, 1.19, 1690), + }, + ("N2O4", "A-50"): { + 1.5: (3000, 21.0, 1.24, 1680), + 1.8: (3120, 21.5, 1.22, 1710), + 2.0: (3180, 22.0, 1.21, 1720), + 2.2: (3210, 22.5, 1.20, 1720), + 2.6: (3180, 23.0, 1.19, 1700), + }, + ("N2O", "Ethanol"): { + 3.0: (2800, 24.0, 1.22, 1550), + 4.0: (2950, 25.0, 1.20, 1580), + 5.0: (3000, 26.0, 1.19, 1570), + 6.0: (2980, 27.0, 1.18, 1540), + }, + ("H2O2", "RP1"): { + 6.0: (2700, 22.5, 1.21, 1580), + 7.0: (2750, 23.0, 1.20, 1590), + 7.5: (2760, 23.5, 1.19, 1580), + 8.0: (2750, 24.0, 1.19, 1570), + }, +} + + +def _interpolate_database( + oxidizer: str, fuel: str, mixture_ratio: float +) -> tuple[float, float, float, float] | None: + """Interpolate propellant database for given mixture ratio.""" + key = (oxidizer, fuel) + if key not in _PROPELLANT_DATABASE: + return None + + data = _PROPELLANT_DATABASE[key] + mrs = sorted(data.keys()) + + # Clamp to available range + if mixture_ratio <= mrs[0]: + return data[mrs[0]] + if mixture_ratio >= mrs[-1]: + return data[mrs[-1]] + + # Find bracketing values + for i in range(len(mrs) - 1): + if mrs[i] <= mixture_ratio <= mrs[i + 1]: + mr_low, mr_high = mrs[i], mrs[i + 1] + break + else: + return data[mrs[-1]] + + # Linear interpolation + t = (mixture_ratio - mr_low) / (mr_high - mr_low) + low = data[mr_low] + high = data[mr_high] + + return ( + low[0] + t * (high[0] - low[0]), # Tc + low[1] + t * (high[1] - low[1]), # MW + low[2] + t * (high[2] - low[2]), # gamma + low[3] + t * (high[3] - low[3]), # cstar + ) + + +# ============================================================================= +# RocketCEA Integration +# ============================================================================= + + +def _get_properties_from_cea( + oxidizer: str, + fuel: str, + mixture_ratio: float, + chamber_pressure_pa: float, +) -> CombustionProperties: + """Get combustion properties using RocketCEA.""" + + # Convert pressure to psia (RocketCEA default) + pc_psia = chamber_pressure_pa / 6894.76 + + # Normalize propellant names + ox_name = _normalize_propellant_name(oxidizer, is_oxidizer=True) + fuel_name = _normalize_propellant_name(fuel, is_oxidizer=False) + + # Create CEA object + cea = CEA_Obj(oxName=ox_name, fuelName=fuel_name) + + # Get chamber properties + # Note: RocketCEA returns (Mw, gamma) from get_Chamber_MolWt_gamma + Tc = cea.get_Tcomb(Pc=pc_psia, MR=mixture_ratio) # Chamber temp in R + Tc_K = Tc * 5 / 9 # Convert Rankine to Kelvin + + mw_gamma = cea.get_Chamber_MolWt_gamma(Pc=pc_psia, MR=mixture_ratio, eps=1.0) + MW = mw_gamma[0] # Molecular weight + gamma = mw_gamma[1] # Gamma + + # Get c* in ft/s, convert to m/s + cstar_fts = cea.get_Cstar(Pc=pc_psia, MR=mixture_ratio) + cstar_ms = cstar_fts * 0.3048 + + # Calculate Cp from gamma and MW + R_universal = 8314.46 # J/(kmol·K) + R_specific = R_universal / MW # J/(kg·K) + Cp = gamma * R_specific / (gamma - 1) + + return CombustionProperties( + chamber_temp_k=Tc_K, + molecular_weight=MW, + gamma=gamma, + specific_heat_cp=Cp, + characteristic_velocity=cstar_ms, + oxidizer=oxidizer, + fuel=fuel, + mixture_ratio=mixture_ratio, + chamber_pressure_pa=chamber_pressure_pa, + source="rocketcea", + ) + + +def _get_properties_from_database( + oxidizer: str, + fuel: str, + mixture_ratio: float, + chamber_pressure_pa: float, +) -> CombustionProperties: + """Get combustion properties from built-in database.""" + # Normalize names + ox_name = _normalize_propellant_name(oxidizer, is_oxidizer=True) + fuel_name = _normalize_propellant_name(fuel, is_oxidizer=False) + + result = _interpolate_database(ox_name, fuel_name, mixture_ratio) + + if result is None: + available = list(_PROPELLANT_DATABASE.keys()) + raise ValueError( + f"Propellant combination ({ox_name}, {fuel_name}) not in database. " + f"Available combinations: {available}. " + f"Install RocketCEA for arbitrary propellant combinations: pip install rocketcea" + ) + + Tc_K, MW, gamma, cstar = result + + # Calculate Cp + R_universal = 8314.46 + R_specific = R_universal / MW + Cp = gamma * R_specific / (gamma - 1) + + return CombustionProperties( + chamber_temp_k=Tc_K, + molecular_weight=MW, + gamma=gamma, + specific_heat_cp=Cp, + characteristic_velocity=cstar, + oxidizer=oxidizer, + fuel=fuel, + mixture_ratio=mixture_ratio, + chamber_pressure_pa=chamber_pressure_pa, + source="database", + ) + + +# ============================================================================= +# Public API +# ============================================================================= + + +@beartype +def get_combustion_properties( + oxidizer: str, + fuel: str, + mixture_ratio: float, + chamber_pressure_pa: float, + use_cea: bool = True, +) -> CombustionProperties: + """Get combustion thermochemistry properties for a propellant combination. + + This function returns the thermochemical properties needed for rocket engine + performance calculations. When RocketCEA is installed and use_cea=True, + it uses NASA CEA for accurate equilibrium calculations. Otherwise, it falls + back to a built-in database of common propellant combinations. + + Args: + oxidizer: Oxidizer name (e.g., "LOX", "N2O4", "N2O", "H2O2") + fuel: Fuel name (e.g., "RP1", "LH2", "CH4", "Ethanol", "MMH") + mixture_ratio: Oxidizer-to-fuel mass ratio (O/F) + chamber_pressure_pa: Chamber pressure in Pascals + use_cea: If True and RocketCEA is installed, use CEA. Otherwise use database. + + Returns: + CombustionProperties containing Tc, MW, gamma, Cp, c* + + Raises: + ValueError: If propellant combination is not available in database + and RocketCEA is not installed + + Example: + >>> props = get_combustion_properties( + ... oxidizer="LOX", + ... fuel="RP1", + ... mixture_ratio=2.7, + ... chamber_pressure_pa=7e6, + ... ) + >>> print(f"Tc = {props.chamber_temp_k:.0f} K, gamma = {props.gamma:.3f}") + """ + if use_cea: + return _get_properties_from_cea( + oxidizer, fuel, mixture_ratio, chamber_pressure_pa + ) + else: + return _get_properties_from_database( + oxidizer, fuel, mixture_ratio, chamber_pressure_pa + ) + + +@beartype +def is_cea_available() -> bool: + """Check if RocketCEA is installed and available. + + Returns: + Always True (RocketCEA is a required dependency) + """ + return True + + +@beartype +def list_database_propellants() -> list[tuple[str, str]]: + """List propellant combinations available in the built-in database. + + Returns: + List of (oxidizer, fuel) tuples available without RocketCEA + """ + return list(_PROPELLANT_DATABASE.keys()) + + +@beartype +def get_optimal_mixture_ratio( + oxidizer: str, + fuel: str, + chamber_pressure_pa: float, + expansion_ratio: float = 40.0, + metric: Literal["isp", "cstar", "density_isp"] = "isp", +) -> tuple[float, float]: + """Find the optimal mixture ratio for maximum performance. + + Searches for the mixture ratio that maximizes the specified metric. + Requires RocketCEA for accurate optimization. + + Args: + oxidizer: Oxidizer name + fuel: Fuel name + chamber_pressure_pa: Chamber pressure in Pascals + expansion_ratio: Nozzle expansion ratio for Isp calculation + metric: Optimization target: + - "isp": Maximize specific impulse + - "cstar": Maximize characteristic velocity + - "density_isp": Maximize density * Isp (important for volume-limited vehicles) + + Returns: + Tuple of (optimal_mixture_ratio, maximum_metric_value) + + Raises: + RuntimeError: If RocketCEA is not installed + """ + pc_psia = chamber_pressure_pa / 6894.76 + ox_name = _normalize_propellant_name(oxidizer, is_oxidizer=True) + fuel_name = _normalize_propellant_name(fuel, is_oxidizer=False) + + cea = CEA_Obj(oxName=ox_name, fuelName=fuel_name) + + # Search over mixture ratios + best_mr = 1.0 + best_value = 0.0 + + # Determine search range based on propellant type + if ox_name == "LOX" and fuel_name == "LH2": + mr_range = [x / 10 for x in range(30, 90, 2)] # 3.0 to 9.0 + elif ox_name == "LOX": + mr_range = [x / 10 for x in range(15, 40, 2)] # 1.5 to 4.0 + else: + mr_range = [x / 10 for x in range(10, 50, 2)] # 1.0 to 5.0 + + for mr in mr_range: + try: + if metric == "isp": + value = cea.get_Isp(Pc=pc_psia, MR=mr, eps=expansion_ratio) + elif metric == "cstar": + value = cea.get_Cstar(Pc=pc_psia, MR=mr) + elif metric == "density_isp": + isp = cea.get_Isp(Pc=pc_psia, MR=mr, eps=expansion_ratio) + # Approximate density Isp (would need propellant densities for accuracy) + value = isp # Simplified - use Isp as proxy + + if value > best_value: + best_value = value + best_mr = mr + except Exception: + continue + + return best_mr, best_value + diff --git a/pyproject.toml b/pyproject.toml index f3da031..bed610e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,6 +26,7 @@ dependencies = [ "beartype>=0.18", "numba>=0.60", "matplotlib>=3.9", + "rocketcea>=1.2.1", ] [project.optional-dependencies] @@ -34,9 +35,6 @@ dev = [ "pytest-cov>=4.0", "ruff>=0.4", ] -cea = [ - "rocketcea>=1.2", -] [project.urls] Homepage = "https://github.com/openrocketengine/openrocketengine" diff --git a/tests/test_propellants.py b/tests/test_propellants.py new file mode 100644 index 0000000..cb6fa6f --- /dev/null +++ b/tests/test_propellants.py @@ -0,0 +1,219 @@ +"""Tests for the propellants module.""" + +import pytest + +from openrocketengine.propellants import ( + CombustionProperties, + get_combustion_properties, + is_cea_available, + list_database_propellants, +) + + +class TestCombustionProperties: + """Test CombustionProperties dataclass.""" + + def test_create_properties(self) -> None: + """Test creating combustion properties.""" + props = CombustionProperties( + chamber_temp_k=3500.0, + molecular_weight=22.0, + gamma=1.2, + specific_heat_cp=2000.0, + characteristic_velocity=1800.0, + oxidizer="LOX", + fuel="RP1", + mixture_ratio=2.7, + chamber_pressure_pa=7e6, + source="test", + ) + assert props.chamber_temp_k == 3500.0 + assert props.gamma == 1.2 + assert props.source == "test" + + +class TestDatabasePropellants: + """Test propellant database functionality.""" + + def test_list_database_propellants(self) -> None: + """Test listing available propellants.""" + propellants = list_database_propellants() + + assert len(propellants) > 0 + assert ("LOX", "RP1") in propellants + assert ("LOX", "LH2") in propellants + assert ("LOX", "CH4") in propellants + + +class TestCEAIntegration: + """Test RocketCEA integration.""" + + def test_is_cea_available(self) -> None: + """Test that CEA is available (required dependency).""" + assert is_cea_available() is True + + def test_get_lox_rp1_properties(self) -> None: + """Test getting LOX/RP1 properties from CEA.""" + props = get_combustion_properties( + oxidizer="LOX", + fuel="RP1", + mixture_ratio=2.7, + chamber_pressure_pa=7e6, + ) + + # Verify reasonable values for LOX/RP1 + assert 3400 < props.chamber_temp_k < 3800 + assert 20 < props.molecular_weight < 26 + assert 1.1 < props.gamma < 1.25 + assert 1700 < props.characteristic_velocity < 1900 + assert props.source == "rocketcea" + + def test_get_lox_lh2_properties(self) -> None: + """Test getting LOX/LH2 properties from CEA.""" + props = get_combustion_properties( + oxidizer="LOX", + fuel="LH2", + mixture_ratio=6.0, + chamber_pressure_pa=10e6, + ) + + # LOX/LH2 has higher Isp, lower MW + assert 3300 < props.chamber_temp_k < 3700 + assert 12 < props.molecular_weight < 18 + assert 1.1 < props.gamma < 1.25 + assert 2300 < props.characteristic_velocity < 2500 + + def test_get_lox_ch4_properties(self) -> None: + """Test getting LOX/CH4 properties from CEA.""" + props = get_combustion_properties( + oxidizer="LOX", + fuel="CH4", + mixture_ratio=3.2, + chamber_pressure_pa=7e6, + ) + + assert 3400 < props.chamber_temp_k < 3700 + assert 18 < props.molecular_weight < 24 + assert props.source == "rocketcea" + + def test_get_n2o4_mmh_properties(self) -> None: + """Test getting N2O4/MMH (hypergolic) properties.""" + props = get_combustion_properties( + oxidizer="N2O4", + fuel="MMH", + mixture_ratio=2.0, + chamber_pressure_pa=1e6, + ) + + assert 3000 < props.chamber_temp_k < 3500 + assert 20 < props.molecular_weight < 25 + + def test_mixture_ratio_variation(self) -> None: + """Test that different mixture ratios give different results.""" + props_low = get_combustion_properties( + oxidizer="LOX", fuel="RP1", mixture_ratio=2.0, chamber_pressure_pa=7e6 + ) + props_mid = get_combustion_properties( + oxidizer="LOX", fuel="RP1", mixture_ratio=2.5, chamber_pressure_pa=7e6 + ) + props_high = get_combustion_properties( + oxidizer="LOX", fuel="RP1", mixture_ratio=3.0, chamber_pressure_pa=7e6 + ) + + # Molecular weight should increase with more oxidizer + assert props_low.molecular_weight < props_high.molecular_weight + # All should be valid + assert props_mid.chamber_temp_k > 3000 + + def test_name_normalization(self) -> None: + """Test that propellant names are normalized.""" + # These should all give similar results + props1 = get_combustion_properties("LOX", "RP1", 2.7, 7e6) + props2 = get_combustion_properties("LO2", "RP-1", 2.7, 7e6) + props3 = get_combustion_properties("OXYGEN", "KEROSENE", 2.7, 7e6) + + # All should map to the same propellant combination + assert props1.chamber_temp_k == pytest.approx(props2.chamber_temp_k, rel=0.01) + assert props1.chamber_temp_k == pytest.approx(props3.chamber_temp_k, rel=0.01) + + +class TestEngineInputsFromPropellants: + """Test EngineInputs.from_propellants() factory method.""" + + def test_from_propellants_basic(self) -> None: + """Test basic from_propellants usage.""" + from openrocketengine.engine import EngineInputs + from openrocketengine.units import kilonewtons, megapascals + + inputs = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="RP1", + thrust=kilonewtons(100), + chamber_pressure=megapascals(7), + mixture_ratio=2.7, + ) + + assert inputs.thrust.to("kN").value == pytest.approx(100) + assert inputs.chamber_pressure.to("MPa").value == pytest.approx(7) + assert inputs.mixture_ratio == pytest.approx(2.7) + # Chamber temp should be set from CEA + assert 3400 < inputs.chamber_temp.to("K").value < 3800 + assert inputs.name == "LOX/RP1 Engine" + + def test_from_propellants_with_defaults(self) -> None: + """Test from_propellants with default parameters.""" + from openrocketengine.engine import EngineInputs + from openrocketengine.units import newtons, pascals + + inputs = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="LH2", + thrust=newtons(50000), + chamber_pressure=pascals(5e6), + mixture_ratio=6.0, + ) + + # Check defaults were applied + assert inputs.exit_pressure.to("Pa").value == pytest.approx(101325) + assert inputs.lstar.to("m").value == pytest.approx(1.0) + assert inputs.contraction_ratio == pytest.approx(4.0) + assert inputs.bell_fraction == pytest.approx(0.8) + + def test_from_propellants_full_workflow(self) -> None: + """Test complete workflow from propellants to geometry.""" + from openrocketengine.engine import EngineInputs, design_engine + from openrocketengine.units import kilonewtons, megapascals + + inputs = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="CH4", + thrust=kilonewtons(50), + chamber_pressure=megapascals(5), + mixture_ratio=3.2, + name="Methane Test Engine", + ) + + # Should be able to compute performance and geometry + performance, geometry = design_engine(inputs) + + # Verify reasonable results + assert 280 < performance.isp.value < 360 # LOX/CH4 Isp range + assert performance.mdot.value > 0 + assert geometry.throat_diameter.value > 0 + assert geometry.expansion_ratio > 1 + + def test_from_propellants_custom_name(self) -> None: + """Test custom engine name.""" + from openrocketengine.engine import EngineInputs + from openrocketengine.units import megapascals, newtons + + inputs = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="RP1", + thrust=newtons(10000), + chamber_pressure=megapascals(2), + mixture_ratio=2.5, + name="My Custom Engine", + ) + + assert inputs.name == "My Custom Engine" diff --git a/uv.lock b/uv.lock index de6c608..8f3d3cf 100644 --- a/uv.lock +++ b/uv.lock @@ -555,12 +555,10 @@ dependencies = [ { name = "matplotlib" }, { name = "numba" }, { name = "numpy" }, + { name = "rocketcea" }, ] [package.optional-dependencies] -cea = [ - { name = "rocketcea" }, -] dev = [ { name = "pytest" }, { name = "pytest-cov" }, @@ -575,10 +573,10 @@ requires-dist = [ { name = "numpy", specifier = ">=2.0" }, { name = "pytest", marker = "extra == 'dev'", specifier = ">=8.0" }, { name = "pytest-cov", marker = "extra == 'dev'", specifier = ">=4.0" }, - { name = "rocketcea", marker = "extra == 'cea'", specifier = ">=1.2" }, + { name = "rocketcea", specifier = ">=1.2.1" }, { name = "ruff", marker = "extra == 'dev'", specifier = ">=0.4" }, ] -provides-extras = ["dev", "cea"] +provides-extras = ["dev"] [[package]] name = "packaging"