Added python implementation MVP

HRAP - Python/hrap_python/.gitignore

@@ -0,0 +1,2 @@
+__pycache__
+**/__pycache__

HRAP - Python/hrap_python/README.md

@@ -0,0 +1,65 @@
+# HRAP Python
+
+## Purpose
+This package provides a Python port of the HRAP utility, allowing users to simulate hybrid rocket engines without the need for a MATLAB license. The package exposes a core API which can be used to generate custom scripts, as well as a GUI for ease of use.
+
+## Getting Started
+To run a simulation, you'll first need to generate a `SimulationConfig` and `InitialConditions`, which are both pretty self-explanatory. A `SimulationConfig` is composed of several child components corresponding to the different parts of a hybrid engine: `PropellantConfig`, `HardwareConfig`, and `NozzleConfig`. `PropellantConfig` requires the creation of a `PropellantData` object, which loads tabulated propellant data from one of the `.npz` files in the `propellants` folder. All told, that looks like: 
+
+```python
+from api.core import sim_loop
+from api.models import *
+from api.sim import shift_OF
+
+material = MaterialData('propellants/Paraffin.npz')
+prop_config = PropellantConfig(ID=LengthValue(2.0, LengthUnit.INCHES),
+                               OD=LengthValue(4.0, LengthUnit.INCHES),
+                               length=LengthValue(12.0, LengthUnit.INCHES),
+                               cstar_eff=1.0,
+                               material=material,
+                               regression_model=shift_OF)
+
+hardware_config = HardwareConfig(tank_volume=VolumeValue(7000.0, VolumeUnit.CU_CENTIMETERS),
+                                 chamber_volume=None,
+                                 injector_Cd=0.6,
+                                 injector_D=LengthValue(0.375, LengthUnit.INCHES),
+                                 injector_N=1,
+                                 vent_state=VentConfig.EXTERNAL,
+                                 vent_Cd=0.6,
+                                 vent_D=LengthValue(0.028, LengthUnit.INCHES))
+
+nozzle_config = NozzleConfig(Cd=1.0,
+                            throat=LengthValue(1.2, LengthUnit.INCHES),
+                            exp_ratio=4.0,
+                            efficiency=1.0)
+
+sim_config = SimulationConfig(timestep=0.001, 
+                              max_run_time=10.0,
+                              max_burn_time=10.0,
+                              ambient_pressure=PressureValue(1.0, PressureUnit.ATM),
+                              hardware=hardware_config,
+                              prop=prop_config,
+                              nozzle=nozzle_config)
+```
+
+Initial conditions are generated in much the same way:
+```python
+init_conditions = InitialConditions(chamber_pressure=PressureValue(1.0, PressureUnit.ATM),
+                                    mass_ox=MassValue(7.0, MassUnit.KILOGRAM),
+                                    tank_temp=TemperatureValue(63.0, TemperatureUnit.FAHRENHEIT))
+```
+
+Finally, the simulation can be run by providing the constructed `SimulationConfig` and `InitialConditions` objects to `sim_loop`, which returns a `SimulationOutput` object containing all relevant simulation data.
+```python
+output = sim_loop(sim_config, init_conditions)
+```
+
+Outputs, much like inputs, will be given as a `*Value`. To convert between a list of `*Value`s to a list of conventional floats, simply use `get_as` in a list comprehension like so:
+```python
+thrust_data = [t.get_as(MassUnit.NEWTONS) for t in output.thrust]
+```
+
+The full code, along with a plotting example, is available in `demo.py`
+
+## Testing
+Unit tests are provided to ensure that behavior is well-defined and doesn't change version-to-version. Tests are located in the `test` module, and can be run from this directory with `python3 -m unittest discover -s test -v`. Some 'magic number' values in the unit tests are copied from the MATLAB output of the program for reference.

HRAP - Python/hrap_python/api/core.py

@@ -0,0 +1,83 @@
+"""
+HRAP Simulation Core (`hrap_python.api.core`)
+
+Core HRAP simulation functions
+"""
+from .models import *
+from .sim import *
+
+def sim_loop(config: SimulationConfig, init_conditions: InitialConditions) -> SimulationOutput:
+    """
+        Runs a simulation.
+
+        Parameters
+        ----------
+        config : SimulationConfig
+            The simulation configuration parameters
+
+        init_conditions : InitialConditions
+            The initial conditions for the simulation state
+
+        Returns
+        -------
+        SimulationOutput
+    """
+    t = 0.0
+    i = 1
+    dt = config.timestep
+
+    state = SimulationState(config, init_conditions)
+    output = SimulationOutput(config)
+
+    while True:
+        t = (i-1)*dt
+        i = i+1
+
+        sim_iteration(config, state, output)
+
+        if state.chamber.grain_ID >= config.prop.OD:
+            output.end_condition = 'Fuel Depleted'
+            break
+        elif state.tank.mass_ox <= 0:
+            output.end_condition = 'Oxidizer Depleted'
+            break
+        elif t >= config.max_run_time:
+            output.end_condition = 'Max Simulation Time Reached'
+            break
+        elif state.chamber.pressure <= config.ambient_pressure:
+            output.end_condition = 'Burn Complete'
+            break
+
+    return output
+
+def sim_iteration(config: SimulationConfig, state: SimulationState, output: SimulationOutput) -> None:
+    """
+        A single iteration of a simulation. Contains all simulation state transitions for a given timestep.
+
+        Parameters
+        ----------
+        config : SimulationConfig
+            The simulation configuration parameters
+
+        state : SimulationState
+            The current simulation state
+
+        output : SimulationOutput
+            The tentative simulation output
+
+        Returns
+        -------
+        None
+    """
+    tank(config, state)
+    config.prop.regression_model(config, state)
+
+    combustion(config, state)
+    chamber(config, state)          
+    nozzle(config, state)
+
+    if config.mass_properties is not None:
+        mass(config, state)
+
+    output.append_state(state)
+    state.time += config.timestep

HRAP - Python/hrap_python/api/models.py

@@ -0,0 +1,242 @@
+"""
+HRAP Data Models (`hrap_python.api.models`)
+
+Configuration and state models for HRAP simulations
+"""
+from enum import Enum
+from util.units import *
+from util.nox import *
+from scipy.optimize import fsolve
+from typing import List
+import numpy as np
+
+###
+### CONFIGURATION MODELS
+###
+class VentConfig(Enum):
+    NONE = 1
+    EXTERNAL = 2
+    INTERNAL = 3
+
+class MaterialData:
+    # TODO: Better way of loading file?
+    def __init__(self, filename: str):
+        container = np.load(filename) 
+        for label in ['metadata', 'regression_coeff', 'OF', 'Pc', 'k', 'M', 'T']:
+            assert label in container, f'Data for "{label}" not found in saved combustion data file {filename}!'
+        # Metadata can be extended to include other mixed data - although it's unclear whether NoneTypes are supported
+        self.name   = container['metadata'][0]
+        self.rho    = float(container['metadata'][1])
+        # All coefficients will be negative (typically -1) if a shifting-OF model is not supported
+        # 0 is not reliable due to floating-point precision errors (may become positive or negative)
+        self.regression_coeff = container['regression_coeff']
+        self.OF     = container['OF']
+        self.Pc     = container['Pc']
+        self.k      = container['k']
+        self.M      = container['M']
+        self.T      = container['T']
+
+class PropellantConfig:
+    def __init__(self, ID: LengthValue, OD: LengthValue, length: LengthValue, cstar_eff: float, 
+                 material: MaterialData, regression_model, const_OF_ratio: float | None = None):
+        self.ID     = ID.base_value         # Grain ID (m)
+        self.OD     = OD.base_value         # Grain OD (m)
+        self.length = length.base_value     # Grain Length (m)
+
+        self.cstar_eff  = cstar_eff # C-star efficiency
+        self.material   = material  # Material Data - imported from file
+
+        assert ((material.regression_coeff[0] > 0 and const_OF_ratio is None)
+                or (const_OF_ratio is not None and const_OF_ratio > 0)), "Material must support shifting OF ratio, or constant OF ratio must be provided!"
+
+        self.regression_model = regression_model
+        self.const_OF_ratio = const_OF_ratio # Constant OF ratio (optional)
+
+class HardwareConfig:
+    def __init__(self, tank_volume: VolumeValue, chamber_volume: VolumeValue | None, injector_Cd: float, injector_D: LengthValue,
+                 injector_N: int, vent_state: VentConfig, vent_Cd: float | None = None, vent_D: LengthValue | None = None):
+        self.tank_volume = tank_volume.base_value       # Tank Volume (m^3)
+        if chamber_volume is None:
+            self.chamber_volume = None
+        else:
+            self.chamber_volume = chamber_volume.base_value # Chamber Volume (m^3) - if None, computed from grain dimensions
+
+        self.injector_cda_N = injector_Cd * d_to_a(injector_D.base_value) * injector_N
+
+        assert vent_state == VentConfig.NONE or (vent_Cd is not None and vent_D is not None), "Must input valid vent configuration!"
+        self.vent_state = vent_state
+        self.vent_cda = vent_Cd * d_to_a(vent_D.base_value)
+
+class NozzleConfig:
+    def __init__(self, Cd: float, throat: LengthValue, exp_ratio: float, efficiency: float):
+        self.Cd = Cd                    # Coefficient of dispersion
+        self.throat = throat.base_value # Throat diameter (m)
+        self.exp_ratio = exp_ratio      # Expansion Ratio (unitless)
+        self.efficiency = efficiency    # Efficiency (>0, <1.0)
+
+class MassProperties:
+    def __init__(self, motor_mass: MassValue, motor_cg: LengthValue, tank_location: LengthValue, grain_location: LengthValue):
+        self.motor_mass = motor_mass.base_value         # Motor Mass (kg) 
+        self.motor_cg = motor_cg.base_value             # Motor center of gravity (m)
+        self.tank_location = tank_location.base_value   # Tank center of gravity (m)
+        self.grain_location = grain_location.base_value # Grain center of gravity (m)
+
+class SimulationConfig:
+    def __init__(self, timestep: float, max_run_time: float, max_burn_time: float, ambient_pressure: PressureValue,
+                 hardware: HardwareConfig, prop: PropellantConfig, nozzle: NozzleConfig, mass_properties: MassProperties | None = None):
+        self.timestep = timestep
+        self.max_run_time = max_run_time
+        self.max_burn_time = max_burn_time
+        self.ambient_pressure = ambient_pressure.base_value
+
+        grain_volume = d_to_a(prop.OD) * prop.length
+        if hardware.chamber_volume is None:
+            hardware.chamber_volume = grain_volume
+
+        assert hardware.chamber_volume >= grain_volume, "Chammber cannot be smaller than fuel grain!"
+
+        self.hardware = hardware
+        self.prop = prop
+        self.nozzle = nozzle
+        self.mass_properties = mass_properties
+
+###
+### STATE MODELS
+###
+class InitialConditions:
+    def __init__(self, chamber_pressure: PressureValue, mass_ox: MassValue, tank_pressure: PressureValue):
+        self.chamber_pressure = chamber_pressure.base_value
+        self.mass_ox = mass_ox.base_value
+        self.tank_temp = None
+        self.tank_pressure = tank_pressure.base_value
+
+    def __init__(self, chamber_pressure: PressureValue, mass_ox: MassValue, tank_temp: TemperatureValue):
+        self.chamber_pressure = chamber_pressure.base_value
+        self.mass_ox = mass_ox.base_value
+        self.tank_temp = tank_temp.base_value
+        self.tank_pressure = None
+
+class TankState:
+    def __init__(self, config: SimulationConfig, init_conditions: InitialConditions):
+        if init_conditions.tank_temp is None:
+            self.temp = fsolve(lambda T: vapor_pressure(T) - init_conditions.tank_pressure, 273.15)[0] # TODO: Error handling (see RootResults), root method selection
+        else:
+            self.temp = init_conditions.tank_temp
+        self.ox_props = get_NOX_properties(self.temp)
+        self.mass_ox = init_conditions.mass_ox
+        self.mass_ox_old = self.mass_ox
+        self.mdot_ox = 0.0
+        self.mdot_vapor = 0.0
+        self.pressure = self.ox_props.Pv
+
+        self.mLiq_new = TankState.liquid_mass(config.hardware.tank_volume, self.mass_ox, self.ox_props)
+        self.mLiq_old = self.mLiq_new + 1.0
+
+    def liquid_mass(tank_vol: float, mass_ox: float, ox_props: NOXProperties) -> float:
+        # This can take in standard floats because unit normalization should have already been done... should have.
+        return ((tank_vol - (mass_ox / ox_props.rho_v)) / 
+                ((1 / ox_props.rho_l) - (1 / ox_props.rho_v)))
+
+class ChamberState:
+    AIR_RHO = 1.225 # Density of air (kg/m^3)
+
+    def __init__(self, config: SimulationConfig, init_conditions: InitialConditions):
+        self.pressure = init_conditions.chamber_pressure
+
+        grain_volume    = (d_to_a(config.prop.OD) - d_to_a(config.prop.ID)) * config.prop.length
+        self.mass_fuel  = grain_volume * config.prop.material.rho
+        self.mass_gas   = ChamberState.AIR_RHO * (config.hardware.chamber_volume - grain_volume) 
+
+        self.vdot_fuel = 0.0
+        self.mdot_fuel = 0.0
+        self.mdot_exit = 0.0
+        self.rdot = 0.0
+
+        self.grain_ID = config.prop.ID
+        self.OF_ratio = config.prop.const_OF_ratio if config.prop.const_OF_ratio is not None else 0.0
+
+class MassState:
+    def __init__(self, config: SimulationConfig, init_conditions: InitialConditions):
+        self.cg = 0.0
+        self.mass_total = 0.0
+
+class CombustionState:
+    def __init__(self):
+        self.k = 0.0
+        self.T = 0.0
+        self.M = 0.0
+        self.R = 0.0
+        self.rho = 0.0
+        self.cstar = 0.0
+
+class SimulationState:
+    def __init__(self, config: SimulationConfig, init_conditions: InitialConditions):
+        self.time = 0.0
+
+        self.tank = TankState(config, init_conditions)
+        self.chamber = ChamberState(config, init_conditions)
+        self.combustion = CombustionState()
+        self.mass = MassState(config, init_conditions)
+        self.thrust = 0.0
+        self.pressure_drop = 0.0
+
+        # For interpolating oxidizer pressure curve
+        self.total_vapor_dP = 0.0
+        self.interpolation_timesteps = 1e-10 # Avoids potential divide-by-zero
+
+###
+### SIMULATION OUTPUT
+###
+class SimulationOutput:
+    def __init__(self, config: SimulationConfig):
+        self.time : List[float] = []
+
+        self.mass_ox : List[MassValue] = []
+        self.mdot_ox : List[MassValue] = []
+        self.tank_pressure : List[PressureValue] = []
+
+        self.mass_fuel : List[MassValue] = []
+        self.mdot_fuel : List[MassValue] = []
+        self.mdot_exit : List[MassValue] = []
+        self.rdot : List[LengthValue] = []
+        self.chamber_pressure : List[PressureValue] = []
+        self.OF_ratio : List[float] = []
+        self.grain_ID : List[LengthValue] = []
+
+        self.pressure_drop : List[PressureValue] = []
+        self.thrust : List[MassValue] = []
+
+        if config.mass_properties is not None:
+            self.total_mass : List[MassValue] = []
+            self.cg : List[LengthValue] = []
+            self.include_mass_props = True
+        else:
+            self.include_mass_props = False
+
+        self.end_condition : str = "Simulation Unfinished"
+
+    def append_state(self, state: SimulationState):
+        self.time.append(state.time)
+
+        # TODO: Mass flux / length 'flux' unit? Probably not necessary - denominator is always 'seconds'
+        self.mass_ox.append(MassValue(state.tank.mass_ox, MassUnit.KILOGRAMS))
+        self.mdot_ox.append(MassValue(state.tank.mdot_ox, MassUnit.KILOGRAMS))
+        self.tank_pressure.append(PressureValue(state.tank.pressure, PressureUnit.PA))
+
+        self.mass_fuel.append(MassValue(state.chamber.mass_fuel, MassUnit.KILOGRAMS))
+        self.mdot_fuel.append(MassValue(state.chamber.mdot_fuel, MassUnit.KILOGRAMS))
+        self.mdot_exit.append(MassValue(state.chamber.mdot_exit, MassUnit.KILOGRAMS))
+        self.rdot.append(LengthValue(state.chamber.rdot, LengthUnit.MILLIMETERS)) # Regression is in mm/s by default
+        self.chamber_pressure.append(PressureValue(state.chamber.pressure, PressureUnit.PA))
+        self.OF_ratio.append(state.chamber.OF_ratio)
+        self.grain_ID.append(LengthValue(state.chamber.grain_ID, LengthUnit.METERS))
+
+        self.pressure_drop.append(PressureValue(state.pressure_drop, PressureUnit.PA))
+        self.thrust.append(MassValue(state.thrust, MassUnit.NEWTONS))
+
+        if self.include_mass_props:
+            self.total_mass.append(MassValue(state.mass.mass_total, MassUnit.KILOGRAMS))
+            self.cg.append(LengthUnit(state.mass.cg, LengthUnit.METERS))
+
+    def max_thrust(self):
+        return MassValue(max([t.base_value for t in self.thrust], MassUnit.NEWTONS))