Skip to content

states.internal_ballistics

machwave.states.internal_ballistics

LiquidEngineState

Bases: MotorState

State for a Liquid Rocket Engine.

The variable names correspond to what they are commonly referred to in books and papers related to Rocket Propulsion. Therefore, PEP8's snake_case will not be followed rigorously.

Source code in machwave/states/internal_ballistics/liquid_engine.py 
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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
class LiquidEngineState(MotorState):
    """
    State for a Liquid Rocket Engine.

    The variable names correspond to what they are commonly referred to in books and papers related to
    Rocket Propulsion. Therefore, PEP8's snake_case will not be followed rigorously.
    """

    motor: LiquidEngine

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

        self.oxidizer_mass = np.array([motor.feed_system.oxidizer_tank.fluid_mass])
        self.fuel_mass = np.array([motor.feed_system.fuel_tank.fluid_mass])
        self.n_cf = np.array([1.0])  # thrust coefficient correction factor

        self.fuel_tank_pressure = np.array([motor.feed_system.fuel_tank.get_pressure()])
        self.oxidizer_tank_pressure = np.array(
            [motor.feed_system.oxidizer_tank.get_pressure()]
        )

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

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

        self._append_time(d_t)
        m_dot_fuel, m_dot_ox = self._compute_nominal_mass_flows()
        m_dot_fuel, m_dot_ox = self._clamp_mass_flows(m_dot_fuel, m_dot_ox, d_t)
        m_dot_fuel, m_dot_ox = self._adjust_flows_for_stoichiometry(
            m_dot_fuel, m_dot_ox, d_t
        )
        self._update_propellant_properties()

        new_P = self._compute_chamber_pressure(d_t, m_dot_fuel, m_dot_ox)
        self._append_chamber_pressure(new_P)
        P_exit = self._compute_exit_pressure()
        self._append_exit_pressure(P_exit)
        self._append_cf_correction(1.0)
        cf, cf_ideal = self._compute_thrust_coefficients(P_exit, P_ext)
        self._append_thrust(cf, cf_ideal)

        self._update_propellant_masses(m_dot_fuel, m_dot_ox, d_t)
        self._update_tank_pressures()
        self._check_burn_end()
        self._check_thrust_end(P_ext)

    def _append_time(self, d_t: float) -> None:
        self.t = np.append(self.t, self.t[-1] + d_t)

    def _update_propellant_properties(self) -> None:
        self.motor.propellant.properties = self.motor.propellant.evaluate(
            chamber_pressure=self.P_0[-1],
            expansion_ratio=self.motor.thrust_chamber.nozzle.expansion_ratio,
        )

    def _compute_nominal_mass_flows(self) -> tuple[float, float]:
        if self.m_prop[-1] > 0:
            fuel_flow = self.motor.feed_system.get_mass_flow_fuel(
                chamber_pressure=self.P_0[-1],
                discharge_coefficient=self.motor.thrust_chamber.injector.discharge_coefficient_fuel,
                injector_area=self.motor.thrust_chamber.injector.area_fuel,
            )
            ox_flow = self.motor.feed_system.get_mass_flow_ox(
                chamber_pressure=self.P_0[-1],
                discharge_coefficient=self.motor.thrust_chamber.injector.discharge_coefficient_oxidizer,
                injector_area=self.motor.thrust_chamber.injector.area_ox,
            )
        else:
            fuel_flow = ox_flow = 0.0
        return fuel_flow, ox_flow

    def _clamp_mass_flows(
        self, m_dot_fuel: float, m_dot_ox: float, d_t: float
    ) -> tuple[float, float]:
        max_fuel_rate = self.fuel_mass[-1] / d_t
        max_ox_rate = self.oxidizer_mass[-1] / d_t
        return min(m_dot_fuel, max_fuel_rate), min(m_dot_ox, max_ox_rate)

    def _adjust_flows_for_stoichiometry(
        self, m_dot_fuel: float, m_dot_ox: float, d_t: float
    ) -> tuple[float, float]:
        of = self.motor.propellant.of_ratio
        fuel_last = self.fuel_mass[-1]
        ox_last = self.oxidizer_mass[-1]
        # convert to consumed mass this step
        cons_f = m_dot_fuel * d_t
        cons_o = m_dot_ox * d_t
        # limiting reagent
        if cons_f >= fuel_last:
            cons_f = fuel_last
            cons_o = of * cons_f
            self.end_burn = True
            self.burn_time = self.t[-1]
        elif cons_o >= ox_last:
            cons_o = ox_last
            cons_f = cons_o / of
            self.end_burn = True
            self.burn_time = self.t[-1]
        # back to rates
        return cons_f / d_t, cons_o / d_t

    def _compute_chamber_pressure(
        self,
        d_t: float,
        m_dot_fuel: float,
        m_dot_ox: float,
    ) -> float:
        props = self.motor.propellant.properties
        assert props is not None
        return rk4th_ode_solver(
            variables={"P0": self.P_0[-1]},
            equation=compute_chamber_pressure_mass_balance_lre,
            d_t=d_t,
            R=props.R_chamber,
            T0=props.adiabatic_flame_temperature,
            V0=self.motor.thrust_chamber.combustion_chamber.internal_volume,
            At=self.motor.thrust_chamber.nozzle.get_throat_area(),
            k=props.gamma_chamber,
            m_dot_ox=m_dot_ox,
            m_dot_fuel=m_dot_fuel,
        )[0]

    def _append_chamber_pressure(self, pressure: float) -> None:
        self.P_0 = np.append(self.P_0, pressure)

    def _compute_exit_pressure(self) -> float:
        assert self.motor.propellant.properties is not None
        return get_exit_pressure(
            self.motor.propellant.properties.gamma_exhaust,
            self.motor.thrust_chamber.nozzle.expansion_ratio,
            self.P_0[-1],
        )

    def _append_exit_pressure(self, pressure: float) -> None:
        self.P_exit = np.append(self.P_exit, pressure)

    def _append_cf_correction(self, value: float) -> None:
        self.n_cf = np.append(self.n_cf, value)

    def _compute_thrust_coefficients(
        self,
        exit_pressure: float,
        P_ext: float,
    ) -> tuple[float, float]:
        assert self.motor.propellant.properties is not None
        cf_ideal = get_ideal_thrust_coefficient(
            self.P_0[-1],
            exit_pressure,
            P_ext,
            self.motor.thrust_chamber.nozzle.expansion_ratio,
            self.motor.propellant.properties.gamma_exhaust,
        )
        cf = apply_thrust_coefficient_correction(cf_ideal, self.n_cf[-1])
        return cf, cf_ideal

    def _append_thrust(self, cf: float, cf_ideal: float) -> None:
        thrust_val = get_thrust_from_thrust_coefficient(
            cf,
            self.P_0[-1],
            self.motor.thrust_chamber.nozzle.get_throat_area(),
        )
        self.C_f = np.append(self.C_f, cf)
        self.C_f_ideal = np.append(self.C_f_ideal, cf_ideal)
        self.thrust = np.append(self.thrust, thrust_val)

    def _update_propellant_masses(
        self, m_dot_fuel: float, m_dot_ox: float, d_t: float
    ) -> None:
        consumed_f = m_dot_fuel * d_t
        consumed_o = m_dot_ox * d_t

        self.motor.feed_system.fuel_tank.remove_propellant(consumed_f)
        self.motor.feed_system.oxidizer_tank.remove_propellant(consumed_o)

        new_fuel = self.fuel_mass[-1] - consumed_f
        new_ox = self.oxidizer_mass[-1] - consumed_o
        self.fuel_mass = np.append(self.fuel_mass, new_fuel)
        self.oxidizer_mass = np.append(self.oxidizer_mass, new_ox)
        self.m_prop = np.append(self.m_prop, new_fuel + new_ox)

    def _update_tank_pressures(self) -> None:
        new_fuel_tank_pressure = self.motor.feed_system.get_fuel_tank_pressure()
        new_oxidizer_tank_pressure = self.motor.feed_system.get_oxidizer_tank_pressure()

        self.fuel_tank_pressure = np.append(
            self.fuel_tank_pressure, new_fuel_tank_pressure
        )
        self.oxidizer_tank_pressure = np.append(
            self.oxidizer_tank_pressure, new_oxidizer_tank_pressure
        )

    def _check_burn_end(self) -> None:
        if self.end_burn:
            return
        if self.m_prop[-1] <= 0:
            self.end_burn = True
            self.burn_time = self.t[-1]

    def _check_thrust_end(self, P_ext: float) -> None:
        assert self.motor.propellant.properties is not None
        if not is_flow_choked(
            self.P_0[-1],
            P_ext,
            get_critical_pressure_ratio(self.motor.propellant.properties.gamma_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.m_prop[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.P_0) * 1e-6:.4f}")
        print(f"  Mean: {np.mean(self.P_0) * 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 = np.trapz(self.thrust, self.t)
        isp = impulse / (self.m_prop[0] * 9.81)
        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)

Initial parameters for a LRE operation.

Source code in machwave/states/internal_ballistics/liquid_engine.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
def __init__(
    self,
    motor: LiquidEngine,
    initial_pressure: float,
    initial_atmospheric_pressure: float,
) -> None:
    """
    Initial parameters for a LRE operation.
    """
    super().__init__(
        motor=motor,
        initial_pressure=initial_pressure,
        initial_atmospheric_pressure=initial_atmospheric_pressure,
    )

    self.oxidizer_mass = np.array([motor.feed_system.oxidizer_tank.fluid_mass])
    self.fuel_mass = np.array([motor.feed_system.fuel_tank.fluid_mass])
    self.n_cf = np.array([1.0])  # thrust coefficient correction factor

    self.fuel_tank_pressure = np.array([motor.feed_system.fuel_tank.get_pressure()])
    self.oxidizer_tank_pressure = np.array(
        [motor.feed_system.oxidizer_tank.get_pressure()]
    )

print_results()

Prints the results obtained during the Liquid Rocket Engine operation.

Source code in machwave/states/internal_ballistics/liquid_engine.py 
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
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.m_prop[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.P_0) * 1e-6:.4f}")
    print(f"  Mean: {np.mean(self.P_0) * 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 = np.trapz(self.thrust, self.t)
    isp = impulse / (self.m_prop[0] * 9.81)
    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, P_ext)

Advance simulation by time step d_t under external pressure P_ext.

Parameters:

Name Type Description Default
d_t float

Time step.

 required
P_ext float

External pressure.

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

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

    self._append_time(d_t)
    m_dot_fuel, m_dot_ox = self._compute_nominal_mass_flows()
    m_dot_fuel, m_dot_ox = self._clamp_mass_flows(m_dot_fuel, m_dot_ox, d_t)
    m_dot_fuel, m_dot_ox = self._adjust_flows_for_stoichiometry(
        m_dot_fuel, m_dot_ox, d_t
    )
    self._update_propellant_properties()

    new_P = self._compute_chamber_pressure(d_t, m_dot_fuel, m_dot_ox)
    self._append_chamber_pressure(new_P)
    P_exit = self._compute_exit_pressure()
    self._append_exit_pressure(P_exit)
    self._append_cf_correction(1.0)
    cf, cf_ideal = self._compute_thrust_coefficients(P_exit, P_ext)
    self._append_thrust(cf, cf_ideal)

    self._update_propellant_masses(m_dot_fuel, m_dot_ox, d_t)
    self._update_tank_pressures()
    self._check_burn_end()
    self._check_thrust_end(P_ext)

MotorState

Bases: State

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

Source code in machwave/states/internal_ballistics/base.py 
 9
10
11
12
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
class MotorState(State):
    """
    Defines a particular motor operation. Stores and processes all attributes
    obtained from the simulation.
    """

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

        self.t = np.array([0])  # time vector

        self.m_prop = np.array([motor.initial_propellant_mass])  # propellant mass
        self.P_0 = np.array([initial_pressure])  # chamber stagnation pressure
        self.P_exit = np.array([initial_atmospheric_pressure])  # exit pressure

        # Thrust coefficients and thrust:
        self.C_f = np.array([0])  # thrust coefficient
        self.C_f_ideal = np.array([0])  # ideal thrust coefficient
        self.thrust = np.array([0])  # thrust force (N)

        # Thrust time:
        self._thrust_time = None

        # If the propellant mass is non zero, 'end_thrust' must be False,
        # since there is still thrust being produced.
        # After the propellant has finished burning and the thrust chamber has
        # stopped producing supersonic flow, 'end_thrust' is changed to True
        # value and the internal ballistics section of the while loop below
        # stops running.
        self.end_thrust = False
        self.end_burn = False

    @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)

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

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

    self.t = np.array([0])  # time vector

    self.m_prop = np.array([motor.initial_propellant_mass])  # propellant mass
    self.P_0 = np.array([initial_pressure])  # chamber stagnation pressure
    self.P_exit = np.array([initial_atmospheric_pressure])  # exit pressure

    # Thrust coefficients and thrust:
    self.C_f = np.array([0])  # thrust coefficient
    self.C_f_ideal = np.array([0])  # ideal thrust coefficient
    self.thrust = np.array([0])  # thrust force (N)

    # Thrust time:
    self._thrust_time = None

    # If the propellant mass is non zero, 'end_thrust' must be False,
    # since there is still thrust being produced.
    # After the propellant has finished burning and the thrust chamber has
    # stopped producing supersonic flow, 'end_thrust' is changed to True
    # value and the internal ballistics section of the while loop below
    # stops running.
    self.end_thrust = False
    self.end_burn = False

print_results() abstractmethod

Prints results obtained during simulation/operation.

Source code in machwave/states/internal_ballistics/base.py 
65
66
67
68
69
70
@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/internal_ballistics/base.py 
50
51
52
53
54
55
56
57
58
59
60
61
62
63
@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.

The variable names correspond to what they are commonly referred to in books and papers related to Solid Rocket Propulsion.

Therefore, PEP8's snake_case will not be followed rigorously.

Source code in machwave/states/internal_ballistics/solid_motor.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
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
class SolidMotorState(ib_base.MotorState):
    """
    State for a Solid Rocket Motor.

    The variable names correspond to what they are commonly referred to in
    books and papers related to Solid Rocket Propulsion.

    Therefore, PEP8's snake_case will not be followed rigorously.
    """

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

        self.motor: motors.SolidMotor = motor

        # Grain and propellant parameters:
        self.V_0 = np.array(
            [motor.thrust_chamber.combustion_chamber.internal_volume]
        )  # empty chamber volume
        self.web = np.array([0])  # instant web thickness
        self.burn_area = np.array([self.motor.grain.get_burn_area(self.web[0])])
        self.propellant_volume = np.array(
            [self.motor.grain.get_propellant_volume(self.web[0])]
        )
        self.burn_rate = np.array([0])  # burn rate

        # Center of gravity and moment of inertia:
        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,
        )
        # Store as lists of 3D arrays (COG) and 3x3 arrays (MOI)
        self.propellant_cog = [initial_cog]  # center of gravity [x, y, z] in meters
        self.propellant_moi = [initial_moi]  # moment of inertia tensor 3x3 in kg-m²

        # Correction factors:
        self.eta_div = np.array([0])  # divergent nozzle correction factor
        self.eta_kin = np.array([0])  # kinetics correction factor
        self.eta_bl = np.array([0])  # boundary layer correction factor
        self.eta_2p = np.array([0])  # two-phase flow correction factor
        self.nozzle_efficiency = np.array([0])  # overall nozzle efficiency
        self.overall_efficiency = np.array([0])  # overall efficiency

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

        Args:
            d_t: Time increment [s].
            P_ext: External pressure [Pa].
        """
        if self.end_thrust:
            return

        self._append_time(d_t)
        self._update_grain_geometry()
        self._update_chamber_volume_and_mass()
        self._update_cog_and_moi()
        self._compute_pressure(d_t, P_ext)
        self._compute_flow(P_ext)
        self._check_burn_end()
        self._check_thrust_end(P_ext)

    def _append_time(self, d_t: float) -> None:
        self.t = np.append(self.t, self.t[-1] + d_t)

    def _update_grain_geometry(self) -> None:
        web_last = self.web[-1]
        area = self.motor.grain.get_burn_area(web_last)
        vol = self.motor.grain.get_propellant_volume(web_last)
        self.burn_area = np.append(self.burn_area, area)
        self.propellant_volume = np.append(self.propellant_volume, vol)
        self.burn_rate = np.append(
            self.burn_rate,
            self.motor.propellant.get_burn_rate(self.P_0[-1]),
        )
        dx = self.burn_rate[-1] * (self.t[-1] - self.t[-2])
        self.web = np.append(self.web, web_last + dx)

    def _update_chamber_volume_and_mass(self) -> None:
        free_vol = self.motor.get_free_chamber_volume(self.propellant_volume[-1])
        self.V_0 = np.append(self.V_0, free_vol)
        m_prop = self.motor.grain.get_propellant_mass(
            web_distance=self.web[-1],
            ideal_density=self.motor.propellant.ideal_density,
        )
        self.m_prop = np.append(self.m_prop, m_prop)

    def _update_cog_and_moi(self) -> None:
        # Update center of gravity and moment of inertia
        cog = self.motor.grain.get_center_of_gravity(
            web_distance=self.web[-1],
        )
        moi = self.motor.grain.get_moment_of_inertia(
            ideal_density=self.motor.propellant.ideal_density,
            web_distance=self.web[-1],
        )
        self.propellant_cog.append(cog)
        self.propellant_moi.append(moi)

    def _compute_pressure(self, d_t: float, P_ext: float) -> None:
        props = self.motor.propellant.properties
        assert props is not None
        new_P = rk4.rk4th_ode_solver(
            variables={"P0": self.P_0[-1]},
            equation=des.compute_chamber_pressure_mass_balance_srm,
            d_t=d_t,
            Pe=P_ext,
            Ab=self.burn_area[-1],
            V0=self.V_0[-1],
            At=self.motor.thrust_chamber.nozzle.get_throat_area(),
            pp=self.motor.grain.get_real_density(
                web_distance=self.web[-1],
                ideal_density=self.motor.propellant.ideal_density,
            ),
            k=props.gamma_chamber,
            R=props.R_chamber,
            T0=props.adiabatic_flame_temperature,
            r=self.burn_rate[-1],
        )[0]
        self.P_0 = np.append(self.P_0, new_P)
        exit_P = isentropic.get_exit_pressure(
            props.gamma_exhaust,
            self.motor.thrust_chamber.nozzle.expansion_ratio,
            new_P,
        )
        self.P_exit = np.append(self.P_exit, exit_P)

    def _compute_flow(self, P_ext: float) -> None:
        P0 = self.P_0[-1]
        chamber_pressure_psi = conversions.convert_pa_to_psi(P0)
        throat_diameter_inch = conversions.convert_meter_to_inch(
            self.motor.thrust_chamber.nozzle.throat_diameter
        )

        eta_div = losses.get_nozzle_divergent_percentage_loss(
            divergent_angle=self.motor.thrust_chamber.nozzle.divergent_angle,
        )
        props = self.motor.propellant.properties
        assert props is not None
        eta_kin = losses.get_kinetics_percentage_loss(
            i_sp_th_frozen=props.i_sp_frozen,
            i_sp_th_shifting=props.i_sp_shifting,
            chamber_pressure_psi=chamber_pressure_psi,
        )
        eta_bl = losses.get_boundary_layer_percentage_loss(
            chamber_pressure_psi=chamber_pressure_psi,
            throat_diameter_inch=throat_diameter_inch,
            expansion_ratio=self.motor.thrust_chamber.nozzle.expansion_ratio,
            time=self.t[-1],
            c_1=self.motor.thrust_chamber.nozzle.material.c_1,
            c_2=self.motor.thrust_chamber.nozzle.material.c_2,
        )
        eta_2p = losses.get_two_phase_flow_percentage_loss(
            chamber_pressure_psi=chamber_pressure_psi,
            mole_fraction_of_condensed_phase=props.qsi_chamber,
            expansion_ratio=self.motor.thrust_chamber.nozzle.expansion_ratio,
            throat_diameter_inch=throat_diameter_inch,
            characteristic_length_inch=conversions.convert_meter_to_inch(
                self.V_0[-1] / self.motor.thrust_chamber.nozzle.get_throat_area()
            ),
        )
        nozzle_efficiency = losses.get_overall_nozzle_efficiency(
            eta_div, eta_kin, eta_bl, eta_2p, other_losses=self.motor.other_losses
        )
        overall_efficiency = (
            nozzle_efficiency * self.motor.propellant.combustion_efficiency
        )

        self.eta_div = np.append(self.eta_div, eta_div)
        self.eta_kin = np.append(self.eta_kin, eta_kin)
        self.eta_bl = np.append(self.eta_bl, eta_bl)
        self.eta_2p = np.append(self.eta_2p, eta_2p)
        self.nozzle_efficiency = np.append(self.nozzle_efficiency, nozzle_efficiency)
        self.overall_efficiency = np.append(self.overall_efficiency, overall_efficiency)

        cf_ideal = nozzle.get_ideal_thrust_coefficient(
            P0,
            self.P_exit[-1],
            P_ext,
            self.motor.thrust_chamber.nozzle.expansion_ratio,
            props.gamma_exhaust,
        )
        cf = nozzle.apply_thrust_coefficient_correction(cf_ideal, overall_efficiency)
        self.C_f = np.append(self.C_f, cf)
        self.C_f_ideal = np.append(self.C_f_ideal, cf_ideal)
        thrust = nozzle.get_thrust_from_thrust_coefficient(
            cf, P0, self.motor.thrust_chamber.nozzle.get_throat_area()
        )
        self.thrust = np.append(self.thrust, thrust)

    def _check_burn_end(self) -> None:
        if self.m_prop[-1] <= 0 and not self.end_burn:
            self.burn_time = self.t[-1]
            self.end_burn = True

    def _check_thrust_end(self, P_ext: float) -> None:
        assert self.motor.propellant.properties is not None
        if not isentropic.is_flow_choked(
            self.P_0[-1],
            P_ext,
            isentropic.get_critical_pressure_ratio(
                self.motor.propellant.properties.gamma_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.m_prop[0] > 1:
            print(f" Propellant initial mass {self.m_prop[0]:.3f} kg")
        else:
            print(f" Propellant initial mass {self.m_prop[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.P_0) * 1e-6:.3f}, "
            f"{np.mean(self.P_0) * 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 correction factor: {np.mean(self.eta_div):.3f}%")
        print(f" Average kinetics correction factor: {np.mean(self.eta_kin):.3f}%")
        print(f" Average boundary layer correction factor: {np.mean(self.eta_bl):.3f}%")
        print(f" Average two-phase flow correction factor: {np.mean(self.eta_2p):.3f}%")

    @property
    def klemmung(self) -> np.ndarray:
        """Get the klemmung values."""
        return (
            self.burn_area[self.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, deviancy: float = 0.02) -> str:
        """Get the burn profile.

        Args:
            deviancy: Deviancy threshold for determining burn profile.
                Defaults to 0.02.

        Returns:
            Burn profile: "regressive", "progressive", or "neutral".
        """
        burn_area = self.burn_area[self.burn_area > 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(
            self.burn_rate,
            self.motor.propellant.ideal_density,
            self.web,
        )

    @property
    def total_impulse(self) -> float:
        """Get the total impulse [N-s]."""
        return np.mean(self.thrust) * self.t[-1]

    @property
    def specific_impulse(self) -> float:
        """Get the specific impulse [s]."""
        return self.total_impulse / self.m_prop[0] / 9.81

burn_profile property

Get the burn profile.

Parameters:

Name Type Description Default
deviancy

Deviancy threshold for determining burn profile. Defaults to 0.02.

 required

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)

Initial parameters for a SRM operation.

Source code in machwave/states/internal_ballistics/solid_motor.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
58
59
60
61
62
63
64
65
66
67
68
69
def __init__(
    self,
    motor: motors.SolidMotor,
    initial_pressure: float,
    initial_atmospheric_pressure: float,
) -> None:
    """
    Initial parameters for a SRM operation.
    """
    super().__init__(
        motor=motor,
        initial_pressure=initial_pressure,
        initial_atmospheric_pressure=initial_atmospheric_pressure,
    )

    self.motor: motors.SolidMotor = motor

    # Grain and propellant parameters:
    self.V_0 = np.array(
        [motor.thrust_chamber.combustion_chamber.internal_volume]
    )  # empty chamber volume
    self.web = np.array([0])  # instant web thickness
    self.burn_area = np.array([self.motor.grain.get_burn_area(self.web[0])])
    self.propellant_volume = np.array(
        [self.motor.grain.get_propellant_volume(self.web[0])]
    )
    self.burn_rate = np.array([0])  # burn rate

    # Center of gravity and moment of inertia:
    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,
    )
    # Store as lists of 3D arrays (COG) and 3x3 arrays (MOI)
    self.propellant_cog = [initial_cog]  # center of gravity [x, y, z] in meters
    self.propellant_moi = [initial_moi]  # moment of inertia tensor 3x3 in kg-m²

    # Correction factors:
    self.eta_div = np.array([0])  # divergent nozzle correction factor
    self.eta_kin = np.array([0])  # kinetics correction factor
    self.eta_bl = np.array([0])  # boundary layer correction factor
    self.eta_2p = np.array([0])  # two-phase flow correction factor
    self.nozzle_efficiency = np.array([0])  # overall nozzle efficiency
    self.overall_efficiency = np.array([0])  # overall efficiency

print_results()

Prints the results obtained during the SRM operation.

Source code in machwave/states/internal_ballistics/solid_motor.py 
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
def print_results(self) -> None:
    """
    Prints the results obtained during the SRM operation.
    """
    print("\nBURN REGRESSION")
    if self.m_prop[0] > 1:
        print(f" Propellant initial mass {self.m_prop[0]:.3f} kg")
    else:
        print(f" Propellant initial mass {self.m_prop[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.P_0) * 1e-6:.3f}, "
        f"{np.mean(self.P_0) * 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 correction factor: {np.mean(self.eta_div):.3f}%")
    print(f" Average kinetics correction factor: {np.mean(self.eta_kin):.3f}%")
    print(f" Average boundary layer correction factor: {np.mean(self.eta_bl):.3f}%")
    print(f" Average two-phase flow correction factor: {np.mean(self.eta_2p):.3f}%")

run_timestep(d_t, P_ext)

Iterate the motor operation by calculating operational parameters.

Parameters:

Name Type Description Default
d_t float

Time increment [s].

 required
P_ext float

External pressure [Pa].

 required
Source code in machwave/states/internal_ballistics/solid_motor.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def run_timestep(
    self,
    d_t: float,
    P_ext: float,
) -> None:
    """Iterate the motor operation by calculating operational parameters.

    Args:
        d_t: Time increment [s].
        P_ext: External pressure [Pa].
    """
    if self.end_thrust:
        return

    self._append_time(d_t)
    self._update_grain_geometry()
    self._update_chamber_volume_and_mass()
    self._update_cog_and_moi()
    self._compute_pressure(d_t, P_ext)
    self._compute_flow(P_ext)
    self._check_burn_end()
    self._check_thrust_end(P_ext)