Skip to content

States

Motor operation state tracking. Each motor type has a corresponding state class (SolidMotorState, LiquidEngineState) that accumulates time-series data during simulation: chamber pressure, exit pressure, thrust, Cf, propellant mass, and efficiency breakdown (divergent, kinetics, boundary layer, two-phase losses).

Solid motor states also track grain regression (web distance, burn area, volume), propellant CoG, and moment-of-inertia tensors over time. Liquid engine states additionally store tank pressures and propellant masses per tank.

machwave.states

LiquidEngineState

Bases: MotorState

State for a Liquid Rocket Engine.

Source code in machwave/states/liquid_engine.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
 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
class LiquidEngineState(MotorState):
    """
    State for a Liquid Rocket Engine.
    """

    motor: LiquidEngine

    SIMULATION_ARRAY_ATTRIBUTE_NAMES = MotorState.SIMULATION_ARRAY_ATTRIBUTE_NAMES + (
        "oxidizer_mass",
        "fuel_mass",
        "nozzle_correction_factor",
        "fuel_tank_pressure",
        "oxidizer_tank_pressure",
    )

    def __init__(
        self,
        motor: LiquidEngine,
        initial_pressure: float,
        initial_atmospheric_pressure: float,
        other_losses: float,
    ) -> None:
        """
        Initial parameters for a LRE operation.
        """
        super().__init__(
            motor=motor,
            initial_pressure=initial_pressure,
            initial_atmospheric_pressure=initial_atmospheric_pressure,
            other_losses=other_losses,
        )

        self.oxidizer_mass: SimulationArray = [
            motor.feed_system.oxidizer_tank.fluid_mass
        ]
        self.fuel_mass: SimulationArray = [motor.feed_system.fuel_tank.fluid_mass]
        self.nozzle_correction_factor: SimulationArray = [1.0]

        self.fuel_tank_pressure: SimulationArray = [
            motor.feed_system.get_fuel_tank_pressure()
        ]
        self.oxidizer_tank_pressure: SimulationArray = [
            motor.feed_system.get_oxidizer_tank_pressure()
        ]

    def get_m_dot_in(self) -> float:
        return self._m_dot_fuel + self._m_dot_ox

    def run_timestep(
        self,
        d_t: float,
        external_pressure: float,
    ) -> None:
        """
        Advance simulation by time step d_t under external pressure.

        Args:
            d_t: Time step.
            external_pressure: External pressure.
        """
        if self.end_thrust:
            return

        nz = self.motor.thrust_chamber.nozzle
        injector = self.motor.thrust_chamber.injector
        feed = self.motor.feed_system

        self.t.append(self.t[-1] + d_t)

        if self.propellant_mass[-1] > 0:
            m_dot_fuel = feed.get_mass_flow_fuel(
                chamber_pressure=self.chamber_pressure[-1],
                discharge_coefficient=injector.discharge_coefficient_fuel,
                injector_area=injector.area_fuel,
            )
            m_dot_ox = feed.get_mass_flow_ox(
                chamber_pressure=self.chamber_pressure[-1],
                discharge_coefficient=injector.discharge_coefficient_oxidizer,
                injector_area=injector.area_ox,
            )
        else:
            m_dot_fuel = m_dot_ox = 0.0

        m_dot_fuel = min(m_dot_fuel, self.fuel_mass[-1] / d_t)
        m_dot_ox = min(m_dot_ox, self.oxidizer_mass[-1] / d_t)

        of = self.motor.propellant.of_ratio
        assert of is not None
        cons_f = m_dot_fuel * d_t
        cons_o = m_dot_ox * d_t
        if cons_f >= self.fuel_mass[-1]:
            cons_f = self.fuel_mass[-1]
            cons_o = of * cons_f
            self.end_burn = True
            self.burn_time = self.t[-1]
        elif cons_o >= self.oxidizer_mass[-1]:
            cons_o = self.oxidizer_mass[-1]
            cons_f = cons_o / of
            self.end_burn = True
            self.burn_time = self.t[-1]
        m_dot_fuel = cons_f / d_t
        m_dot_ox = cons_o / d_t
        self._m_dot_fuel = m_dot_fuel
        self._m_dot_ox = m_dot_ox

        self.motor.propellant.properties = self.motor.propellant.evaluate(
            chamber_pressure=self.chamber_pressure[-1],
            expansion_ratio=nz.expansion_ratio,
        )
        props = self.motor.propellant.properties
        assert props is not None

        new_chamber_pressure = rk4th_ode_solver(
            variables={"chamber_pressure": self.chamber_pressure[-1]},
            equation=compute_chamber_pressure_mass_balance,
            d_t=d_t,
            external_pressure=external_pressure,
            mass_flow_in=self.get_m_dot_in(),
            free_chamber_volume=self.motor.thrust_chamber.combustion_chamber.internal_volume,
            throat_area=nz.get_throat_area(),
            k=props.k_chamber,
            R=props.R_chamber,
            flame_temperature=props.adiabatic_flame_temperature,
            discharge_coefficient=nz.discharge_coefficient,
        )[0]
        self.chamber_pressure.append(new_chamber_pressure)
        exit_pressure = get_exit_pressure(
            props.k_exhaust, nz.expansion_ratio, new_chamber_pressure
        )
        self.exit_pressure.append(exit_pressure)

        chamber_pressure_psi = conversions.convert_pa_to_psi(new_chamber_pressure)
        divergent_loss = get_nozzle_divergent_loss_fraction(
            divergent_angle=nz.divergent_angle,
        )
        kinetics_loss = get_kinetics_loss_fraction(
            i_sp_th_frozen=props.i_sp_frozen,
            i_sp_th_shifting=props.i_sp_shifting,
            chamber_pressure_psi=chamber_pressure_psi,
        )
        nozzle_efficiency = get_overall_nozzle_efficiency(
            divergent_loss, kinetics_loss, 0.0, 0.0, other_losses=self.other_losses
        )
        nozzle_correction_factor = (
            nozzle_efficiency * self.motor.propellant.combustion_efficiency
        )
        self.nozzle_correction_factor.append(nozzle_correction_factor)

        thrust_coefficient_ideal = get_ideal_thrust_coefficient(
            new_chamber_pressure,
            exit_pressure,
            external_pressure,
            nz.expansion_ratio,
            props.k_exhaust,
        )
        thrust_coefficient = apply_thrust_coefficient_correction(
            thrust_coefficient_ideal, nozzle_correction_factor
        )
        self.thrust_coefficient.append(thrust_coefficient)
        self.thrust_coefficient_ideal.append(thrust_coefficient_ideal)
        self.thrust.append(
            get_thrust_from_thrust_coefficient(
                thrust_coefficient, new_chamber_pressure, nz.get_throat_area()
            )
        )

        feed.fuel_tank.remove_propellant(cons_f)
        feed.oxidizer_tank.remove_propellant(cons_o)
        new_fuel = self.fuel_mass[-1] - cons_f
        new_ox = self.oxidizer_mass[-1] - cons_o
        self.fuel_mass.append(new_fuel)
        self.oxidizer_mass.append(new_ox)
        self.propellant_mass.append(new_fuel + new_ox)

        self.fuel_tank_pressure.append(feed.get_fuel_tank_pressure())
        self.oxidizer_tank_pressure.append(feed.get_oxidizer_tank_pressure())

        if self.propellant_mass[-1] <= 0 and not self.end_burn:
            self.end_burn = True
            self.burn_time = self.t[-1]
        if not is_flow_choked(
            new_chamber_pressure,
            external_pressure,
            get_critical_pressure_ratio(props.k_chamber),
        ):
            self._thrust_time = self.t[-1]
            self.end_thrust = True

    def print_results(self) -> None:
        """
        Prints the results obtained during the Liquid Rocket Engine operation.
        """
        # HEADER
        print("\nLIQUID ENGINE OPERATION RESULTS")

        # Propellant summary
        print(f"Initial propellant mass: {self.propellant_mass[0]:.4f} kg")
        try:
            print(f"Burnout time: {self.burn_time:.4f} s")
        except AttributeError:
            print("Burnout time: (not reached)")
        print(f"Thrust time: {self.thrust_time:.4f} s")

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

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

        # Impulse
        impulse = performance.get_total_impulse(
            np.asarray(self.thrust), np.asarray(self.t)
        )
        isp = performance.get_specific_impulse(impulse, self.propellant_mass[0])
        print("\nIMPULSE AND I_SP")
        print(f"  Total impulse: {impulse:.4f} N·s")
        print(f"  Specific impulse: {isp:.4f} s")

        # Remaining propellant masses
        print("\nPROPELLANT REMAINING (kg)")
        print(f"  Oxidizer: {self.oxidizer_mass[-1]:.4f}")
        print(f"  Fuel:     {self.fuel_mass[-1]:.4f}")

__init__(motor, initial_pressure, initial_atmospheric_pressure, other_losses)

Initial parameters for a LRE operation.

Source code in machwave/states/liquid_engine.py 
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
def __init__(
    self,
    motor: LiquidEngine,
    initial_pressure: float,
    initial_atmospheric_pressure: float,
    other_losses: float,
) -> None:
    """
    Initial parameters for a LRE operation.
    """
    super().__init__(
        motor=motor,
        initial_pressure=initial_pressure,
        initial_atmospheric_pressure=initial_atmospheric_pressure,
        other_losses=other_losses,
    )

    self.oxidizer_mass: SimulationArray = [
        motor.feed_system.oxidizer_tank.fluid_mass
    ]
    self.fuel_mass: SimulationArray = [motor.feed_system.fuel_tank.fluid_mass]
    self.nozzle_correction_factor: SimulationArray = [1.0]

    self.fuel_tank_pressure: SimulationArray = [
        motor.feed_system.get_fuel_tank_pressure()
    ]
    self.oxidizer_tank_pressure: SimulationArray = [
        motor.feed_system.get_oxidizer_tank_pressure()
    ]

print_results()

Prints the results obtained during the Liquid Rocket Engine operation.

Source code in machwave/states/liquid_engine.py 
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
def print_results(self) -> None:
    """
    Prints the results obtained during the Liquid Rocket Engine operation.
    """
    # HEADER
    print("\nLIQUID ENGINE OPERATION RESULTS")

    # Propellant summary
    print(f"Initial propellant mass: {self.propellant_mass[0]:.4f} kg")
    try:
        print(f"Burnout time: {self.burn_time:.4f} s")
    except AttributeError:
        print("Burnout time: (not reached)")
    print(f"Thrust time: {self.thrust_time:.4f} s")

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

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

    # Impulse
    impulse = performance.get_total_impulse(
        np.asarray(self.thrust), np.asarray(self.t)
    )
    isp = performance.get_specific_impulse(impulse, self.propellant_mass[0])
    print("\nIMPULSE AND I_SP")
    print(f"  Total impulse: {impulse:.4f} N·s")
    print(f"  Specific impulse: {isp:.4f} s")

    # Remaining propellant masses
    print("\nPROPELLANT REMAINING (kg)")
    print(f"  Oxidizer: {self.oxidizer_mass[-1]:.4f}")
    print(f"  Fuel:     {self.fuel_mass[-1]:.4f}")

run_timestep(d_t, external_pressure)

Advance simulation by time step d_t under external pressure.

Parameters:

Name Type Description Default
d_t float

Time step.

 required
external_pressure float

External pressure.

 required
Source code in machwave/states/liquid_engine.py 
 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
def run_timestep(
    self,
    d_t: float,
    external_pressure: float,
) -> None:
    """
    Advance simulation by time step d_t under external pressure.

    Args:
        d_t: Time step.
        external_pressure: External pressure.
    """
    if self.end_thrust:
        return

    nz = self.motor.thrust_chamber.nozzle
    injector = self.motor.thrust_chamber.injector
    feed = self.motor.feed_system

    self.t.append(self.t[-1] + d_t)

    if self.propellant_mass[-1] > 0:
        m_dot_fuel = feed.get_mass_flow_fuel(
            chamber_pressure=self.chamber_pressure[-1],
            discharge_coefficient=injector.discharge_coefficient_fuel,
            injector_area=injector.area_fuel,
        )
        m_dot_ox = feed.get_mass_flow_ox(
            chamber_pressure=self.chamber_pressure[-1],
            discharge_coefficient=injector.discharge_coefficient_oxidizer,
            injector_area=injector.area_ox,
        )
    else:
        m_dot_fuel = m_dot_ox = 0.0

    m_dot_fuel = min(m_dot_fuel, self.fuel_mass[-1] / d_t)
    m_dot_ox = min(m_dot_ox, self.oxidizer_mass[-1] / d_t)

    of = self.motor.propellant.of_ratio
    assert of is not None
    cons_f = m_dot_fuel * d_t
    cons_o = m_dot_ox * d_t
    if cons_f >= self.fuel_mass[-1]:
        cons_f = self.fuel_mass[-1]
        cons_o = of * cons_f
        self.end_burn = True
        self.burn_time = self.t[-1]
    elif cons_o >= self.oxidizer_mass[-1]:
        cons_o = self.oxidizer_mass[-1]
        cons_f = cons_o / of
        self.end_burn = True
        self.burn_time = self.t[-1]
    m_dot_fuel = cons_f / d_t
    m_dot_ox = cons_o / d_t
    self._m_dot_fuel = m_dot_fuel
    self._m_dot_ox = m_dot_ox

    self.motor.propellant.properties = self.motor.propellant.evaluate(
        chamber_pressure=self.chamber_pressure[-1],
        expansion_ratio=nz.expansion_ratio,
    )
    props = self.motor.propellant.properties
    assert props is not None

    new_chamber_pressure = rk4th_ode_solver(
        variables={"chamber_pressure": self.chamber_pressure[-1]},
        equation=compute_chamber_pressure_mass_balance,
        d_t=d_t,
        external_pressure=external_pressure,
        mass_flow_in=self.get_m_dot_in(),
        free_chamber_volume=self.motor.thrust_chamber.combustion_chamber.internal_volume,
        throat_area=nz.get_throat_area(),
        k=props.k_chamber,
        R=props.R_chamber,
        flame_temperature=props.adiabatic_flame_temperature,
        discharge_coefficient=nz.discharge_coefficient,
    )[0]
    self.chamber_pressure.append(new_chamber_pressure)
    exit_pressure = get_exit_pressure(
        props.k_exhaust, nz.expansion_ratio, new_chamber_pressure
    )
    self.exit_pressure.append(exit_pressure)

    chamber_pressure_psi = conversions.convert_pa_to_psi(new_chamber_pressure)
    divergent_loss = get_nozzle_divergent_loss_fraction(
        divergent_angle=nz.divergent_angle,
    )
    kinetics_loss = get_kinetics_loss_fraction(
        i_sp_th_frozen=props.i_sp_frozen,
        i_sp_th_shifting=props.i_sp_shifting,
        chamber_pressure_psi=chamber_pressure_psi,
    )
    nozzle_efficiency = get_overall_nozzle_efficiency(
        divergent_loss, kinetics_loss, 0.0, 0.0, other_losses=self.other_losses
    )
    nozzle_correction_factor = (
        nozzle_efficiency * self.motor.propellant.combustion_efficiency
    )
    self.nozzle_correction_factor.append(nozzle_correction_factor)

    thrust_coefficient_ideal = get_ideal_thrust_coefficient(
        new_chamber_pressure,
        exit_pressure,
        external_pressure,
        nz.expansion_ratio,
        props.k_exhaust,
    )
    thrust_coefficient = apply_thrust_coefficient_correction(
        thrust_coefficient_ideal, nozzle_correction_factor
    )
    self.thrust_coefficient.append(thrust_coefficient)
    self.thrust_coefficient_ideal.append(thrust_coefficient_ideal)
    self.thrust.append(
        get_thrust_from_thrust_coefficient(
            thrust_coefficient, new_chamber_pressure, nz.get_throat_area()
        )
    )

    feed.fuel_tank.remove_propellant(cons_f)
    feed.oxidizer_tank.remove_propellant(cons_o)
    new_fuel = self.fuel_mass[-1] - cons_f
    new_ox = self.oxidizer_mass[-1] - cons_o
    self.fuel_mass.append(new_fuel)
    self.oxidizer_mass.append(new_ox)
    self.propellant_mass.append(new_fuel + new_ox)

    self.fuel_tank_pressure.append(feed.get_fuel_tank_pressure())
    self.oxidizer_tank_pressure.append(feed.get_oxidizer_tank_pressure())

    if self.propellant_mass[-1] <= 0 and not self.end_burn:
        self.end_burn = True
        self.burn_time = self.t[-1]
    if not is_flow_choked(
        new_chamber_pressure,
        external_pressure,
        get_critical_pressure_ratio(props.k_chamber),
    ):
        self._thrust_time = self.t[-1]
        self.end_thrust = True

MotorState

Bases: ABC

Defines a particular motor operation. Stores and processes all attributes obtained from the simulation.

Source code in machwave/states/base.py 
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 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
class MotorState(ABC):
    """
    Defines a particular motor operation. Stores and processes all attributes
    obtained from the simulation.
    """

    SIMULATION_ARRAY_ATTRIBUTE_NAMES: tuple[str, ...] = (
        "t",
        "propellant_mass",
        "chamber_pressure",
        "exit_pressure",
        "thrust_coefficient",
        "thrust_coefficient_ideal",
        "thrust",
    )

    def __init__(
        self,
        motor: Motor,
        initial_pressure: float,
        initial_atmospheric_pressure: float,
        other_losses: float,
    ) -> None:
        """
        Initializes attributes for the motor operation.
        Each motor category will contain a particular set of attributes.
        """
        self.motor = motor
        self.other_losses = other_losses

        self.t: SimulationArray = [0.0]

        self.propellant_mass: SimulationArray = [motor.initial_propellant_mass]
        self.chamber_pressure: SimulationArray = [initial_pressure]
        self.exit_pressure: SimulationArray = [initial_atmospheric_pressure]

        self.thrust_coefficient: SimulationArray = [0.0]
        self.thrust_coefficient_ideal: SimulationArray = [0.0]
        self.thrust: SimulationArray = [0.0]

        self._thrust_time = None

        self.end_thrust = False
        self.end_burn = False

    def convert_simulation_arrays_to_numpy(self) -> None:
        """Convert accumulated lists into ndarrays for downstream consumers."""
        for name in self.SIMULATION_ARRAY_ATTRIBUTE_NAMES:
            value = getattr(self, name)
            if not isinstance(value, np.ndarray):
                setattr(self, name, np.asarray(value))

    @abstractmethod
    def get_m_dot_in(self) -> float:
        """Mass flow rate into the combustion chamber [kg/s]."""
        pass

    @abstractmethod
    def run_timestep(self, *args, **kwargs) -> None:
        """
        Calculates and stores operational parameters in the corresponding
        vectors.

        This method will depend on the motor category. While a SRM will have
        to use/store operational parameters such as burn area and propellant
        volume, a HRE or LRE would not have to.

        When executed, the method must increment the necessary attributes
        according to a differential property (time, distance or other).
        """
        pass

    @abstractmethod
    def print_results(self) -> None:
        """
        Prints results obtained during simulation/operation.
        """
        pass

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

    @property
    def thrust_time(self) -> float:
        """Total time of thrust production [s]."""
        if self._thrust_time is None:
            raise ValueError("Thrust time has not been set, run the simulation.")

        return self._thrust_time

initial_propellant_mass property

Get the initial propellant mass [kg].

thrust_time property

Total time of thrust production [s].

__init__(motor, initial_pressure, initial_atmospheric_pressure, other_losses)

Initializes attributes for the motor operation. Each motor category will contain a particular set of attributes.

Source code in machwave/states/base.py 
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 __init__(
    self,
    motor: Motor,
    initial_pressure: float,
    initial_atmospheric_pressure: float,
    other_losses: float,
) -> None:
    """
    Initializes attributes for the motor operation.
    Each motor category will contain a particular set of attributes.
    """
    self.motor = motor
    self.other_losses = other_losses

    self.t: SimulationArray = [0.0]

    self.propellant_mass: SimulationArray = [motor.initial_propellant_mass]
    self.chamber_pressure: SimulationArray = [initial_pressure]
    self.exit_pressure: SimulationArray = [initial_atmospheric_pressure]

    self.thrust_coefficient: SimulationArray = [0.0]
    self.thrust_coefficient_ideal: SimulationArray = [0.0]
    self.thrust: SimulationArray = [0.0]

    self._thrust_time = None

    self.end_thrust = False
    self.end_burn = False

convert_simulation_arrays_to_numpy()

Convert accumulated lists into ndarrays for downstream consumers.

Source code in machwave/states/base.py 
58
59
60
61
62
63
def convert_simulation_arrays_to_numpy(self) -> None:
    """Convert accumulated lists into ndarrays for downstream consumers."""
    for name in self.SIMULATION_ARRAY_ATTRIBUTE_NAMES:
        value = getattr(self, name)
        if not isinstance(value, np.ndarray):
            setattr(self, name, np.asarray(value))

get_m_dot_in() abstractmethod

Mass flow rate into the combustion chamber [kg/s].

Source code in machwave/states/base.py 
65
66
67
68
@abstractmethod
def get_m_dot_in(self) -> float:
    """Mass flow rate into the combustion chamber [kg/s]."""
    pass

print_results() abstractmethod

Prints results obtained during simulation/operation.

Source code in machwave/states/base.py 
85
86
87
88
89
90
@abstractmethod
def print_results(self) -> None:
    """
    Prints results obtained during simulation/operation.
    """
    pass

run_timestep(*args, **kwargs) abstractmethod

Calculates and stores operational parameters in the corresponding vectors.

This method will depend on the motor category. While a SRM will have to use/store operational parameters such as burn area and propellant volume, a HRE or LRE would not have to.

When executed, the method must increment the necessary attributes according to a differential property (time, distance or other).

Source code in machwave/states/base.py 
70
71
72
73
74
75
76
77
78
79
80
81
82
83
@abstractmethod
def run_timestep(self, *args, **kwargs) -> None:
    """
    Calculates and stores operational parameters in the corresponding
    vectors.

    This method will depend on the motor category. While a SRM will have
    to use/store operational parameters such as burn area and propellant
    volume, a HRE or LRE would not have to.

    When executed, the method must increment the necessary attributes
    according to a differential property (time, distance or other).
    """
    pass

SolidMotorState

Bases: MotorState

State for a Solid Rocket Motor.

Source code in machwave/states/solid_motor.py 
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
class SolidMotorState(states_base.MotorState):
    """
    State for a Solid Rocket Motor.
    """

    SIMULATION_ARRAY_ATTRIBUTE_NAMES = (
        states_base.MotorState.SIMULATION_ARRAY_ATTRIBUTE_NAMES
        + (
            "free_chamber_volume",
            "web",
            "burn_area",
            "propellant_volume",
            "burn_rate",
            "divergent_loss",
            "kinetics_loss",
            "boundary_layer_loss",
            "two_phase_loss",
            "nozzle_efficiency",
            "overall_efficiency",
        )
    )

    def __init__(
        self,
        motor: motors.SolidMotor,
        initial_pressure: float,
        initial_atmospheric_pressure: float,
        other_losses: float,
    ) -> None:
        """
        Initial parameters for a SRM operation.
        """
        super().__init__(
            motor=motor,
            initial_pressure=initial_pressure,
            initial_atmospheric_pressure=initial_atmospheric_pressure,
            other_losses=other_losses,
        )

        self.motor: motors.SolidMotor = motor

        self.free_chamber_volume: states_base.SimulationArray = [
            motor.thrust_chamber.combustion_chamber.internal_volume
        ]
        self.web: states_base.SimulationArray = [0.0]
        self.burn_area: states_base.SimulationArray = [
            self.motor.grain.get_burn_area(0.0)
        ]
        self.propellant_volume: states_base.SimulationArray = [
            self.motor.grain.get_propellant_volume(0.0)
        ]
        self.burn_rate: states_base.SimulationArray = [0.0]

        initial_cog = motor.grain.get_center_of_gravity(
            web_distance=0.0,
        )
        initial_moi = motor.grain.get_moment_of_inertia(
            ideal_density=motor.propellant.ideal_density,
            web_distance=0.0,
        )
        self.propellant_cog = [initial_cog]  # [x, y, z] in meters
        self.propellant_moi = [initial_moi]  # 3x3 tensor in kg-m²

        self.divergent_loss: states_base.SimulationArray = [0.0]
        self.kinetics_loss: states_base.SimulationArray = [0.0]
        self.boundary_layer_loss: states_base.SimulationArray = [0.0]
        self.two_phase_loss: states_base.SimulationArray = [0.0]
        self.nozzle_efficiency: states_base.SimulationArray = [0.0]
        self.overall_efficiency: states_base.SimulationArray = [0.0]

    def get_m_dot_in(self) -> float:
        propellant_density = self.motor.grain.get_real_density(
            web_distance=self.web[-1],
            ideal_density=self.motor.propellant.ideal_density,
        )
        return propellant_density * self.burn_rate[-1] * self.burn_area[-1]

    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].
        """
        if self.end_thrust:
            return

        propellant_properties = self.motor.propellant.properties
        if propellant_properties is None:
            raise ValueError(
                "Propellant properties must be defined to run the simulation."
            )

        nozzle = self.motor.thrust_chamber.nozzle
        ideal_density = self.motor.propellant.ideal_density

        self.t.append(self.t[-1] + d_t)

        web_last = self.web[-1]
        self.burn_area.append(self.motor.grain.get_burn_area(web_last))
        self.propellant_volume.append(self.motor.grain.get_propellant_volume(web_last))
        burn_rate = self.motor.propellant.get_burn_rate(self.chamber_pressure[-1])
        self.burn_rate.append(burn_rate)
        self.web.append(web_last + burn_rate * (self.t[-1] - self.t[-2]))

        self.free_chamber_volume.append(
            self.motor.get_free_chamber_volume(self.propellant_volume[-1])
        )
        self.propellant_mass.append(
            self.motor.grain.get_propellant_mass(
                web_distance=self.web[-1],
                ideal_density=ideal_density,
            )
        )

        self.propellant_cog.append(
            self.motor.grain.get_center_of_gravity(web_distance=self.web[-1])
        )
        self.propellant_moi.append(
            self.motor.grain.get_moment_of_inertia(
                ideal_density=ideal_density,
                web_distance=self.web[-1],
            )
        )

        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=self.get_m_dot_in(),
            free_chamber_volume=self.free_chamber_volume[-1],
            throat_area=nozzle.get_throat_area(),
            k=propellant_properties.k_chamber,
            R=propellant_properties.R_chamber,
            flame_temperature=propellant_properties.adiabatic_flame_temperature,
            discharge_coefficient=nozzle.discharge_coefficient,
        )[0]
        self.chamber_pressure.append(new_chamber_pressure)
        self.exit_pressure.append(
            isentropic.get_exit_pressure(
                propellant_properties.k_exhaust,
                nozzle.expansion_ratio,
                new_chamber_pressure,
            )
        )

        chamber_pressure_psi = conversions.convert_pa_to_psi(new_chamber_pressure)
        throat_diameter_inch = conversions.convert_meter_to_inch(nozzle.throat_diameter)
        divergent_loss = losses.get_nozzle_divergent_loss_fraction(
            divergent_angle=nozzle.divergent_angle,
        )
        kinetics_loss = losses.get_kinetics_loss_fraction(
            i_sp_th_frozen=propellant_properties.i_sp_frozen,
            i_sp_th_shifting=propellant_properties.i_sp_shifting,
            chamber_pressure_psi=chamber_pressure_psi,
        )
        boundary_layer_loss = losses.get_boundary_layer_loss_fraction(
            chamber_pressure_psi=chamber_pressure_psi,
            throat_diameter_inch=throat_diameter_inch,
            expansion_ratio=nozzle.expansion_ratio,
            time=self.t[-1],
            c_1=nozzle.c_1,
            c_2=nozzle.c_2,
        )
        two_phase_loss = losses.get_two_phase_flow_loss_fraction(
            chamber_pressure_psi=chamber_pressure_psi,
            mass_fraction_of_condensed_phase=propellant_properties.qsi_chamber,
            expansion_ratio=nozzle.expansion_ratio,
            throat_diameter_inch=throat_diameter_inch,
            characteristic_length_inch=conversions.convert_meter_to_inch(
                self.free_chamber_volume[-1] / nozzle.get_throat_area()
            ),
        )
        nozzle_efficiency = losses.get_overall_nozzle_efficiency(
            divergent_loss,
            kinetics_loss,
            boundary_layer_loss,
            two_phase_loss,
            other_losses=self.other_losses,
        )
        overall_efficiency = (
            nozzle_efficiency * self.motor.propellant.combustion_efficiency
        )
        self.divergent_loss.append(divergent_loss)
        self.kinetics_loss.append(kinetics_loss)
        self.boundary_layer_loss.append(boundary_layer_loss)
        self.two_phase_loss.append(two_phase_loss)
        self.nozzle_efficiency.append(nozzle_efficiency)
        self.overall_efficiency.append(overall_efficiency)

        thrust_coefficient_ideal = nozzle_core.get_ideal_thrust_coefficient(
            new_chamber_pressure,
            self.exit_pressure[-1],
            external_pressure,
            nozzle.expansion_ratio,
            propellant_properties.k_exhaust,
        )
        thrust_coefficient = nozzle_core.apply_thrust_coefficient_correction(
            thrust_coefficient_ideal, overall_efficiency
        )
        self.thrust_coefficient.append(thrust_coefficient)
        self.thrust_coefficient_ideal.append(thrust_coefficient_ideal)
        self.thrust.append(
            nozzle_core.get_thrust_from_thrust_coefficient(
                thrust_coefficient, new_chamber_pressure, nozzle.get_throat_area()
            )
        )

        if self.propellant_mass[-1] <= 0 and not self.end_burn:
            self.burn_time = self.t[-1]
            self.end_burn = True
        if not isentropic.is_flow_choked(
            new_chamber_pressure,
            external_pressure,
            isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
        ):
            self._thrust_time = self.t[-1]
            self.end_thrust = True

    def print_results(self) -> None:
        """
        Prints the results obtained during the SRM operation.
        """
        print("\nBURN REGRESSION")
        if self.propellant_mass[0] > 1:
            print(f" Propellant initial mass {self.propellant_mass[0]:.3f} kg")
        else:
            print(f" Propellant initial mass {self.propellant_mass[0] * 1e3:.3f} g")
        print(" Mean Kn: %.2f" % np.mean(self.klemmung))
        print(" Max Kn: %.2f" % np.max(self.klemmung))
        print(f" Initial to final Kn ratio: {self.initial_to_final_klemmung_ratio:.3f}")
        print(f" Volumetric efficiency: {self.volumetric_efficiency:.3%}")
        print(" Burn profile: " + self.burn_profile)
        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"
        )

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

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

        print("\nNOZZLE DESIGN")
        print(f" Average nozzle efficiency: {np.mean(self.nozzle_efficiency):.3%}")
        print(f" Average overall efficiency: {np.mean(self.overall_efficiency):.3%}")
        print(f" Divergent nozzle loss fraction: {np.mean(self.divergent_loss):.3%}")
        print(f" Average kinetics loss fraction: {np.mean(self.kinetics_loss):.3%}")
        print(
            f" Average boundary layer loss fraction: {np.mean(self.boundary_layer_loss):.3%}"
        )
        print(
            f" Average two-phase flow loss fraction: {np.mean(self.two_phase_loss):.3%}"
        )

    @property
    def klemmung(self) -> np.ndarray:
        """Get the klemmung values."""
        burn_area = np.asarray(self.burn_area)
        return (
            burn_area[burn_area > 0]
            / self.motor.thrust_chamber.nozzle.get_throat_area()
        )

    @property
    def initial_to_final_klemmung_ratio(self) -> float:
        """Get the ratio of the initial to final klemmung."""
        return self.klemmung[0] / self.klemmung[-1]

    @property
    def volumetric_efficiency(self) -> float:
        """Get the volumetric efficiency."""
        return (
            self.propellant_volume[0]
            / self.motor.thrust_chamber.combustion_chamber.internal_volume
        )

    @property
    def burn_profile(self) -> str:
        """Get the burn profile.

        Returns:
            Burn profile: "regressive", "progressive", or "neutral".
        """
        deviancy = 0.02
        burn_area_arr = np.asarray(self.burn_area)
        burn_area = burn_area_arr[burn_area_arr > 0]

        if burn_area[0] / burn_area[-1] > 1 + deviancy:
            return "regressive"
        elif burn_area[0] / burn_area[-1] < 1 - deviancy:
            return "progressive"
        else:
            return "neutral"

    @property
    def max_mass_flux(self) -> float:
        """Get the maximum mass flux."""
        return np.max(self.grain_mass_flux)

    @property
    def grain_mass_flux(self) -> np.ndarray:
        """Get the grain mass flux."""
        return self.motor.grain.get_mass_flux_per_segment(
            np.asarray(self.burn_rate),
            self.motor.propellant.ideal_density,
            np.asarray(self.web),
        )

    @property
    def total_impulse(self) -> float:
        """Get the total impulse [N-s]."""
        return performance.get_total_impulse(
            np.asarray(self.thrust), np.asarray(self.t)
        )

    @property
    def specific_impulse(self) -> float:
        """Get the specific impulse [s]."""
        return performance.get_specific_impulse(
            total_impulse=self.total_impulse,
            initial_propellant_mass=self.propellant_mass[0],
        )

burn_profile property

Get the burn profile.

Returns:

Type Description
str

Burn profile: "regressive", "progressive", or "neutral".

grain_mass_flux property

Get the grain mass flux.

initial_to_final_klemmung_ratio property

Get the ratio of the initial to final klemmung.

klemmung property

Get the klemmung values.

max_mass_flux property

Get the maximum mass flux.

specific_impulse property

Get the specific impulse [s].

total_impulse property

Get the total impulse [N-s].

volumetric_efficiency property

Get the volumetric efficiency.

__init__(motor, initial_pressure, initial_atmospheric_pressure, other_losses)

Initial parameters for a SRM operation.

Source code in machwave/states/solid_motor.py 
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
def __init__(
    self,
    motor: motors.SolidMotor,
    initial_pressure: float,
    initial_atmospheric_pressure: float,
    other_losses: float,
) -> None:
    """
    Initial parameters for a SRM operation.
    """
    super().__init__(
        motor=motor,
        initial_pressure=initial_pressure,
        initial_atmospheric_pressure=initial_atmospheric_pressure,
        other_losses=other_losses,
    )

    self.motor: motors.SolidMotor = motor

    self.free_chamber_volume: states_base.SimulationArray = [
        motor.thrust_chamber.combustion_chamber.internal_volume
    ]
    self.web: states_base.SimulationArray = [0.0]
    self.burn_area: states_base.SimulationArray = [
        self.motor.grain.get_burn_area(0.0)
    ]
    self.propellant_volume: states_base.SimulationArray = [
        self.motor.grain.get_propellant_volume(0.0)
    ]
    self.burn_rate: states_base.SimulationArray = [0.0]

    initial_cog = motor.grain.get_center_of_gravity(
        web_distance=0.0,
    )
    initial_moi = motor.grain.get_moment_of_inertia(
        ideal_density=motor.propellant.ideal_density,
        web_distance=0.0,
    )
    self.propellant_cog = [initial_cog]  # [x, y, z] in meters
    self.propellant_moi = [initial_moi]  # 3x3 tensor in kg-m²

    self.divergent_loss: states_base.SimulationArray = [0.0]
    self.kinetics_loss: states_base.SimulationArray = [0.0]
    self.boundary_layer_loss: states_base.SimulationArray = [0.0]
    self.two_phase_loss: states_base.SimulationArray = [0.0]
    self.nozzle_efficiency: states_base.SimulationArray = [0.0]
    self.overall_efficiency: states_base.SimulationArray = [0.0]

print_results()

Prints the results obtained during the SRM operation.

Source code in machwave/states/solid_motor.py 
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
def print_results(self) -> None:
    """
    Prints the results obtained during the SRM operation.
    """
    print("\nBURN REGRESSION")
    if self.propellant_mass[0] > 1:
        print(f" Propellant initial mass {self.propellant_mass[0]:.3f} kg")
    else:
        print(f" Propellant initial mass {self.propellant_mass[0] * 1e3:.3f} g")
    print(" Mean Kn: %.2f" % np.mean(self.klemmung))
    print(" Max Kn: %.2f" % np.max(self.klemmung))
    print(f" Initial to final Kn ratio: {self.initial_to_final_klemmung_ratio:.3f}")
    print(f" Volumetric efficiency: {self.volumetric_efficiency:.3%}")
    print(" Burn profile: " + self.burn_profile)
    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"
    )

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

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

    print("\nNOZZLE DESIGN")
    print(f" Average nozzle efficiency: {np.mean(self.nozzle_efficiency):.3%}")
    print(f" Average overall efficiency: {np.mean(self.overall_efficiency):.3%}")
    print(f" Divergent nozzle loss fraction: {np.mean(self.divergent_loss):.3%}")
    print(f" Average kinetics loss fraction: {np.mean(self.kinetics_loss):.3%}")
    print(
        f" Average boundary layer loss fraction: {np.mean(self.boundary_layer_loss):.3%}"
    )
    print(
        f" Average two-phase flow loss fraction: {np.mean(self.two_phase_loss):.3%}"
    )

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/states/solid_motor.py 
 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
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].
    """
    if self.end_thrust:
        return

    propellant_properties = self.motor.propellant.properties
    if propellant_properties is None:
        raise ValueError(
            "Propellant properties must be defined to run the simulation."
        )

    nozzle = self.motor.thrust_chamber.nozzle
    ideal_density = self.motor.propellant.ideal_density

    self.t.append(self.t[-1] + d_t)

    web_last = self.web[-1]
    self.burn_area.append(self.motor.grain.get_burn_area(web_last))
    self.propellant_volume.append(self.motor.grain.get_propellant_volume(web_last))
    burn_rate = self.motor.propellant.get_burn_rate(self.chamber_pressure[-1])
    self.burn_rate.append(burn_rate)
    self.web.append(web_last + burn_rate * (self.t[-1] - self.t[-2]))

    self.free_chamber_volume.append(
        self.motor.get_free_chamber_volume(self.propellant_volume[-1])
    )
    self.propellant_mass.append(
        self.motor.grain.get_propellant_mass(
            web_distance=self.web[-1],
            ideal_density=ideal_density,
        )
    )

    self.propellant_cog.append(
        self.motor.grain.get_center_of_gravity(web_distance=self.web[-1])
    )
    self.propellant_moi.append(
        self.motor.grain.get_moment_of_inertia(
            ideal_density=ideal_density,
            web_distance=self.web[-1],
        )
    )

    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=self.get_m_dot_in(),
        free_chamber_volume=self.free_chamber_volume[-1],
        throat_area=nozzle.get_throat_area(),
        k=propellant_properties.k_chamber,
        R=propellant_properties.R_chamber,
        flame_temperature=propellant_properties.adiabatic_flame_temperature,
        discharge_coefficient=nozzle.discharge_coefficient,
    )[0]
    self.chamber_pressure.append(new_chamber_pressure)
    self.exit_pressure.append(
        isentropic.get_exit_pressure(
            propellant_properties.k_exhaust,
            nozzle.expansion_ratio,
            new_chamber_pressure,
        )
    )

    chamber_pressure_psi = conversions.convert_pa_to_psi(new_chamber_pressure)
    throat_diameter_inch = conversions.convert_meter_to_inch(nozzle.throat_diameter)
    divergent_loss = losses.get_nozzle_divergent_loss_fraction(
        divergent_angle=nozzle.divergent_angle,
    )
    kinetics_loss = losses.get_kinetics_loss_fraction(
        i_sp_th_frozen=propellant_properties.i_sp_frozen,
        i_sp_th_shifting=propellant_properties.i_sp_shifting,
        chamber_pressure_psi=chamber_pressure_psi,
    )
    boundary_layer_loss = losses.get_boundary_layer_loss_fraction(
        chamber_pressure_psi=chamber_pressure_psi,
        throat_diameter_inch=throat_diameter_inch,
        expansion_ratio=nozzle.expansion_ratio,
        time=self.t[-1],
        c_1=nozzle.c_1,
        c_2=nozzle.c_2,
    )
    two_phase_loss = losses.get_two_phase_flow_loss_fraction(
        chamber_pressure_psi=chamber_pressure_psi,
        mass_fraction_of_condensed_phase=propellant_properties.qsi_chamber,
        expansion_ratio=nozzle.expansion_ratio,
        throat_diameter_inch=throat_diameter_inch,
        characteristic_length_inch=conversions.convert_meter_to_inch(
            self.free_chamber_volume[-1] / nozzle.get_throat_area()
        ),
    )
    nozzle_efficiency = losses.get_overall_nozzle_efficiency(
        divergent_loss,
        kinetics_loss,
        boundary_layer_loss,
        two_phase_loss,
        other_losses=self.other_losses,
    )
    overall_efficiency = (
        nozzle_efficiency * self.motor.propellant.combustion_efficiency
    )
    self.divergent_loss.append(divergent_loss)
    self.kinetics_loss.append(kinetics_loss)
    self.boundary_layer_loss.append(boundary_layer_loss)
    self.two_phase_loss.append(two_phase_loss)
    self.nozzle_efficiency.append(nozzle_efficiency)
    self.overall_efficiency.append(overall_efficiency)

    thrust_coefficient_ideal = nozzle_core.get_ideal_thrust_coefficient(
        new_chamber_pressure,
        self.exit_pressure[-1],
        external_pressure,
        nozzle.expansion_ratio,
        propellant_properties.k_exhaust,
    )
    thrust_coefficient = nozzle_core.apply_thrust_coefficient_correction(
        thrust_coefficient_ideal, overall_efficiency
    )
    self.thrust_coefficient.append(thrust_coefficient)
    self.thrust_coefficient_ideal.append(thrust_coefficient_ideal)
    self.thrust.append(
        nozzle_core.get_thrust_from_thrust_coefficient(
            thrust_coefficient, new_chamber_pressure, nozzle.get_throat_area()
        )
    )

    if self.propellant_mass[-1] <= 0 and not self.end_burn:
        self.burn_time = self.t[-1]
        self.end_burn = True
    if not isentropic.is_flow_choked(
        new_chamber_pressure,
        external_pressure,
        isentropic.get_critical_pressure_ratio(propellant_properties.k_chamber),
    ):
        self._thrust_time = self.t[-1]
        self.end_thrust = True