Source code for rocketpy.motors.solid_motor

from functools import cached_property

import numpy as np
from scipy import integrate

from ..mathutils.function import Function, funcify_method, reset_funcified_methods
from ..plots.solid_motor_plots import _SolidMotorPlots
from ..prints.solid_motor_prints import _SolidMotorPrints
from .motor import Motor


[docs] class SolidMotor(Motor): """Class to specify characteristics and useful operations for solid motors. Inherits from the abstract class rocketpy.Motor. See Also -------- Motor Attributes ---------- SolidMotor.coordinate_system_orientation : str Orientation of the motor's coordinate system. The coordinate system is defined by the motor's axis of symmetry. The origin of the coordinate system may be placed anywhere along such axis, such as at the nozzle area, and must be kept the same for all other positions specified. Options are "nozzle_to_combustion_chamber" and "combustion_chamber_to_nozzle". SolidMotor.nozzle_radius : float Radius of motor nozzle outlet in meters. SolidMotor.nozzle_position : float Motor's nozzle outlet position in meters, specified in the motor's coordinate system. See :doc:`Positions and Coordinate Systems </user/positions>` for more information. SolidMotor.throat_radius : float Radius of motor nozzle throat in meters. SolidMotor.grain_number : int Number of solid grains. SolidMotor.grains_center_of_mass_position : float Position of the center of mass of the grains in meters, specified in the motor's coordinate system. See :doc:`Positions and Coordinate Systems </user/positions>` for more information. SolidMotor.grain_separation : float Distance between two grains in meters. SolidMotor.grain_density : float Density of each grain in kg/meters cubed. SolidMotor.grain_outer_radius : float Outer radius of each grain in meters. SolidMotor.grain_initial_inner_radius : float Initial inner radius of each grain in meters. SolidMotor.grain_initial_height : float Initial height of each grain in meters. SolidMotor.grain_initial_volume : float Initial volume of each grain in meters cubed. SolidMotor.grain_inner_radius : Function Inner radius of each grain in meters as a function of time. SolidMotor.grain_height : Function Height of each grain in meters as a function of time. SolidMotor.grain_initial_mass : float Initial mass of each grain in kg. SolidMotor.dry_mass : float Same as in Motor class. See the :class:`Motor <rocketpy.Motor>` docs. SolidMotor.propellant_initial_mass : float Total propellant initial mass in kg. SolidMotor.total_mass : Function Total motor mass in kg as a function of time, defined as the sum of propellant and dry mass. SolidMotor.propellant_mass : Function Total propellant mass in kg as a function of time. SolidMotor.total_mass_flow_rate : Function Time derivative of propellant total mass in kg/s as a function of time as obtained by the thrust source. SolidMotor.center_of_mass : Function Position of the motor center of mass in meters as a function of time, with respect to the motor's coordinate system. See :doc:`Positions and Coordinate Systems </user/positions>` for more information regarding the motor's coordinate system. SolidMotor.center_of_propellant_mass : Function Position of the motor propellant center of mass in meters as a function of time. See :doc:`Positions and Coordinate Systems </user/positions>` for more information regarding the motor's coordinate system. SolidMotor.I_11 : Function Component of the motor's inertia tensor relative to the e_1 axis in kg*m^2, as a function of time. The e_1 axis is the direction perpendicular to the motor body axis of symmetry, centered at the instantaneous motor center of mass. SolidMotor.I_22 : Function Component of the motor's inertia tensor relative to the e_2 axis in kg*m^2, as a function of time. The e_2 axis is the direction perpendicular to the motor body axis of symmetry, centered at the instantaneous motor center of mass. Numerically equivalent to I_11 due to symmetry. SolidMotor.I_33 : Function Component of the motor's inertia tensor relative to the e_3 axis in kg*m^2, as a function of time. The e_3 axis is the direction of the motor body axis of symmetry, centered at the instantaneous motor center of mass. SolidMotor.I_12 : Function Component of the motor's inertia tensor relative to the e_1 and e_2 axes in kg*m^2, as a function of time. See SolidMotor.I_11 and SolidMotor.I_22 for more information. SolidMotor.I_13 : Function Component of the motor's inertia tensor relative to the e_1 and e_3 axes in kg*m^2, as a function of time. See SolidMotor.I_11 and SolidMotor.I_33 for more information. SolidMotor.I_23 : Function Component of the motor's inertia tensor relative to the e_2 and e_3 axes in kg*m^2, as a function of time. See SolidMotor.I_22 and SolidMotor.I_33 for more information. SolidMotor.propellant_I_11 : Function Component of the propellant inertia tensor relative to the e_1 axis in kg*m^2, as a function of time. The e_1 axis is the direction perpendicular to the motor body axis of symmetry, centered at the instantaneous propellant center of mass. SolidMotor.propellant_I_22 : Function Component of the propellant inertia tensor relative to the e_2 axis in kg*m^2, as a function of time. The e_2 axis is the direction perpendicular to the motor body axis of symmetry, centered at the instantaneous propellant center of mass. Numerically equivalent to propellant_I_11 due to symmetry. SolidMotor.propellant_I_33 : Function Component of the propellant inertia tensor relative to the e_3 axis in kg*m^2, as a function of time. The e_3 axis is the direction of the motor body axis of symmetry, centered at the instantaneous propellant center of mass. SolidMotor.propellant_I_12 : Function Component of the propellant inertia tensor relative to the e_1 and e_2 axes in kg*m^2, as a function of time. See SolidMotor.propellant_I_11 and SolidMotor.propellant_I_22 for more information. SolidMotor.propellant_I_13 : Function Component of the propellant inertia tensor relative to the e_1 and e_3 axes in kg*m^2, as a function of time. See SolidMotor.propellant_I_11 and SolidMotor.propellant_I_33 for more information. SolidMotor.propellant_I_23 : Function Component of the propellant inertia tensor relative to the e_2 and e_3 axes in kg*m^2, as a function of time. See SolidMotor.propellant_I_22 and SolidMotor.propellant_I_33 for more information. SolidMotor.thrust : Function Motor thrust force, in Newtons, as a function of time. SolidMotor.total_impulse : float Total impulse of the thrust curve in N*s. SolidMotor.max_thrust : float Maximum thrust value of the given thrust curve, in N. SolidMotor.max_thrust_time : float Time, in seconds, in which the maximum thrust value is achieved. SolidMotor.average_thrust : float Average thrust of the motor, given in N. SolidMotor.burn_time : tuple of float Tuple containing the initial and final time of the motor's burn time in seconds. SolidMotor.burn_start_time : float Motor burn start time, in seconds. SolidMotor.burn_out_time : float Motor burn out time, in seconds. SolidMotor.burn_duration : float Total motor burn duration, in seconds. It is the difference between the ``burn_out_time`` and the ``burn_start_time``. SolidMotor.exhaust_velocity : Function Propulsion gases exhaust velocity, assumed constant, in m/s. SolidMotor.burn_area : Function Total burn area considering all grains, made out of inner cylindrical burn area and grain top and bottom faces. Expressed in meters squared as a function of time. SolidMotor.Kn : Function Motor Kn as a function of time. Defined as burn_area divided by nozzle throat cross sectional area. Has no units. SolidMotor.burn_rate : Function Propellant burn rate in meter/second as a function of time. SolidMotor.interpolate : string Method of interpolation used in case thrust curve is given by data set in .csv or .eng, or as an array. Options are 'spline' 'akima' and 'linear'. Default is "linear". """
[docs] def __init__( self, thrust_source, dry_mass, dry_inertia, nozzle_radius, grain_number, grain_density, grain_outer_radius, grain_initial_inner_radius, grain_initial_height, grain_separation, grains_center_of_mass_position, center_of_dry_mass_position, nozzle_position=0.0, burn_time=None, throat_radius=0.01, reshape_thrust_curve=False, interpolation_method="linear", coordinate_system_orientation="nozzle_to_combustion_chamber", ): """Initialize Motor class, process thrust curve and geometrical parameters and store results. Parameters ---------- thrust_source : int, float, callable, string, array, Function Motor's thrust curve. Can be given as an int or float, in which case the thrust will be considered constant in time. It can also be given as a callable function, whose argument is time in seconds and returns the thrust supplied by the motor in the instant. If a string is given, it must point to a .csv or .eng file. The .csv file can contain a single line header and the first column must specify time in seconds, while the second column specifies thrust. Arrays may also be specified, following rules set by the class Function. Thrust units are Newtons. .. seealso:: :doc:`Thrust Source Details </user/motors/thrust>` nozzle_radius : int, float Motor's nozzle outlet radius in meters. dry_mass : int, float Same as in Motor class. See the :class:`Motor <rocketpy.Motor>` docs dry_inertia : tuple, list Tuple or list containing the motor's dry mass inertia tensor components, in kg*m^2. This inertia is defined with respect to the the `center_of_dry_mass_position` position. Assuming e_3 is the rocket's axis of symmetry, e_1 and e_2 are orthogonal and form a plane perpendicular to e_3, the dry mass inertia tensor components must be given in the following order: (I_11, I_22, I_33, I_12, I_13, I_23), where I_ij is the component of the inertia tensor in the direction of e_i x e_j. Alternatively, the inertia tensor can be given as (I_11, I_22, I_33), where I_12 = I_13 = I_23 = 0. grain_number : int Number of solid grains grain_density : int, float Solid grain density in kg/m3. grain_outer_radius : int, float Solid grain outer radius in meters. grain_initial_inner_radius : int, float Solid grain initial inner radius in meters. grain_initial_height : int, float Solid grain initial height in meters. grain_separation : int, float Distance between grains, in meters. grains_center_of_mass_position : float Position of the center of mass of the grains in meters. More specifically, the coordinate of the center of mass specified in the motor's coordinate system. See :doc:`Positions and Coordinate Systems </user/positions>` for more information. center_of_dry_mass_position : int, float The position, in meters, of the motor's center of mass with respect to the motor's coordinate system when it is devoid of propellant. See :doc:`Positions and Coordinate Systems </user/positions>`. nozzle_position : int, float, optional Motor's nozzle outlet position in meters, in the motor's coordinate system. See :doc:`Positions and Coordinate Systems </user/positions>` for details. Default is 0, in which case the origin of the coordinate system is placed at the motor's nozzle outlet. burn_time: float, tuple of float, optional Motor's burn time. If a float is given, the burn time is assumed to be between 0 and the given float, in seconds. If a tuple of float is given, the burn time is assumed to be between the first and second elements of the tuple, in seconds. If not specified, automatically sourced as the range between the first- and last-time step of the motor's thrust curve. This can only be used if the motor's thrust is defined by a list of points, such as a .csv file, a .eng file or a Function instance whose source is a list. throat_radius : int, float, optional Motor's nozzle throat radius in meters. Used to calculate Kn curve. Optional if the Kn curve is not interesting. Its value does not impact trajectory simulation. reshape_thrust_curve : boolean, tuple, optional If False, the original thrust curve supplied is not altered. If a tuple is given, whose first parameter is a new burn out time and whose second parameter is a new total impulse in Ns, the thrust curve is reshaped to match the new specifications. May be useful for motors whose thrust curve shape is expected to remain similar in case the impulse and burn time varies slightly. Default is False. interpolation_method : string, optional Method of interpolation to be used in case thrust curve is given by data set in .csv or .eng, or as an array. Options are 'spline' 'akima' and 'linear'. Default is "linear". coordinate_system_orientation : string, optional Orientation of the motor's coordinate system. The coordinate system is defined by the motor's axis of symmetry. The origin of the coordinate system may be placed anywhere along such axis, such as at the nozzle area, and must be kept the same for all other positions specified. Options are "nozzle_to_combustion_chamber" and "combustion_chamber_to_nozzle". Default is "nozzle_to_combustion_chamber". Returns ------- None """ super().__init__( thrust_source, dry_mass, dry_inertia, nozzle_radius, center_of_dry_mass_position, nozzle_position, burn_time, reshape_thrust_curve, interpolation_method, coordinate_system_orientation, ) # Nozzle parameters self.throat_radius = throat_radius self.throat_area = np.pi * throat_radius**2 # Grain parameters self.grains_center_of_mass_position = grains_center_of_mass_position self.grain_number = grain_number self.grain_separation = grain_separation self.grain_density = grain_density self.grain_outer_radius = grain_outer_radius self.grain_initial_inner_radius = grain_initial_inner_radius self.grain_initial_height = grain_initial_height # Grains initial geometrical parameters self.grain_initial_volume = ( self.grain_initial_height * np.pi * (self.grain_outer_radius**2 - self.grain_initial_inner_radius**2) ) self.grain_initial_mass = self.grain_density * self.grain_initial_volume self.evaluate_geometry() # Initialize plots and prints object self.prints = _SolidMotorPrints(self) self.plots = _SolidMotorPlots(self) return None
@funcify_method("Time (s)", "Mass (kg)") def propellant_mass(self): """Evaluates the total propellant mass as a function of time. Returns ------- Function Mass of the motor, in kg. """ return self.grain_volume * self.grain_density * self.grain_number @funcify_method("Time (s)", "Grain volume (m³)") def grain_volume(self): """Evaluates the total propellant volume as a function of time. The propellant is assumed to be a cylindrical Bates grain under uniform burn. Returns ------- Function Propellant volume as a function of time. """ cross_section_area = np.pi * ( self.grain_outer_radius**2 - self.grain_inner_radius**2 ) return cross_section_area * self.grain_height @funcify_method("Time (s)", "Exhaust velocity (m/s)") def exhaust_velocity(self): """Exhaust velocity by assuming it as a constant. The formula used is total impulse/propellant initial mass. Returns ------- self.exhaust_velocity : Function Gas exhaust velocity of the motor. """ return Function( self.total_impulse / self.propellant_initial_mass ).set_discrete_based_on_model(self.thrust) @property def propellant_initial_mass(self): """Returns the initial propellant mass. Returns ------- float Initial propellant mass in kg. """ return self.grain_number * self.grain_initial_mass @property def mass_flow_rate(self): """Time derivative of propellant mass. Assumes constant exhaust velocity. The formula used is the opposite of thrust divided by exhaust velocity. Returns ------- self.mass_flow_rate : Function Time derivative of total propellant mass as a function of time. See Also -------- Motor.total_mass_flow_rate : Calculates the total mass flow rate of the motor assuming constant exhaust velocity. """ try: return self._mass_flow_rate except AttributeError: self._mass_flow_rate = self.total_mass_flow_rate return self._mass_flow_rate @mass_flow_rate.setter def mass_flow_rate(self, value): """Sets the mass flow rate of the motor. This includes all the grains burning all at once. Parameters ---------- value : Function Mass flow rate in kg/s. Returns ------- None """ self._mass_flow_rate = value.reset("Time (s)", "Grain mass flow rate (kg/s)") self.evaluate_geometry() @funcify_method("Time (s)", "Center of Propellant Mass (m)", "linear") def center_of_propellant_mass(self): """Position of the propellant center of mass as a function of time. The position is specified as a scalar, relative to the motor's coordinate system. Returns ------- Function Position of the propellant center of mass as a function of time. """ time_source = self.grain_inner_radius.x_array center_of_mass = np.full_like(time_source, self.grains_center_of_mass_position) return np.column_stack((time_source, center_of_mass))
[docs] def evaluate_geometry(self): """Calculates grain inner radius and grain height as a function of time by assuming that every propellant mass burnt is exhausted. In order to do that, a system of differential equations is solved using scipy.integrate.solve_ivp. Returns ------- None """ # Define initial conditions for integration y0 = [self.grain_initial_inner_radius, self.grain_initial_height] # Define time mesh t = self.thrust.source[:, 0] t_span = t[0], t[-1] density = self.grain_density rO = self.grain_outer_radius n_grain = self.grain_number # Define system of differential equations def geometry_dot(t, y): # Store physical parameters volume_diff = self.mass_flow_rate(t) / (n_grain * density) # Compute state vector derivative rI, h = y burn_area = 2 * np.pi * (rO**2 - rI**2 + rI * h) rI_dot = -volume_diff / burn_area h_dot = -2 * rI_dot return [rI_dot, h_dot] # Define jacobian of the system of differential equations def geometry_jacobian(t, y): # Store physical parameters volume_diff = self.mass_flow_rate(t) / (n_grain * density) # Compute jacobian rI, h = y factor = volume_diff / (2 * np.pi * (rO**2 - rI**2 + rI * h) ** 2) drI_dot_drI = factor * (h - 2 * rI) drI_dot_dh = factor * rI dh_dot_drI = -2 * drI_dot_drI dh_dot_dh = -2 * drI_dot_dh return [[drI_dot_drI, drI_dot_dh], [dh_dot_drI, dh_dot_dh]] def terminate_burn(t, y): end_function = (self.grain_outer_radius - y[0]) * y[1] return end_function terminate_burn.terminal = True terminate_burn.direction = -1 # Solve the system of differential equations sol = integrate.solve_ivp( geometry_dot, t_span, y0, jac=geometry_jacobian, events=terminate_burn, atol=1e-12, rtol=1e-11, method="LSODA", ) self.grain_burn_out = sol.t[-1] # Write down functions for innerRadius and height self.grain_inner_radius = Function( np.concatenate(([sol.t], [sol.y[0]])).transpose().tolist(), "Time (s)", "Grain Inner Radius (m)", self.interpolate, "constant", ) self.grain_height = Function( np.concatenate(([sol.t], [sol.y[1]])).transpose().tolist(), "Time (s)", "Grain Height (m)", self.interpolate, "constant", ) reset_funcified_methods(self) return None
@funcify_method("Time (s)", "burn area (m²)") def burn_area(self): """Calculates the BurnArea of the grain for each time. Assuming that the grains are cylindrical BATES grains. Returns ------- burn_area : Function Function representing the burn area progression with the time. """ burn_area = ( 2 * np.pi * ( self.grain_outer_radius**2 - self.grain_inner_radius**2 + self.grain_inner_radius * self.grain_height ) * self.grain_number ) return burn_area @funcify_method("Time (s)", "burn rate (m/s)") def burn_rate(self): """Calculates the burn_rate with respect to time. This evaluation assumes that it was already calculated the mass_dot, burn_area time series. Returns ------- burn_rate : Function Rate of progression of the inner radius during the combustion. """ return -1 * self.mass_flow_rate / (self.burn_area * self.grain_density) @cached_property def Kn(self): """Calculates the motor Kn as a function of time. Defined as burn_area divided by the nozzle throat cross sectional area. Returns ------- Kn : Function Kn as a function of time. """ Kn_source = ( np.concatenate( ( [self.grain_inner_radius.source[:, 1]], [self.burn_area.source[:, 1] / self.throat_area], ) ).transpose() ).tolist() Kn = Function( Kn_source, "Grain Inner Radius (m)", "Kn (m2/m2)", self.interpolate, "constant", ) return Kn @funcify_method("Time (s)", "Inertia I_11 (kg m²)") def propellant_I_11(self): """Inertia tensor 11 component of the propellant, the inertia is relative to the e_1 axis, centered at the instantaneous propellant center of mass. Returns ------- Function Propellant inertia tensor 11 component at time t. Notes ----- The e_1 direction is assumed to be the direction perpendicular to the motor body axis. References ---------- .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ grain_mass = self.propellant_mass / self.grain_number grain_number = self.grain_number grain_inertia11 = grain_mass * ( (1 / 4) * (self.grain_outer_radius**2 + self.grain_inner_radius**2) + (1 / 12) * self.grain_height**2 ) # Calculate each grain's distance d to propellant center of mass # Assuming each grain's COM are evenly spaced initial_value = (grain_number - 1) / 2 d = np.linspace(-initial_value, initial_value, grain_number) d = d * (self.grain_initial_height + self.grain_separation) # Calculate inertia for all grains I_11 = grain_number * grain_inertia11 + grain_mass * np.sum(d**2) return I_11 @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def propellant_I_22(self): """Inertia tensor 22 component of the propellant, the inertia is relative to the e_2 axis, centered at the instantaneous propellant center of mass. Returns ------- Function Propellant inertia tensor 22 component at time t. Notes ----- The e_2 direction is assumed to be the direction perpendicular to the motor body axis, and perpendicular to e_1. References ---------- .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ return self.propellant_I_11 @funcify_method("Time (s)", "Inertia I_33 (kg m²)") def propellant_I_33(self): """Inertia tensor 33 component of the propellant, the inertia is relative to the e_3 axis, centered at the instantaneous propellant center of mass. Returns ------- Function Propellant inertia tensor 33 component at time t. Notes ----- The e_3 direction is assumed to be the axial direction of the rocket motor. References ---------- .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ I_33 = ( (1 / 2.0) * self.propellant_mass * (self.grain_outer_radius**2 + self.grain_inner_radius**2) ) return I_33 @funcify_method("Time (s)", "Inertia I_12 (kg m²)") def propellant_I_12(self): return 0 @funcify_method("Time (s)", "Inertia I_13 (kg m²)") def propellant_I_13(self): return 0 @funcify_method("Time (s)", "Inertia I_23 (kg m²)") def propellant_I_23(self): return 0
[docs] def draw(self): """Draw a representation of the SolidMotor.""" self.plots.draw()
[docs] def info(self): """Prints out basic data about the SolidMotor.""" self.prints.all() self.plots.thrust() return None
[docs] def all_info(self): """Prints out all data and graphs available about the SolidMotor.""" self.prints.all() self.plots.all() return None