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.

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 
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
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 
40
41
42
43
44
45
46
47
48
49
50
51
52
53
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 
67
68
69
70
71
72
73
74
75
76
77
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.
igniter_pressure float

Igniter pressure.
external_pressure float

External pressure.
Source code in machwave/simulation/base.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
@dataclass
class InternalBallisticsSimulationParams:
    """
    Parameters for an internal ballistics simulation.

    Attributes:
        d_t: Time step.
        igniter_pressure: Igniter pressure.
        external_pressure: External pressure.
    """

    d_t: float
    igniter_pressure: float
    external_pressure: float

MotorState

Bases: ABC

Defines the states and iteration step for a motor operation.

Source code in machwave/simulation/states.py 
 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
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.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)

    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:
        """
        Return the burn time [s].

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

burn_time property

Return the burn time [s].

Raises:

Type Description
ValueError

If the simulation has not yet completed.

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 
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
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.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 
147
148
149
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 
59
60
61
@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
@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]
    burn_time: float
    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 summary(self) -> dict[str, float]:
        """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 
41
42
43
44
45
46
47
@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 
86
87
88
89
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 
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def summary(self) -> dict[str, float]:
    """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 
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
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 
40
41
42
43
44
45
46
47
48
49
50
51
52
53
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 
67
68
69
70
71
72
73
74
75
76
77
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.
igniter_pressure float

Igniter pressure.
external_pressure float

External pressure.
Source code in machwave/simulation/base.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
@dataclass
class InternalBallisticsSimulationParams:
    """
    Parameters for an internal ballistics simulation.

    Attributes:
        d_t: Time step.
        igniter_pressure: Igniter pressure.
        external_pressure: External pressure.
    """

    d_t: float
    igniter_pressure: float
    external_pressure: float

biliquid

BiliquidEngineState

Bases: MotorState

State for a biliquid rocket engine.

Source code in machwave/simulation/biliquid/states.py 
 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
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,
        )

        self.oxidizer_mass: simulation_states.SimulationStateArray = [
            motor.feed_system.oxidizer_tank.initial_fluid_mass
        ]
        self.fuel_mass: simulation_states.SimulationStateArray = [
            motor.feed_system.fuel_tank.initial_fluid_mass
        ]

        self.fuel_mass_flow_rate: simulation_states.SimulationStateArray = []
        self.oxidizer_mass_flow_rate: simulation_states.SimulationStateArray = []
        self.oxidizer_to_fuel_ratio: simulation_states.SimulationStateArray = []
        self.fuel_tank_pressure: simulation_states.SimulationStateArray = []
        self.oxidizer_tank_pressure: 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 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]
        fuel_mass = self.fuel_mass[-1]
        oxidizer_mass = self.oxidizer_mass[-1]

        propellant_mass = fuel_mass + oxidizer_mass
        self.propellant_mass.append(propellant_mass)

        fuel_tank_pressure = feed_system.get_fuel_tank_pressure(
            oxidizer_mass=oxidizer_mass, fuel_mass=fuel_mass
        )
        self.fuel_tank_pressure.append(fuel_tank_pressure)
        oxidizer_tank_pressure = feed_system.get_oxidizer_tank_pressure(
            oxidizer_mass=oxidizer_mass
        )
        self.oxidizer_tank_pressure.append(oxidizer_tank_pressure)

        is_feeding = (
            propellant_mass > 0
            and fuel_tank_pressure > chamber_pressure
            and oxidizer_tank_pressure > chamber_pressure
        )
        injector_flows = functools.partial(
            get_injector_mass_flows,
            feed_system=feed_system,
            injector=self.motor.thrust_chamber.injector,
            fuel_mass=fuel_mass,
            oxidizer_mass=oxidizer_mass,
            fuel_tank_pressure=fuel_tank_pressure,
            oxidizer_tank_pressure=oxidizer_tank_pressure,
            is_feeding=is_feeding,
            d_t=d_t,
        )
        m_dot_fuel, m_dot_ox = injector_flows(chamber_pressure)
        self.fuel_mass_flow_rate.append(m_dot_fuel)
        self.oxidizer_mass_flow_rate.append(m_dot_ox)
        fuel_consumed = m_dot_fuel * d_t
        oxidizer_consumed = m_dot_ox * d_t

        if m_dot_fuel > 0.0:
            oxidizer_to_fuel_ratio = m_dot_ox / m_dot_fuel
        else:
            design_ratio = self.motor.propellant.oxidizer_to_fuel_ratio
            assert design_ratio is not None
            oxidizer_to_fuel_ratio = design_ratio
        self.oxidizer_to_fuel_ratio.append(oxidizer_to_fuel_ratio)

        propellant_properties = self._evaluate_propellant_properties(
            chamber_pressure=chamber_pressure,
            mixture_ratio=oxidizer_to_fuel_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=m_dot_fuel + m_dot_ox,
            nozzle=nozzle,
            propellant_properties=propellant_properties,
            fuel_mass=fuel_mass,
            oxidizer_mass=oxidizer_mass,
            fuel_mass_flow_rate=m_dot_fuel,
            oxidizer_mass_flow_rate=m_dot_ox,
            oxidizer_to_fuel_ratio=oxidizer_to_fuel_ratio,
            fuel_tank_pressure=fuel_tank_pressure,
            oxidizer_tank_pressure=oxidizer_tank_pressure,
        )
        self._apply_nozzle_losses(
            ideal_momentum_term,
            ideal_pressure_term,
            timestep_conditions,
            chamber_pressure,
        )

        if (
            not is_feeding
            or fuel_consumed >= fuel_mass
            or oxidizer_consumed >= oxidizer_mass
        ) and not self.end_burn:
            self.end_burn = True
            self._burn_time = time + d_t

        if not isentropic.is_flow_choked(
            chamber_pressure,
            external_pressure,
            isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
        ):
            if self._burn_time is None:
                self._burn_time = time
            self._thrust_time = time
            self.end_thrust = True
            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)
        self.fuel_mass.append(fuel_mass - fuel_consumed)
        self.oxidizer_mass.append(oxidizer_mass - oxidizer_consumed)
__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 
 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
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,
    )

    self.oxidizer_mass: simulation_states.SimulationStateArray = [
        motor.feed_system.oxidizer_tank.initial_fluid_mass
    ]
    self.fuel_mass: simulation_states.SimulationStateArray = [
        motor.feed_system.fuel_tank.initial_fluid_mass
    ]

    self.fuel_mass_flow_rate: simulation_states.SimulationStateArray = []
    self.oxidizer_mass_flow_rate: simulation_states.SimulationStateArray = []
    self.oxidizer_to_fuel_ratio: simulation_states.SimulationStateArray = []
    self.fuel_tank_pressure: simulation_states.SimulationStateArray = []
    self.oxidizer_tank_pressure: 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 
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
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]
    fuel_mass = self.fuel_mass[-1]
    oxidizer_mass = self.oxidizer_mass[-1]

    propellant_mass = fuel_mass + oxidizer_mass
    self.propellant_mass.append(propellant_mass)

    fuel_tank_pressure = feed_system.get_fuel_tank_pressure(
        oxidizer_mass=oxidizer_mass, fuel_mass=fuel_mass
    )
    self.fuel_tank_pressure.append(fuel_tank_pressure)
    oxidizer_tank_pressure = feed_system.get_oxidizer_tank_pressure(
        oxidizer_mass=oxidizer_mass
    )
    self.oxidizer_tank_pressure.append(oxidizer_tank_pressure)

    is_feeding = (
        propellant_mass > 0
        and fuel_tank_pressure > chamber_pressure
        and oxidizer_tank_pressure > chamber_pressure
    )
    injector_flows = functools.partial(
        get_injector_mass_flows,
        feed_system=feed_system,
        injector=self.motor.thrust_chamber.injector,
        fuel_mass=fuel_mass,
        oxidizer_mass=oxidizer_mass,
        fuel_tank_pressure=fuel_tank_pressure,
        oxidizer_tank_pressure=oxidizer_tank_pressure,
        is_feeding=is_feeding,
        d_t=d_t,
    )
    m_dot_fuel, m_dot_ox = injector_flows(chamber_pressure)
    self.fuel_mass_flow_rate.append(m_dot_fuel)
    self.oxidizer_mass_flow_rate.append(m_dot_ox)
    fuel_consumed = m_dot_fuel * d_t
    oxidizer_consumed = m_dot_ox * d_t

    if m_dot_fuel > 0.0:
        oxidizer_to_fuel_ratio = m_dot_ox / m_dot_fuel
    else:
        design_ratio = self.motor.propellant.oxidizer_to_fuel_ratio
        assert design_ratio is not None
        oxidizer_to_fuel_ratio = design_ratio
    self.oxidizer_to_fuel_ratio.append(oxidizer_to_fuel_ratio)

    propellant_properties = self._evaluate_propellant_properties(
        chamber_pressure=chamber_pressure,
        mixture_ratio=oxidizer_to_fuel_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=m_dot_fuel + m_dot_ox,
        nozzle=nozzle,
        propellant_properties=propellant_properties,
        fuel_mass=fuel_mass,
        oxidizer_mass=oxidizer_mass,
        fuel_mass_flow_rate=m_dot_fuel,
        oxidizer_mass_flow_rate=m_dot_ox,
        oxidizer_to_fuel_ratio=oxidizer_to_fuel_ratio,
        fuel_tank_pressure=fuel_tank_pressure,
        oxidizer_tank_pressure=oxidizer_tank_pressure,
    )
    self._apply_nozzle_losses(
        ideal_momentum_term,
        ideal_pressure_term,
        timestep_conditions,
        chamber_pressure,
    )

    if (
        not is_feeding
        or fuel_consumed >= fuel_mass
        or oxidizer_consumed >= oxidizer_mass
    ) and not self.end_burn:
        self.end_burn = True
        self._burn_time = time + d_t

    if not isentropic.is_flow_choked(
        chamber_pressure,
        external_pressure,
        isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
    ):
        if self._burn_time is None:
            self._burn_time = time
        self._thrust_time = time
        self.end_thrust = True
        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)
    self.fuel_mass.append(fuel_mass - fuel_consumed)
    self.oxidizer_mass.append(oxidizer_mass - oxidizer_consumed)

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
@dataclass(frozen=True, kw_only=True)
class BiliquidSimulationResult(
    simulation_results.SimulationResult["biliquid_states.BiliquidEngineState"]
):
    """Simulation result for a biliquid engine run."""

    oxidizer_mass: simulation_results.SimulationResultArray
    fuel_mass: simulation_results.SimulationResultArray
    fuel_tank_pressure: simulation_results.SimulationResultArray
    oxidizer_tank_pressure: simulation_results.SimulationResultArray
    final_oxidizer_mass: float
    final_fuel_mass: float

    @classmethod
    def _collect_extra_fields(
        cls, state: "biliquid_states.BiliquidEngineState"
    ) -> dict[str, Any]:
        oxidizer_mass = np.asarray(state.oxidizer_mass)
        fuel_mass = np.asarray(state.fuel_mass)
        return {
            "oxidizer_mass": oxidizer_mass,
            "fuel_mass": fuel_mass,
            "fuel_tank_pressure": np.asarray(state.fuel_tank_pressure),
            "oxidizer_tank_pressure": np.asarray(state.oxidizer_tank_pressure),
            "final_oxidizer_mass": float(oxidizer_mass[-1]),
            "final_fuel_mass": float(fuel_mass[-1]),
        }

    def _extra_summary(self) -> dict[str, float]:
        return {
            "final_oxidizer_mass": self.final_oxidizer_mass,
            "final_fuel_mass": self.final_fuel_mass,
        }

    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.burn_time:.4f} s", 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)
        print(f"  Oxidizer: {self.final_oxidizer_mass:.4f}", file=file)
        print(f"  Fuel:     {self.final_fuel_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
@dataclass(frozen=True, kw_only=True)
class BiliquidSimulationResult(
    simulation_results.SimulationResult["biliquid_states.BiliquidEngineState"]
):
    """Simulation result for a biliquid engine run."""

    oxidizer_mass: simulation_results.SimulationResultArray
    fuel_mass: simulation_results.SimulationResultArray
    fuel_tank_pressure: simulation_results.SimulationResultArray
    oxidizer_tank_pressure: simulation_results.SimulationResultArray
    final_oxidizer_mass: float
    final_fuel_mass: float

    @classmethod
    def _collect_extra_fields(
        cls, state: "biliquid_states.BiliquidEngineState"
    ) -> dict[str, Any]:
        oxidizer_mass = np.asarray(state.oxidizer_mass)
        fuel_mass = np.asarray(state.fuel_mass)
        return {
            "oxidizer_mass": oxidizer_mass,
            "fuel_mass": fuel_mass,
            "fuel_tank_pressure": np.asarray(state.fuel_tank_pressure),
            "oxidizer_tank_pressure": np.asarray(state.oxidizer_tank_pressure),
            "final_oxidizer_mass": float(oxidizer_mass[-1]),
            "final_fuel_mass": float(fuel_mass[-1]),
        }

    def _extra_summary(self) -> dict[str, float]:
        return {
            "final_oxidizer_mass": self.final_oxidizer_mass,
            "final_fuel_mass": self.final_fuel_mass,
        }

    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.burn_time:.4f} s", 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)
        print(f"  Oxidizer: {self.final_oxidizer_mass:.4f}", file=file)
        print(f"  Fuel:     {self.final_fuel_mass:.4f}", file=file)

states

BiliquidEngineState

Bases: MotorState

State for a biliquid rocket engine.

Source code in machwave/simulation/biliquid/states.py 
 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
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,
        )

        self.oxidizer_mass: simulation_states.SimulationStateArray = [
            motor.feed_system.oxidizer_tank.initial_fluid_mass
        ]
        self.fuel_mass: simulation_states.SimulationStateArray = [
            motor.feed_system.fuel_tank.initial_fluid_mass
        ]

        self.fuel_mass_flow_rate: simulation_states.SimulationStateArray = []
        self.oxidizer_mass_flow_rate: simulation_states.SimulationStateArray = []
        self.oxidizer_to_fuel_ratio: simulation_states.SimulationStateArray = []
        self.fuel_tank_pressure: simulation_states.SimulationStateArray = []
        self.oxidizer_tank_pressure: 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 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]
        fuel_mass = self.fuel_mass[-1]
        oxidizer_mass = self.oxidizer_mass[-1]

        propellant_mass = fuel_mass + oxidizer_mass
        self.propellant_mass.append(propellant_mass)

        fuel_tank_pressure = feed_system.get_fuel_tank_pressure(
            oxidizer_mass=oxidizer_mass, fuel_mass=fuel_mass
        )
        self.fuel_tank_pressure.append(fuel_tank_pressure)
        oxidizer_tank_pressure = feed_system.get_oxidizer_tank_pressure(
            oxidizer_mass=oxidizer_mass
        )
        self.oxidizer_tank_pressure.append(oxidizer_tank_pressure)

        is_feeding = (
            propellant_mass > 0
            and fuel_tank_pressure > chamber_pressure
            and oxidizer_tank_pressure > chamber_pressure
        )
        injector_flows = functools.partial(
            get_injector_mass_flows,
            feed_system=feed_system,
            injector=self.motor.thrust_chamber.injector,
            fuel_mass=fuel_mass,
            oxidizer_mass=oxidizer_mass,
            fuel_tank_pressure=fuel_tank_pressure,
            oxidizer_tank_pressure=oxidizer_tank_pressure,
            is_feeding=is_feeding,
            d_t=d_t,
        )
        m_dot_fuel, m_dot_ox = injector_flows(chamber_pressure)
        self.fuel_mass_flow_rate.append(m_dot_fuel)
        self.oxidizer_mass_flow_rate.append(m_dot_ox)
        fuel_consumed = m_dot_fuel * d_t
        oxidizer_consumed = m_dot_ox * d_t

        if m_dot_fuel > 0.0:
            oxidizer_to_fuel_ratio = m_dot_ox / m_dot_fuel
        else:
            design_ratio = self.motor.propellant.oxidizer_to_fuel_ratio
            assert design_ratio is not None
            oxidizer_to_fuel_ratio = design_ratio
        self.oxidizer_to_fuel_ratio.append(oxidizer_to_fuel_ratio)

        propellant_properties = self._evaluate_propellant_properties(
            chamber_pressure=chamber_pressure,
            mixture_ratio=oxidizer_to_fuel_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=m_dot_fuel + m_dot_ox,
            nozzle=nozzle,
            propellant_properties=propellant_properties,
            fuel_mass=fuel_mass,
            oxidizer_mass=oxidizer_mass,
            fuel_mass_flow_rate=m_dot_fuel,
            oxidizer_mass_flow_rate=m_dot_ox,
            oxidizer_to_fuel_ratio=oxidizer_to_fuel_ratio,
            fuel_tank_pressure=fuel_tank_pressure,
            oxidizer_tank_pressure=oxidizer_tank_pressure,
        )
        self._apply_nozzle_losses(
            ideal_momentum_term,
            ideal_pressure_term,
            timestep_conditions,
            chamber_pressure,
        )

        if (
            not is_feeding
            or fuel_consumed >= fuel_mass
            or oxidizer_consumed >= oxidizer_mass
        ) and not self.end_burn:
            self.end_burn = True
            self._burn_time = time + d_t

        if not isentropic.is_flow_choked(
            chamber_pressure,
            external_pressure,
            isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
        ):
            if self._burn_time is None:
                self._burn_time = time
            self._thrust_time = time
            self.end_thrust = True
            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)
        self.fuel_mass.append(fuel_mass - fuel_consumed)
        self.oxidizer_mass.append(oxidizer_mass - oxidizer_consumed)
__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 
 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
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,
    )

    self.oxidizer_mass: simulation_states.SimulationStateArray = [
        motor.feed_system.oxidizer_tank.initial_fluid_mass
    ]
    self.fuel_mass: simulation_states.SimulationStateArray = [
        motor.feed_system.fuel_tank.initial_fluid_mass
    ]

    self.fuel_mass_flow_rate: simulation_states.SimulationStateArray = []
    self.oxidizer_mass_flow_rate: simulation_states.SimulationStateArray = []
    self.oxidizer_to_fuel_ratio: simulation_states.SimulationStateArray = []
    self.fuel_tank_pressure: simulation_states.SimulationStateArray = []
    self.oxidizer_tank_pressure: 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 
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
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]
    fuel_mass = self.fuel_mass[-1]
    oxidizer_mass = self.oxidizer_mass[-1]

    propellant_mass = fuel_mass + oxidizer_mass
    self.propellant_mass.append(propellant_mass)

    fuel_tank_pressure = feed_system.get_fuel_tank_pressure(
        oxidizer_mass=oxidizer_mass, fuel_mass=fuel_mass
    )
    self.fuel_tank_pressure.append(fuel_tank_pressure)
    oxidizer_tank_pressure = feed_system.get_oxidizer_tank_pressure(
        oxidizer_mass=oxidizer_mass
    )
    self.oxidizer_tank_pressure.append(oxidizer_tank_pressure)

    is_feeding = (
        propellant_mass > 0
        and fuel_tank_pressure > chamber_pressure
        and oxidizer_tank_pressure > chamber_pressure
    )
    injector_flows = functools.partial(
        get_injector_mass_flows,
        feed_system=feed_system,
        injector=self.motor.thrust_chamber.injector,
        fuel_mass=fuel_mass,
        oxidizer_mass=oxidizer_mass,
        fuel_tank_pressure=fuel_tank_pressure,
        oxidizer_tank_pressure=oxidizer_tank_pressure,
        is_feeding=is_feeding,
        d_t=d_t,
    )
    m_dot_fuel, m_dot_ox = injector_flows(chamber_pressure)
    self.fuel_mass_flow_rate.append(m_dot_fuel)
    self.oxidizer_mass_flow_rate.append(m_dot_ox)
    fuel_consumed = m_dot_fuel * d_t
    oxidizer_consumed = m_dot_ox * d_t

    if m_dot_fuel > 0.0:
        oxidizer_to_fuel_ratio = m_dot_ox / m_dot_fuel
    else:
        design_ratio = self.motor.propellant.oxidizer_to_fuel_ratio
        assert design_ratio is not None
        oxidizer_to_fuel_ratio = design_ratio
    self.oxidizer_to_fuel_ratio.append(oxidizer_to_fuel_ratio)

    propellant_properties = self._evaluate_propellant_properties(
        chamber_pressure=chamber_pressure,
        mixture_ratio=oxidizer_to_fuel_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=m_dot_fuel + m_dot_ox,
        nozzle=nozzle,
        propellant_properties=propellant_properties,
        fuel_mass=fuel_mass,
        oxidizer_mass=oxidizer_mass,
        fuel_mass_flow_rate=m_dot_fuel,
        oxidizer_mass_flow_rate=m_dot_ox,
        oxidizer_to_fuel_ratio=oxidizer_to_fuel_ratio,
        fuel_tank_pressure=fuel_tank_pressure,
        oxidizer_tank_pressure=oxidizer_tank_pressure,
    )
    self._apply_nozzle_losses(
        ideal_momentum_term,
        ideal_pressure_term,
        timestep_conditions,
        chamber_pressure,
    )

    if (
        not is_feeding
        or fuel_consumed >= fuel_mass
        or oxidizer_consumed >= oxidizer_mass
    ) and not self.end_burn:
        self.end_burn = True
        self._burn_time = time + d_t

    if not isentropic.is_flow_choked(
        chamber_pressure,
        external_pressure,
        isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
    ):
        if self._burn_time is None:
            self._burn_time = time
        self._thrust_time = time
        self.end_thrust = True
        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)
    self.fuel_mass.append(fuel_mass - fuel_consumed)
    self.oxidizer_mass.append(oxidizer_mass - oxidizer_consumed)
BiliquidTimestepConditions dataclass

Bases: TimestepConditions

Timestep conditions for a biliquid engine.

Source code in machwave/simulation/biliquid/states.py 
68
69
70
71
72
73
74
75
76
77
78
@dataclasses.dataclass(frozen=True, kw_only=True)
class BiliquidTimestepConditions(simulation_states.TimestepConditions):
    """Timestep conditions for a biliquid engine."""

    fuel_mass: float
    oxidizer_mass: float
    fuel_mass_flow_rate: float
    oxidizer_mass_flow_rate: float
    oxidizer_to_fuel_ratio: float
    fuel_tank_pressure: float
    oxidizer_tank_pressure: float
get_injector_mass_flows(chamber_pressure, *, feed_system, injector, fuel_mass, oxidizer_mass, fuel_tank_pressure, oxidizer_tank_pressure, is_feeding, d_t)

Fuel and oxidizer injector flows at the given chamber pressure [kg/s].

Source code in machwave/simulation/biliquid/states.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
def get_injector_mass_flows(
    chamber_pressure: float,
    *,
    feed_system: feed_systems.FeedSystem,
    injector: injector_models.BipropellantInjector,
    fuel_mass: float,
    oxidizer_mass: float,
    fuel_tank_pressure: float,
    oxidizer_tank_pressure: float,
    is_feeding: bool,
    d_t: float,
) -> tuple[float, float]:
    """Fuel and oxidizer injector flows at the given chamber pressure [kg/s]."""
    if not is_feeding:
        return 0.0, 0.0
    fuel_flow = (
        feed_system.get_mass_flow_fuel(
            chamber_pressure=chamber_pressure,
            injector=injector,
            fuel_mass=fuel_mass,
            oxidizer_mass=oxidizer_mass,
        )
        if fuel_tank_pressure > chamber_pressure
        else 0.0
    )
    oxidizer_flow = (
        feed_system.get_mass_flow_ox(
            chamber_pressure=chamber_pressure,
            injector=injector,
            oxidizer_mass=oxidizer_mass,
        )
        if oxidizer_tank_pressure > chamber_pressure
        else 0.0
    )
    return (
        min(fuel_flow, fuel_mass / d_t),
        min(oxidizer_flow, oxidizer_mass / d_t),
    )
get_total_injector_mass_flow(chamber_pressure, *, injector_flows)

Total injector mass flow (fuel + oxidizer) at the given pressure [kg/s].

Source code in machwave/simulation/biliquid/states.py 
59
60
61
62
63
64
65
def get_total_injector_mass_flow(
    chamber_pressure: float,
    *,
    injector_flows: Callable[[float], tuple[float, float]],
) -> float:
    """Total injector mass flow (fuel + oxidizer) at the given pressure [kg/s]."""
    return sum(injector_flows(chamber_pressure))

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
@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]
    burn_time: float
    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 summary(self) -> dict[str, float]:
        """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 
41
42
43
44
45
46
47
@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 
86
87
88
89
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 
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def summary(self) -> dict[str, float]:
    """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
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].

        Raises:
            ValueError: If the motor's propellant has no thermochemical
                properties. Solid propellant properties are fixed for the
                entire run, so a missing value is a configuration error and
                the simulation refuses to start.
        """
        propellant_properties = motor.propellant.properties
        if propellant_properties is None:
            raise ValueError(
                "Propellant properties must be defined to run the simulation."
            )

        super().__init__(
            motor=motor,
            igniter_pressure=igniter_pressure,
            external_pressure=external_pressure,
        )

        self.propellant_properties = propellant_properties
        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 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].
        """
        propellant_properties = self.propellant_properties
        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]

        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
            )
            propellant_moi = self.motor.grain.get_moment_of_inertia(
                ideal_density=ideal_propellant_density, web_distance=web_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 not isentropic.is_flow_choked(
            chamber_pressure,
            external_pressure,
            isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
        ):
            if self._burn_time is None:
                self._burn_time = time
            self._thrust_time = time
            self.end_thrust = True
            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

Raises:

Type Description
ValueError

If the motor's propellant has no thermochemical properties. Solid propellant properties are fixed for the entire run, so a missing value is a configuration error and the simulation refuses to start.
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
108
109
110
111
112
113
114
115
116
117
118
119
120
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].

    Raises:
        ValueError: If the motor's propellant has no thermochemical
            properties. Solid propellant properties are fixed for the
            entire run, so a missing value is a configuration error and
            the simulation refuses to start.
    """
    propellant_properties = motor.propellant.properties
    if propellant_properties is None:
        raise ValueError(
            "Propellant properties must be defined to run the simulation."
        )

    super().__init__(
        motor=motor,
        igniter_pressure=igniter_pressure,
        external_pressure=external_pressure,
    )

    self.propellant_properties = propellant_properties
    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 
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
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].
    """
    propellant_properties = self.propellant_properties
    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]

    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
        )
        propellant_moi = self.motor.grain.get_moment_of_inertia(
            ideal_density=ideal_propellant_density, web_distance=web_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 not isentropic.is_flow_choked(
        chamber_pressure,
        external_pressure,
        isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
    ):
        if self._burn_time is None:
            self._burn_time = time
        self._thrust_time = time
        self.end_thrust = True
        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
@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: 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, thrust time: {self.burn_time:.3f}, {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
@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: 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, thrust time: {self.burn_time:.3f}, {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
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].

        Raises:
            ValueError: If the motor's propellant has no thermochemical
                properties. Solid propellant properties are fixed for the
                entire run, so a missing value is a configuration error and
                the simulation refuses to start.
        """
        propellant_properties = motor.propellant.properties
        if propellant_properties is None:
            raise ValueError(
                "Propellant properties must be defined to run the simulation."
            )

        super().__init__(
            motor=motor,
            igniter_pressure=igniter_pressure,
            external_pressure=external_pressure,
        )

        self.propellant_properties = propellant_properties
        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 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].
        """
        propellant_properties = self.propellant_properties
        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]

        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
            )
            propellant_moi = self.motor.grain.get_moment_of_inertia(
                ideal_density=ideal_propellant_density, web_distance=web_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 not isentropic.is_flow_choked(
            chamber_pressure,
            external_pressure,
            isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
        ):
            if self._burn_time is None:
                self._burn_time = time
            self._thrust_time = time
            self.end_thrust = True
            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

Raises:

Type Description
ValueError

If the motor's propellant has no thermochemical properties. Solid propellant properties are fixed for the entire run, so a missing value is a configuration error and the simulation refuses to start.
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
108
109
110
111
112
113
114
115
116
117
118
119
120
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].

    Raises:
        ValueError: If the motor's propellant has no thermochemical
            properties. Solid propellant properties are fixed for the
            entire run, so a missing value is a configuration error and
            the simulation refuses to start.
    """
    propellant_properties = motor.propellant.properties
    if propellant_properties is None:
        raise ValueError(
            "Propellant properties must be defined to run the simulation."
        )

    super().__init__(
        motor=motor,
        igniter_pressure=igniter_pressure,
        external_pressure=external_pressure,
    )

    self.propellant_properties = propellant_properties
    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 
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
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].
    """
    propellant_properties = self.propellant_properties
    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]

    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
        )
        propellant_moi = self.motor.grain.get_moment_of_inertia(
            ideal_density=ideal_propellant_density, web_distance=web_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 not isentropic.is_flow_choked(
        chamber_pressure,
        external_pressure,
        isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
    ):
        if self._burn_time is None:
            self._burn_time = time
        self._thrust_time = time
        self.end_thrust = True
        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 
 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
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.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)

    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:
        """
        Return the burn time [s].

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

Return the burn time [s].

Raises:

Type Description
ValueError

If the simulation has not yet completed.
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 
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
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.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 
147
148
149
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 
59
60
61
@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 
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
@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