Skip to content

Simulation

Main entry point for running internal ballistics simulations. InternalBallisticsSimulation takes a motor model and simulation parameters (time step, igniter pressure, external pressure), then marches through time using an RK4 solver.

The run ends on whichever of these two conditions fires first:

Condition Fires when Ends
Loss of choking Chamber pressure falls to the critical ratio of the ambient pressure, roughly 1.7 to 1.8 times external_pressure Runs against an atmosphere
Tail-off The motor has burnt out and thrust has decayed to TAIL_OFF_THRUST_FRACTION (0.1%) of its peak Vacuum and upper stage runs, where the nozzle stays choked all the way down

run() returns a frozen SimulationResult (subclassed per motor type — SolidSimulationResult, BiliquidSimulationResult) carrying the full time-series data (thrust, chamber pressure, propellant mass, efficiency losses, …) along with derived scalars (total impulse, specific impulse, burn time). Each result class provides a report() method to print a human-readable summary and a summary() method that returns the scalar metrics as a dict.

The per-step accumulator state used by the integrator (MotorState and its subclasses SolidMotorState, BiliquidEngineState) lives in this same package — it is rarely consumed directly outside the simulation loop. Alongside it, the abstract TimestepConditions and its per-engine subclasses SolidTimestepConditions, BiliquidTimestepConditions snapshot the operating quantities of one timestep and are handed to the nozzle loss model.

machwave.simulation

InternalBallisticsSimulation

Internal ballistics simulation class.

Attributes:

Name Type Description
motor Motor

Motor object.
params InternalBallisticsSimulationParams

Simulation parameters.
Source code in machwave/simulation/base.py 
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
class InternalBallisticsSimulation:
    """
    Internal ballistics simulation class.

    Attributes:
        motor: Motor object.
        params: Simulation parameters.
    """

    _STATE_CLASS_BY_MOTOR_TYPE = (
        (motors.SolidMotor, solid_states.SolidMotorState),
        (motors.BiliquidEngine, biliquid_states.BiliquidEngineState),
    )

    def __init__(
        self,
        motor: motors.Motor,
        params: InternalBallisticsSimulationParams,
    ) -> None:
        """
        Initialize an internal ballistics simulation.

        Args:
            motor: Motor model to simulate.
            params: Simulation parameters.
        """
        self.motor: motors.Motor = motor
        self.params: InternalBallisticsSimulationParams = params

    def _build_motor_state(self) -> simulation_states.MotorState:
        """Build the motor state matching the configured motor type."""
        state_kwargs = {
            "motor": self.motor,
            "igniter_pressure": self.params.igniter_pressure,
            "external_pressure": self.params.external_pressure,
        }
        for motor_type, state_class in self._STATE_CLASS_BY_MOTOR_TYPE:
            if isinstance(self.motor, motor_type):
                return state_class(**state_kwargs)
        raise ValueError(f"Unsupported motor type: {type(self.motor).__name__}.")

    def run(self) -> simulation_results.SimulationResult:
        """Run the simulation to thrust termination and return its result."""
        motor_state = self._build_motor_state()

        d_t = self.params.d_t
        external_pressure = self.params.external_pressure

        while not motor_state.end_thrust:
            motor_state.run_timestep(d_t, external_pressure)

        return motor_state.build_result()

__init__(motor, params)

Initialize an internal ballistics simulation.

Parameters:

Name Type Description Default
motor Motor

Motor model to simulate.

 required
params InternalBallisticsSimulationParams

Simulation parameters.

 required
Source code in machwave/simulation/base.py 
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def __init__(
    self,
    motor: motors.Motor,
    params: InternalBallisticsSimulationParams,
) -> None:
    """
    Initialize an internal ballistics simulation.

    Args:
        motor: Motor model to simulate.
        params: Simulation parameters.
    """
    self.motor: motors.Motor = motor
    self.params: InternalBallisticsSimulationParams = params

run()

Run the simulation to thrust termination and return its result.

Source code in machwave/simulation/base.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def run(self) -> simulation_results.SimulationResult:
    """Run the simulation to thrust termination and return its result."""
    motor_state = self._build_motor_state()

    d_t = self.params.d_t
    external_pressure = self.params.external_pressure

    while not motor_state.end_thrust:
        motor_state.run_timestep(d_t, external_pressure)

    return motor_state.build_result()

InternalBallisticsSimulationParams dataclass

Parameters for an internal ballistics simulation.

Attributes:

Name Type Description
d_t float

Time step [s].
igniter_pressure float

Igniter pressure [Pa].
external_pressure float

External pressure [Pa].
Source code in machwave/simulation/base.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
@dataclass
class InternalBallisticsSimulationParams:
    """
    Parameters for an internal ballistics simulation.

    Attributes:
        d_t: Time step [s].
        igniter_pressure: Igniter pressure [Pa].
        external_pressure: External pressure [Pa].
    """

    d_t: float
    igniter_pressure: float
    external_pressure: float

    def __post_init__(self) -> None:
        """
        Validate the simulation parameters.

        Raises:
            ValueError: If any field is outside its valid physical range.
        """
        if not 0.0 < self.d_t <= MAX_TIME_STEP:
            raise ValueError(f"d_t must be in (0, {MAX_TIME_STEP}] s, got {self.d_t}")

        if self.igniter_pressure < SEA_LEVEL_PRESSURE:
            raise ValueError(
                f"igniter_pressure must be at least {SEA_LEVEL_PRESSURE} Pa, got "
                f"{self.igniter_pressure}"
            )

        if self.external_pressure < 0.0:
            raise ValueError(
                f"external_pressure must be non-negative, got {self.external_pressure}"
            )

__post_init__()

Validate the simulation parameters.

Raises:

Type Description
ValueError

If any field is outside its valid physical range.
Source code in machwave/simulation/base.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
def __post_init__(self) -> None:
    """
    Validate the simulation parameters.

    Raises:
        ValueError: If any field is outside its valid physical range.
    """
    if not 0.0 < self.d_t <= MAX_TIME_STEP:
        raise ValueError(f"d_t must be in (0, {MAX_TIME_STEP}] s, got {self.d_t}")

    if self.igniter_pressure < SEA_LEVEL_PRESSURE:
        raise ValueError(
            f"igniter_pressure must be at least {SEA_LEVEL_PRESSURE} Pa, got "
            f"{self.igniter_pressure}"
        )

    if self.external_pressure < 0.0:
        raise ValueError(
            f"external_pressure must be non-negative, got {self.external_pressure}"
        )

MotorState

Bases: ABC

Defines the states and iteration step for a motor operation.

Source code in machwave/simulation/states.py 
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
class MotorState(ABC):
    """Defines the states and iteration step for a motor operation."""

    result_class: ClassVar[type["simulation_results.SimulationResult"]]

    def __init__(
        self,
        motor: motors.Motor,
        igniter_pressure: float,
        external_pressure: float,
    ) -> None:
        """
        Initialize a motor state.

        Args:
            motor: Motor to track.
            igniter_pressure: Initial chamber pressure from the igniter [Pa].
            external_pressure: Ambient pressure [Pa].
        """
        self.motor = motor
        self.external_pressure = external_pressure

        self.time: SimulationStateArray = [0.0]
        self.chamber_pressure: SimulationStateArray = [igniter_pressure]

        self.propellant_mass: SimulationStateArray = []
        self.exit_pressure: SimulationStateArray = []
        self.ideal_thrust_coefficient: SimulationStateArray = []
        self.thrust_coefficient: SimulationStateArray = []
        self.thrust: SimulationStateArray = []
        self.peak_thrust: float = 0.0
        self.nozzle_efficiency: SimulationStateArray = []
        self.loss_fractions: dict[str, SimulationStateArray] = {
            name: [] for name in motor.nozzle_loss_model.component_names
        }

        self._thrust_time: float | None = None
        self._burn_time: float | None = None

        self.end_thrust: bool = False
        self.end_burn: bool = False

    @abstractmethod
    def run_timestep(self, *args, **kwargs) -> None:
        """Advance the per-step accumulators by one time increment."""

    def _ideal_thrust_coefficient_terms(
        self,
        k_exhaust: float,
        chamber_pressure: float,
        external_pressure: float,
    ) -> tuple[float, float, float, float]:
        """
        Resolve the separated exit conditions and ideal thrust coefficient terms.

        Appends the effective exit pressure and the ideal thrust coefficient for the
        timestep.

        Args:
            k_exhaust: Isentropic exponent at the nozzle exit.
            chamber_pressure: Chamber pressure [Pa].
            external_pressure: Ambient pressure [Pa].

        Returns:
            The effective expansion ratio, effective exit pressure [Pa], and the
            momentum and pressure terms of the ideal thrust coefficient.
        """
        nozzle = self.motor.thrust_chamber.nozzle
        effective_expansion_ratio, exit_pressure = (
            nozzle_core.get_separated_exit_conditions(
                k_exhaust,
                nozzle.expansion_ratio,
                chamber_pressure,
                external_pressure,
                nozzle.separation_pressure_ratio,
            )
        )
        self.exit_pressure.append(exit_pressure)

        ideal_momentum_term, ideal_pressure_term = (
            nozzle_core.get_ideal_thrust_coefficient_terms(
                chamber_pressure,
                exit_pressure,
                external_pressure,
                effective_expansion_ratio,
                k_exhaust,
            )
        )
        self.ideal_thrust_coefficient.append(ideal_momentum_term + ideal_pressure_term)
        return (
            effective_expansion_ratio,
            exit_pressure,
            ideal_momentum_term,
            ideal_pressure_term,
        )

    def _apply_nozzle_losses(
        self,
        ideal_momentum_term: float,
        ideal_pressure_term: float,
        timestep_conditions: TimestepConditions,
        chamber_pressure: float,
    ) -> None:
        """
        Derate the ideal thrust coefficient terms and record the loss outputs.

        Appends the realized nozzle efficiency, each component loss fraction, the
        corrected thrust coefficient, and the thrust for the timestep.

        Args:
            ideal_momentum_term: Momentum term of the ideal thrust coefficient.
            ideal_pressure_term: Pressure term of the ideal thrust coefficient.
            timestep_conditions: Engine conditions at a point in time.
            chamber_pressure: Chamber pressure [Pa].
        """
        nozzle = self.motor.thrust_chamber.nozzle
        loss_result = self.motor.nozzle_loss_model.evaluate(
            ideal_momentum_term, ideal_pressure_term, timestep_conditions
        )
        self.nozzle_efficiency.append(loss_result.nozzle_efficiency)
        for name, fraction in loss_result.loss_fractions.items():
            self.loss_fractions[name].append(fraction)

        thrust_coefficient = loss_result.momentum_term + loss_result.pressure_term
        self.thrust_coefficient.append(thrust_coefficient)
        thrust = nozzle_core.get_thrust_from_thrust_coefficient(
            thrust_coefficient, chamber_pressure, nozzle.get_throat_area()
        )
        self.thrust.append(thrust)
        self.peak_thrust = max(self.peak_thrust, thrust)

    def _update_thrust_termination(
        self,
        time: float,
        chamber_pressure: float,
        external_pressure: float,
        k_chamber: float,
    ) -> bool:
        """
        Flag thrust termination once the nozzle un-chokes or the tail-off ends.

        Loss of choking never fires against a vacuum, where the nozzle stays choked
        for any chamber pressure, so a burnt out motor also terminates once its
        thrust decays to ``TAIL_OFF_THRUST_FRACTION`` of the peak.

        Args:
            time: Time at the start of the timestep [s].
            chamber_pressure: Chamber pressure [Pa].
            external_pressure: Ambient pressure [Pa].
            k_chamber: Isentropic exponent in the chamber.

        Returns:
            True if thrust has terminated on this timestep.
        """
        is_choked = isentropic.is_flow_choked(
            chamber_pressure,
            external_pressure,
            isentropic.get_critical_pressure_ratio(k_chamber),
        )
        has_tailed_off = (
            self.end_burn
            and self.thrust[-1] < TAIL_OFF_THRUST_FRACTION * self.peak_thrust
        )
        if is_choked and not has_tailed_off:
            return False

        self._thrust_time = time
        self.end_thrust = True
        return True

    def build_result(self) -> "simulation_results.SimulationResult":
        """Return a frozen ``SimulationResult`` snapshot of this state."""
        return self.result_class.from_state(self)

    @property
    def initial_propellant_mass(self) -> float:
        """Return the initial propellant mass [kg]."""
        return self.motor.initial_propellant_mass

    @property
    def thrust_time(self) -> float:
        """
        Return the thrust time [s].

        Raises:
            ValueError: If the simulation has not yet completed.
        """
        if self._thrust_time is None:
            raise ValueError("Thrust time has not been set, run the simulation.")
        return self._thrust_time

    @property
    def burn_time(self) -> float | None:
        """
        Return the burn time [s], or None if burnout was never reached.

        Only propellant depletion defines a burn time: a run whose thrust
        terminates with propellant remaining leaves it undefined.
        """
        return self._burn_time

burn_time property

Return the burn time [s], or None if burnout was never reached.

Only propellant depletion defines a burn time: a run whose thrust terminates with propellant remaining leaves it undefined.

initial_propellant_mass property

Return the initial propellant mass [kg].

thrust_time property

Return the thrust time [s].

Raises:

Type Description
ValueError

If the simulation has not yet completed.

__init__(motor, igniter_pressure, external_pressure)

Initialize a motor state.

Parameters:

Name Type Description Default
motor Motor

Motor to track.

 required
igniter_pressure float

Initial chamber pressure from the igniter [Pa].

 required
external_pressure float

Ambient pressure [Pa].

 required
Source code in machwave/simulation/states.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def __init__(
    self,
    motor: motors.Motor,
    igniter_pressure: float,
    external_pressure: float,
) -> None:
    """
    Initialize a motor state.

    Args:
        motor: Motor to track.
        igniter_pressure: Initial chamber pressure from the igniter [Pa].
        external_pressure: Ambient pressure [Pa].
    """
    self.motor = motor
    self.external_pressure = external_pressure

    self.time: SimulationStateArray = [0.0]
    self.chamber_pressure: SimulationStateArray = [igniter_pressure]

    self.propellant_mass: SimulationStateArray = []
    self.exit_pressure: SimulationStateArray = []
    self.ideal_thrust_coefficient: SimulationStateArray = []
    self.thrust_coefficient: SimulationStateArray = []
    self.thrust: SimulationStateArray = []
    self.peak_thrust: float = 0.0
    self.nozzle_efficiency: SimulationStateArray = []
    self.loss_fractions: dict[str, SimulationStateArray] = {
        name: [] for name in motor.nozzle_loss_model.component_names
    }

    self._thrust_time: float | None = None
    self._burn_time: float | None = None

    self.end_thrust: bool = False
    self.end_burn: bool = False

build_result()

Return a frozen SimulationResult snapshot of this state.

Source code in machwave/simulation/states.py 
191
192
193
def build_result(self) -> "simulation_results.SimulationResult":
    """Return a frozen ``SimulationResult`` snapshot of this state."""
    return self.result_class.from_state(self)

run_timestep(*args, **kwargs) abstractmethod

Advance the per-step accumulators by one time increment.

Source code in machwave/simulation/states.py 
63
64
65
@abstractmethod
def run_timestep(self, *args, **kwargs) -> None:
    """Advance the per-step accumulators by one time increment."""

SimulationResult dataclass

Bases: ABC, Generic[StateT]

Results of a finished internal ballistics simulation.

Source code in machwave/simulation/results.py 
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
@dataclass(frozen=True, kw_only=True)
class SimulationResult(ABC, Generic[StateT]):
    """Results of a finished internal ballistics simulation."""

    time: SimulationResultArray
    propellant_mass: SimulationResultArray
    chamber_pressure: SimulationResultArray
    exit_pressure: SimulationResultArray
    thrust_coefficient: SimulationResultArray
    ideal_thrust_coefficient: SimulationResultArray
    thrust: SimulationResultArray
    nozzle_efficiency: SimulationResultArray
    loss_fractions: dict[str, SimulationResultArray]
    loss_labels: dict[str, str]
    # None when the run terminated with propellant remaining.
    burn_time: float | None
    thrust_time: float
    end_thrust: bool
    end_burn: bool
    initial_propellant_mass: float
    total_impulse: float
    specific_impulse: float

    @classmethod
    def from_state(cls, state: StateT) -> "SimulationResult":
        """Build a `SimulationResult` from a finished motor state."""
        return cls(
            **cls._collect_base_fields(state),
            **cls._collect_extra_fields(state),
        )

    @classmethod
    def _collect_base_fields(cls, state: StateT) -> dict[str, Any]:
        time = np.asarray(state.time)
        thrust = np.asarray(state.thrust)
        total_impulse = performance.get_total_impulse(thrust, time)
        initial_propellant_mass = state.motor.initial_propellant_mass
        return {
            "time": time,
            "propellant_mass": np.asarray(state.propellant_mass),
            "chamber_pressure": np.asarray(state.chamber_pressure),
            "exit_pressure": np.asarray(state.exit_pressure),
            "thrust_coefficient": np.asarray(state.thrust_coefficient),
            "ideal_thrust_coefficient": np.asarray(state.ideal_thrust_coefficient),
            "thrust": thrust,
            "nozzle_efficiency": np.asarray(state.nozzle_efficiency),
            "loss_fractions": {
                name: np.asarray(series)
                for name, series in state.loss_fractions.items()
            },
            "loss_labels": dict(state.motor.nozzle_loss_model.component_labels),
            "burn_time": state.burn_time,
            "thrust_time": state.thrust_time,
            "end_thrust": state.end_thrust,
            "end_burn": state.end_burn,
            "initial_propellant_mass": initial_propellant_mass,
            "total_impulse": total_impulse,
            "specific_impulse": performance.get_specific_impulse(
                total_impulse=total_impulse,
                initial_propellant_mass=initial_propellant_mass,
            ),
        }

    @classmethod
    @abstractmethod
    def _collect_extra_fields(cls, state: StateT) -> dict[str, Any]:
        """Build constructor kwargs for fields specific to the subclass."""

    def report(self, file: IO = sys.stdout) -> None:
        """Print a human-readable report of the simulation result."""
        print("\nINTERNAL BALLISTICS SIMULATION RESULTS", file=file)
        self._report_body(file)

    @abstractmethod
    def _report_body(self, file: IO) -> None:
        """Print the subclass-specific portion of the report."""

    def _report_nozzle_losses(self, file: IO) -> None:
        """Print the nozzle efficiency and each loss component's mean fraction."""
        print("\nNOZZLE", file=file)
        print(
            f"  Average nozzle efficiency: {np.mean(self.nozzle_efficiency):.3%}",
            file=file,
        )
        for name, series in self.loss_fractions.items():
            label = self.loss_labels[name]
            print(f"  Average {label} fraction: {np.mean(series):.3%}", file=file)

    def _format_burn_time(self, decimals: int = 3) -> str:
        """Format the burn time for a report, or flag it as never reached."""
        if self.burn_time is None:
            return "not reached"
        return f"{self.burn_time:.{decimals}f} s"

    def summary(self) -> dict[str, float | None]:
        """Return a mapping of headline scalar metrics for this result."""
        return {
            "burn_time": self.burn_time,
            "thrust_time": self.thrust_time,
            "initial_propellant_mass": self.initial_propellant_mass,
            "total_impulse": self.total_impulse,
            "specific_impulse": self.specific_impulse,
            "peak_chamber_pressure": float(np.max(self.chamber_pressure)),
            "mean_chamber_pressure": float(np.mean(self.chamber_pressure)),
            "peak_thrust": float(np.max(self.thrust)),
            "mean_thrust": float(np.mean(self.thrust)),
            **self._extra_summary(),
        }

    def _extra_summary(self) -> dict[str, float]:
        """Return subclass-specific scalar metrics to merge into `summary()`."""
        return {}

from_state(state) classmethod

Build a SimulationResult from a finished motor state.

Source code in machwave/simulation/results.py 
42
43
44
45
46
47
48
@classmethod
def from_state(cls, state: StateT) -> "SimulationResult":
    """Build a `SimulationResult` from a finished motor state."""
    return cls(
        **cls._collect_base_fields(state),
        **cls._collect_extra_fields(state),
    )

report(file=sys.stdout)

Print a human-readable report of the simulation result.

Source code in machwave/simulation/results.py 
87
88
89
90
def report(self, file: IO = sys.stdout) -> None:
    """Print a human-readable report of the simulation result."""
    print("\nINTERNAL BALLISTICS SIMULATION RESULTS", file=file)
    self._report_body(file)

summary()

Return a mapping of headline scalar metrics for this result.

Source code in machwave/simulation/results.py 
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def summary(self) -> dict[str, float | None]:
    """Return a mapping of headline scalar metrics for this result."""
    return {
        "burn_time": self.burn_time,
        "thrust_time": self.thrust_time,
        "initial_propellant_mass": self.initial_propellant_mass,
        "total_impulse": self.total_impulse,
        "specific_impulse": self.specific_impulse,
        "peak_chamber_pressure": float(np.max(self.chamber_pressure)),
        "mean_chamber_pressure": float(np.mean(self.chamber_pressure)),
        "peak_thrust": float(np.max(self.thrust)),
        "mean_thrust": float(np.mean(self.thrust)),
        **self._extra_summary(),
    }

base

InternalBallisticsSimulation

Internal ballistics simulation class.

Attributes:

Name Type Description
motor Motor

Motor object.
params InternalBallisticsSimulationParams

Simulation parameters.
Source code in machwave/simulation/base.py 
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
class InternalBallisticsSimulation:
    """
    Internal ballistics simulation class.

    Attributes:
        motor: Motor object.
        params: Simulation parameters.
    """

    _STATE_CLASS_BY_MOTOR_TYPE = (
        (motors.SolidMotor, solid_states.SolidMotorState),
        (motors.BiliquidEngine, biliquid_states.BiliquidEngineState),
    )

    def __init__(
        self,
        motor: motors.Motor,
        params: InternalBallisticsSimulationParams,
    ) -> None:
        """
        Initialize an internal ballistics simulation.

        Args:
            motor: Motor model to simulate.
            params: Simulation parameters.
        """
        self.motor: motors.Motor = motor
        self.params: InternalBallisticsSimulationParams = params

    def _build_motor_state(self) -> simulation_states.MotorState:
        """Build the motor state matching the configured motor type."""
        state_kwargs = {
            "motor": self.motor,
            "igniter_pressure": self.params.igniter_pressure,
            "external_pressure": self.params.external_pressure,
        }
        for motor_type, state_class in self._STATE_CLASS_BY_MOTOR_TYPE:
            if isinstance(self.motor, motor_type):
                return state_class(**state_kwargs)
        raise ValueError(f"Unsupported motor type: {type(self.motor).__name__}.")

    def run(self) -> simulation_results.SimulationResult:
        """Run the simulation to thrust termination and return its result."""
        motor_state = self._build_motor_state()

        d_t = self.params.d_t
        external_pressure = self.params.external_pressure

        while not motor_state.end_thrust:
            motor_state.run_timestep(d_t, external_pressure)

        return motor_state.build_result()
__init__(motor, params)

Initialize an internal ballistics simulation.

Parameters:

Name Type Description Default
motor Motor

Motor model to simulate.

 required
params InternalBallisticsSimulationParams

Simulation parameters.

 required
Source code in machwave/simulation/base.py 
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def __init__(
    self,
    motor: motors.Motor,
    params: InternalBallisticsSimulationParams,
) -> None:
    """
    Initialize an internal ballistics simulation.

    Args:
        motor: Motor model to simulate.
        params: Simulation parameters.
    """
    self.motor: motors.Motor = motor
    self.params: InternalBallisticsSimulationParams = params
run()

Run the simulation to thrust termination and return its result.

Source code in machwave/simulation/base.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def run(self) -> simulation_results.SimulationResult:
    """Run the simulation to thrust termination and return its result."""
    motor_state = self._build_motor_state()

    d_t = self.params.d_t
    external_pressure = self.params.external_pressure

    while not motor_state.end_thrust:
        motor_state.run_timestep(d_t, external_pressure)

    return motor_state.build_result()

InternalBallisticsSimulationParams dataclass

Parameters for an internal ballistics simulation.

Attributes:

Name Type Description
d_t float

Time step [s].
igniter_pressure float

Igniter pressure [Pa].
external_pressure float

External pressure [Pa].
Source code in machwave/simulation/base.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
@dataclass
class InternalBallisticsSimulationParams:
    """
    Parameters for an internal ballistics simulation.

    Attributes:
        d_t: Time step [s].
        igniter_pressure: Igniter pressure [Pa].
        external_pressure: External pressure [Pa].
    """

    d_t: float
    igniter_pressure: float
    external_pressure: float

    def __post_init__(self) -> None:
        """
        Validate the simulation parameters.

        Raises:
            ValueError: If any field is outside its valid physical range.
        """
        if not 0.0 < self.d_t <= MAX_TIME_STEP:
            raise ValueError(f"d_t must be in (0, {MAX_TIME_STEP}] s, got {self.d_t}")

        if self.igniter_pressure < SEA_LEVEL_PRESSURE:
            raise ValueError(
                f"igniter_pressure must be at least {SEA_LEVEL_PRESSURE} Pa, got "
                f"{self.igniter_pressure}"
            )

        if self.external_pressure < 0.0:
            raise ValueError(
                f"external_pressure must be non-negative, got {self.external_pressure}"
            )
__post_init__()

Validate the simulation parameters.

Raises:

Type Description
ValueError

If any field is outside its valid physical range.
Source code in machwave/simulation/base.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
def __post_init__(self) -> None:
    """
    Validate the simulation parameters.

    Raises:
        ValueError: If any field is outside its valid physical range.
    """
    if not 0.0 < self.d_t <= MAX_TIME_STEP:
        raise ValueError(f"d_t must be in (0, {MAX_TIME_STEP}] s, got {self.d_t}")

    if self.igniter_pressure < SEA_LEVEL_PRESSURE:
        raise ValueError(
            f"igniter_pressure must be at least {SEA_LEVEL_PRESSURE} Pa, got "
            f"{self.igniter_pressure}"
        )

    if self.external_pressure < 0.0:
        raise ValueError(
            f"external_pressure must be non-negative, got {self.external_pressure}"
        )

biliquid

BiliquidEngineState

Bases: MotorState

State for a biliquid rocket engine.

Source code in machwave/simulation/biliquid/states.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
class BiliquidEngineState(simulation_states.MotorState):
    """State for a biliquid rocket engine."""

    motor: motors.BiliquidEngine
    result_class = biliquid_results.BiliquidSimulationResult

    def __init__(
        self,
        motor: motors.BiliquidEngine,
        igniter_pressure: float,
        external_pressure: float,
    ) -> None:
        """
        Initialize a biliquid engine state.

        Args:
            motor: Biliquid engine to track.
            igniter_pressure: Initial chamber pressure from the igniter [Pa].
            external_pressure: Ambient pressure [Pa].
        """
        super().__init__(
            motor=motor,
            igniter_pressure=igniter_pressure,
            external_pressure=external_pressure,
        )

        feed_system = motor.feed_system
        self.lines: dict[str, line_models.PropellantLine] = feed_system.lines
        self.oxidizer_line_names = tuple(
            line.name
            for line in feed_system.get_lines_with_role(
                propellant_components.ComponentRole.OXIDIZER
            )
        )
        self.fuel_line_names = tuple(
            line.name
            for line in feed_system.get_lines_with_role(
                propellant_components.ComponentRole.FUEL
            )
        )

        initial_states = {name: line.initial_state for name, line in self.lines.items()}
        self.fluid_mass_per_line: dict[str, simulation_states.SimulationStateArray] = {
            name: [state.fluid_mass] for name, state in initial_states.items()
        }

        # Second state variable of a tank running an energy balance, which the
        # integrator carries beside the mass. An isothermal tank has none.
        self.internal_energy_per_line: dict[str, float | None] = {
            name: state.internal_energy for name, state in initial_states.items()
        }

        self.mass_flow_rate_per_line: dict[
            str, simulation_states.SimulationStateArray
        ] = {name: [] for name in self.lines}
        self.tank_pressure_per_line: dict[
            str, simulation_states.SimulationStateArray
        ] = {name: [] for name in self.lines}
        self.tank_temperature_per_line: dict[
            str, simulation_states.SimulationStateArray
        ] = {name: [] for name in self.lines}
        self.oxidizer_to_fuel_ratio: simulation_states.SimulationStateArray = []

    def _evaluate_propellant_properties(
        self,
        chamber_pressure: float,
        mixture_ratio: float | None,
    ) -> propellant_properties_models.ThermochemicalProperties:
        """Evaluate the propellant at the chamber pressure and mixture ratio."""
        return self.motor.propellant.evaluate(
            chamber_pressure=chamber_pressure,
            expansion_ratio=self.motor.thrust_chamber.nozzle.expansion_ratio,
            mixture_ratio=mixture_ratio,
        )

    def _get_mixture_ratio(self, mass_flows: Mapping[str, float]) -> float | None:
        """
        Total oxidizer flow over total fuel flow.

        None where the ratio has no meaning: a monoliquid, or an engine whose
        lines have stopped flowing.
        """
        oxidizer_flow = sum(mass_flows[name] for name in self.oxidizer_line_names)
        fuel_flow = sum(mass_flows[name] for name in self.fuel_line_names)

        if oxidizer_flow <= 0.0 or fuel_flow <= 0.0:
            return None

        return oxidizer_flow / fuel_flow

    def run_timestep(
        self,
        d_t: float,
        external_pressure: float,
    ) -> None:
        """
        Iterate the engine operation by calculating operational parameters.

        Args:
            d_t: Time increment [s].
            external_pressure: External pressure [Pa].
        """
        nozzle = self.motor.thrust_chamber.nozzle
        feed_system = self.motor.feed_system

        time = self.time[-1]
        chamber_pressure = self.chamber_pressure[-1]
        fluid_masses = {
            name: series[-1] for name, series in self.fluid_mass_per_line.items()
        }

        propellant_mass = sum(fluid_masses.values())
        self.propellant_mass.append(propellant_mass)

        inlet_states = feed_system.get_inlet_states(
            {
                name: line_models.LineState(
                    fluid_mass=fluid_masses[name],
                    internal_energy=self.internal_energy_per_line[name],
                )
                for name in self.lines
            }
        )

        tank_pressures = {name: inlet.pressure for name, inlet in inlet_states.items()}
        tank_temperatures = {
            name: inlet.temperature for name, inlet in inlet_states.items()
        }
        for name in self.lines:
            self.tank_pressure_per_line[name].append(tank_pressures[name])
            self.tank_temperature_per_line[name].append(tank_temperatures[name])

        is_feeding = (
            not self.end_burn
            and all(fluid_mass > 0 for fluid_mass in fluid_masses.values())
            and all(pressure > chamber_pressure for pressure in tank_pressures.values())
        )
        injector_flows = functools.partial(
            get_injector_mass_flows,
            injector=self.motor.thrust_chamber.injector,
            inlet_states=inlet_states,
            line_masses=fluid_masses,
            is_feeding=is_feeding,
            d_t=d_t,
        )
        mass_flows = injector_flows(chamber_pressure)
        for name, mass_flow in mass_flows.items():
            self.mass_flow_rate_per_line[name].append(mass_flow)
        mass_consumed = {name: flow * d_t for name, flow in mass_flows.items()}

        mixture_ratio = self._get_mixture_ratio(mass_flows)
        oxidizer_to_fuel_ratio = math.nan if mixture_ratio is None else mixture_ratio
        self.oxidizer_to_fuel_ratio.append(oxidizer_to_fuel_ratio)

        propellant_properties = self._evaluate_propellant_properties(
            chamber_pressure=chamber_pressure,
            mixture_ratio=mixture_ratio,
        )

        (
            effective_expansion_ratio,
            exit_pressure,
            ideal_momentum_term,
            ideal_pressure_term,
        ) = self._ideal_thrust_coefficient_terms(
            propellant_properties.k_exhaust, chamber_pressure, external_pressure
        )

        timestep_conditions = BiliquidTimestepConditions(
            time=time,
            chamber_pressure=chamber_pressure,
            external_pressure=external_pressure,
            exit_pressure=exit_pressure,
            effective_expansion_ratio=effective_expansion_ratio,
            free_chamber_volume=(
                self.motor.thrust_chamber.combustion_chamber.internal_volume
            ),
            propellant_mass=propellant_mass,
            propellant_mass_flow_rate=sum(mass_flows.values()),
            nozzle=nozzle,
            propellant_properties=propellant_properties,
            fluid_mass_per_line=fluid_masses,
            mass_flow_rate_per_line=mass_flows,
            oxidizer_to_fuel_ratio=oxidizer_to_fuel_ratio,
            tank_pressure_per_line=tank_pressures,
            tank_temperature_per_line=tank_temperatures,
        )
        self._apply_nozzle_losses(
            ideal_momentum_term,
            ideal_pressure_term,
            timestep_conditions,
            chamber_pressure,
        )

        if (
            not is_feeding
            or any(mass_consumed[name] >= fluid_masses[name] for name in self.lines)
        ) and not self.end_burn:
            self.end_burn = True
            self._burn_time = time + d_t

        if self._update_thrust_termination(
            time,
            chamber_pressure,
            external_pressure,
            propellant_properties.k_chamber,
        ):
            return

        new_time = time + d_t
        self.time.append(new_time)
        effective_flame_temperature = performance.get_effective_flame_temperature(
            adiabatic_flame_temperature=propellant_properties.adiabatic_flame_temperature,
            combustion_efficiency=self.motor.combustion_efficiency,
        )
        new_chamber_pressure = rk4.rk4th_ode_solver(
            variables={"chamber_pressure": chamber_pressure},
            equation=mass_balance.compute_chamber_pressure_mass_balance,
            d_t=d_t,
            external_pressure=external_pressure,
            mass_flow_in=functools.partial(
                get_total_injector_mass_flow, injector_flows=injector_flows
            ),
            free_chamber_volume=self.motor.thrust_chamber.combustion_chamber.internal_volume,
            throat_area=nozzle.get_throat_area(),
            k=propellant_properties.k_chamber,
            R=propellant_properties.R_chamber,
            flame_temperature=effective_flame_temperature,
            nozzle_discharge_coefficient=nozzle.discharge_coefficient,
        )[0]
        self.chamber_pressure.append(new_chamber_pressure)

        for name, line in self.lines.items():
            self.fluid_mass_per_line[name].append(
                fluid_masses[name] - mass_consumed[name]
            )
            self.internal_energy_per_line[name] = self._drain_internal_energy(
                tank=line.tank,
                internal_energy=self.internal_energy_per_line[name],
                fluid_mass=fluid_masses[name],
                mass_drained=mass_consumed[name],
            )

    @staticmethod
    def _drain_internal_energy(
        *,
        tank: tank_models.Tank,
        internal_energy: float | None,
        fluid_mass: float,
        mass_drained: float,
    ) -> float | None:
        """
        Take the enthalpy the drained fluid carries out of the tank [J].

        The tank is adiabatic and does no work on anything but the fluid it
        pushes out, so its internal energy falls by the enthalpy of what left.
        The fluid behind boils to refill the ullage and cools doing it, which
        is what walks the saturation pressure down over the burn.

        Returns:
            The internal energy left in the tank [J], or None for an
            isothermal tank, which runs no energy balance.
        """
        if internal_energy is None or mass_drained <= 0.0:
            return internal_energy

        return internal_energy - mass_drained * tank.get_outflow_specific_enthalpy(
            fluid_mass, internal_energy
        )
__init__(motor, igniter_pressure, external_pressure)

Initialize a biliquid engine state.

Parameters:

Name Type Description Default
motor BiliquidEngine

Biliquid engine to track.

 required
igniter_pressure float

Initial chamber pressure from the igniter [Pa].

 required
external_pressure float

Ambient pressure [Pa].

 required
Source code in machwave/simulation/biliquid/states.py 
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
def __init__(
    self,
    motor: motors.BiliquidEngine,
    igniter_pressure: float,
    external_pressure: float,
) -> None:
    """
    Initialize a biliquid engine state.

    Args:
        motor: Biliquid engine to track.
        igniter_pressure: Initial chamber pressure from the igniter [Pa].
        external_pressure: Ambient pressure [Pa].
    """
    super().__init__(
        motor=motor,
        igniter_pressure=igniter_pressure,
        external_pressure=external_pressure,
    )

    feed_system = motor.feed_system
    self.lines: dict[str, line_models.PropellantLine] = feed_system.lines
    self.oxidizer_line_names = tuple(
        line.name
        for line in feed_system.get_lines_with_role(
            propellant_components.ComponentRole.OXIDIZER
        )
    )
    self.fuel_line_names = tuple(
        line.name
        for line in feed_system.get_lines_with_role(
            propellant_components.ComponentRole.FUEL
        )
    )

    initial_states = {name: line.initial_state for name, line in self.lines.items()}
    self.fluid_mass_per_line: dict[str, simulation_states.SimulationStateArray] = {
        name: [state.fluid_mass] for name, state in initial_states.items()
    }

    # Second state variable of a tank running an energy balance, which the
    # integrator carries beside the mass. An isothermal tank has none.
    self.internal_energy_per_line: dict[str, float | None] = {
        name: state.internal_energy for name, state in initial_states.items()
    }

    self.mass_flow_rate_per_line: dict[
        str, simulation_states.SimulationStateArray
    ] = {name: [] for name in self.lines}
    self.tank_pressure_per_line: dict[
        str, simulation_states.SimulationStateArray
    ] = {name: [] for name in self.lines}
    self.tank_temperature_per_line: dict[
        str, simulation_states.SimulationStateArray
    ] = {name: [] for name in self.lines}
    self.oxidizer_to_fuel_ratio: simulation_states.SimulationStateArray = []
run_timestep(d_t, external_pressure)

Iterate the engine operation by calculating operational parameters.

Parameters:

Name Type Description Default
d_t float

Time increment [s].

 required
external_pressure float

External pressure [Pa].

 required
Source code in machwave/simulation/biliquid/states.py 
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
def run_timestep(
    self,
    d_t: float,
    external_pressure: float,
) -> None:
    """
    Iterate the engine operation by calculating operational parameters.

    Args:
        d_t: Time increment [s].
        external_pressure: External pressure [Pa].
    """
    nozzle = self.motor.thrust_chamber.nozzle
    feed_system = self.motor.feed_system

    time = self.time[-1]
    chamber_pressure = self.chamber_pressure[-1]
    fluid_masses = {
        name: series[-1] for name, series in self.fluid_mass_per_line.items()
    }

    propellant_mass = sum(fluid_masses.values())
    self.propellant_mass.append(propellant_mass)

    inlet_states = feed_system.get_inlet_states(
        {
            name: line_models.LineState(
                fluid_mass=fluid_masses[name],
                internal_energy=self.internal_energy_per_line[name],
            )
            for name in self.lines
        }
    )

    tank_pressures = {name: inlet.pressure for name, inlet in inlet_states.items()}
    tank_temperatures = {
        name: inlet.temperature for name, inlet in inlet_states.items()
    }
    for name in self.lines:
        self.tank_pressure_per_line[name].append(tank_pressures[name])
        self.tank_temperature_per_line[name].append(tank_temperatures[name])

    is_feeding = (
        not self.end_burn
        and all(fluid_mass > 0 for fluid_mass in fluid_masses.values())
        and all(pressure > chamber_pressure for pressure in tank_pressures.values())
    )
    injector_flows = functools.partial(
        get_injector_mass_flows,
        injector=self.motor.thrust_chamber.injector,
        inlet_states=inlet_states,
        line_masses=fluid_masses,
        is_feeding=is_feeding,
        d_t=d_t,
    )
    mass_flows = injector_flows(chamber_pressure)
    for name, mass_flow in mass_flows.items():
        self.mass_flow_rate_per_line[name].append(mass_flow)
    mass_consumed = {name: flow * d_t for name, flow in mass_flows.items()}

    mixture_ratio = self._get_mixture_ratio(mass_flows)
    oxidizer_to_fuel_ratio = math.nan if mixture_ratio is None else mixture_ratio
    self.oxidizer_to_fuel_ratio.append(oxidizer_to_fuel_ratio)

    propellant_properties = self._evaluate_propellant_properties(
        chamber_pressure=chamber_pressure,
        mixture_ratio=mixture_ratio,
    )

    (
        effective_expansion_ratio,
        exit_pressure,
        ideal_momentum_term,
        ideal_pressure_term,
    ) = self._ideal_thrust_coefficient_terms(
        propellant_properties.k_exhaust, chamber_pressure, external_pressure
    )

    timestep_conditions = BiliquidTimestepConditions(
        time=time,
        chamber_pressure=chamber_pressure,
        external_pressure=external_pressure,
        exit_pressure=exit_pressure,
        effective_expansion_ratio=effective_expansion_ratio,
        free_chamber_volume=(
            self.motor.thrust_chamber.combustion_chamber.internal_volume
        ),
        propellant_mass=propellant_mass,
        propellant_mass_flow_rate=sum(mass_flows.values()),
        nozzle=nozzle,
        propellant_properties=propellant_properties,
        fluid_mass_per_line=fluid_masses,
        mass_flow_rate_per_line=mass_flows,
        oxidizer_to_fuel_ratio=oxidizer_to_fuel_ratio,
        tank_pressure_per_line=tank_pressures,
        tank_temperature_per_line=tank_temperatures,
    )
    self._apply_nozzle_losses(
        ideal_momentum_term,
        ideal_pressure_term,
        timestep_conditions,
        chamber_pressure,
    )

    if (
        not is_feeding
        or any(mass_consumed[name] >= fluid_masses[name] for name in self.lines)
    ) and not self.end_burn:
        self.end_burn = True
        self._burn_time = time + d_t

    if self._update_thrust_termination(
        time,
        chamber_pressure,
        external_pressure,
        propellant_properties.k_chamber,
    ):
        return

    new_time = time + d_t
    self.time.append(new_time)
    effective_flame_temperature = performance.get_effective_flame_temperature(
        adiabatic_flame_temperature=propellant_properties.adiabatic_flame_temperature,
        combustion_efficiency=self.motor.combustion_efficiency,
    )
    new_chamber_pressure = rk4.rk4th_ode_solver(
        variables={"chamber_pressure": chamber_pressure},
        equation=mass_balance.compute_chamber_pressure_mass_balance,
        d_t=d_t,
        external_pressure=external_pressure,
        mass_flow_in=functools.partial(
            get_total_injector_mass_flow, injector_flows=injector_flows
        ),
        free_chamber_volume=self.motor.thrust_chamber.combustion_chamber.internal_volume,
        throat_area=nozzle.get_throat_area(),
        k=propellant_properties.k_chamber,
        R=propellant_properties.R_chamber,
        flame_temperature=effective_flame_temperature,
        nozzle_discharge_coefficient=nozzle.discharge_coefficient,
    )[0]
    self.chamber_pressure.append(new_chamber_pressure)

    for name, line in self.lines.items():
        self.fluid_mass_per_line[name].append(
            fluid_masses[name] - mass_consumed[name]
        )
        self.internal_energy_per_line[name] = self._drain_internal_energy(
            tank=line.tank,
            internal_energy=self.internal_energy_per_line[name],
            fluid_mass=fluid_masses[name],
            mass_drained=mass_consumed[name],
        )

BiliquidSimulationResult dataclass

Bases: SimulationResult['biliquid_states.BiliquidEngineState']

Simulation result for a biliquid engine run.

Source code in machwave/simulation/biliquid/results.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
@dataclass(frozen=True, kw_only=True)
class BiliquidSimulationResult(
    simulation_results.SimulationResult["biliquid_states.BiliquidEngineState"]
):
    """Simulation result for a biliquid engine run."""

    fluid_mass_per_line: dict[str, simulation_results.SimulationResultArray]
    tank_pressure_per_line: dict[str, simulation_results.SimulationResultArray]
    # Flat for an isothermal tank, decaying for one running an energy balance.
    tank_temperature_per_line: dict[str, simulation_results.SimulationResultArray]
    final_fluid_mass_per_line: dict[str, float]

    @classmethod
    def _collect_extra_fields(
        cls, state: "biliquid_states.BiliquidEngineState"
    ) -> dict[str, Any]:
        fluid_mass_per_line = {
            name: np.asarray(series)
            for name, series in state.fluid_mass_per_line.items()
        }
        return {
            "fluid_mass_per_line": fluid_mass_per_line,
            "tank_pressure_per_line": {
                name: np.asarray(series)
                for name, series in state.tank_pressure_per_line.items()
            },
            "tank_temperature_per_line": {
                name: np.asarray(series)
                for name, series in state.tank_temperature_per_line.items()
            },
            "final_fluid_mass_per_line": {
                name: float(series[-1]) for name, series in fluid_mass_per_line.items()
            },
        }

    def _extra_summary(self) -> dict[str, float]:
        return {
            f"final_{name}_mass": mass
            for name, mass in self.final_fluid_mass_per_line.items()
        }

    def _report_body(self, file: IO) -> None:
        print("\nBILIQUID ENGINE OPERATION RESULTS", file=file)

        print(f"Initial propellant mass: {self.propellant_mass[0]:.4f} kg", file=file)
        print(f"Burnout time: {self._format_burn_time(decimals=4)}", file=file)
        print(f"Thrust time: {self.thrust_time:.4f} s", file=file)

        print("\nCHAMBER PRESSURE (MPa)", file=file)
        print(f"  Max: {np.max(self.chamber_pressure) * 1e-6:.4f}", file=file)
        print(f"  Mean: {np.mean(self.chamber_pressure) * 1e-6:.4f}", file=file)

        print("\nTHRUST (N)", file=file)
        print(f"  Max: {np.max(self.thrust):.4f}", file=file)
        print(f"  Mean: {np.mean(self.thrust):.4f}", file=file)

        print("\nIMPULSE AND I_SP", file=file)
        print(f"  Total impulse: {self.total_impulse:.4f} N-s", file=file)
        print(f"  Specific impulse: {self.specific_impulse:.4f} s", file=file)

        self._report_nozzle_losses(file)

        print("\nPROPELLANT REMAINING (kg)", file=file)
        label_width = max(len(name) for name in self.final_fluid_mass_per_line)
        for name, mass in self.final_fluid_mass_per_line.items():
            print(
                f"  {name.capitalize() + ':':<{label_width + 1}} {mass:.4f}", file=file
            )

results

BiliquidSimulationResult dataclass

Bases: SimulationResult['biliquid_states.BiliquidEngineState']

Simulation result for a biliquid engine run.

Source code in machwave/simulation/biliquid/results.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
@dataclass(frozen=True, kw_only=True)
class BiliquidSimulationResult(
    simulation_results.SimulationResult["biliquid_states.BiliquidEngineState"]
):
    """Simulation result for a biliquid engine run."""

    fluid_mass_per_line: dict[str, simulation_results.SimulationResultArray]
    tank_pressure_per_line: dict[str, simulation_results.SimulationResultArray]
    # Flat for an isothermal tank, decaying for one running an energy balance.
    tank_temperature_per_line: dict[str, simulation_results.SimulationResultArray]
    final_fluid_mass_per_line: dict[str, float]

    @classmethod
    def _collect_extra_fields(
        cls, state: "biliquid_states.BiliquidEngineState"
    ) -> dict[str, Any]:
        fluid_mass_per_line = {
            name: np.asarray(series)
            for name, series in state.fluid_mass_per_line.items()
        }
        return {
            "fluid_mass_per_line": fluid_mass_per_line,
            "tank_pressure_per_line": {
                name: np.asarray(series)
                for name, series in state.tank_pressure_per_line.items()
            },
            "tank_temperature_per_line": {
                name: np.asarray(series)
                for name, series in state.tank_temperature_per_line.items()
            },
            "final_fluid_mass_per_line": {
                name: float(series[-1]) for name, series in fluid_mass_per_line.items()
            },
        }

    def _extra_summary(self) -> dict[str, float]:
        return {
            f"final_{name}_mass": mass
            for name, mass in self.final_fluid_mass_per_line.items()
        }

    def _report_body(self, file: IO) -> None:
        print("\nBILIQUID ENGINE OPERATION RESULTS", file=file)

        print(f"Initial propellant mass: {self.propellant_mass[0]:.4f} kg", file=file)
        print(f"Burnout time: {self._format_burn_time(decimals=4)}", file=file)
        print(f"Thrust time: {self.thrust_time:.4f} s", file=file)

        print("\nCHAMBER PRESSURE (MPa)", file=file)
        print(f"  Max: {np.max(self.chamber_pressure) * 1e-6:.4f}", file=file)
        print(f"  Mean: {np.mean(self.chamber_pressure) * 1e-6:.4f}", file=file)

        print("\nTHRUST (N)", file=file)
        print(f"  Max: {np.max(self.thrust):.4f}", file=file)
        print(f"  Mean: {np.mean(self.thrust):.4f}", file=file)

        print("\nIMPULSE AND I_SP", file=file)
        print(f"  Total impulse: {self.total_impulse:.4f} N-s", file=file)
        print(f"  Specific impulse: {self.specific_impulse:.4f} s", file=file)

        self._report_nozzle_losses(file)

        print("\nPROPELLANT REMAINING (kg)", file=file)
        label_width = max(len(name) for name in self.final_fluid_mass_per_line)
        for name, mass in self.final_fluid_mass_per_line.items():
            print(
                f"  {name.capitalize() + ':':<{label_width + 1}} {mass:.4f}", file=file
            )

states

BiliquidEngineState

Bases: MotorState

State for a biliquid rocket engine.

Source code in machwave/simulation/biliquid/states.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
class BiliquidEngineState(simulation_states.MotorState):
    """State for a biliquid rocket engine."""

    motor: motors.BiliquidEngine
    result_class = biliquid_results.BiliquidSimulationResult

    def __init__(
        self,
        motor: motors.BiliquidEngine,
        igniter_pressure: float,
        external_pressure: float,
    ) -> None:
        """
        Initialize a biliquid engine state.

        Args:
            motor: Biliquid engine to track.
            igniter_pressure: Initial chamber pressure from the igniter [Pa].
            external_pressure: Ambient pressure [Pa].
        """
        super().__init__(
            motor=motor,
            igniter_pressure=igniter_pressure,
            external_pressure=external_pressure,
        )

        feed_system = motor.feed_system
        self.lines: dict[str, line_models.PropellantLine] = feed_system.lines
        self.oxidizer_line_names = tuple(
            line.name
            for line in feed_system.get_lines_with_role(
                propellant_components.ComponentRole.OXIDIZER
            )
        )
        self.fuel_line_names = tuple(
            line.name
            for line in feed_system.get_lines_with_role(
                propellant_components.ComponentRole.FUEL
            )
        )

        initial_states = {name: line.initial_state for name, line in self.lines.items()}
        self.fluid_mass_per_line: dict[str, simulation_states.SimulationStateArray] = {
            name: [state.fluid_mass] for name, state in initial_states.items()
        }

        # Second state variable of a tank running an energy balance, which the
        # integrator carries beside the mass. An isothermal tank has none.
        self.internal_energy_per_line: dict[str, float | None] = {
            name: state.internal_energy for name, state in initial_states.items()
        }

        self.mass_flow_rate_per_line: dict[
            str, simulation_states.SimulationStateArray
        ] = {name: [] for name in self.lines}
        self.tank_pressure_per_line: dict[
            str, simulation_states.SimulationStateArray
        ] = {name: [] for name in self.lines}
        self.tank_temperature_per_line: dict[
            str, simulation_states.SimulationStateArray
        ] = {name: [] for name in self.lines}
        self.oxidizer_to_fuel_ratio: simulation_states.SimulationStateArray = []

    def _evaluate_propellant_properties(
        self,
        chamber_pressure: float,
        mixture_ratio: float | None,
    ) -> propellant_properties_models.ThermochemicalProperties:
        """Evaluate the propellant at the chamber pressure and mixture ratio."""
        return self.motor.propellant.evaluate(
            chamber_pressure=chamber_pressure,
            expansion_ratio=self.motor.thrust_chamber.nozzle.expansion_ratio,
            mixture_ratio=mixture_ratio,
        )

    def _get_mixture_ratio(self, mass_flows: Mapping[str, float]) -> float | None:
        """
        Total oxidizer flow over total fuel flow.

        None where the ratio has no meaning: a monoliquid, or an engine whose
        lines have stopped flowing.
        """
        oxidizer_flow = sum(mass_flows[name] for name in self.oxidizer_line_names)
        fuel_flow = sum(mass_flows[name] for name in self.fuel_line_names)

        if oxidizer_flow <= 0.0 or fuel_flow <= 0.0:
            return None

        return oxidizer_flow / fuel_flow

    def run_timestep(
        self,
        d_t: float,
        external_pressure: float,
    ) -> None:
        """
        Iterate the engine operation by calculating operational parameters.

        Args:
            d_t: Time increment [s].
            external_pressure: External pressure [Pa].
        """
        nozzle = self.motor.thrust_chamber.nozzle
        feed_system = self.motor.feed_system

        time = self.time[-1]
        chamber_pressure = self.chamber_pressure[-1]
        fluid_masses = {
            name: series[-1] for name, series in self.fluid_mass_per_line.items()
        }

        propellant_mass = sum(fluid_masses.values())
        self.propellant_mass.append(propellant_mass)

        inlet_states = feed_system.get_inlet_states(
            {
                name: line_models.LineState(
                    fluid_mass=fluid_masses[name],
                    internal_energy=self.internal_energy_per_line[name],
                )
                for name in self.lines
            }
        )

        tank_pressures = {name: inlet.pressure for name, inlet in inlet_states.items()}
        tank_temperatures = {
            name: inlet.temperature for name, inlet in inlet_states.items()
        }
        for name in self.lines:
            self.tank_pressure_per_line[name].append(tank_pressures[name])
            self.tank_temperature_per_line[name].append(tank_temperatures[name])

        is_feeding = (
            not self.end_burn
            and all(fluid_mass > 0 for fluid_mass in fluid_masses.values())
            and all(pressure > chamber_pressure for pressure in tank_pressures.values())
        )
        injector_flows = functools.partial(
            get_injector_mass_flows,
            injector=self.motor.thrust_chamber.injector,
            inlet_states=inlet_states,
            line_masses=fluid_masses,
            is_feeding=is_feeding,
            d_t=d_t,
        )
        mass_flows = injector_flows(chamber_pressure)
        for name, mass_flow in mass_flows.items():
            self.mass_flow_rate_per_line[name].append(mass_flow)
        mass_consumed = {name: flow * d_t for name, flow in mass_flows.items()}

        mixture_ratio = self._get_mixture_ratio(mass_flows)
        oxidizer_to_fuel_ratio = math.nan if mixture_ratio is None else mixture_ratio
        self.oxidizer_to_fuel_ratio.append(oxidizer_to_fuel_ratio)

        propellant_properties = self._evaluate_propellant_properties(
            chamber_pressure=chamber_pressure,
            mixture_ratio=mixture_ratio,
        )

        (
            effective_expansion_ratio,
            exit_pressure,
            ideal_momentum_term,
            ideal_pressure_term,
        ) = self._ideal_thrust_coefficient_terms(
            propellant_properties.k_exhaust, chamber_pressure, external_pressure
        )

        timestep_conditions = BiliquidTimestepConditions(
            time=time,
            chamber_pressure=chamber_pressure,
            external_pressure=external_pressure,
            exit_pressure=exit_pressure,
            effective_expansion_ratio=effective_expansion_ratio,
            free_chamber_volume=(
                self.motor.thrust_chamber.combustion_chamber.internal_volume
            ),
            propellant_mass=propellant_mass,
            propellant_mass_flow_rate=sum(mass_flows.values()),
            nozzle=nozzle,
            propellant_properties=propellant_properties,
            fluid_mass_per_line=fluid_masses,
            mass_flow_rate_per_line=mass_flows,
            oxidizer_to_fuel_ratio=oxidizer_to_fuel_ratio,
            tank_pressure_per_line=tank_pressures,
            tank_temperature_per_line=tank_temperatures,
        )
        self._apply_nozzle_losses(
            ideal_momentum_term,
            ideal_pressure_term,
            timestep_conditions,
            chamber_pressure,
        )

        if (
            not is_feeding
            or any(mass_consumed[name] >= fluid_masses[name] for name in self.lines)
        ) and not self.end_burn:
            self.end_burn = True
            self._burn_time = time + d_t

        if self._update_thrust_termination(
            time,
            chamber_pressure,
            external_pressure,
            propellant_properties.k_chamber,
        ):
            return

        new_time = time + d_t
        self.time.append(new_time)
        effective_flame_temperature = performance.get_effective_flame_temperature(
            adiabatic_flame_temperature=propellant_properties.adiabatic_flame_temperature,
            combustion_efficiency=self.motor.combustion_efficiency,
        )
        new_chamber_pressure = rk4.rk4th_ode_solver(
            variables={"chamber_pressure": chamber_pressure},
            equation=mass_balance.compute_chamber_pressure_mass_balance,
            d_t=d_t,
            external_pressure=external_pressure,
            mass_flow_in=functools.partial(
                get_total_injector_mass_flow, injector_flows=injector_flows
            ),
            free_chamber_volume=self.motor.thrust_chamber.combustion_chamber.internal_volume,
            throat_area=nozzle.get_throat_area(),
            k=propellant_properties.k_chamber,
            R=propellant_properties.R_chamber,
            flame_temperature=effective_flame_temperature,
            nozzle_discharge_coefficient=nozzle.discharge_coefficient,
        )[0]
        self.chamber_pressure.append(new_chamber_pressure)

        for name, line in self.lines.items():
            self.fluid_mass_per_line[name].append(
                fluid_masses[name] - mass_consumed[name]
            )
            self.internal_energy_per_line[name] = self._drain_internal_energy(
                tank=line.tank,
                internal_energy=self.internal_energy_per_line[name],
                fluid_mass=fluid_masses[name],
                mass_drained=mass_consumed[name],
            )

    @staticmethod
    def _drain_internal_energy(
        *,
        tank: tank_models.Tank,
        internal_energy: float | None,
        fluid_mass: float,
        mass_drained: float,
    ) -> float | None:
        """
        Take the enthalpy the drained fluid carries out of the tank [J].

        The tank is adiabatic and does no work on anything but the fluid it
        pushes out, so its internal energy falls by the enthalpy of what left.
        The fluid behind boils to refill the ullage and cools doing it, which
        is what walks the saturation pressure down over the burn.

        Returns:
            The internal energy left in the tank [J], or None for an
            isothermal tank, which runs no energy balance.
        """
        if internal_energy is None or mass_drained <= 0.0:
            return internal_energy

        return internal_energy - mass_drained * tank.get_outflow_specific_enthalpy(
            fluid_mass, internal_energy
        )
__init__(motor, igniter_pressure, external_pressure)

Initialize a biliquid engine state.

Parameters:

Name Type Description Default
motor BiliquidEngine

Biliquid engine to track.

 required
igniter_pressure float

Initial chamber pressure from the igniter [Pa].

 required
external_pressure float

Ambient pressure [Pa].

 required
Source code in machwave/simulation/biliquid/states.py 
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
def __init__(
    self,
    motor: motors.BiliquidEngine,
    igniter_pressure: float,
    external_pressure: float,
) -> None:
    """
    Initialize a biliquid engine state.

    Args:
        motor: Biliquid engine to track.
        igniter_pressure: Initial chamber pressure from the igniter [Pa].
        external_pressure: Ambient pressure [Pa].
    """
    super().__init__(
        motor=motor,
        igniter_pressure=igniter_pressure,
        external_pressure=external_pressure,
    )

    feed_system = motor.feed_system
    self.lines: dict[str, line_models.PropellantLine] = feed_system.lines
    self.oxidizer_line_names = tuple(
        line.name
        for line in feed_system.get_lines_with_role(
            propellant_components.ComponentRole.OXIDIZER
        )
    )
    self.fuel_line_names = tuple(
        line.name
        for line in feed_system.get_lines_with_role(
            propellant_components.ComponentRole.FUEL
        )
    )

    initial_states = {name: line.initial_state for name, line in self.lines.items()}
    self.fluid_mass_per_line: dict[str, simulation_states.SimulationStateArray] = {
        name: [state.fluid_mass] for name, state in initial_states.items()
    }

    # Second state variable of a tank running an energy balance, which the
    # integrator carries beside the mass. An isothermal tank has none.
    self.internal_energy_per_line: dict[str, float | None] = {
        name: state.internal_energy for name, state in initial_states.items()
    }

    self.mass_flow_rate_per_line: dict[
        str, simulation_states.SimulationStateArray
    ] = {name: [] for name in self.lines}
    self.tank_pressure_per_line: dict[
        str, simulation_states.SimulationStateArray
    ] = {name: [] for name in self.lines}
    self.tank_temperature_per_line: dict[
        str, simulation_states.SimulationStateArray
    ] = {name: [] for name in self.lines}
    self.oxidizer_to_fuel_ratio: simulation_states.SimulationStateArray = []
run_timestep(d_t, external_pressure)

Iterate the engine operation by calculating operational parameters.

Parameters:

Name Type Description Default
d_t float

Time increment [s].

 required
external_pressure float

External pressure [Pa].

 required
Source code in machwave/simulation/biliquid/states.py 
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
def run_timestep(
    self,
    d_t: float,
    external_pressure: float,
) -> None:
    """
    Iterate the engine operation by calculating operational parameters.

    Args:
        d_t: Time increment [s].
        external_pressure: External pressure [Pa].
    """
    nozzle = self.motor.thrust_chamber.nozzle
    feed_system = self.motor.feed_system

    time = self.time[-1]
    chamber_pressure = self.chamber_pressure[-1]
    fluid_masses = {
        name: series[-1] for name, series in self.fluid_mass_per_line.items()
    }

    propellant_mass = sum(fluid_masses.values())
    self.propellant_mass.append(propellant_mass)

    inlet_states = feed_system.get_inlet_states(
        {
            name: line_models.LineState(
                fluid_mass=fluid_masses[name],
                internal_energy=self.internal_energy_per_line[name],
            )
            for name in self.lines
        }
    )

    tank_pressures = {name: inlet.pressure for name, inlet in inlet_states.items()}
    tank_temperatures = {
        name: inlet.temperature for name, inlet in inlet_states.items()
    }
    for name in self.lines:
        self.tank_pressure_per_line[name].append(tank_pressures[name])
        self.tank_temperature_per_line[name].append(tank_temperatures[name])

    is_feeding = (
        not self.end_burn
        and all(fluid_mass > 0 for fluid_mass in fluid_masses.values())
        and all(pressure > chamber_pressure for pressure in tank_pressures.values())
    )
    injector_flows = functools.partial(
        get_injector_mass_flows,
        injector=self.motor.thrust_chamber.injector,
        inlet_states=inlet_states,
        line_masses=fluid_masses,
        is_feeding=is_feeding,
        d_t=d_t,
    )
    mass_flows = injector_flows(chamber_pressure)
    for name, mass_flow in mass_flows.items():
        self.mass_flow_rate_per_line[name].append(mass_flow)
    mass_consumed = {name: flow * d_t for name, flow in mass_flows.items()}

    mixture_ratio = self._get_mixture_ratio(mass_flows)
    oxidizer_to_fuel_ratio = math.nan if mixture_ratio is None else mixture_ratio
    self.oxidizer_to_fuel_ratio.append(oxidizer_to_fuel_ratio)

    propellant_properties = self._evaluate_propellant_properties(
        chamber_pressure=chamber_pressure,
        mixture_ratio=mixture_ratio,
    )

    (
        effective_expansion_ratio,
        exit_pressure,
        ideal_momentum_term,
        ideal_pressure_term,
    ) = self._ideal_thrust_coefficient_terms(
        propellant_properties.k_exhaust, chamber_pressure, external_pressure
    )

    timestep_conditions = BiliquidTimestepConditions(
        time=time,
        chamber_pressure=chamber_pressure,
        external_pressure=external_pressure,
        exit_pressure=exit_pressure,
        effective_expansion_ratio=effective_expansion_ratio,
        free_chamber_volume=(
            self.motor.thrust_chamber.combustion_chamber.internal_volume
        ),
        propellant_mass=propellant_mass,
        propellant_mass_flow_rate=sum(mass_flows.values()),
        nozzle=nozzle,
        propellant_properties=propellant_properties,
        fluid_mass_per_line=fluid_masses,
        mass_flow_rate_per_line=mass_flows,
        oxidizer_to_fuel_ratio=oxidizer_to_fuel_ratio,
        tank_pressure_per_line=tank_pressures,
        tank_temperature_per_line=tank_temperatures,
    )
    self._apply_nozzle_losses(
        ideal_momentum_term,
        ideal_pressure_term,
        timestep_conditions,
        chamber_pressure,
    )

    if (
        not is_feeding
        or any(mass_consumed[name] >= fluid_masses[name] for name in self.lines)
    ) and not self.end_burn:
        self.end_burn = True
        self._burn_time = time + d_t

    if self._update_thrust_termination(
        time,
        chamber_pressure,
        external_pressure,
        propellant_properties.k_chamber,
    ):
        return

    new_time = time + d_t
    self.time.append(new_time)
    effective_flame_temperature = performance.get_effective_flame_temperature(
        adiabatic_flame_temperature=propellant_properties.adiabatic_flame_temperature,
        combustion_efficiency=self.motor.combustion_efficiency,
    )
    new_chamber_pressure = rk4.rk4th_ode_solver(
        variables={"chamber_pressure": chamber_pressure},
        equation=mass_balance.compute_chamber_pressure_mass_balance,
        d_t=d_t,
        external_pressure=external_pressure,
        mass_flow_in=functools.partial(
            get_total_injector_mass_flow, injector_flows=injector_flows
        ),
        free_chamber_volume=self.motor.thrust_chamber.combustion_chamber.internal_volume,
        throat_area=nozzle.get_throat_area(),
        k=propellant_properties.k_chamber,
        R=propellant_properties.R_chamber,
        flame_temperature=effective_flame_temperature,
        nozzle_discharge_coefficient=nozzle.discharge_coefficient,
    )[0]
    self.chamber_pressure.append(new_chamber_pressure)

    for name, line in self.lines.items():
        self.fluid_mass_per_line[name].append(
            fluid_masses[name] - mass_consumed[name]
        )
        self.internal_energy_per_line[name] = self._drain_internal_energy(
            tank=line.tank,
            internal_energy=self.internal_energy_per_line[name],
            fluid_mass=fluid_masses[name],
            mass_drained=mass_consumed[name],
        )
BiliquidTimestepConditions dataclass

Bases: TimestepConditions

Timestep conditions for a biliquid engine.

Source code in machwave/simulation/biliquid/states.py 
55
56
57
58
59
60
61
62
63
@dataclasses.dataclass(frozen=True, kw_only=True)
class BiliquidTimestepConditions(simulation_states.TimestepConditions):
    """Timestep conditions for a biliquid engine."""

    fluid_mass_per_line: dict[str, float]
    mass_flow_rate_per_line: dict[str, float]
    oxidizer_to_fuel_ratio: float
    tank_pressure_per_line: dict[str, float]
    tank_temperature_per_line: dict[str, float]
get_injector_mass_flows(chamber_pressure, *, injector, inlet_states, line_masses, is_feeding, d_t)

Injector flow on every line at the given chamber pressure [kg/s].

No line delivers more than the mass it has left over the timestep.

Source code in machwave/simulation/biliquid/states.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def get_injector_mass_flows(
    chamber_pressure: float,
    *,
    injector: injector_models.Injector,
    inlet_states: Mapping[str, fluid_state_models.FluidState],
    line_masses: Mapping[str, float],
    is_feeding: bool,
    d_t: float,
) -> dict[str, float]:
    """
    Injector flow on every line at the given chamber pressure [kg/s].

    No line delivers more than the mass it has left over the timestep.
    """
    if not is_feeding:
        return {name: 0.0 for name in inlet_states}

    flows = injector.get_mass_flows(
        inlet_states=inlet_states, chamber_pressure=chamber_pressure
    )
    return {name: min(flow, line_masses[name] / d_t) for name, flow in flows.items()}
get_total_injector_mass_flow(chamber_pressure, *, injector_flows)

Total injector mass flow over every line at the given pressure [kg/s].

Source code in machwave/simulation/biliquid/states.py 
46
47
48
49
50
51
52
def get_total_injector_mass_flow(
    chamber_pressure: float,
    *,
    injector_flows: Callable[[float], Mapping[str, float]],
) -> float:
    """Total injector mass flow over every line at the given pressure [kg/s]."""
    return sum(injector_flows(chamber_pressure).values())

results

SimulationResult dataclass

Bases: ABC, Generic[StateT]

Results of a finished internal ballistics simulation.

Source code in machwave/simulation/results.py 
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
@dataclass(frozen=True, kw_only=True)
class SimulationResult(ABC, Generic[StateT]):
    """Results of a finished internal ballistics simulation."""

    time: SimulationResultArray
    propellant_mass: SimulationResultArray
    chamber_pressure: SimulationResultArray
    exit_pressure: SimulationResultArray
    thrust_coefficient: SimulationResultArray
    ideal_thrust_coefficient: SimulationResultArray
    thrust: SimulationResultArray
    nozzle_efficiency: SimulationResultArray
    loss_fractions: dict[str, SimulationResultArray]
    loss_labels: dict[str, str]
    # None when the run terminated with propellant remaining.
    burn_time: float | None
    thrust_time: float
    end_thrust: bool
    end_burn: bool
    initial_propellant_mass: float
    total_impulse: float
    specific_impulse: float

    @classmethod
    def from_state(cls, state: StateT) -> "SimulationResult":
        """Build a `SimulationResult` from a finished motor state."""
        return cls(
            **cls._collect_base_fields(state),
            **cls._collect_extra_fields(state),
        )

    @classmethod
    def _collect_base_fields(cls, state: StateT) -> dict[str, Any]:
        time = np.asarray(state.time)
        thrust = np.asarray(state.thrust)
        total_impulse = performance.get_total_impulse(thrust, time)
        initial_propellant_mass = state.motor.initial_propellant_mass
        return {
            "time": time,
            "propellant_mass": np.asarray(state.propellant_mass),
            "chamber_pressure": np.asarray(state.chamber_pressure),
            "exit_pressure": np.asarray(state.exit_pressure),
            "thrust_coefficient": np.asarray(state.thrust_coefficient),
            "ideal_thrust_coefficient": np.asarray(state.ideal_thrust_coefficient),
            "thrust": thrust,
            "nozzle_efficiency": np.asarray(state.nozzle_efficiency),
            "loss_fractions": {
                name: np.asarray(series)
                for name, series in state.loss_fractions.items()
            },
            "loss_labels": dict(state.motor.nozzle_loss_model.component_labels),
            "burn_time": state.burn_time,
            "thrust_time": state.thrust_time,
            "end_thrust": state.end_thrust,
            "end_burn": state.end_burn,
            "initial_propellant_mass": initial_propellant_mass,
            "total_impulse": total_impulse,
            "specific_impulse": performance.get_specific_impulse(
                total_impulse=total_impulse,
                initial_propellant_mass=initial_propellant_mass,
            ),
        }

    @classmethod
    @abstractmethod
    def _collect_extra_fields(cls, state: StateT) -> dict[str, Any]:
        """Build constructor kwargs for fields specific to the subclass."""

    def report(self, file: IO = sys.stdout) -> None:
        """Print a human-readable report of the simulation result."""
        print("\nINTERNAL BALLISTICS SIMULATION RESULTS", file=file)
        self._report_body(file)

    @abstractmethod
    def _report_body(self, file: IO) -> None:
        """Print the subclass-specific portion of the report."""

    def _report_nozzle_losses(self, file: IO) -> None:
        """Print the nozzle efficiency and each loss component's mean fraction."""
        print("\nNOZZLE", file=file)
        print(
            f"  Average nozzle efficiency: {np.mean(self.nozzle_efficiency):.3%}",
            file=file,
        )
        for name, series in self.loss_fractions.items():
            label = self.loss_labels[name]
            print(f"  Average {label} fraction: {np.mean(series):.3%}", file=file)

    def _format_burn_time(self, decimals: int = 3) -> str:
        """Format the burn time for a report, or flag it as never reached."""
        if self.burn_time is None:
            return "not reached"
        return f"{self.burn_time:.{decimals}f} s"

    def summary(self) -> dict[str, float | None]:
        """Return a mapping of headline scalar metrics for this result."""
        return {
            "burn_time": self.burn_time,
            "thrust_time": self.thrust_time,
            "initial_propellant_mass": self.initial_propellant_mass,
            "total_impulse": self.total_impulse,
            "specific_impulse": self.specific_impulse,
            "peak_chamber_pressure": float(np.max(self.chamber_pressure)),
            "mean_chamber_pressure": float(np.mean(self.chamber_pressure)),
            "peak_thrust": float(np.max(self.thrust)),
            "mean_thrust": float(np.mean(self.thrust)),
            **self._extra_summary(),
        }

    def _extra_summary(self) -> dict[str, float]:
        """Return subclass-specific scalar metrics to merge into `summary()`."""
        return {}
from_state(state) classmethod

Build a SimulationResult from a finished motor state.

Source code in machwave/simulation/results.py 
42
43
44
45
46
47
48
@classmethod
def from_state(cls, state: StateT) -> "SimulationResult":
    """Build a `SimulationResult` from a finished motor state."""
    return cls(
        **cls._collect_base_fields(state),
        **cls._collect_extra_fields(state),
    )
report(file=sys.stdout)

Print a human-readable report of the simulation result.

Source code in machwave/simulation/results.py 
87
88
89
90
def report(self, file: IO = sys.stdout) -> None:
    """Print a human-readable report of the simulation result."""
    print("\nINTERNAL BALLISTICS SIMULATION RESULTS", file=file)
    self._report_body(file)
summary()

Return a mapping of headline scalar metrics for this result.

Source code in machwave/simulation/results.py 
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def summary(self) -> dict[str, float | None]:
    """Return a mapping of headline scalar metrics for this result."""
    return {
        "burn_time": self.burn_time,
        "thrust_time": self.thrust_time,
        "initial_propellant_mass": self.initial_propellant_mass,
        "total_impulse": self.total_impulse,
        "specific_impulse": self.specific_impulse,
        "peak_chamber_pressure": float(np.max(self.chamber_pressure)),
        "mean_chamber_pressure": float(np.mean(self.chamber_pressure)),
        "peak_thrust": float(np.max(self.thrust)),
        "mean_thrust": float(np.mean(self.thrust)),
        **self._extra_summary(),
    }

solid

SolidMotorState

Bases: MotorState

State for a Solid Rocket Motor.

Source code in machwave/simulation/solid/states.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
class SolidMotorState(simulation_states.MotorState):
    """State for a Solid Rocket Motor."""

    motor: motors.SolidMotor
    result_class = solid_results.SolidSimulationResult

    def __init__(
        self,
        motor: motors.SolidMotor,
        igniter_pressure: float,
        external_pressure: float,
    ) -> None:
        """
        Initialize a solid motor state.

        Args:
            motor: Solid motor to track.
            igniter_pressure: Initial chamber pressure from the igniter [Pa].
            external_pressure: Ambient pressure [Pa].
        """
        super().__init__(
            motor=motor,
            igniter_pressure=igniter_pressure,
            external_pressure=external_pressure,
        )

        self.segment_density_ratios = motor.grain.get_density_ratio_per_segment()

        self.web: simulation_states.SimulationStateArray = [0.0]

        self.burn_area: simulation_states.SimulationStateArray = []
        self.propellant_volume: simulation_states.SimulationStateArray = []
        self.burn_area_per_segment: list[npt.NDArray[np.float64]] = []
        self.propellant_volume_per_segment: list[npt.NDArray[np.float64]] = []
        self.propellant_mass_per_segment: list[npt.NDArray[np.float64]] = []
        self.burn_rate: simulation_states.SimulationStateArray = []
        self.free_chamber_volume: simulation_states.SimulationStateArray = []
        self.free_chamber_volume_rate: simulation_states.SimulationStateArray = []
        self.grain_segment_mass_flow: list[npt.NDArray[np.float64]] = []

        self.propellant_cog: list[npt.NDArray[np.float64]] = []
        self.propellant_moi: list[npt.NDArray[np.float64]] = []

    def _evaluate_propellant_properties(
        self,
        chamber_pressure: float,
    ) -> propellant_properties_models.ThermochemicalProperties:
        """Evaluate the propellant at the chamber pressure."""
        return self.motor.propellant.evaluate(
            chamber_pressure=chamber_pressure,
            expansion_ratio=self.motor.thrust_chamber.nozzle.expansion_ratio,
        )

    def run_timestep(
        self,
        d_t: float,
        external_pressure: float,
    ) -> None:
        """
        Iterate the motor operation by calculating operational parameters.

        Args:
            d_t: Time increment [s].
            external_pressure: External pressure [Pa].
        """
        nozzle = self.motor.thrust_chamber.nozzle
        ideal_propellant_density = self.motor.propellant.ideal_density

        time = self.time[-1]
        web_distance = self.web[-1]
        chamber_pressure = self.chamber_pressure[-1]

        propellant_properties = self._evaluate_propellant_properties(
            chamber_pressure=chamber_pressure
        )

        burn_area_per_segment = self.motor.grain.get_burn_area_per_segment(web_distance)
        self.burn_area_per_segment.append(burn_area_per_segment)
        burn_area = float(np.sum(burn_area_per_segment))
        self.burn_area.append(burn_area)

        # Pressure-dependent inflow: recorded at start-of-step, re-evaluated per stage.
        mass_flow_per_segment = functools.partial(
            get_grain_mass_flow_per_segment,
            propellant=self.motor.propellant,
            burn_area_per_segment=burn_area_per_segment,
            segment_density_ratios=self.segment_density_ratios,
        )
        free_chamber_volume_rate = functools.partial(
            get_free_chamber_volume_rate,
            propellant=self.motor.propellant,
            burn_area=burn_area,
        )

        propellant_volume_per_segment = (
            self.motor.grain.get_propellant_volume_per_segment(web_distance)
        )
        self.propellant_volume_per_segment.append(propellant_volume_per_segment)
        propellant_volume = float(np.sum(propellant_volume_per_segment))
        self.propellant_volume.append(propellant_volume)

        burn_rate = self.motor.propellant.get_burn_rate(chamber_pressure)
        self.burn_rate.append(burn_rate)
        web_consumed = burn_rate * d_t

        free_chamber_volume = self.motor.get_free_chamber_volume(propellant_volume)
        self.free_chamber_volume.append(free_chamber_volume)
        self.free_chamber_volume_rate.append(free_chamber_volume_rate(chamber_pressure))
        propellant_mass_per_segment = (
            propellant_volume_per_segment
            * self.segment_density_ratios
            * ideal_propellant_density
        )
        self.propellant_mass_per_segment.append(propellant_mass_per_segment)
        propellant_mass = float(np.sum(propellant_mass_per_segment))
        self.propellant_mass.append(propellant_mass)

        if propellant_mass > 0.0:
            propellant_cog = self.motor.grain.get_center_of_gravity(
                web_distance=web_distance,
                volume_per_segment=propellant_volume_per_segment,
            )
            propellant_moi = self.motor.grain.get_moment_of_inertia(
                ideal_density=ideal_propellant_density,
                web_distance=web_distance,
                volume_per_segment=propellant_volume_per_segment,
                center_of_gravity=propellant_cog,
            )
            # Shift after the parallel axis step, which needs the grain frame.
            propellant_cog[0] += (
                self.motor.thrust_chamber.nozzle_exit_to_grain_port_distance
            )
        else:  # no propellant left: CoG is undefined and inertia is zero
            propellant_cog = np.full(3, np.nan)
            propellant_moi = np.zeros((3, 3))
        self.propellant_cog.append(propellant_cog)
        self.propellant_moi.append(propellant_moi)

        self.grain_segment_mass_flow.append(mass_flow_per_segment(chamber_pressure))

        (
            effective_expansion_ratio,
            exit_pressure,
            ideal_momentum_term,
            ideal_pressure_term,
        ) = self._ideal_thrust_coefficient_terms(
            propellant_properties.k_exhaust, chamber_pressure, external_pressure
        )

        timestep_conditions = SolidTimestepConditions(
            time=time,
            chamber_pressure=chamber_pressure,
            external_pressure=external_pressure,
            exit_pressure=exit_pressure,
            effective_expansion_ratio=effective_expansion_ratio,
            free_chamber_volume=free_chamber_volume,
            propellant_mass=propellant_mass,
            propellant_mass_flow_rate=float(np.sum(self.grain_segment_mass_flow[-1])),
            nozzle=nozzle,
            propellant_properties=propellant_properties,
            burn_area=burn_area,
            burn_rate=burn_rate,
            propellant_volume=propellant_volume,
            web_distance=web_distance,
            free_chamber_volume_rate=self.free_chamber_volume_rate[-1],
        )
        self._apply_nozzle_losses(
            ideal_momentum_term,
            ideal_pressure_term,
            timestep_conditions,
            chamber_pressure,
        )

        if propellant_mass <= 0 and not self.end_burn:
            self._burn_time = time
            self.end_burn = True

        if self._update_thrust_termination(
            time,
            chamber_pressure,
            external_pressure,
            propellant_properties.k_chamber,
        ):
            return

        new_time = time + d_t
        self.time.append(new_time)
        new_web_distance = web_distance + web_consumed
        self.web.append(new_web_distance)
        effective_flame_temperature = performance.get_effective_flame_temperature(
            adiabatic_flame_temperature=propellant_properties.adiabatic_flame_temperature,
            combustion_efficiency=self.motor.combustion_efficiency,
        )

        new_chamber_pressure = rk4.rk4th_ode_solver(
            variables={"chamber_pressure": self.chamber_pressure[-1]},
            equation=mass_balance.compute_chamber_pressure_mass_balance,
            d_t=d_t,
            external_pressure=external_pressure,
            mass_flow_in=functools.partial(
                get_grain_mass_flow, mass_flow_per_segment=mass_flow_per_segment
            ),
            free_chamber_volume=free_chamber_volume,
            throat_area=nozzle.get_throat_area(),
            k=propellant_properties.k_chamber,
            R=propellant_properties.R_chamber,
            flame_temperature=effective_flame_temperature,
            nozzle_discharge_coefficient=nozzle.discharge_coefficient,
            free_chamber_volume_rate=free_chamber_volume_rate,
        )[0]
        self.chamber_pressure.append(new_chamber_pressure)
__init__(motor, igniter_pressure, external_pressure)

Initialize a solid motor state.

Parameters:

Name Type Description Default
motor SolidMotor

Solid motor to track.

 required
igniter_pressure float

Initial chamber pressure from the igniter [Pa].

 required
external_pressure float

Ambient pressure [Pa].

 required
Source code in machwave/simulation/solid/states.py 
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def __init__(
    self,
    motor: motors.SolidMotor,
    igniter_pressure: float,
    external_pressure: float,
) -> None:
    """
    Initialize a solid motor state.

    Args:
        motor: Solid motor to track.
        igniter_pressure: Initial chamber pressure from the igniter [Pa].
        external_pressure: Ambient pressure [Pa].
    """
    super().__init__(
        motor=motor,
        igniter_pressure=igniter_pressure,
        external_pressure=external_pressure,
    )

    self.segment_density_ratios = motor.grain.get_density_ratio_per_segment()

    self.web: simulation_states.SimulationStateArray = [0.0]

    self.burn_area: simulation_states.SimulationStateArray = []
    self.propellant_volume: simulation_states.SimulationStateArray = []
    self.burn_area_per_segment: list[npt.NDArray[np.float64]] = []
    self.propellant_volume_per_segment: list[npt.NDArray[np.float64]] = []
    self.propellant_mass_per_segment: list[npt.NDArray[np.float64]] = []
    self.burn_rate: simulation_states.SimulationStateArray = []
    self.free_chamber_volume: simulation_states.SimulationStateArray = []
    self.free_chamber_volume_rate: simulation_states.SimulationStateArray = []
    self.grain_segment_mass_flow: list[npt.NDArray[np.float64]] = []

    self.propellant_cog: list[npt.NDArray[np.float64]] = []
    self.propellant_moi: list[npt.NDArray[np.float64]] = []
run_timestep(d_t, external_pressure)

Iterate the motor operation by calculating operational parameters.

Parameters:

Name Type Description Default
d_t float

Time increment [s].

 required
external_pressure float

External pressure [Pa].

 required
Source code in machwave/simulation/solid/states.py 
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
def run_timestep(
    self,
    d_t: float,
    external_pressure: float,
) -> None:
    """
    Iterate the motor operation by calculating operational parameters.

    Args:
        d_t: Time increment [s].
        external_pressure: External pressure [Pa].
    """
    nozzle = self.motor.thrust_chamber.nozzle
    ideal_propellant_density = self.motor.propellant.ideal_density

    time = self.time[-1]
    web_distance = self.web[-1]
    chamber_pressure = self.chamber_pressure[-1]

    propellant_properties = self._evaluate_propellant_properties(
        chamber_pressure=chamber_pressure
    )

    burn_area_per_segment = self.motor.grain.get_burn_area_per_segment(web_distance)
    self.burn_area_per_segment.append(burn_area_per_segment)
    burn_area = float(np.sum(burn_area_per_segment))
    self.burn_area.append(burn_area)

    # Pressure-dependent inflow: recorded at start-of-step, re-evaluated per stage.
    mass_flow_per_segment = functools.partial(
        get_grain_mass_flow_per_segment,
        propellant=self.motor.propellant,
        burn_area_per_segment=burn_area_per_segment,
        segment_density_ratios=self.segment_density_ratios,
    )
    free_chamber_volume_rate = functools.partial(
        get_free_chamber_volume_rate,
        propellant=self.motor.propellant,
        burn_area=burn_area,
    )

    propellant_volume_per_segment = (
        self.motor.grain.get_propellant_volume_per_segment(web_distance)
    )
    self.propellant_volume_per_segment.append(propellant_volume_per_segment)
    propellant_volume = float(np.sum(propellant_volume_per_segment))
    self.propellant_volume.append(propellant_volume)

    burn_rate = self.motor.propellant.get_burn_rate(chamber_pressure)
    self.burn_rate.append(burn_rate)
    web_consumed = burn_rate * d_t

    free_chamber_volume = self.motor.get_free_chamber_volume(propellant_volume)
    self.free_chamber_volume.append(free_chamber_volume)
    self.free_chamber_volume_rate.append(free_chamber_volume_rate(chamber_pressure))
    propellant_mass_per_segment = (
        propellant_volume_per_segment
        * self.segment_density_ratios
        * ideal_propellant_density
    )
    self.propellant_mass_per_segment.append(propellant_mass_per_segment)
    propellant_mass = float(np.sum(propellant_mass_per_segment))
    self.propellant_mass.append(propellant_mass)

    if propellant_mass > 0.0:
        propellant_cog = self.motor.grain.get_center_of_gravity(
            web_distance=web_distance,
            volume_per_segment=propellant_volume_per_segment,
        )
        propellant_moi = self.motor.grain.get_moment_of_inertia(
            ideal_density=ideal_propellant_density,
            web_distance=web_distance,
            volume_per_segment=propellant_volume_per_segment,
            center_of_gravity=propellant_cog,
        )
        # Shift after the parallel axis step, which needs the grain frame.
        propellant_cog[0] += (
            self.motor.thrust_chamber.nozzle_exit_to_grain_port_distance
        )
    else:  # no propellant left: CoG is undefined and inertia is zero
        propellant_cog = np.full(3, np.nan)
        propellant_moi = np.zeros((3, 3))
    self.propellant_cog.append(propellant_cog)
    self.propellant_moi.append(propellant_moi)

    self.grain_segment_mass_flow.append(mass_flow_per_segment(chamber_pressure))

    (
        effective_expansion_ratio,
        exit_pressure,
        ideal_momentum_term,
        ideal_pressure_term,
    ) = self._ideal_thrust_coefficient_terms(
        propellant_properties.k_exhaust, chamber_pressure, external_pressure
    )

    timestep_conditions = SolidTimestepConditions(
        time=time,
        chamber_pressure=chamber_pressure,
        external_pressure=external_pressure,
        exit_pressure=exit_pressure,
        effective_expansion_ratio=effective_expansion_ratio,
        free_chamber_volume=free_chamber_volume,
        propellant_mass=propellant_mass,
        propellant_mass_flow_rate=float(np.sum(self.grain_segment_mass_flow[-1])),
        nozzle=nozzle,
        propellant_properties=propellant_properties,
        burn_area=burn_area,
        burn_rate=burn_rate,
        propellant_volume=propellant_volume,
        web_distance=web_distance,
        free_chamber_volume_rate=self.free_chamber_volume_rate[-1],
    )
    self._apply_nozzle_losses(
        ideal_momentum_term,
        ideal_pressure_term,
        timestep_conditions,
        chamber_pressure,
    )

    if propellant_mass <= 0 and not self.end_burn:
        self._burn_time = time
        self.end_burn = True

    if self._update_thrust_termination(
        time,
        chamber_pressure,
        external_pressure,
        propellant_properties.k_chamber,
    ):
        return

    new_time = time + d_t
    self.time.append(new_time)
    new_web_distance = web_distance + web_consumed
    self.web.append(new_web_distance)
    effective_flame_temperature = performance.get_effective_flame_temperature(
        adiabatic_flame_temperature=propellant_properties.adiabatic_flame_temperature,
        combustion_efficiency=self.motor.combustion_efficiency,
    )

    new_chamber_pressure = rk4.rk4th_ode_solver(
        variables={"chamber_pressure": self.chamber_pressure[-1]},
        equation=mass_balance.compute_chamber_pressure_mass_balance,
        d_t=d_t,
        external_pressure=external_pressure,
        mass_flow_in=functools.partial(
            get_grain_mass_flow, mass_flow_per_segment=mass_flow_per_segment
        ),
        free_chamber_volume=free_chamber_volume,
        throat_area=nozzle.get_throat_area(),
        k=propellant_properties.k_chamber,
        R=propellant_properties.R_chamber,
        flame_temperature=effective_flame_temperature,
        nozzle_discharge_coefficient=nozzle.discharge_coefficient,
        free_chamber_volume_rate=free_chamber_volume_rate,
    )[0]
    self.chamber_pressure.append(new_chamber_pressure)

SolidSimulationResult dataclass

Bases: SimulationResult['solid_states.SolidMotorState']

Simulation result for a solid motor run.

Source code in machwave/simulation/solid/results.py 
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
@dataclass(frozen=True, kw_only=True)
class SolidSimulationResult(
    simulation_results.SimulationResult["solid_states.SolidMotorState"]
):
    """Simulation result for a solid motor run."""

    free_chamber_volume: simulation_results.SimulationResultArray
    free_chamber_volume_rate: simulation_results.SimulationResultArray
    web: simulation_results.SimulationResultArray
    burn_area: simulation_results.SimulationResultArray
    propellant_volume: simulation_results.SimulationResultArray
    burn_area_per_segment: simulation_results.SimulationResultArray
    propellant_volume_per_segment: simulation_results.SimulationResultArray
    propellant_mass_per_segment: simulation_results.SimulationResultArray
    burn_rate: simulation_results.SimulationResultArray
    # propellant_cog is measured from the nozzle exit, positive toward the bulkhead.
    propellant_cog: simulation_results.SimulationResultArray
    propellant_moi: simulation_results.SimulationResultArray
    # klemmung is filtered to burn_area > 0, so its length is < len(time).
    klemmung: simulation_results.SimulationResultArray = field(
        metadata={"non_aligned": True}
    )
    # grain_mass_flux is shaped [segment_count, time_count]; axis 0 is segments.
    grain_mass_flux: simulation_results.SimulationResultArray = field(
        metadata={"non_aligned": True}
    )
    initial_to_final_klemmung_ratio: float
    volumetric_efficiency: float
    burn_profile: str
    max_mass_flux: float

    @classmethod
    def _collect_extra_fields(
        cls, state: "solid_states.SolidMotorState"
    ) -> dict[str, Any]:
        motor = state.motor
        nozzle = motor.thrust_chamber.nozzle
        chamber_volume = motor.thrust_chamber.combustion_chamber.internal_volume

        burn_area = np.asarray(state.burn_area)
        propellant_volume = np.asarray(state.propellant_volume)
        burn_rate = np.asarray(state.burn_rate)
        web = np.asarray(state.web)

        klemmung = cls._get_klemmung(burn_area, nozzle.get_throat_area())
        grain_mass_flux = motor.grain.get_mass_flux_per_segment(
            burn_rate, motor.propellant.ideal_density, web
        )

        return {
            "free_chamber_volume": np.asarray(state.free_chamber_volume),
            "free_chamber_volume_rate": np.asarray(state.free_chamber_volume_rate),
            "web": web,
            "burn_area": burn_area,
            "propellant_volume": propellant_volume,
            "burn_area_per_segment": np.stack(state.burn_area_per_segment),
            "propellant_volume_per_segment": np.stack(
                state.propellant_volume_per_segment
            ),
            "propellant_mass_per_segment": np.stack(state.propellant_mass_per_segment),
            "burn_rate": burn_rate,
            "propellant_cog": np.stack(
                [np.asarray(cog) for cog in state.propellant_cog]
            ),
            "propellant_moi": np.stack(
                [np.asarray(moi) for moi in state.propellant_moi]
            ),
            "klemmung": klemmung,
            "grain_mass_flux": grain_mass_flux,
            "initial_to_final_klemmung_ratio": float(klemmung[0] / klemmung[-1]),
            "volumetric_efficiency": cls._get_volumetric_efficiency(
                propellant_volume[0], chamber_volume
            ),
            "burn_profile": cls._classify_burn_profile(klemmung),
            "max_mass_flux": float(np.max(grain_mass_flux)),
        }

    @staticmethod
    def _get_klemmung(
        burn_area: simulation_results.SimulationResultArray, throat_area: float
    ) -> simulation_results.SimulationResultArray:
        """
        Return Klemmung (Kn) over non-zero burn-area samples.

        Returns Kn values where burn area is positive; length is at most
        `len(burn_area)`.
        """
        return burn_area[burn_area > 0] / throat_area

    @staticmethod
    def _classify_burn_profile(
        klemmung: simulation_results.SimulationResultArray, deviancy: float = 0.02
    ) -> str:
        """
        Classify a burn profile as "regressive", "progressive", or "neutral".

        Args:
            klemmung: Klemmung samples.
            deviancy: Fractional threshold around 1.0 for the initial-to-final
                ratio that delimits the neutral band.
        """
        ratio = float(klemmung[0] / klemmung[-1])
        if ratio > 1 + deviancy:
            return "regressive"
        if ratio < 1 - deviancy:
            return "progressive"
        return "neutral"

    @staticmethod
    def _get_volumetric_efficiency(
        initial_propellant_volume: float, internal_chamber_volume: float
    ) -> float:
        """Fraction of the chamber initially occupied by propellant."""
        return float(initial_propellant_volume / internal_chamber_volume)

    def _extra_summary(self) -> dict[str, float]:
        return {
            "peak_klemmung": float(np.max(self.klemmung)),
            "mean_klemmung": float(np.mean(self.klemmung)),
            "initial_to_final_klemmung_ratio": self.initial_to_final_klemmung_ratio,
            "volumetric_efficiency": self.volumetric_efficiency,
            "max_mass_flux": self.max_mass_flux,
        }

    def _report_body(self, file: IO) -> None:
        print("\nBURN REGRESSION", file=file)
        if self.propellant_mass[0] > 1:
            print(
                f" Propellant initial mass {self.propellant_mass[0]:.3f} kg", file=file
            )
        else:
            print(
                f" Propellant initial mass {self.propellant_mass[0] * 1e3:.3f} g",
                file=file,
            )
        print(
            f" Initial propellant volume: {self.propellant_volume[0] * 1e6:.3f} cm^3",
            file=file,
        )
        print(f" Initial burn area: {self.burn_area[0] * 1e4:.3f} cm^2", file=file)
        print(
            f" Peak burn area: {float(np.max(self.burn_area)) * 1e4:.3f} cm^2",
            file=file,
        )
        print(" Mean Kn: %.2f" % np.mean(self.klemmung), file=file)
        print(" Max Kn: %.2f" % np.max(self.klemmung), file=file)
        print(
            f" Initial to final Kn ratio: {self.initial_to_final_klemmung_ratio:.3f}",
            file=file,
        )
        print(f" Volumetric efficiency: {self.volumetric_efficiency:.3%}", file=file)
        print(" Burn profile: " + self.burn_profile, file=file)
        print(
            f" Max initial mass flux: {self.max_mass_flux:.3f} kg/s-m-m or "
            f"{conversions.convert_mass_flux_metric_to_imperial(self.max_mass_flux):.3f} "
            "lb/s-in-in",
            file=file,
        )

        print("\nCHAMBER PRESSURE", file=file)
        print(
            f" Maximum, average chamber pressure: {np.max(self.chamber_pressure) * 1e-6:.3f}, "
            f"{np.mean(self.chamber_pressure) * 1e-6:.3f} MPa",
            file=file,
        )

        print("\nTHRUST AND IMPULSE", file=file)
        print(
            f" Maximum, average thrust: {np.max(self.thrust):.3f}, {np.mean(self.thrust):.3f} N",
            file=file,
        )
        print(
            f" Total, specific impulses: {self.total_impulse:.3f} N-s, {self.specific_impulse:.3f} s",
            file=file,
        )
        print(
            f" Burnout time: {self._format_burn_time()}, "
            f"thrust time: {self.thrust_time:.3f} s",
            file=file,
        )

        self._report_nozzle_losses(file)

results

SolidSimulationResult dataclass

Bases: SimulationResult['solid_states.SolidMotorState']

Simulation result for a solid motor run.

Source code in machwave/simulation/solid/results.py 
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
@dataclass(frozen=True, kw_only=True)
class SolidSimulationResult(
    simulation_results.SimulationResult["solid_states.SolidMotorState"]
):
    """Simulation result for a solid motor run."""

    free_chamber_volume: simulation_results.SimulationResultArray
    free_chamber_volume_rate: simulation_results.SimulationResultArray
    web: simulation_results.SimulationResultArray
    burn_area: simulation_results.SimulationResultArray
    propellant_volume: simulation_results.SimulationResultArray
    burn_area_per_segment: simulation_results.SimulationResultArray
    propellant_volume_per_segment: simulation_results.SimulationResultArray
    propellant_mass_per_segment: simulation_results.SimulationResultArray
    burn_rate: simulation_results.SimulationResultArray
    # propellant_cog is measured from the nozzle exit, positive toward the bulkhead.
    propellant_cog: simulation_results.SimulationResultArray
    propellant_moi: simulation_results.SimulationResultArray
    # klemmung is filtered to burn_area > 0, so its length is < len(time).
    klemmung: simulation_results.SimulationResultArray = field(
        metadata={"non_aligned": True}
    )
    # grain_mass_flux is shaped [segment_count, time_count]; axis 0 is segments.
    grain_mass_flux: simulation_results.SimulationResultArray = field(
        metadata={"non_aligned": True}
    )
    initial_to_final_klemmung_ratio: float
    volumetric_efficiency: float
    burn_profile: str
    max_mass_flux: float

    @classmethod
    def _collect_extra_fields(
        cls, state: "solid_states.SolidMotorState"
    ) -> dict[str, Any]:
        motor = state.motor
        nozzle = motor.thrust_chamber.nozzle
        chamber_volume = motor.thrust_chamber.combustion_chamber.internal_volume

        burn_area = np.asarray(state.burn_area)
        propellant_volume = np.asarray(state.propellant_volume)
        burn_rate = np.asarray(state.burn_rate)
        web = np.asarray(state.web)

        klemmung = cls._get_klemmung(burn_area, nozzle.get_throat_area())
        grain_mass_flux = motor.grain.get_mass_flux_per_segment(
            burn_rate, motor.propellant.ideal_density, web
        )

        return {
            "free_chamber_volume": np.asarray(state.free_chamber_volume),
            "free_chamber_volume_rate": np.asarray(state.free_chamber_volume_rate),
            "web": web,
            "burn_area": burn_area,
            "propellant_volume": propellant_volume,
            "burn_area_per_segment": np.stack(state.burn_area_per_segment),
            "propellant_volume_per_segment": np.stack(
                state.propellant_volume_per_segment
            ),
            "propellant_mass_per_segment": np.stack(state.propellant_mass_per_segment),
            "burn_rate": burn_rate,
            "propellant_cog": np.stack(
                [np.asarray(cog) for cog in state.propellant_cog]
            ),
            "propellant_moi": np.stack(
                [np.asarray(moi) for moi in state.propellant_moi]
            ),
            "klemmung": klemmung,
            "grain_mass_flux": grain_mass_flux,
            "initial_to_final_klemmung_ratio": float(klemmung[0] / klemmung[-1]),
            "volumetric_efficiency": cls._get_volumetric_efficiency(
                propellant_volume[0], chamber_volume
            ),
            "burn_profile": cls._classify_burn_profile(klemmung),
            "max_mass_flux": float(np.max(grain_mass_flux)),
        }

    @staticmethod
    def _get_klemmung(
        burn_area: simulation_results.SimulationResultArray, throat_area: float
    ) -> simulation_results.SimulationResultArray:
        """
        Return Klemmung (Kn) over non-zero burn-area samples.

        Returns Kn values where burn area is positive; length is at most
        `len(burn_area)`.
        """
        return burn_area[burn_area > 0] / throat_area

    @staticmethod
    def _classify_burn_profile(
        klemmung: simulation_results.SimulationResultArray, deviancy: float = 0.02
    ) -> str:
        """
        Classify a burn profile as "regressive", "progressive", or "neutral".

        Args:
            klemmung: Klemmung samples.
            deviancy: Fractional threshold around 1.0 for the initial-to-final
                ratio that delimits the neutral band.
        """
        ratio = float(klemmung[0] / klemmung[-1])
        if ratio > 1 + deviancy:
            return "regressive"
        if ratio < 1 - deviancy:
            return "progressive"
        return "neutral"

    @staticmethod
    def _get_volumetric_efficiency(
        initial_propellant_volume: float, internal_chamber_volume: float
    ) -> float:
        """Fraction of the chamber initially occupied by propellant."""
        return float(initial_propellant_volume / internal_chamber_volume)

    def _extra_summary(self) -> dict[str, float]:
        return {
            "peak_klemmung": float(np.max(self.klemmung)),
            "mean_klemmung": float(np.mean(self.klemmung)),
            "initial_to_final_klemmung_ratio": self.initial_to_final_klemmung_ratio,
            "volumetric_efficiency": self.volumetric_efficiency,
            "max_mass_flux": self.max_mass_flux,
        }

    def _report_body(self, file: IO) -> None:
        print("\nBURN REGRESSION", file=file)
        if self.propellant_mass[0] > 1:
            print(
                f" Propellant initial mass {self.propellant_mass[0]:.3f} kg", file=file
            )
        else:
            print(
                f" Propellant initial mass {self.propellant_mass[0] * 1e3:.3f} g",
                file=file,
            )
        print(
            f" Initial propellant volume: {self.propellant_volume[0] * 1e6:.3f} cm^3",
            file=file,
        )
        print(f" Initial burn area: {self.burn_area[0] * 1e4:.3f} cm^2", file=file)
        print(
            f" Peak burn area: {float(np.max(self.burn_area)) * 1e4:.3f} cm^2",
            file=file,
        )
        print(" Mean Kn: %.2f" % np.mean(self.klemmung), file=file)
        print(" Max Kn: %.2f" % np.max(self.klemmung), file=file)
        print(
            f" Initial to final Kn ratio: {self.initial_to_final_klemmung_ratio:.3f}",
            file=file,
        )
        print(f" Volumetric efficiency: {self.volumetric_efficiency:.3%}", file=file)
        print(" Burn profile: " + self.burn_profile, file=file)
        print(
            f" Max initial mass flux: {self.max_mass_flux:.3f} kg/s-m-m or "
            f"{conversions.convert_mass_flux_metric_to_imperial(self.max_mass_flux):.3f} "
            "lb/s-in-in",
            file=file,
        )

        print("\nCHAMBER PRESSURE", file=file)
        print(
            f" Maximum, average chamber pressure: {np.max(self.chamber_pressure) * 1e-6:.3f}, "
            f"{np.mean(self.chamber_pressure) * 1e-6:.3f} MPa",
            file=file,
        )

        print("\nTHRUST AND IMPULSE", file=file)
        print(
            f" Maximum, average thrust: {np.max(self.thrust):.3f}, {np.mean(self.thrust):.3f} N",
            file=file,
        )
        print(
            f" Total, specific impulses: {self.total_impulse:.3f} N-s, {self.specific_impulse:.3f} s",
            file=file,
        )
        print(
            f" Burnout time: {self._format_burn_time()}, "
            f"thrust time: {self.thrust_time:.3f} s",
            file=file,
        )

        self._report_nozzle_losses(file)

states

SolidMotorState

Bases: MotorState

State for a Solid Rocket Motor.

Source code in machwave/simulation/solid/states.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
class SolidMotorState(simulation_states.MotorState):
    """State for a Solid Rocket Motor."""

    motor: motors.SolidMotor
    result_class = solid_results.SolidSimulationResult

    def __init__(
        self,
        motor: motors.SolidMotor,
        igniter_pressure: float,
        external_pressure: float,
    ) -> None:
        """
        Initialize a solid motor state.

        Args:
            motor: Solid motor to track.
            igniter_pressure: Initial chamber pressure from the igniter [Pa].
            external_pressure: Ambient pressure [Pa].
        """
        super().__init__(
            motor=motor,
            igniter_pressure=igniter_pressure,
            external_pressure=external_pressure,
        )

        self.segment_density_ratios = motor.grain.get_density_ratio_per_segment()

        self.web: simulation_states.SimulationStateArray = [0.0]

        self.burn_area: simulation_states.SimulationStateArray = []
        self.propellant_volume: simulation_states.SimulationStateArray = []
        self.burn_area_per_segment: list[npt.NDArray[np.float64]] = []
        self.propellant_volume_per_segment: list[npt.NDArray[np.float64]] = []
        self.propellant_mass_per_segment: list[npt.NDArray[np.float64]] = []
        self.burn_rate: simulation_states.SimulationStateArray = []
        self.free_chamber_volume: simulation_states.SimulationStateArray = []
        self.free_chamber_volume_rate: simulation_states.SimulationStateArray = []
        self.grain_segment_mass_flow: list[npt.NDArray[np.float64]] = []

        self.propellant_cog: list[npt.NDArray[np.float64]] = []
        self.propellant_moi: list[npt.NDArray[np.float64]] = []

    def _evaluate_propellant_properties(
        self,
        chamber_pressure: float,
    ) -> propellant_properties_models.ThermochemicalProperties:
        """Evaluate the propellant at the chamber pressure."""
        return self.motor.propellant.evaluate(
            chamber_pressure=chamber_pressure,
            expansion_ratio=self.motor.thrust_chamber.nozzle.expansion_ratio,
        )

    def run_timestep(
        self,
        d_t: float,
        external_pressure: float,
    ) -> None:
        """
        Iterate the motor operation by calculating operational parameters.

        Args:
            d_t: Time increment [s].
            external_pressure: External pressure [Pa].
        """
        nozzle = self.motor.thrust_chamber.nozzle
        ideal_propellant_density = self.motor.propellant.ideal_density

        time = self.time[-1]
        web_distance = self.web[-1]
        chamber_pressure = self.chamber_pressure[-1]

        propellant_properties = self._evaluate_propellant_properties(
            chamber_pressure=chamber_pressure
        )

        burn_area_per_segment = self.motor.grain.get_burn_area_per_segment(web_distance)
        self.burn_area_per_segment.append(burn_area_per_segment)
        burn_area = float(np.sum(burn_area_per_segment))
        self.burn_area.append(burn_area)

        # Pressure-dependent inflow: recorded at start-of-step, re-evaluated per stage.
        mass_flow_per_segment = functools.partial(
            get_grain_mass_flow_per_segment,
            propellant=self.motor.propellant,
            burn_area_per_segment=burn_area_per_segment,
            segment_density_ratios=self.segment_density_ratios,
        )
        free_chamber_volume_rate = functools.partial(
            get_free_chamber_volume_rate,
            propellant=self.motor.propellant,
            burn_area=burn_area,
        )

        propellant_volume_per_segment = (
            self.motor.grain.get_propellant_volume_per_segment(web_distance)
        )
        self.propellant_volume_per_segment.append(propellant_volume_per_segment)
        propellant_volume = float(np.sum(propellant_volume_per_segment))
        self.propellant_volume.append(propellant_volume)

        burn_rate = self.motor.propellant.get_burn_rate(chamber_pressure)
        self.burn_rate.append(burn_rate)
        web_consumed = burn_rate * d_t

        free_chamber_volume = self.motor.get_free_chamber_volume(propellant_volume)
        self.free_chamber_volume.append(free_chamber_volume)
        self.free_chamber_volume_rate.append(free_chamber_volume_rate(chamber_pressure))
        propellant_mass_per_segment = (
            propellant_volume_per_segment
            * self.segment_density_ratios
            * ideal_propellant_density
        )
        self.propellant_mass_per_segment.append(propellant_mass_per_segment)
        propellant_mass = float(np.sum(propellant_mass_per_segment))
        self.propellant_mass.append(propellant_mass)

        if propellant_mass > 0.0:
            propellant_cog = self.motor.grain.get_center_of_gravity(
                web_distance=web_distance,
                volume_per_segment=propellant_volume_per_segment,
            )
            propellant_moi = self.motor.grain.get_moment_of_inertia(
                ideal_density=ideal_propellant_density,
                web_distance=web_distance,
                volume_per_segment=propellant_volume_per_segment,
                center_of_gravity=propellant_cog,
            )
            # Shift after the parallel axis step, which needs the grain frame.
            propellant_cog[0] += (
                self.motor.thrust_chamber.nozzle_exit_to_grain_port_distance
            )
        else:  # no propellant left: CoG is undefined and inertia is zero
            propellant_cog = np.full(3, np.nan)
            propellant_moi = np.zeros((3, 3))
        self.propellant_cog.append(propellant_cog)
        self.propellant_moi.append(propellant_moi)

        self.grain_segment_mass_flow.append(mass_flow_per_segment(chamber_pressure))

        (
            effective_expansion_ratio,
            exit_pressure,
            ideal_momentum_term,
            ideal_pressure_term,
        ) = self._ideal_thrust_coefficient_terms(
            propellant_properties.k_exhaust, chamber_pressure, external_pressure
        )

        timestep_conditions = SolidTimestepConditions(
            time=time,
            chamber_pressure=chamber_pressure,
            external_pressure=external_pressure,
            exit_pressure=exit_pressure,
            effective_expansion_ratio=effective_expansion_ratio,
            free_chamber_volume=free_chamber_volume,
            propellant_mass=propellant_mass,
            propellant_mass_flow_rate=float(np.sum(self.grain_segment_mass_flow[-1])),
            nozzle=nozzle,
            propellant_properties=propellant_properties,
            burn_area=burn_area,
            burn_rate=burn_rate,
            propellant_volume=propellant_volume,
            web_distance=web_distance,
            free_chamber_volume_rate=self.free_chamber_volume_rate[-1],
        )
        self._apply_nozzle_losses(
            ideal_momentum_term,
            ideal_pressure_term,
            timestep_conditions,
            chamber_pressure,
        )

        if propellant_mass <= 0 and not self.end_burn:
            self._burn_time = time
            self.end_burn = True

        if self._update_thrust_termination(
            time,
            chamber_pressure,
            external_pressure,
            propellant_properties.k_chamber,
        ):
            return

        new_time = time + d_t
        self.time.append(new_time)
        new_web_distance = web_distance + web_consumed
        self.web.append(new_web_distance)
        effective_flame_temperature = performance.get_effective_flame_temperature(
            adiabatic_flame_temperature=propellant_properties.adiabatic_flame_temperature,
            combustion_efficiency=self.motor.combustion_efficiency,
        )

        new_chamber_pressure = rk4.rk4th_ode_solver(
            variables={"chamber_pressure": self.chamber_pressure[-1]},
            equation=mass_balance.compute_chamber_pressure_mass_balance,
            d_t=d_t,
            external_pressure=external_pressure,
            mass_flow_in=functools.partial(
                get_grain_mass_flow, mass_flow_per_segment=mass_flow_per_segment
            ),
            free_chamber_volume=free_chamber_volume,
            throat_area=nozzle.get_throat_area(),
            k=propellant_properties.k_chamber,
            R=propellant_properties.R_chamber,
            flame_temperature=effective_flame_temperature,
            nozzle_discharge_coefficient=nozzle.discharge_coefficient,
            free_chamber_volume_rate=free_chamber_volume_rate,
        )[0]
        self.chamber_pressure.append(new_chamber_pressure)
__init__(motor, igniter_pressure, external_pressure)

Initialize a solid motor state.

Parameters:

Name Type Description Default
motor SolidMotor

Solid motor to track.

 required
igniter_pressure float

Initial chamber pressure from the igniter [Pa].

 required
external_pressure float

Ambient pressure [Pa].

 required
Source code in machwave/simulation/solid/states.py 
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def __init__(
    self,
    motor: motors.SolidMotor,
    igniter_pressure: float,
    external_pressure: float,
) -> None:
    """
    Initialize a solid motor state.

    Args:
        motor: Solid motor to track.
        igniter_pressure: Initial chamber pressure from the igniter [Pa].
        external_pressure: Ambient pressure [Pa].
    """
    super().__init__(
        motor=motor,
        igniter_pressure=igniter_pressure,
        external_pressure=external_pressure,
    )

    self.segment_density_ratios = motor.grain.get_density_ratio_per_segment()

    self.web: simulation_states.SimulationStateArray = [0.0]

    self.burn_area: simulation_states.SimulationStateArray = []
    self.propellant_volume: simulation_states.SimulationStateArray = []
    self.burn_area_per_segment: list[npt.NDArray[np.float64]] = []
    self.propellant_volume_per_segment: list[npt.NDArray[np.float64]] = []
    self.propellant_mass_per_segment: list[npt.NDArray[np.float64]] = []
    self.burn_rate: simulation_states.SimulationStateArray = []
    self.free_chamber_volume: simulation_states.SimulationStateArray = []
    self.free_chamber_volume_rate: simulation_states.SimulationStateArray = []
    self.grain_segment_mass_flow: list[npt.NDArray[np.float64]] = []

    self.propellant_cog: list[npt.NDArray[np.float64]] = []
    self.propellant_moi: list[npt.NDArray[np.float64]] = []
run_timestep(d_t, external_pressure)

Iterate the motor operation by calculating operational parameters.

Parameters:

Name Type Description Default
d_t float

Time increment [s].

 required
external_pressure float

External pressure [Pa].

 required
Source code in machwave/simulation/solid/states.py 
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
def run_timestep(
    self,
    d_t: float,
    external_pressure: float,
) -> None:
    """
    Iterate the motor operation by calculating operational parameters.

    Args:
        d_t: Time increment [s].
        external_pressure: External pressure [Pa].
    """
    nozzle = self.motor.thrust_chamber.nozzle
    ideal_propellant_density = self.motor.propellant.ideal_density

    time = self.time[-1]
    web_distance = self.web[-1]
    chamber_pressure = self.chamber_pressure[-1]

    propellant_properties = self._evaluate_propellant_properties(
        chamber_pressure=chamber_pressure
    )

    burn_area_per_segment = self.motor.grain.get_burn_area_per_segment(web_distance)
    self.burn_area_per_segment.append(burn_area_per_segment)
    burn_area = float(np.sum(burn_area_per_segment))
    self.burn_area.append(burn_area)

    # Pressure-dependent inflow: recorded at start-of-step, re-evaluated per stage.
    mass_flow_per_segment = functools.partial(
        get_grain_mass_flow_per_segment,
        propellant=self.motor.propellant,
        burn_area_per_segment=burn_area_per_segment,
        segment_density_ratios=self.segment_density_ratios,
    )
    free_chamber_volume_rate = functools.partial(
        get_free_chamber_volume_rate,
        propellant=self.motor.propellant,
        burn_area=burn_area,
    )

    propellant_volume_per_segment = (
        self.motor.grain.get_propellant_volume_per_segment(web_distance)
    )
    self.propellant_volume_per_segment.append(propellant_volume_per_segment)
    propellant_volume = float(np.sum(propellant_volume_per_segment))
    self.propellant_volume.append(propellant_volume)

    burn_rate = self.motor.propellant.get_burn_rate(chamber_pressure)
    self.burn_rate.append(burn_rate)
    web_consumed = burn_rate * d_t

    free_chamber_volume = self.motor.get_free_chamber_volume(propellant_volume)
    self.free_chamber_volume.append(free_chamber_volume)
    self.free_chamber_volume_rate.append(free_chamber_volume_rate(chamber_pressure))
    propellant_mass_per_segment = (
        propellant_volume_per_segment
        * self.segment_density_ratios
        * ideal_propellant_density
    )
    self.propellant_mass_per_segment.append(propellant_mass_per_segment)
    propellant_mass = float(np.sum(propellant_mass_per_segment))
    self.propellant_mass.append(propellant_mass)

    if propellant_mass > 0.0:
        propellant_cog = self.motor.grain.get_center_of_gravity(
            web_distance=web_distance,
            volume_per_segment=propellant_volume_per_segment,
        )
        propellant_moi = self.motor.grain.get_moment_of_inertia(
            ideal_density=ideal_propellant_density,
            web_distance=web_distance,
            volume_per_segment=propellant_volume_per_segment,
            center_of_gravity=propellant_cog,
        )
        # Shift after the parallel axis step, which needs the grain frame.
        propellant_cog[0] += (
            self.motor.thrust_chamber.nozzle_exit_to_grain_port_distance
        )
    else:  # no propellant left: CoG is undefined and inertia is zero
        propellant_cog = np.full(3, np.nan)
        propellant_moi = np.zeros((3, 3))
    self.propellant_cog.append(propellant_cog)
    self.propellant_moi.append(propellant_moi)

    self.grain_segment_mass_flow.append(mass_flow_per_segment(chamber_pressure))

    (
        effective_expansion_ratio,
        exit_pressure,
        ideal_momentum_term,
        ideal_pressure_term,
    ) = self._ideal_thrust_coefficient_terms(
        propellant_properties.k_exhaust, chamber_pressure, external_pressure
    )

    timestep_conditions = SolidTimestepConditions(
        time=time,
        chamber_pressure=chamber_pressure,
        external_pressure=external_pressure,
        exit_pressure=exit_pressure,
        effective_expansion_ratio=effective_expansion_ratio,
        free_chamber_volume=free_chamber_volume,
        propellant_mass=propellant_mass,
        propellant_mass_flow_rate=float(np.sum(self.grain_segment_mass_flow[-1])),
        nozzle=nozzle,
        propellant_properties=propellant_properties,
        burn_area=burn_area,
        burn_rate=burn_rate,
        propellant_volume=propellant_volume,
        web_distance=web_distance,
        free_chamber_volume_rate=self.free_chamber_volume_rate[-1],
    )
    self._apply_nozzle_losses(
        ideal_momentum_term,
        ideal_pressure_term,
        timestep_conditions,
        chamber_pressure,
    )

    if propellant_mass <= 0 and not self.end_burn:
        self._burn_time = time
        self.end_burn = True

    if self._update_thrust_termination(
        time,
        chamber_pressure,
        external_pressure,
        propellant_properties.k_chamber,
    ):
        return

    new_time = time + d_t
    self.time.append(new_time)
    new_web_distance = web_distance + web_consumed
    self.web.append(new_web_distance)
    effective_flame_temperature = performance.get_effective_flame_temperature(
        adiabatic_flame_temperature=propellant_properties.adiabatic_flame_temperature,
        combustion_efficiency=self.motor.combustion_efficiency,
    )

    new_chamber_pressure = rk4.rk4th_ode_solver(
        variables={"chamber_pressure": self.chamber_pressure[-1]},
        equation=mass_balance.compute_chamber_pressure_mass_balance,
        d_t=d_t,
        external_pressure=external_pressure,
        mass_flow_in=functools.partial(
            get_grain_mass_flow, mass_flow_per_segment=mass_flow_per_segment
        ),
        free_chamber_volume=free_chamber_volume,
        throat_area=nozzle.get_throat_area(),
        k=propellant_properties.k_chamber,
        R=propellant_properties.R_chamber,
        flame_temperature=effective_flame_temperature,
        nozzle_discharge_coefficient=nozzle.discharge_coefficient,
        free_chamber_volume_rate=free_chamber_volume_rate,
    )[0]
    self.chamber_pressure.append(new_chamber_pressure)
SolidTimestepConditions dataclass

Bases: TimestepConditions

Timestep conditions for a solid motor.

Source code in machwave/simulation/solid/states.py 
55
56
57
58
59
60
61
62
63
@dataclasses.dataclass(frozen=True, kw_only=True)
class SolidTimestepConditions(simulation_states.TimestepConditions):
    """Timestep conditions for a solid motor."""

    burn_area: float
    burn_rate: float
    propellant_volume: float
    web_distance: float
    free_chamber_volume_rate: float
get_free_chamber_volume_rate(chamber_pressure, *, propellant, burn_area)

Free chamber volume growth rate at the given chamber pressure [m^3/s].

Source code in machwave/simulation/solid/states.py 
45
46
47
48
49
50
51
52
def get_free_chamber_volume_rate(
    chamber_pressure: float,
    *,
    propellant: propellants.SolidPropellant,
    burn_area: float,
) -> float:
    """Free chamber volume growth rate at the given chamber pressure [m^3/s]."""
    return propellant.get_burn_rate(chamber_pressure) * burn_area
get_grain_mass_flow(chamber_pressure, *, mass_flow_per_segment)

Total grain mass generation rate at the given chamber pressure [kg/s].

Source code in machwave/simulation/solid/states.py 
36
37
38
39
40
41
42
def get_grain_mass_flow(
    chamber_pressure: float,
    *,
    mass_flow_per_segment: Callable[[float], npt.NDArray[np.float64]],
) -> float:
    """Total grain mass generation rate at the given chamber pressure [kg/s]."""
    return float(np.sum(mass_flow_per_segment(chamber_pressure)))
get_grain_mass_flow_per_segment(chamber_pressure, *, propellant, burn_area_per_segment, segment_density_ratios)

Per-segment grain mass generation rate at the given chamber pressure [kg/s].

Source code in machwave/simulation/solid/states.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def get_grain_mass_flow_per_segment(
    chamber_pressure: float,
    *,
    propellant: propellants.SolidPropellant,
    burn_area_per_segment: npt.NDArray[np.float64],
    segment_density_ratios: npt.NDArray[np.float64],
) -> npt.NDArray[np.float64]:
    """Per-segment grain mass generation rate at the given chamber pressure [kg/s]."""
    return (
        propellant.ideal_density
        * propellant.get_burn_rate(chamber_pressure)
        * burn_area_per_segment
        * segment_density_ratios
    )

states

MotorState

Bases: ABC

Defines the states and iteration step for a motor operation.

Source code in machwave/simulation/states.py 
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
class MotorState(ABC):
    """Defines the states and iteration step for a motor operation."""

    result_class: ClassVar[type["simulation_results.SimulationResult"]]

    def __init__(
        self,
        motor: motors.Motor,
        igniter_pressure: float,
        external_pressure: float,
    ) -> None:
        """
        Initialize a motor state.

        Args:
            motor: Motor to track.
            igniter_pressure: Initial chamber pressure from the igniter [Pa].
            external_pressure: Ambient pressure [Pa].
        """
        self.motor = motor
        self.external_pressure = external_pressure

        self.time: SimulationStateArray = [0.0]
        self.chamber_pressure: SimulationStateArray = [igniter_pressure]

        self.propellant_mass: SimulationStateArray = []
        self.exit_pressure: SimulationStateArray = []
        self.ideal_thrust_coefficient: SimulationStateArray = []
        self.thrust_coefficient: SimulationStateArray = []
        self.thrust: SimulationStateArray = []
        self.peak_thrust: float = 0.0
        self.nozzle_efficiency: SimulationStateArray = []
        self.loss_fractions: dict[str, SimulationStateArray] = {
            name: [] for name in motor.nozzle_loss_model.component_names
        }

        self._thrust_time: float | None = None
        self._burn_time: float | None = None

        self.end_thrust: bool = False
        self.end_burn: bool = False

    @abstractmethod
    def run_timestep(self, *args, **kwargs) -> None:
        """Advance the per-step accumulators by one time increment."""

    def _ideal_thrust_coefficient_terms(
        self,
        k_exhaust: float,
        chamber_pressure: float,
        external_pressure: float,
    ) -> tuple[float, float, float, float]:
        """
        Resolve the separated exit conditions and ideal thrust coefficient terms.

        Appends the effective exit pressure and the ideal thrust coefficient for the
        timestep.

        Args:
            k_exhaust: Isentropic exponent at the nozzle exit.
            chamber_pressure: Chamber pressure [Pa].
            external_pressure: Ambient pressure [Pa].

        Returns:
            The effective expansion ratio, effective exit pressure [Pa], and the
            momentum and pressure terms of the ideal thrust coefficient.
        """
        nozzle = self.motor.thrust_chamber.nozzle
        effective_expansion_ratio, exit_pressure = (
            nozzle_core.get_separated_exit_conditions(
                k_exhaust,
                nozzle.expansion_ratio,
                chamber_pressure,
                external_pressure,
                nozzle.separation_pressure_ratio,
            )
        )
        self.exit_pressure.append(exit_pressure)

        ideal_momentum_term, ideal_pressure_term = (
            nozzle_core.get_ideal_thrust_coefficient_terms(
                chamber_pressure,
                exit_pressure,
                external_pressure,
                effective_expansion_ratio,
                k_exhaust,
            )
        )
        self.ideal_thrust_coefficient.append(ideal_momentum_term + ideal_pressure_term)
        return (
            effective_expansion_ratio,
            exit_pressure,
            ideal_momentum_term,
            ideal_pressure_term,
        )

    def _apply_nozzle_losses(
        self,
        ideal_momentum_term: float,
        ideal_pressure_term: float,
        timestep_conditions: TimestepConditions,
        chamber_pressure: float,
    ) -> None:
        """
        Derate the ideal thrust coefficient terms and record the loss outputs.

        Appends the realized nozzle efficiency, each component loss fraction, the
        corrected thrust coefficient, and the thrust for the timestep.

        Args:
            ideal_momentum_term: Momentum term of the ideal thrust coefficient.
            ideal_pressure_term: Pressure term of the ideal thrust coefficient.
            timestep_conditions: Engine conditions at a point in time.
            chamber_pressure: Chamber pressure [Pa].
        """
        nozzle = self.motor.thrust_chamber.nozzle
        loss_result = self.motor.nozzle_loss_model.evaluate(
            ideal_momentum_term, ideal_pressure_term, timestep_conditions
        )
        self.nozzle_efficiency.append(loss_result.nozzle_efficiency)
        for name, fraction in loss_result.loss_fractions.items():
            self.loss_fractions[name].append(fraction)

        thrust_coefficient = loss_result.momentum_term + loss_result.pressure_term
        self.thrust_coefficient.append(thrust_coefficient)
        thrust = nozzle_core.get_thrust_from_thrust_coefficient(
            thrust_coefficient, chamber_pressure, nozzle.get_throat_area()
        )
        self.thrust.append(thrust)
        self.peak_thrust = max(self.peak_thrust, thrust)

    def _update_thrust_termination(
        self,
        time: float,
        chamber_pressure: float,
        external_pressure: float,
        k_chamber: float,
    ) -> bool:
        """
        Flag thrust termination once the nozzle un-chokes or the tail-off ends.

        Loss of choking never fires against a vacuum, where the nozzle stays choked
        for any chamber pressure, so a burnt out motor also terminates once its
        thrust decays to ``TAIL_OFF_THRUST_FRACTION`` of the peak.

        Args:
            time: Time at the start of the timestep [s].
            chamber_pressure: Chamber pressure [Pa].
            external_pressure: Ambient pressure [Pa].
            k_chamber: Isentropic exponent in the chamber.

        Returns:
            True if thrust has terminated on this timestep.
        """
        is_choked = isentropic.is_flow_choked(
            chamber_pressure,
            external_pressure,
            isentropic.get_critical_pressure_ratio(k_chamber),
        )
        has_tailed_off = (
            self.end_burn
            and self.thrust[-1] < TAIL_OFF_THRUST_FRACTION * self.peak_thrust
        )
        if is_choked and not has_tailed_off:
            return False

        self._thrust_time = time
        self.end_thrust = True
        return True

    def build_result(self) -> "simulation_results.SimulationResult":
        """Return a frozen ``SimulationResult`` snapshot of this state."""
        return self.result_class.from_state(self)

    @property
    def initial_propellant_mass(self) -> float:
        """Return the initial propellant mass [kg]."""
        return self.motor.initial_propellant_mass

    @property
    def thrust_time(self) -> float:
        """
        Return the thrust time [s].

        Raises:
            ValueError: If the simulation has not yet completed.
        """
        if self._thrust_time is None:
            raise ValueError("Thrust time has not been set, run the simulation.")
        return self._thrust_time

    @property
    def burn_time(self) -> float | None:
        """
        Return the burn time [s], or None if burnout was never reached.

        Only propellant depletion defines a burn time: a run whose thrust
        terminates with propellant remaining leaves it undefined.
        """
        return self._burn_time
burn_time property

Return the burn time [s], or None if burnout was never reached.

Only propellant depletion defines a burn time: a run whose thrust terminates with propellant remaining leaves it undefined.

initial_propellant_mass property

Return the initial propellant mass [kg].

thrust_time property

Return the thrust time [s].

Raises:

Type Description
ValueError

If the simulation has not yet completed.
__init__(motor, igniter_pressure, external_pressure)

Initialize a motor state.

Parameters:

Name Type Description Default
motor Motor

Motor to track.

 required
igniter_pressure float

Initial chamber pressure from the igniter [Pa].

 required
external_pressure float

Ambient pressure [Pa].

 required
Source code in machwave/simulation/states.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def __init__(
    self,
    motor: motors.Motor,
    igniter_pressure: float,
    external_pressure: float,
) -> None:
    """
    Initialize a motor state.

    Args:
        motor: Motor to track.
        igniter_pressure: Initial chamber pressure from the igniter [Pa].
        external_pressure: Ambient pressure [Pa].
    """
    self.motor = motor
    self.external_pressure = external_pressure

    self.time: SimulationStateArray = [0.0]
    self.chamber_pressure: SimulationStateArray = [igniter_pressure]

    self.propellant_mass: SimulationStateArray = []
    self.exit_pressure: SimulationStateArray = []
    self.ideal_thrust_coefficient: SimulationStateArray = []
    self.thrust_coefficient: SimulationStateArray = []
    self.thrust: SimulationStateArray = []
    self.peak_thrust: float = 0.0
    self.nozzle_efficiency: SimulationStateArray = []
    self.loss_fractions: dict[str, SimulationStateArray] = {
        name: [] for name in motor.nozzle_loss_model.component_names
    }

    self._thrust_time: float | None = None
    self._burn_time: float | None = None

    self.end_thrust: bool = False
    self.end_burn: bool = False
build_result()

Return a frozen SimulationResult snapshot of this state.

Source code in machwave/simulation/states.py 
191
192
193
def build_result(self) -> "simulation_results.SimulationResult":
    """Return a frozen ``SimulationResult`` snapshot of this state."""
    return self.result_class.from_state(self)
run_timestep(*args, **kwargs) abstractmethod

Advance the per-step accumulators by one time increment.

Source code in machwave/simulation/states.py 
63
64
65
@abstractmethod
def run_timestep(self, *args, **kwargs) -> None:
    """Advance the per-step accumulators by one time increment."""

TimestepConditions dataclass

Engine conditions at one simulation timestep, in SI units.

Holds the scalar operating quantities every engine/motor type computes for the step, except the performance-related ones (thrust, thrust coefficient, nozzle efficiency).

Source code in machwave/simulation/states.py 
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
@dataclasses.dataclass(frozen=True, kw_only=True)
class TimestepConditions:
    """
    Engine conditions at one simulation timestep, in SI units.

    Holds the scalar operating quantities every engine/motor type computes for the step,
    except the performance-related ones (thrust, thrust coefficient, nozzle efficiency).
    """

    time: float
    chamber_pressure: float
    external_pressure: float
    exit_pressure: float
    effective_expansion_ratio: float
    free_chamber_volume: float
    propellant_mass: float
    propellant_mass_flow_rate: float
    nozzle: thrust_chamber_models.Nozzle
    propellant_properties: propellant_properties_models.ThermochemicalProperties