From 0765dcb8ad15af1276b8e6e76cbf9c42c98e2ba2 Mon Sep 17 00:00:00 2001 From: ayoubdsp Date: Thu, 6 Nov 2025 15:19:00 +0100 Subject: [PATCH] Implement ClusterMotor class for motor aggregation - The solver now calls `motor.get_total_thrust_vector(t)` and `motor.get_total_moment(t, rocket_cm)` to get the full 3D forces and torques from the propulsion system. - These vectors are used to solve the full rigid body equations of motion (Euler's equations) for angular acceleration. - This enables correct 6-DOF simulation of motor clusters, vectored thrust, and thrust misalignments. Update Rocket class for 3D inertia and ClusterMotor Refactors the `Rocket` class to integrate the new `ClusterMotor` and handle 3D centers of mass and dynamic inertia. - `add_motor()` is updated to detect `ClusterMotor` objects. - `add_cluster_motor()` helper method added. - `evaluate_center_of_mass` and `evaluate_center_of_dry_mass` now perform weighted 3D vector averaging for CoM. - `evaluate_dry_inertias` and `evaluate_inertias` are updated to use the new dynamic `parallel_axis_theorem` functions. - Internal attributes like `center_of_dry_mass_position` are now 3D `Vector`s instead of Z-axis scalars. Update HybridMotor to new dynamic inertia contract Refactors `HybridMotor.__init__` to comply with the new base class inertia contract. - The class now calculates the individual inertias of the liquid (tanks) and solid (grain) components relative to their own CoMs. - It then uses the new dynamic PAT functions to aggregate these inertias relative to the *total propellant center of mass*. - The result is stored in `self.propellant_I_xx_from_propellant_CM` for the `Motor` base class to consume. Update LiquidMotor to new dynamic inertia contract Refactors `LiquidMotor` to comply with the new base class inertia contract. - The class now correctly calculates the aggregated inertia of all tanks relative to the total propellant center of mass (logic handled by base class and PAT tools). - Corrects the `propellant_I_11`, `_I_22`, etc. methods to return the pre-calculated `propellant_I_xx_from_propellant_CM` attribute, preventing an incorrect double application of the Parallel Axis Theorem. SolidMotor inertia logic for dynamic contract 1. **`.eng`/`.rse` File Pre-parsing:** * The logic to read `.eng` and `.rse` files (using `Motor.import_eng` and `Motor.import_rse`) has been moved from the base class into the `SolidMotor.__init__` method. * This ensures that the `thrust_source` data array is loaded and available *before* `super().__init__` is called. 2. **Inertia Re-calculation:** * The original `SolidMotor` relied on the `motor.py` base class to handle all inertia calculations within `super().__init__`. * This was problematic because `evaluate_geometry()` (which defines the grain properties needed for inertia) was called *after* `super().__init__`. * This commit fixes this by adding a new block of code at the **end** of `SolidMotor.__init__` (after `evaluate_geometry()` has run). * This new block explicitly: 1. Assigns the now-valid propellant inertia methods (e.g., `self.propellant_I_11`) to the new "contract" attributes (e.g., `self.propellant_I_11_from_propellant_CM`). 2. Imports the dynamic `parallel_axis_theorem` tools. 3. **Re-calculates** the propellant inertia relative to the motor origin (using the PAT logic from `motor.py`). 4. **Re-calculates** the final total motor inertia (`self.I_11`, `self.I_33`, etc.) by summing the dry and (now correct) propellant inertias. This overwrites the incorrect values that were set during the initial `super().__init__` call. Refactor base Motor class for dynamic 6-DOF inertia 1. **Dynamic Inertia Calculation (Major Change):** * The `__init__` method now imports the new dynamic Parallel Axis Theorem (PAT) functions from `tools.py` (e.g., `parallel_axis_theorem_I11`, `_I12`, etc.). * A new "Inertia Contract" is established: `__init__` defines new abstract attributes (e.g., `self.propellant_I_11_from_propellant_CM = Function(0)`). * Subclasses (like `SolidMotor`, `HybridMotor`) are now *required* to provide their propellant inertia relative to their *own propellant CoM* by overriding these `_from_propellant_CM` attributes. * `__init__` (lines 280-333) immediately uses the dynamic PAT functions to calculate the propellant inertia tensor (e.g., `self.propellant_I_11`) relative to the **motor's origin**, based on these "contract" attributes. * `__init__` (lines 336-351) then calculates the **total motor inertia** (`self.I_11`, `self.I_33`, etc.) relative to the **motor's origin** by adding the (constant) dry inertia. *Note: This differs completely from the GitHub version, where inertia was calculated *later* inside the `@funcify_method` definitions (`I_11`, `I_33`, etc.) and relative to the *instantaneous center of mass*.* 2. **`.eng` File Parsing Fix:** * The `__init__` method now includes logic (lines 235-240) to detect if `thrust_source` is a `.eng` file. * If true, it sets the `delimiter` to a space (" ") and `comments` to a semicolon (";"), correcting the parsing failure that occurs in the original `Function` constructor which defaults to commas. 3. **`reference_pressure` Added:** * The `__init__` signature (line 229) now accepts `reference_pressure=None`. * This value is stored as `self.reference_pressure` (line 257) to be used in vacuum thrust calculations. 4. **Modified Inertia Methods (`I_11`, `I_33`, etc.):** * The instance methods (e.g., `@funcify_method def I_11(self)`) starting at line 641 have been modified. * While the GitHub version used these methods to calculate inertia relative to the instantaneous CoM, this version's logic is updated, although it appears redundant given that the primary inertia calculation (relative to the origin) is now finalized in `__init__`. Update Rocket.draw() to visualize ClusterMotor configurations 1. **Import `ClusterMotor`:** * The file now imports `ClusterMotor` (Line 5) to check the motor type. 2. **Refactored `_draw_motor` Method (Line 234):** * This method is completely refactored. It now acts as a dispatcher, calling the new `_generate_motor_patches` helper function. * It checks `isinstance(self.rocket.motor, ClusterMotor)`. * If it's a cluster, it adds all patches generated by the helper. * If it's a simple motor, it uses the original logic (drawing the nozzle and chamber centered at y=0). * It also correctly calls `_draw_nozzle_tube` to connect the airframe to the start of the cluster or single motor. 3. **New `_generate_motor_patches` Method (Line 259):** * This is an **entirely new** helper function. * It contains the core logic for drawing clusters. * It iterates through `cluster.motors` and `cluster.positions`. * For each sub-motor, it correctly calculates the 2D plot offset (using `sub_pos[0]` for the 'xz' plane or `sub_pos[1]` for 'yz') and the longitudinal position (`sub_pos[2]`). * It re-uses the plotting logic of the individual sub-motors (e.g., `_generate_combustion_chamber`, `_generate_grains`) but applies the correct `translate=(pos_z, offset)` transform. * This allows multiple motors to be drawn side-by-side. 4. **Fix `_draw_center_of_mass_and_pressure` (Line 389):** * This method is updated to handle the new 3D `Vector` object returned by `self.rocket.center_of_mass(0)`. * It now accesses the Z-coordinate correctly using `cm_z = float(cm_vector.z.real)` instead of assuming the return value is a simple float. Update function.py Clear comments and cleaned code Updated ClusterMotor Clean comments and cleared code Hybrid Motor correction Correction of the comments and cleared the code Correction of LiquidMotor Clean comments and cleared code Correction of Motor Clean comments and clear code Correction of SolidMotor Clean comments and cleared code Correction of Rocketplots Clean comments and cleared code from copilot Correction of Rocket Clean comments and cleared code with copilot review Correction of Flight Clean comments and clear code Correction of tools Clean comments and clear code ClusterMotor initialization and add docstring Addresses pull request feedback on the new ClusterMotor class. The monolithic __init__ method was large and difficult to maintain. This commit refactors the constructor by breaking its logic into several smaller, private helper methods: - `_validate_inputs`: Handles all input validation and normalization. - `_initialize_basic_properties`: Calculates scalar values (mass, impulse, burn time). - `_initialize_thrust_and_mass`: Sets up thrust, mass flow rate, and propellant mass Functions. - `_initialize_center_of_mass`: Sets up CoM and CoPM Functions. - `_initialize_inertia_properties`: Calculates dry inertia and sets up propellant inertia Functions. This improves readability and separates concerns. Additionally: - Adds a NumPy-style docstring to the `__init__` method. - Cleans up inline comments, retaining all essential docstrings. Refactor imports in solid_motor.py for clarity Simplify return statement for inertia tensor calculation Removing the # .real generated comment Removing the generated unuseful comment Implement deprecation warning decorator Added a deprecation decorator to warn users about obsolete functions and their alternatives. make format make lint partially fix tests few adjusts --- rocketpy/_encoders.py | 29 +- rocketpy/mathutils/inertia.py | 217 ++++++++++ rocketpy/motors/__init__.py | 1 + rocketpy/motors/cluster_motor.py | 410 +++++++++++++++++++ rocketpy/motors/hybrid_motor.py | 652 ++++++++++++++++++------------- rocketpy/motors/liquid_motor.py | 130 +----- rocketpy/motors/motor.py | 491 +++++++++++------------ rocketpy/motors/solid_motor.py | 112 +++++- rocketpy/plots/plot_helpers.py | 2 +- rocketpy/plots/rocket_plots.py | 223 +++++++---- rocketpy/rocket/rocket.py | 610 ++++++++++++----------------- rocketpy/simulation/flight.py | 271 +++++++------ rocketpy/tools.py | 1 + 13 files changed, 1903 insertions(+), 1246 deletions(-) create mode 100644 rocketpy/mathutils/inertia.py create mode 100644 rocketpy/motors/cluster_motor.py diff --git a/rocketpy/_encoders.py b/rocketpy/_encoders.py index 58eeae809..e7f93135a 100644 --- a/rocketpy/_encoders.py +++ b/rocketpy/_encoders.py @@ -3,6 +3,7 @@ import json from datetime import datetime from importlib import import_module +from inspect import Parameter, signature import numpy as np @@ -71,11 +72,12 @@ def default(self, o): encoding["signature"] = get_class_signature(o) return encoding elif hasattr(o, "to_dict"): - encoding = o.to_dict( - include_outputs=self.include_outputs, - discretize=self.discretize, - allow_pickle=self.allow_pickle, - ) + call_kwargs = { + "include_outputs": self.include_outputs, + "discretize": self.discretize, + "allow_pickle": self.allow_pickle, + } + encoding = _call_to_dict_with_supported_kwargs(o, call_kwargs) encoding = remove_circular_references(encoding) encoding["signature"] = get_class_signature(o) @@ -195,6 +197,23 @@ def set_minimal_flight_attributes(flight, obj): flight.t_initial = flight.initial_solution[0] +def _call_to_dict_with_supported_kwargs(obj, candidate_kwargs): + """Call obj.to_dict passing only supported keyword arguments.""" + method = obj.to_dict + sig = signature(method) + params = list(sig.parameters.values()) + + if any(param.kind == Parameter.VAR_KEYWORD for param in params): + return method(**candidate_kwargs) + + supported_kwargs = { + name: candidate_kwargs[name] + for name in sig.parameters + if name != "self" and name in candidate_kwargs + } + return method(**supported_kwargs) + + def get_class_signature(obj): """Returns the signature of a class so it can be identified on decoding. The signature is a dictionary with the module and diff --git a/rocketpy/mathutils/inertia.py b/rocketpy/mathutils/inertia.py new file mode 100644 index 000000000..c3420c42b --- /dev/null +++ b/rocketpy/mathutils/inertia.py @@ -0,0 +1,217 @@ +"""Utilities related to inertia tensor transformations. + +This module centralizes dynamic helpers for applying the parallel axis + theorem (PAT). It lives inside ``rocketpy.mathutils`` so that functionality + depending on :class:`rocketpy.mathutils.function.Function` does not leak into + generic utility modules such as ``rocketpy.tools``. +""" + +from rocketpy.mathutils.function import Function +from rocketpy.mathutils.vector_matrix import Vector + + +def _pat_dynamic_helper(com_inertia_moment, mass, distance_vec_3d, axes_term_lambda): + """Apply the PAT to inertia moments, supporting static and dynamic inputs.""" + + is_dynamic = ( + isinstance(com_inertia_moment, Function) + or isinstance(mass, Function) + or isinstance(distance_vec_3d, Function) + ) + + def get_val(arg, t): + return arg(t) if isinstance(arg, Function) else arg + + if not is_dynamic: + d_vec = Vector(distance_vec_3d) + mass_term = mass * axes_term_lambda(d_vec) + return com_inertia_moment + mass_term + + def new_source(t): + d_vec_t = get_val(distance_vec_3d, t) + mass_t = get_val(mass, t) + inertia_t = get_val(com_inertia_moment, t) + mass_term = mass_t * axes_term_lambda(d_vec_t) + return inertia_t + mass_term + + return Function(new_source, inputs="t", outputs="Inertia (kg*m^2)") + + +def _pat_dynamic_product_helper( + com_inertia_product, mass, distance_vec_3d, product_term_lambda +): + """Apply the PAT to inertia products, supporting static and dynamic inputs.""" + + is_dynamic = ( + isinstance(com_inertia_product, Function) + or isinstance(mass, Function) + or isinstance(distance_vec_3d, Function) + ) + + def get_val(arg, t): + return arg(t) if isinstance(arg, Function) else arg + + if not is_dynamic: + d_vec = Vector(distance_vec_3d) + mass_term = mass * product_term_lambda(d_vec) + return com_inertia_product + mass_term + + def new_source(t): + d_vec_t = get_val(distance_vec_3d, t) + mass_t = get_val(mass, t) + inertia_t = get_val(com_inertia_product, t) + mass_term = mass_t * product_term_lambda(d_vec_t) + return inertia_t + mass_term + + return Function(new_source, inputs="t", outputs="Inertia (kg*m^2)") + + +# --- Public functions for the Parallel Axis Theorem --- + + +def parallel_axis_theorem_I11(com_inertia_moment, mass, distance_vec_3d): + """Apply PAT to the I11 inertia term. + + Parameters + ---------- + com_inertia_moment : float or Function + Inertia moment relative to the component center of mass. + mass : float or Function + Mass of the component. If a Function, it must map time to mass. + distance_vec_3d : array-like or Function + Displacement vector from the component COM to the reference COM. + + Returns + ------- + float or Function + Updated I11 value referenced to the new axis. + """ + + return _pat_dynamic_helper( + com_inertia_moment, mass, distance_vec_3d, lambda d_vec: d_vec.y**2 + d_vec.z**2 + ) + + +def parallel_axis_theorem_I22(com_inertia_moment, mass, distance_vec_3d): + """Apply PAT to the I22 inertia term. + + Parameters + ---------- + com_inertia_moment : float or Function + Inertia moment relative to the component center of mass. + mass : float or Function + Mass of the component. If a Function, it must map time to mass. + distance_vec_3d : array-like or Function + Displacement vector from the component COM to the reference COM. + + Returns + ------- + float or Function + Updated I22 value referenced to the new axis. + """ + + return _pat_dynamic_helper( + com_inertia_moment, mass, distance_vec_3d, lambda d_vec: d_vec.x**2 + d_vec.z**2 + ) + + +def parallel_axis_theorem_I33(com_inertia_moment, mass, distance_vec_3d): + """Apply PAT to the I33 inertia term. + + Parameters + ---------- + com_inertia_moment : float or Function + Inertia moment relative to the component center of mass. + mass : float or Function + Mass of the component. If a Function, it must map time to mass. + distance_vec_3d : array-like or Function + Displacement vector from the component COM to the reference COM. + + Returns + ------- + float or Function + Updated I33 value referenced to the new axis. + """ + + return _pat_dynamic_helper( + com_inertia_moment, mass, distance_vec_3d, lambda d_vec: d_vec.x**2 + d_vec.y**2 + ) + + +def parallel_axis_theorem_I12(com_inertia_product, mass, distance_vec_3d): + """Apply PAT to the I12 inertia product. + + Parameters + ---------- + com_inertia_product : float or Function + Product of inertia relative to the component center of mass. + mass : float or Function + Mass of the component. If a Function, it must map time to mass. + distance_vec_3d : array-like or Function + Displacement vector from the component COM to the reference COM. + + Returns + ------- + float or Function + Updated I12 value referenced to the new axis. + """ + + return _pat_dynamic_product_helper( + com_inertia_product, mass, distance_vec_3d, lambda d_vec: d_vec.x * d_vec.y + ) + + +def parallel_axis_theorem_I13(com_inertia_product, mass, distance_vec_3d): + """Apply PAT to the I13 inertia product. + + Parameters + ---------- + com_inertia_product : float or Function + Product of inertia relative to the component center of mass. + mass : float or Function + Mass of the component. If a Function, it must map time to mass. + distance_vec_3d : array-like or Function + Displacement vector from the component COM to the reference COM. + + Returns + ------- + float or Function + Updated I13 value referenced to the new axis. + """ + + return _pat_dynamic_product_helper( + com_inertia_product, mass, distance_vec_3d, lambda d_vec: d_vec.x * d_vec.z + ) + + +def parallel_axis_theorem_I23(com_inertia_product, mass, distance_vec_3d): + """Apply PAT to the I23 inertia product. + + Parameters + ---------- + com_inertia_product : float or Function + Product of inertia relative to the component center of mass. + mass : float or Function + Mass of the component. If a Function, it must map time to mass. + distance_vec_3d : array-like or Function + Displacement vector from the component COM to the reference COM. + + Returns + ------- + float or Function + Updated I23 value referenced to the new axis. + """ + + return _pat_dynamic_product_helper( + com_inertia_product, mass, distance_vec_3d, lambda d_vec: d_vec.y * d_vec.z + ) + + +__all__ = [ + "parallel_axis_theorem_I11", + "parallel_axis_theorem_I22", + "parallel_axis_theorem_I33", + "parallel_axis_theorem_I12", + "parallel_axis_theorem_I13", + "parallel_axis_theorem_I23", +] diff --git a/rocketpy/motors/__init__.py b/rocketpy/motors/__init__.py index b13ff9392..718548a0c 100644 --- a/rocketpy/motors/__init__.py +++ b/rocketpy/motors/__init__.py @@ -1,3 +1,4 @@ +from .cluster_motor import ClusterMotor from .empty_motor import EmptyMotor from .fluid import Fluid from .hybrid_motor import HybridMotor diff --git a/rocketpy/motors/cluster_motor.py b/rocketpy/motors/cluster_motor.py new file mode 100644 index 000000000..a08f4e422 --- /dev/null +++ b/rocketpy/motors/cluster_motor.py @@ -0,0 +1,410 @@ +import matplotlib.patches as patches +import matplotlib.pyplot as plt +import numpy as np + +from rocketpy.mathutils.function import Function +from rocketpy.mathutils.inertia import ( + parallel_axis_theorem_I11, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I22, + parallel_axis_theorem_I23, + parallel_axis_theorem_I33, +) +from rocketpy.mathutils.vector_matrix import Vector + + +class ClusterMotor: + """ + Manages a cluster of motors. + + This class behaves like a single motor by aggregating the properties + of several individual motors and implementing the full interface + expected by the Rocket class. It handles the calculation of + aggregated thrust, mass, center of mass, and inertia tensor as + functions of time, considering the 3D position and orientation + of each motor. + """ + + def __init__(self, motors, positions, orientations=None): + """Initializes the ClusterMotor, aggregating multiple motors. + + Parameters + ---------- + motors : list[Motor] + A list of instantiated RocketPy Motor objects to be part of the cluster. + positions : list[tuple] or list[list] or list[np.array] + A list of 3D position vectors [x, y, z] for each motor. + The position is relative to the rocket's coordinate system origin, + which is also the cluster's origin. The coordinate system + orientation is assumed to be 'tail_to_nose'. + orientations : list[tuple] or list[list] or list[np.array], optional + A list of 3D unit vectors [x, y, z] specifying the thrust + direction for each motor in the rocket's coordinate system. + If None, all motors are assumed to thrust along the rocket's + positive Z-axis (e.g., [0, 0, 1]). Defaults to None. + """ + self._validate_inputs(motors, positions, orientations) + + self.coordinate_system_orientation = "tail_to_nose" + self._csys = 1 + + self._initialize_basic_properties() + self._initialize_thrust_and_mass() + self._initialize_center_of_mass() + self._initialize_inertia_properties() + + def _validate_inputs(self, motors, positions, orientations): + """Validates and stores the primary inputs for the cluster.""" + if not motors: + raise ValueError("The list of motors cannot be empty.") + + self.motors = motors + self.positions = [np.array(pos) for pos in positions] + + if orientations is None: + self.orientations = [np.array([0, 0, 1.0]) for _ in motors] + else: + self.orientations = [ + np.array(ori) / np.linalg.norm(ori) for ori in orientations + ] + + if not len(self.motors) == len(self.positions) == len(self.orientations): + raise ValueError( + "The 'motors', 'positions', and 'orientations' lists must have the same length." + ) + + def _initialize_basic_properties(self): + """Calculates simple aggregated scalar properties.""" + self.propellant_initial_mass = sum( + m.propellant_initial_mass for m in self.motors + ) + self.dry_mass = sum(m.dry_mass for m in self.motors) + self.total_impulse = sum(m.total_impulse for m in self.motors) + + self.burn_start_time = min(motor.burn_start_time for motor in self.motors) + self.burn_out_time = max(motor.burn_out_time for motor in self.motors) + self.burn_time = (self.burn_start_time, self.burn_out_time) + + self.nozzle_radius = np.sqrt(sum(m.nozzle_radius**2 for m in self.motors)) + self.throat_radius = np.sqrt(sum(m.throat_radius**2 for m in self.motors)) + + self.nozzle_position = np.mean( + [ + pos[2] + motor.nozzle_position * motor._csys + for motor, pos in zip(self.motors, self.positions) + ] + ) + + def _initialize_thrust_and_mass(self): + """Initializes thrust and mass-related Function objects.""" + self.thrust = Function( + self._calculate_total_thrust_scalar, inputs="t", outputs="Thrust (N)" + ) + if self.burn_time[1] > self.burn_time[0]: + self.thrust.set_discrete(self.burn_start_time, self.burn_out_time, 100) + + if self.propellant_initial_mass > 1e-9: + average_exhaust_velocity = self.total_impulse / self.propellant_initial_mass + self.total_mass_flow_rate = self.thrust / -average_exhaust_velocity + else: + self.total_mass_flow_rate = Function(0) + + self.mass_flow_rate = self.total_mass_flow_rate + + self.propellant_mass = Function( + self._calculate_propellant_mass, inputs="t", outputs="Propellant Mass (kg)" + ) + self.total_mass = Function( + self._calculate_total_mass, inputs="t", outputs="Total Mass (kg)" + ) + + def _initialize_center_of_mass(self): + """Initializes center of mass Function objects.""" + self.center_of_mass = Function( + self._calculate_center_of_mass, inputs="t", outputs="Cluster CoM Vector (m)" + ) + self.center_of_propellant_mass = Function( + self._calculate_center_of_propellant_mass, + inputs="t", + outputs="Cluster CoPM Vector (m)", + ) + + com_time = ( + self.burn_out_time + if self.burn_out_time > self.burn_start_time + else self.burn_start_time + ) + self.center_of_dry_mass_position = self._calculate_center_of_mass(com_time) + + def _initialize_inertia_properties(self): + """Initializes dry and propellant inertia properties.""" + ( + self.dry_I_11, + self.dry_I_22, + self.dry_I_33, + self.dry_I_12, + self.dry_I_13, + self.dry_I_23, + ) = self._calculate_total_dry_inertia() + + def propellant_inertia_func(t): + return self._calculate_total_propellant_inertia(t) + + self.propellant_I_11 = Function( + lambda t: propellant_inertia_func(t)[0], inputs="t", outputs="I_11 (kg*m^2)" + ) + self.propellant_I_22 = Function( + lambda t: propellant_inertia_func(t)[1], inputs="t", outputs="I_22 (kg*m^2)" + ) + self.propellant_I_33 = Function( + lambda t: propellant_inertia_func(t)[2], inputs="t", outputs="I_33 (kg*m^2)" + ) + self.propellant_I_12 = Function( + lambda t: propellant_inertia_func(t)[3], inputs="t", outputs="I_12 (kg*m^2)" + ) + self.propellant_I_13 = Function( + lambda t: propellant_inertia_func(t)[4], inputs="t", outputs="I_13 (kg*m^2)" + ) + self.propellant_I_23 = Function( + lambda t: propellant_inertia_func(t)[5], inputs="t", outputs="I_23 (kg*m^2)" + ) + + def _calculate_total_mass(self, t): + """Calculates total cluster mass at time t.""" + return self.dry_mass + self._calculate_propellant_mass(t) + + def _calculate_propellant_mass(self, t): + """ + Calculates the total propellant mass at time t by integrating the + total mass flow rate. This ensures consistency between mass and thrust. + """ + return self.propellant_initial_mass + self.total_mass_flow_rate.integral( + self.burn_start_time, t + ) + + def get_total_thrust_vector(self, t): + """ + Calculates the total thrust vector of the cluster at a given time t. + This vector is the sum of all individual motor thrust vectors, + considering their orientation. + """ + total_thrust = np.zeros(3) + for motor, orientation in zip(self.motors, self.orientations): + scalar_thrust = motor.thrust.get_value_opt(t) + total_thrust += scalar_thrust * orientation + return Vector(total_thrust) + + def _calculate_total_thrust_scalar(self, t): + """ + Calculates the magnitude of the total thrust vector. + This is what is wrapped by the `self.thrust` Function. + """ + return abs(self.get_total_thrust_vector(t)) + + def get_total_moment(self, t, ref_point): + """ + Calculates the total moment (torque) generated by the cluster + about a given reference point (e.g., the rocket's CoM). + """ + total_moment = np.zeros(3, dtype=np.float64) + ref_point_arr = np.array(ref_point, dtype=np.float64) + + for motor, pos, orientation in zip( + self.motors, self.positions, self.orientations + ): + force_magnitude = motor.thrust.get_value_opt(t) + force = force_magnitude * orientation + arm = pos - ref_point_arr + total_moment += np.cross(arm, force) + + if hasattr(motor, "thrust_eccentricity_y") and hasattr( + motor, "thrust_eccentricity_x" + ): + total_moment[0] += motor.thrust_eccentricity_y * force_magnitude + total_moment[1] -= motor.thrust_eccentricity_x * force_magnitude + + return Vector(total_moment) + + def pressure_thrust(self, pressure): + """Calculates the total pressure thrust correction for the cluster.""" + return sum(motor.pressure_thrust(pressure) for motor in self.motors) + + def _calculate_center_of_mass(self, t): + """Calculates the aggregated center of mass of the cluster at time t.""" + total_mass_val = self._calculate_total_mass(t) + if total_mass_val < 1e-9: + return Vector( + np.mean(self.positions, axis=0) if self.positions else np.zeros(3) + ) + + weighted_sum = np.zeros(3, dtype=np.float64) + for motor, pos in zip(self.motors, self.positions): + motor_com_local_z = motor.center_of_mass.get_value_opt(t) + motor_com_global = pos + np.array([0, 0, motor_com_local_z * motor._csys]) + weighted_sum += motor.total_mass.get_value_opt(t) * motor_com_global + + return Vector(weighted_sum / total_mass_val) + + def _calculate_center_of_propellant_mass(self, t): + """ + Calculates the aggregated center of mass of the cluster's propellant. + This calculation is based on the individual motor properties. + """ + total_prop_mass = 0.0 + weighted_sum = np.zeros(3, dtype=np.float64) + + for motor in self.motors: + total_prop_mass += motor.propellant_mass.get_value_opt(t) + + if total_prop_mass < 1e-9: + return self.center_of_dry_mass_position + + for motor, pos in zip(self.motors, self.positions): + prop_mass_t = motor.propellant_mass.get_value_opt(t) + if prop_mass_t > 1e-9: + prop_com_local_z = motor.center_of_propellant_mass.get_value_opt(t) + prop_com_global = pos + np.array([0, 0, prop_com_local_z * motor._csys]) + weighted_sum += prop_mass_t * prop_com_global + + return Vector(weighted_sum / total_prop_mass) + + def _calculate_total_dry_inertia(self): + """ + Calculates the cluster's total dry inertia tensor relative to the + cluster's center of dry mass. + """ + I_11, I_22, I_33 = 0.0, 0.0, 0.0 + I_12, I_13, I_23 = 0.0, 0.0, 0.0 + + ref_point = self.center_of_dry_mass_position + + for motor, pos in zip(self.motors, self.positions): + m = motor.dry_mass + motor_cdm_local_z = motor.center_of_dry_mass_position + motor_cdm_global = pos + np.array([0, 0, motor_cdm_local_z * motor._csys]) + r_vec = Vector(motor_cdm_global) - ref_point + + I_11 += parallel_axis_theorem_I11(motor.dry_I_11, m, r_vec) + I_22 += parallel_axis_theorem_I22(motor.dry_I_22, m, r_vec) + I_33 += parallel_axis_theorem_I33(motor.dry_I_33, m, r_vec) + I_12 += parallel_axis_theorem_I12(motor.dry_I_12, m, r_vec) + I_13 += parallel_axis_theorem_I13(motor.dry_I_13, m, r_vec) + I_23 += parallel_axis_theorem_I23(motor.dry_I_23, m, r_vec) + + return I_11, I_22, I_33, I_12, I_13, I_23 + + def _calculate_total_propellant_inertia(self, t): + """ + Calculates the cluster's total propellant inertia tensor relative to the + cluster's center of propellant mass at time t. + """ + I_11, I_22, I_33 = 0.0, 0.0, 0.0 + I_12, I_13, I_23 = 0.0, 0.0, 0.0 + + ref_point = self._calculate_center_of_propellant_mass(t) + + for motor, pos in zip(self.motors, self.positions): + m = motor.propellant_mass.get_value_opt(t) + if m < 1e-9: + continue + + prop_com_local_z = motor.center_of_propellant_mass.get_value_opt(t) + prop_com_global = pos + np.array([0, 0, prop_com_local_z * motor._csys]) + r_vec = Vector(prop_com_global) - ref_point + + I_11 += parallel_axis_theorem_I11( + motor.propellant_I_11.get_value_opt(t), m, r_vec + ) + I_22 += parallel_axis_theorem_I22( + motor.propellant_I_22.get_value_opt(t), m, r_vec + ) + I_33 += parallel_axis_theorem_I33( + motor.propellant_I_33.get_value_opt(t), m, r_vec + ) + I_12 += parallel_axis_theorem_I12( + motor.propellant_I_12.get_value_opt(t), m, r_vec + ) + I_13 += parallel_axis_theorem_I13( + motor.propellant_I_13.get_value_opt(t), m, r_vec + ) + I_23 += parallel_axis_theorem_I23( + motor.propellant_I_23.get_value_opt(t), m, r_vec + ) + + return I_11, I_22, I_33, I_12, I_13, I_23 + + def draw_rear_view(self, rocket_radius, tail_radius=None, filename=None): + """ + Plots a 2D rear view of the motor cluster, showing the main + rocket body and optionally the tail cone diameter. + + Args: + rocket_radius (float): The main radius of the rocket body. + tail_radius (float, optional): The rocket's radius at the tail (boattail). + If provided, a second circle will be drawn. + filename (str, optional): If provided, saves the plot to a file. + Otherwise, shows the plot. + """ + _, ax = plt.subplots(figsize=(8.5, 8.5)) + + rocket_body = patches.Circle( + (0, 0), + radius=rocket_radius, + facecolor="lightgrey", + edgecolor="black", + linewidth=2, + label=f"Main Body ({rocket_radius * 200:.0f} cm)", + ) + ax.add_patch(rocket_body) + + if tail_radius is not None: + tail_body = patches.Circle( + (0, 0), + radius=tail_radius, + facecolor="silver", + edgecolor="black", + linestyle="--", + linewidth=1.5, + label=f"Shrunk ({tail_radius * 200:.0f} cm)", + ) + ax.add_patch(tail_body) + + for i, (motor, pos) in enumerate(zip(self.motors, self.positions)): + motor_body = patches.Circle( + (pos[0], pos[1]), + radius=getattr(motor, "grain_outer_radius", motor.nozzle_radius / 2), + facecolor="dimgrey", + edgecolor="black", + linewidth=1.5, + ) + ax.add_patch(motor_body) + ax.text( + pos[0], + pos[1], + f"M{i + 1}", + ha="center", + va="center", + color="white", + fontsize=12, + fontweight="bold", + ) + + ax.set_aspect("equal", adjustable="box") + plot_limit = rocket_radius * 1.2 + ax.set_xlim(-plot_limit, plot_limit) + ax.set_ylim(-plot_limit, plot_limit) + ax.set_title("Detailed Rear View of the Cluster", fontsize=16) + ax.set_xlabel("Axis X (m)") + ax.set_ylabel("Axis Y (m)") + ax.grid(True, linestyle="--", alpha=0.6) + ax.axhline(0, color="black", linewidth=0.5) + ax.axvline(0, color="black", linewidth=0.5) + + plt.legend() + + if filename: + plt.savefig(filename) + print(f"Rear view plot saved to: {filename}") + else: + plt.show() diff --git a/rocketpy/motors/hybrid_motor.py b/rocketpy/motors/hybrid_motor.py index efee47867..e32276953 100644 --- a/rocketpy/motors/hybrid_motor.py +++ b/rocketpy/motors/hybrid_motor.py @@ -1,11 +1,16 @@ from functools import cached_property -from rocketpy.tools import parallel_axis_theorem_from_com - from ..mathutils.function import Function, funcify_method, reset_funcified_methods +from ..mathutils.inertia import ( + parallel_axis_theorem_I11, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I22, + parallel_axis_theorem_I23, + parallel_axis_theorem_I33, +) +from ..mathutils.vector_matrix import Vector from ..plots.hybrid_motor_plots import _HybridMotorPlots -from ..prints.hybrid_motor_prints import _HybridMotorPrints -from .liquid_motor import LiquidMotor from .motor import Motor from .solid_motor import SolidMotor @@ -198,200 +203,421 @@ class HybridMotor(Motor): Grain length remains constant throughout the burn. Default is True. """ - # pylint: disable=too-many-arguments def __init__( # pylint: disable=too-many-arguments self, thrust_source, dry_mass, dry_inertia, nozzle_radius, + burn_time, + center_of_dry_mass_position, grain_number, + grain_separation, grain_density, grain_outer_radius, grain_initial_inner_radius, grain_initial_height, - grain_separation, grains_center_of_mass_position, - center_of_dry_mass_position, + tanks_mass, + oxidizer_initial_mass, + oxidizer_mass_flow_rate_curve, + oxidizer_density, + oxidizer_tanks_geometries, + oxidizer_tanks_positions, + oxidizer_initial_inertia=(0, 0, 0), nozzle_position=0, - burn_time=None, - throat_radius=0.01, reshape_thrust_curve=False, interpolation_method="linear", coordinate_system_orientation="nozzle_to_combustion_chamber", reference_pressure=None, only_radial_burn=True, ): - """Initialize Motor class, process thrust curve and geometrical + """Initializes HybridMotor 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 ` + thrust_source : int, float, callable, string, numpy.ndarray, list + Motor's thrust curve. Can be given as Thrust-Time pairs array or + file path (csv or eng format). Is passed to the Function class, see + help(Function) for more information. Thrust units are Newtons, time + units are seconds. dry_mass : int, float - Same as in Motor class. See the :class:`Motor ` docs + Motor's dry mass in kg. This is the mass of the motor without + propellant. 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. + motor's center of dry mass position. Assuming e_3 is the motor'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). + Alternatively, the inertia tensor can be given as (I_11, I_22, I_33), + where I_12 = I_13 = I_23 = 0. nozzle_radius : int, float - Motor's nozzle outlet radius in meters. + Motor's nozzle radius in meters. + burn_time: int, float, tuple of int, float + Motor's burn time. + If a tuple is passed, the first value is the ignition time and the + second value is the end of burn time. If a single number is passed, + the ignition time is assumed to be 0 and the end of burn time is + the number passed. + center_of_dry_mass_position : int, float + Position of the motor's center of dry mass (i.e. center of mass + without propellant) relative to the motor's coordinate system + origin, in meters. See the ``coordinate_system_orientation`` + parameter for details on the coordinate system. grain_number : int - Number of solid grains + Number of solid grains. + grain_separation : int, float + Distance between grains, in meters. grain_density : int, float - Solid grain density in kg/m3. + Solid grain density in kg/m³. 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 ` 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 `. + grains_center_of_mass_position : int, float + Position of the center of mass of the grains relative to the motor's + coordinate system origin in meters. Generally equal to + ``center_of_dry_mass_position``. + grain_initial_inertia : tuple, list, optional + Tuple or list containing the initial inertia tensor components of a + single grain, in kg*m^2. This inertia is defined with respect to the + the grains_center_of_mass_position position. If not specified, the + grain is assumed to be a hollow cylinder with the initial dimensions. + Assuming e_3 is the grain's axis of symmetry, e_1 and e_2 are + orthogonal and form a plane perpendicular to e_3, the initial inertia + tensor components must be given in the following order: + (I_11, I_22, I_33, I_12, I_13, I_23). Alternatively, the inertia + tensor can be given as (I_11, I_22, I_33), where I_12 = I_13 = I_23 = 0. + Default is (0, 0, 0). + tanks_mass : float + Total mass of the oxidizer tanks structures in kg. Includes the mass + of the tanks themselves, valves, pipes, etc. It is assumed constant + over time. + oxidizer_initial_mass : float + Initial mass of the oxidizer, including liquid and gas phases, in kg. + oxidizer_mass_flow_rate_curve : int, float, callable, string, numpy.ndarray, list + Oxidizer mass flow rate curve. Can be given as MassFlowRate-Time + pairs array or file path (csv format). It is used to calculate the + oxidizer mass and center of mass position as a function of time. + If int or float is given, it is assumed constant. Mass flow rate + units are kg/s, time units are seconds. Passed to the Function + class, see help(Function) for more information. + oxidizer_density : float + Density of the oxidizer in kg/m³. It is used to calculate the volume + and height of the oxidizer in the tanks. It is assumed constant over + time. + oxidizer_tanks_geometries : list + List of tuples, where each tuple represents the geometry of an + oxidizer tank. Accepted geometries are: + ('cylinder', (top_radius, bottom_radius, height)) + ('sphere', radius) + ('ullage', volume) + Dimensions should be in meters and volume in cubic meters. + The list must contain at least one tank geometry. Ullage tanks can only be + placed at the top or bottom of the tanks stack. + Example: [('ullage', 0.01), ('cylinder', (0.1, 0.1, 0.5)), ('cylinder', (0.1, 0.05, 0.2))] + oxidizer_tanks_positions : list + List of floats, representing the position of the centroid of each + oxidizer tank's geometry with respect to the motor's coordinate system + origin, in meters. The list must have the same length as + ``oxidizer_tanks_geometries``. + See the ``coordinate_system_orientation`` parameter for details on the coordinate system. + oxidizer_tanks_initial_liquid_level : float, optional + Initial liquid level in the tanks, measured in meters from the bottom + of the tanks stack. If specified, this parameter overrides the initial + oxidizer mass calculation based on ``oxidizer_initial_mass``, allowing + precise control over the starting volume of the liquid oxidizer. If + None, the initial liquid level is derived from ``oxidizer_initial_mass``. + Default is None. + oxidizer_tanks_initial_ullage_mass : float, optional + Initial mass of the ullage gas in kg. If not specified, it is assumed + to be 0. Default is 0. + oxidizer_tanks_initial_ullage_volume : float, optional + Initial volume of the ullage gas in cubic meters. If not specified, it + is automatically calculated based on the tanks geometries and the + initial liquid level. Default is None. + oxidizer_initial_inertia : tuple, list, optional + Tuple or list containing the initial inertia tensor components of the + oxidizer (liquid + gas), in kg*m^2. This inertia is defined with + respect to the initial oxidizer center of mass position. If not + specified, the oxidizer is assumed to be a point mass. Default is (0, 0, 0). + Assuming e_3 is the motor's axis of symmetry, e_1 and e_2 are + orthogonal and form a plane perpendicular to e_3, the initial inertia + tensor components must be given in the following order: + (I_11, I_22, I_33, I_12, I_13, I_23). Alternatively, the inertia + tensor can be given as (I_11, I_22, I_33), where I_12 = I_13 = I_23 = 0. nozzle_position : int, float, optional - Motor's nozzle outlet position in meters, in the motor's coordinate - system. See :doc:`Positions and Coordinate Systems ` - 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. + Motor's nozzle outlet position in meters, specified in the motor's + coordinate system. Default is 0, which corresponds to the motor's + origin. See the ``coordinate_system_orientation`` parameter for + details on the coordinate system. 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. + If False, the original thrust curve supplied is used. If a tuple is + given, the thrust curve is reshaped to match the new grain mass + flow rate and burn time. The tuple should contain the initial grain + mass and the final grain mass, in kg. 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". + by data points. 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". + has its origin at the motor's center of mass and is oriented + according to the following options: + "nozzle_to_combustion_chamber" : the coordinate system is oriented + with the z-axis pointing from the nozzle towards the combustion + chamber. + "combustion_chamber_to_nozzle" : the coordinate system is oriented + with the z-axis pointing from the combustion chamber towards the + nozzle. + Default is "nozzle_to_combustion_chamber". reference_pressure : int, float, optional - Atmospheric pressure in Pa at which the thrust data was recorded. - only_radial_burn : boolean, optional - If True, inhibits the grain from burning axially, only computing - radial burn. If False, allows the grain to also burn - axially. May be useful for axially inhibited grains or hybrid motors. - Default is False. + Reference pressure in Pa used to calculate vacuum thrust and + pressure thrust. This corresponds to the atmospheric pressure + measured during the static test of the motor. Default is None, + which means no pressure correction is applied. Returns ------- None """ - super().__init__( + # Call SolidMotor init to initialize grain parameters + # Note: dry_mass and dry_inertia are temporarily set to 0, they will be + # calculated later considering the tanks mass. + SolidMotor.__init__( + self, thrust_source=thrust_source, - dry_inertia=dry_inertia, + dry_mass=0, + dry_inertia=(0, 0, 0), nozzle_radius=nozzle_radius, + burn_time=burn_time, center_of_dry_mass_position=center_of_dry_mass_position, - dry_mass=dry_mass, + grain_number=grain_number, + grain_separation=grain_separation, + grain_density=grain_density, + grain_outer_radius=grain_outer_radius, + grain_initial_inner_radius=grain_initial_inner_radius, + grain_initial_height=grain_initial_height, + grains_center_of_mass_position=grains_center_of_mass_position, nozzle_position=nozzle_position, - burn_time=burn_time, reshape_thrust_curve=reshape_thrust_curve, interpolation_method=interpolation_method, coordinate_system_orientation=coordinate_system_orientation, reference_pressure=reference_pressure, ) - self.liquid = LiquidMotor( - 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, - reference_pressure, + + # Oxidizer parameters initialization + self.tanks_mass = tanks_mass + self.oxidizer_initial_mass = oxidizer_initial_mass + self.oxidizer_density = oxidizer_density + self.oxidizer_tanks_geometries = oxidizer_tanks_geometries + self.oxidizer_tanks_positions = oxidizer_tanks_positions + self.oxidizer_initial_inertia = ( + (*oxidizer_initial_inertia, 0, 0, 0) + if len(oxidizer_initial_inertia) == 3 + else oxidizer_initial_inertia ) - self.solid = SolidMotor( - 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, - burn_time, - throat_radius, - reshape_thrust_curve, + + # Oxidizer mass flow rate definition and processing + self.oxidizer_mass_flow_rate = Function( + oxidizer_mass_flow_rate_curve, + "Time (s)", + "Oxidizer Mass Flow Rate (kg/s)", interpolation_method, - coordinate_system_orientation, - reference_pressure, - only_radial_burn, + extrapolation="zero", + ) + + # Correct dry mass and dry inertia to include tanks mass + self.dry_mass = dry_mass + tanks_mass + dry_inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia + self.dry_I_11 = dry_inertia[0] + self.dry_I_22 = dry_inertia[1] + self.dry_I_33 = dry_inertia[2] + self.dry_I_12 = dry_inertia[3] + self.dry_I_13 = dry_inertia[4] + self.dry_I_23 = dry_inertia[5] + # TODO: Calculate tanks inertia tensor based on their geometry and mass + """ + # Initialize Tanks object + self.tanks = Tank( + geometries=oxidizer_tanks_geometries, + positions=oxidizer_tanks_positions, + fluid_mass=self.oxidizer_mass, + fluid_density=self.oxidizer_density, + initial_liquid_level=oxidizer_tanks_initial_liquid_level, + initial_ullage_mass=oxidizer_tanks_initial_ullage_mass, + initial_ullage_volume=oxidizer_tanks_initial_ullage_volume, + ) + """ + # Store important functions + self.liquid_propellant_mass = self.tanks.liquid_mass + self.gas_propellant_mass = self.tanks.gas_mass + self.center_of_liquid_propellant_mass = self.tanks.liquid_center_of_mass + self.center_of_gas_propellant_mass = self.tanks.gas_center_of_mass + self.liquid_propellant_I_11 = self.tanks.liquid_I_11 + self.liquid_propellant_I_22 = self.tanks.liquid_I_22 + self.liquid_propellant_I_33 = self.tanks.liquid_I_33 + self.gas_propellant_I_11 = self.tanks.gas_I_11 + self.gas_propellant_I_22 = self.tanks.gas_I_22 + self.gas_propellant_I_33 = self.tanks.gas_I_33 + + # Rename grain attributes for clarity + self.grain_propellant_mass = self.grain_mass + self.center_of_grain_propellant_mass = self.grains_center_of_mass_position + self.grain_propellant_I_11 = self.grains_I_11 + self.grain_propellant_I_22 = self.grains_I_22 + self.grain_propellant_I_33 = self.grains_I_33 + + # Overall propellant inertia tensor components relative to propellant CoM + # We need to recalculate the total propellant CoM function first + # (Assuming self.liquid_propellant_mass and self.grain_propellant_mass exist and are Functions) + # (Assuming self.center_of_liquid_propellant_mass and self.center_of_grain_propellant_mass exist and are Functions returning scalars) + self._propellant_mass = self.liquid_propellant_mass + self.grain_propellant_mass + self._center_of_propellant_mass = ( + self.center_of_liquid_propellant_mass * self.liquid_propellant_mass + + self.center_of_grain_propellant_mass * self.grain_propellant_mass + ) / self._propellant_mass + # Ensure division by zero is handled if needed, although propellant mass shouldn't be zero initially + + # Create Functions returning distance vectors relative to the overall propellant CoM + liquid_com_to_prop_com = ( + self.center_of_liquid_propellant_mass - self._center_of_propellant_mass + ) + grain_com_to_prop_com = ( + self.center_of_grain_propellant_mass - self._center_of_propellant_mass + ) + + # Convert scalar distances to 3D vectors for PAT functions + # The distance is along the Z-axis in the motor's coordinate system + liquid_dist_vec_func = Function( + lambda t: Vector([0, 0, liquid_com_to_prop_com(t)]), inputs="t" + ) + grain_dist_vec_func = Function( + lambda t: Vector([0, 0, grain_com_to_prop_com(t)]), inputs="t" + ) + + # Apply PAT using the new specific functions + # Inertias relative to component CoMs are needed (e.g., self.liquid_propellant_I_11_from_liquid_CM) + # Assuming these exist, otherwise adjust the first argument of the PAT functions + + # --- I_11 --- + # Assuming self.liquid_propellant_I_11 refers to inertia relative to liquid CoM + liquid_I_11_prop_com = parallel_axis_theorem_I11( + self.liquid_propellant_I_11, # Inertia relative to liquid's own CoM + self.liquid_propellant_mass, + liquid_dist_vec_func, # Distance from total prop CoM to liquid CoM + ) + # Assuming self.grain_propellant_I_11 refers to inertia relative to grain CoM + grain_I_11_prop_com = parallel_axis_theorem_I11( + self.grain_propellant_I_11, # Inertia relative to grain's own CoM + self.grain_propellant_mass, + grain_dist_vec_func, # Distance from total prop CoM to grain CoM + ) + self.propellant_I_11_from_propellant_CM = ( + liquid_I_11_prop_com + grain_I_11_prop_com + ) + + # --- I_22 --- + liquid_I_22_prop_com = parallel_axis_theorem_I22( + self.liquid_propellant_I_22, # Inertia relative to liquid's own CoM + self.liquid_propellant_mass, + liquid_dist_vec_func, + ) + grain_I_22_prop_com = parallel_axis_theorem_I22( + self.grain_propellant_I_22, # Inertia relative to grain's own CoM + self.grain_propellant_mass, + grain_dist_vec_func, + ) + self.propellant_I_22_from_propellant_CM = ( + liquid_I_22_prop_com + grain_I_22_prop_com + ) + + # --- I_33 --- + liquid_I_33_prop_com = parallel_axis_theorem_I33( + self.liquid_propellant_I_33, # Inertia relative to liquid's own CoM + self.liquid_propellant_mass, + liquid_dist_vec_func, + ) + grain_I_33_prop_com = parallel_axis_theorem_I33( + self.grain_propellant_I_33, # Inertia relative to grain's own CoM + self.grain_propellant_mass, + grain_dist_vec_func, + ) + self.propellant_I_33_from_propellant_CM = ( + liquid_I_33_prop_com + grain_I_33_prop_com + ) + + # --- Products of Inertia (I_12, I_13, I_23) --- + # Assume components PoI are 0 relative to their own CoM due to axisymmetry + # PAT calculation will correctly handle the axisymmetry (result should be 0) + + # I_12 + liquid_I_12_prop_com = parallel_axis_theorem_I12( + Function(0), self.liquid_propellant_mass, liquid_dist_vec_func + ) + grain_I_12_prop_com = parallel_axis_theorem_I12( + Function(0), self.grain_propellant_mass, grain_dist_vec_func + ) + # Store intermediate result if needed by Motor.__init__ later, prefix with '_' if not part of public API + self._propellant_I_12_from_propellant_CM = ( + liquid_I_12_prop_com + grain_I_12_prop_com + ) + + # I_13 + liquid_I_13_prop_com = parallel_axis_theorem_I13( + Function(0), self.liquid_propellant_mass, liquid_dist_vec_func + ) + grain_I_13_prop_com = parallel_axis_theorem_I13( + Function(0), self.grain_propellant_mass, grain_dist_vec_func + ) + self._propellant_I_13_from_propellant_CM = ( + liquid_I_13_prop_com + grain_I_13_prop_com + ) + + # I_23 + liquid_I_23_prop_com = parallel_axis_theorem_I23( + Function(0), self.liquid_propellant_mass, liquid_dist_vec_func + ) + grain_I_23_prop_com = parallel_axis_theorem_I23( + Function(0), self.grain_propellant_mass, grain_dist_vec_func + ) + self._propellant_I_23_from_propellant_CM = ( + liquid_I_23_prop_com + grain_I_23_prop_com + ) + + # IMPORTANT: Call the parent __init__ AFTER calculating component inertias + # because the parent __init__ uses these calculated values. + super().__init__( + thrust_source=thrust_source, + dry_mass=self.dry_mass, # Use the corrected dry mass + dry_inertia=( + self.dry_I_11, + self.dry_I_22, + self.dry_I_33, + self.dry_I_12, + self.dry_I_13, + self.dry_I_23, + ), # Use corrected dry inertia + nozzle_radius=nozzle_radius, + burn_time=burn_time, + center_of_dry_mass_position=center_of_dry_mass_position, + nozzle_position=nozzle_position, + reshape_thrust_curve=reshape_thrust_curve, + interpolation_method=interpolation_method, + coordinate_system_orientation=coordinate_system_orientation, + reference_pressure=reference_pressure, ) + # The parent __init__ will now correctly use the calculated + # self.propellant_I_xx_from_propellant_CM values. - self.positioned_tanks = self.liquid.positioned_tanks - self.grain_number = grain_number - 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 - self.grain_separation = grain_separation - self.grains_center_of_mass_position = grains_center_of_mass_position - self.throat_radius = throat_radius - - # Initialize plots and prints object - self.prints = _HybridMotorPrints(self) + # Initialize plots object specific to HybridMotor self.plots = _HybridMotorPlots(self) - @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. @@ -469,145 +695,51 @@ def center_of_propellant_mass(self): ) return mass_balance / self.propellant_mass + @funcify_method("Time (s)", "Inertia I_11 (kg m²)") + @property @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 - ---------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ + # Returns the value calculated in __init__ + return self.propellant_I_11_from_propellant_CM - solid_mass = self.solid.propellant_mass - liquid_mass = self.liquid.propellant_mass - - cm = self.center_of_propellant_mass - solid_cm_to_cm = self.solid.center_of_propellant_mass - cm - liquid_cm_to_cm = self.liquid.center_of_propellant_mass - cm - - solid_prop_inertia = self.solid.propellant_I_11 - liquid_prop_inertia = self.liquid.propellant_I_11 - - I_11 = parallel_axis_theorem_from_com( - solid_prop_inertia, solid_mass, solid_cm_to_cm - ) + parallel_axis_theorem_from_com( - liquid_prop_inertia, liquid_mass, liquid_cm_to_cm - ) - - return I_11 - + @property @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. + """Inertia tensor 22 component of the propellant... (Identical to I_11)""" - References - ---------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor - """ - return self.propellant_I_11 + return self.propellant_I_22_from_propellant_CM + @property @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. + """Inertia tensor 33 component of the propellant...""" - 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 - ---------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor - """ - return self.solid.propellant_I_33 + self.liquid.propellant_I_33 + return self.propellant_I_33_from_propellant_CM + @property @funcify_method("Time (s)", "Inertia I_12 (kg m²)") def propellant_I_12(self): - """Inertia tensor 12 component of the propellant, the inertia is - relative to the e_1 and e_2 axes, centered at the instantaneous - propellant center of mass. + """Inertia tensor 12 component of the propellant...""" - Returns - ------- - Function - Propellant inertia tensor 12 component at time t. - - Notes - ----- - This is assumed to be zero due to axial symmetry of the motor. This - could be improved in the future to account for the fact that the - motor is not perfectly symmetric. - """ - return 0 + return self._propellant_I_12_from_propellant_CM + @property @funcify_method("Time (s)", "Inertia I_13 (kg m²)") def propellant_I_13(self): - """Inertia tensor 13 component of the propellant, the inertia is - relative to the e_1 and e_3 axes, centered at the instantaneous - propellant center of mass. + """Inertia tensor 13 component of the propellant...""" - Returns - ------- - Function - Propellant inertia tensor 13 component at time t. - - Notes - ----- - This is assumed to be zero due to axial symmetry of the motor. This - could be improved in the future to account for the fact that the - motor is not perfectly symmetric. - """ - return 0 + return self._propellant_I_13_from_propellant_CM + @property @funcify_method("Time (s)", "Inertia I_23 (kg m²)") def propellant_I_23(self): - """Inertia tensor 23 component of the propellant, the inertia is - relative to the e_2 and e_3 axes, centered at the instantaneous - propellant center of mass. + """Inertia tensor 23 component of the propellant...""" - Returns - ------- - Function - Propellant inertia tensor 23 component at time t. - - Notes - ----- - This is assumed to be zero due to axial symmetry of the motor. This - could be improved in the future to account for the fact that the - motor is not perfectly symmetric. - """ - return 0 + return self._propellant_I_23_from_propellant_CM def add_tank(self, tank, position): """Adds a tank to the motor. @@ -652,8 +784,8 @@ def draw(self, *, filename=None): """ self.plots.draw(filename=filename) - def to_dict(self, **kwargs): - data = super().to_dict(**kwargs) + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs) data.update( { "grain_number": self.grain_number, @@ -671,18 +803,13 @@ def to_dict(self, **kwargs): } ) - if kwargs.get("include_outputs", False): - burn_rate = self.solid.burn_rate - if kwargs.get("discretize", False): - burn_rate = burn_rate.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) + if include_outputs: data.update( { "grain_inner_radius": self.solid.grain_inner_radius, "grain_height": self.solid.grain_height, "burn_area": self.solid.burn_area, - "burn_rate": burn_rate, + "burn_rate": self.solid.burn_rate, } ) @@ -714,11 +841,12 @@ def from_dict(cls, data): grain_separation=data["grain_separation"], grains_center_of_mass_position=data["grains_center_of_mass_position"], nozzle_position=data["nozzle_position"], - throat_radius=data["throat_radius"], + tanks_mass=data["tanks_mass"], + oxidizer_initial_mass=data["oxidizer_initial_mass"], + oxidizer_mass_flow_rate_curve=data["oxidizer_mass_flow_rate_curve"], + oxidizer_density=data["oxidizer_density"], + oxidizer_tanks_geometries=data["oxidizer_tanks_geometries"], + oxidizer_tanks_positions=data["oxidizer_tanks_positions"], reference_pressure=data.get("reference_pressure"), ) - - for tank in data["positioned_tanks"]: - motor.add_tank(tank["tank"], tank["position"]) - return motor diff --git a/rocketpy/motors/liquid_motor.py b/rocketpy/motors/liquid_motor.py index 1cd76a52f..f4a951fdd 100644 --- a/rocketpy/motors/liquid_motor.py +++ b/rocketpy/motors/liquid_motor.py @@ -3,21 +3,17 @@ import numpy as np from rocketpy.mathutils.function import funcify_method, reset_funcified_methods -from rocketpy.tools import parallel_axis_theorem_from_com from ..plots.liquid_motor_plots import _LiquidMotorPlots -from ..prints.liquid_motor_prints import _LiquidMotorPrints from .motor import Motor class LiquidMotor(Motor): """Class to specify characteristics and useful operations for Liquid motors. This class inherits from the Motor class. - See Also -------- Motor - Attributes ---------- LiquidMotor.coordinate_system_orientation : str @@ -157,113 +153,43 @@ class LiquidMotor(Motor): It will allow to obtain the net thrust in the Flight class. """ + # pylint: disable=too-many-locals, too-many-statements, too-many-arguments def __init__( self, thrust_source, dry_mass, dry_inertia, nozzle_radius, + burn_time, center_of_dry_mass_position, nozzle_position=0, - burn_time=None, reshape_thrust_curve=False, interpolation_method="linear", coordinate_system_orientation="nozzle_to_combustion_chamber", reference_pressure=None, + tanks_mass=0, ): - """Initialize LiquidMotor class, process thrust curve and geometrical - parameters and store results. + # Initialize the list of tanks + self.positioned_tanks = [] + + # Correct the dry mass to include the mass of the tanks + dry_mass = dry_mass + tanks_mass + dry_inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia - 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 ` - dry_mass : int, float - Same as in Motor class. See the :class:`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. - nozzle_radius : int, float - Motor's nozzle outlet radius in meters. - 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 ` - nozzle_position : float - Motor's nozzle outlet position in meters, specified in the motor's - coordinate system. See - :doc:`Positions and Coordinate Systems ` for - more information. - 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. - 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". - reference_pressure : int, float, optional - Atmospheric pressure in Pa at which the thrust data was recorded. - """ super().__init__( thrust_source=thrust_source, + dry_mass=dry_mass, dry_inertia=dry_inertia, nozzle_radius=nozzle_radius, + burn_time=burn_time, center_of_dry_mass_position=center_of_dry_mass_position, - dry_mass=dry_mass, nozzle_position=nozzle_position, - burn_time=burn_time, reshape_thrust_curve=reshape_thrust_curve, interpolation_method=interpolation_method, coordinate_system_orientation=coordinate_system_orientation, reference_pressure=reference_pressure, ) - self.positioned_tanks = [] - - # Initialize plots and prints object - self.prints = _LiquidMotorPrints(self) self.plots = _LiquidMotorPlots(self) @funcify_method("Time (s)", "Exhaust Velocity (m/s)") @@ -374,36 +300,16 @@ def center_of_propellant_mass(self): @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. + """Inertia tensor 11 component of the total propellant, the inertia is + relative to the e_1 axis, centered at the instantaneous total propellant + center of mass. Recalculated here relative to the instantaneous CoM. 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 - ---------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor + Total propellant inertia tensor 11 component at time t relative to total propellant CoM. """ - I_11 = 0 - center_of_mass = self.center_of_propellant_mass - - for positioned_tank in self.positioned_tanks: - tank = positioned_tank.get("tank") - tank_position = positioned_tank.get("position") - distance = tank_position + tank.center_of_mass - center_of_mass - I_11 += parallel_axis_theorem_from_com( - tank.inertia, tank.fluid_mass, distance - ) - - return I_11 + return self.propellant_I_11_from_propellant_CM @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def propellant_I_22(self): @@ -497,8 +403,8 @@ def draw(self, *, filename=None): """ self.plots.draw(filename=filename) - def to_dict(self, **kwargs): - data = super().to_dict(**kwargs) + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs) data.update( { "positioned_tanks": [ diff --git a/rocketpy/motors/motor.py b/rocketpy/motors/motor.py index 733653712..0a130c6d0 100644 --- a/rocketpy/motors/motor.py +++ b/rocketpy/motors/motor.py @@ -12,9 +12,18 @@ import requests from ..mathutils.function import Function, funcify_method +from ..mathutils.inertia import ( + parallel_axis_theorem_I11, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I22, + parallel_axis_theorem_I23, + parallel_axis_theorem_I33, +) +from ..mathutils.vector_matrix import Matrix, Vector from ..plots.motor_plots import _MotorPlots from ..prints.motor_prints import _MotorPrints -from ..tools import parallel_axis_theorem_from_com, tuple_handler +from ..tools import tuple_handler # pylint: disable=too-many-public-methods # ThrustCurve API cache @@ -172,177 +181,191 @@ class Motor(ABC): def __init__( self, thrust_source, + dry_mass, dry_inertia, nozzle_radius, + burn_time, center_of_dry_mass_position, - dry_mass=None, nozzle_position=0, - burn_time=None, reshape_thrust_curve=False, interpolation_method="linear", coordinate_system_orientation="nozzle_to_combustion_chamber", reference_pressure=None, ): - """Initialize Motor class, process thrust curve and geometrical + """Initializes 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 ` - - dry_mass : int, float - Same as in Motor class. See the :class:`Motor ` docs - 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 ` - 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. - nozzle_radius : int, float, optional - Motor's nozzle outlet radius in meters. - 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. - nozzle_position : int, float, optional - Motor's nozzle outlet position in meters, in the motor's coordinate - system. See :doc:`Positions and Coordinate Systems ` - for details. Default is 0, in which case the origin of the - coordinate system is placed at the motor's nozzle outlet. - 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. Note that the Motor burn_time parameter must include the new - reshaped burn time. - 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". - reference_pressure : int, float, optional - Atmospheric pressure in Pa at which the thrust data was recorded. - - Returns - ------- - None + [... Docstring inchangée ...] """ + # Motor attributes + self.nozzle_radius = nozzle_radius + self.nozzle_position = nozzle_position + self.dry_mass = dry_mass + self.center_of_dry_mass_position = center_of_dry_mass_position + dry_inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia + self.dry_I_11 = dry_inertia[0] + self.dry_I_22 = dry_inertia[1] + self.dry_I_33 = dry_inertia[2] + self.dry_I_12 = dry_inertia[3] + self.dry_I_13 = dry_inertia[4] + self.dry_I_23 = dry_inertia[5] + # Define coordinate system orientation - self.coordinate_system_orientation = coordinate_system_orientation if coordinate_system_orientation == "nozzle_to_combustion_chamber": self._csys = 1 elif coordinate_system_orientation == "combustion_chamber_to_nozzle": self._csys = -1 - else: # pragma: no cover - raise ValueError( - "Invalid coordinate system orientation. Options are " - "'nozzle_to_combustion_chamber' and 'combustion_chamber_to_nozzle'." + else: + raise TypeError( + "Invalid coordinate system orientation. Please choose between " + + '"nozzle_to_combustion_chamber" and ' + + '"combustion_chamber_to_nozzle".' ) + self.coordinate_system_orientation = coordinate_system_orientation - # Motor parameters - self.interpolate = interpolation_method - self.nozzle_position = nozzle_position - self.nozzle_radius = nozzle_radius - self.nozzle_area = np.pi * nozzle_radius**2 - self.center_of_dry_mass_position = center_of_dry_mass_position - self.reference_pressure = reference_pressure - - # Inertia tensor setup - inertia = (*dry_inertia, 0, 0, 0) if len(dry_inertia) == 3 else dry_inertia - self.dry_I_11 = inertia[0] - self.dry_I_22 = inertia[1] - self.dry_I_33 = inertia[2] - self.dry_I_12 = inertia[3] - self.dry_I_13 = inertia[4] - self.dry_I_23 = inertia[5] - - # Handle .eng or .rse file inputs - self.description_eng_file = None - self.rse_motor_data = None - if isinstance(thrust_source, str): - if ( - path.exists(thrust_source) - and path.splitext(path.basename(thrust_source))[1] == ".eng" - ): - _, self.description_eng_file, points = Motor.import_eng(thrust_source) - thrust_source = points - elif ( - path.exists(thrust_source) - and path.splitext(path.basename(thrust_source))[1] == ".rse" - ): - self.rse_motor_data, points = Motor.import_rse(thrust_source) - thrust_source = points - # Evaluate raw thrust source - self.thrust_source = thrust_source + self.reference_pressure = ( + reference_pressure # <-- CORRECTION 1 : Stockage de l'argument + ) + + # Thrust curve definition and processing self.thrust = Function( - thrust_source, "Time (s)", "Thrust (N)", self.interpolate, "zero" + thrust_source, + "Time (s)", + "Thrust (N)", + interpolation_method, + extrapolation="zero", ) - # Handle dry_mass input - self.dry_mass = dry_mass + self.interpolate = interpolation_method + # Burn time definition + self.burn_time = tuple_handler(burn_time) + self.burn_start_time = self.burn_time[0] + self.burn_out_time = self.burn_time[1] - # Handle burn_time input - self.burn_time = burn_time + # Reshape thrust curve if needed + reshape_thrust_curve_input = reshape_thrust_curve + if reshape_thrust_curve_input: + try: + new_burn_time, desired_impulse = reshape_thrust_curve_input + except (TypeError, ValueError) as exc: + raise TypeError( + "reshape_thrust_curve must be an iterable with two elements:" + " (new_burn_time, total_impulse)." + ) from exc + + self.thrust = self.reshape_thrust_curve( + self.thrust, new_burn_time, desired_impulse + ) - if callable(self.thrust.source): - self.thrust.set_discrete(*self.burn_time, 50, self.interpolate, "zero") + # Basic calculations and attributes + self.burn_duration = self.burn_out_time - self.burn_start_time + self.total_impulse = self.thrust.integral( + self.burn_start_time, self.burn_out_time + ) - # Reshape thrust_source if needed - if reshape_thrust_curve: - # Overwrites burn_time and thrust - self.thrust = Motor.reshape_thrust_curve(self.thrust, *reshape_thrust_curve) - self.burn_time = (self.thrust.x_array[0], self.thrust.x_array[-1]) + # Calculate max_thrust and max_thrust_time + try: + self.max_thrust = self.thrust.max + # Find time of max thrust by evaluating thrust at multiple points + if hasattr(self.thrust, "x_array"): + max_thrust_index = np.argmax(self.thrust.y_array) + self.max_thrust_time = self.thrust.x_array[max_thrust_index] + else: + # For lambda functions, sample over burn time + time_samples = np.linspace( + self.burn_start_time, self.burn_out_time, 1000 + ) + thrust_samples = [self.thrust(t) for t in time_samples] + max_thrust_index = np.argmax(thrust_samples) + self.max_thrust = thrust_samples[max_thrust_index] + self.max_thrust_time = time_samples[max_thrust_index] + except AttributeError: + # If thrust is lambda-based, sample to find max + time_samples = np.linspace(self.burn_start_time, self.burn_out_time, 1000) + thrust_samples = [self.thrust(t) for t in time_samples] + max_thrust_index = np.argmax(thrust_samples) + self.max_thrust = thrust_samples[max_thrust_index] + self.max_thrust_time = time_samples[max_thrust_index] - # Post process thrust - self.thrust = Motor.clip_thrust(self.thrust, self.burn_time) + self.average_thrust = self.total_impulse / self.burn_duration + self.reshape_thrust_curve_request = reshape_thrust_curve_input + + # Abstract methods - must be implemented by subclasses + self._propellant_initial_mass = 0 + self.exhaust_velocity = Function(0) + self.propellant_mass = Function(0) + self.total_mass = Function(0) + self.propellant_I_11_from_propellant_CM = Function(0) + self.propellant_I_22_from_propellant_CM = Function(0) + self.propellant_I_33_from_propellant_CM = Function(0) + self.center_of_propellant_mass = Function(0) + self.center_of_mass = Function(0) + + # --- DÉBUT CORRECTION 3 : Utilisation des nouvelles fonctions PAT --- + + # Inertia tensor for propellant referenced to the MOTOR's origin + # Note: distance_vec_3d is the vector from the motor origin to the propellant CoM + propellant_com_func = self.center_of_propellant_mass + + # Create a Function that returns the distance vector [0, 0, center_of_propellant_mass(t)] + propellant_com_vector_func = Function( + lambda t: Vector([0, 0, propellant_com_func(t)]), + inputs="t", + outputs="Vector (m)", + ) - # Auxiliary quantities - self.burn_start_time = self.burn_time[0] - self.burn_out_time = self.burn_time[1] - self.burn_duration = self.burn_time[1] - self.burn_time[0] + # Use the new specific PAT functions + self.propellant_I_11 = parallel_axis_theorem_I11( + self.propellant_I_11_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_11.set_outputs("Propellant I_11 (kg*m^2)") - # Compute thrust metrics - self.max_thrust = np.amax(self.thrust.y_array) - max_thrust_index = np.argmax(self.thrust.y_array) - self.max_thrust_time = self.thrust.source[max_thrust_index, 0] - self.average_thrust = self.total_impulse / self.burn_duration + self.propellant_I_22 = parallel_axis_theorem_I22( + self.propellant_I_22_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_22.set_outputs("Propellant I_22 (kg*m^2)") + + self.propellant_I_33 = parallel_axis_theorem_I33( + self.propellant_I_33_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_33.set_outputs("Propellant I_33 (kg*m^2)") + + self.propellant_I_12 = parallel_axis_theorem_I12( + Function(0), self.propellant_mass, propellant_com_vector_func + ) + self.propellant_I_12.set_outputs("Propellant I_12 (kg*m^2)") + + self.propellant_I_13 = parallel_axis_theorem_I13( + Function(0), self.propellant_mass, propellant_com_vector_func + ) + self.propellant_I_13.set_outputs("Propellant I_13 (kg*m^2)") + + self.propellant_I_23 = parallel_axis_theorem_I23( + Function(0), self.propellant_mass, propellant_com_vector_func + ) + self.propellant_I_23.set_outputs("Propellant I_23 (kg*m^2)") + + # Calculate total motor inertia relative to motor's origin + self.I_11 = Function(lambda t: self.dry_I_11) + self.propellant_I_11 + self.I_11.set_outputs("Motor I_11 (kg*m^2)") + self.I_22 = Function(lambda t: self.dry_I_22) + self.propellant_I_22 + self.I_22.set_outputs("Motor I_22 (kg*m^2)") + self.I_33 = Function(lambda t: self.dry_I_33) + self.propellant_I_33 + self.I_33.set_outputs("Motor I_33 (kg*m^2)") + + # Calculate total products of inertia relative to motor's origin + self.I_12 = Function(lambda t: self.dry_I_12) + self.propellant_I_12 + self.I_12.set_outputs("Motor I_12 (kg*m^2)") + self.I_13 = Function(lambda t: self.dry_I_13) + self.propellant_I_13 + self.I_13.set_outputs("Motor I_13 (kg*m^2)") + self.I_23 = Function(lambda t: self.dry_I_23) + self.propellant_I_23 + self.I_23.set_outputs("Motor I_23 (kg*m^2)") # Initialize plots and prints object self.prints = _MotorPrints(self) @@ -582,36 +605,30 @@ def center_of_propellant_mass(self): @funcify_method("Time (s)", "Inertia I_11 (kg m²)") def I_11(self): - """Inertia tensor 11 component, which corresponds to the inertia - relative to the e_1 axis, centered at the instantaneous 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. Also, due to symmetry, I_11 = I_22. - - See Also - -------- - https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor + """Builds a Function for the total I_11 inertia component relative + to the instantaneous center of mass of the motor. """ + # Inertias relative to the motor origin (calculated in __init__) + prop_I_11_origin = self.propellant_I_11 + dry_I_11_origin = Function(lambda t: self.dry_I_11) # Dry inertia is constant + + # Distance vectors FROM the instantaneous motor CoM TO the component CoMs + prop_com_to_inst_com = self.center_of_propellant_mass - self.center_of_mass + dry_com_to_inst_com = ( + Function(lambda t: Vector([0, 0, self.center_of_dry_mass_position])) + - self.center_of_mass + ) - prop_I_11 = self.propellant_I_11 - dry_I_11 = self.dry_I_11 - - prop_to_cm = self.center_of_propellant_mass - self.center_of_mass - dry_to_cm = self.center_of_dry_mass_position - self.center_of_mass - - prop_I_11 = parallel_axis_theorem_from_com( - prop_I_11, self.propellant_mass, prop_to_cm + # Apply PAT relative to the instantaneous motor CoM + # Note: We need the negative of the distance vector for PAT formula (origin to point) + prop_I_11_cm = parallel_axis_theorem_I11( + prop_I_11_origin, self.propellant_mass, -prop_com_to_inst_com + ) + dry_I_11_cm = parallel_axis_theorem_I11( + dry_I_11_origin, self.dry_mass, -dry_com_to_inst_com ) - dry_I_11 = parallel_axis_theorem_from_com(dry_I_11, self.dry_mass, dry_to_cm) - return prop_I_11 + dry_I_11 + return prop_I_11_cm + dry_I_11_cm @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def I_22(self): @@ -866,7 +883,7 @@ def propellant_I_13(self): Returns ------- Function - Propellant inertia tensor 13 component at time t. + Propellant inertia tensor 13 components at time t. Notes ----- @@ -1238,9 +1255,16 @@ def get_attr_value(obj, attr_name, multiplier=1): # Write last line file.write(f"{self.thrust.source[-1, 0]:.4f} {0:.3f}\n") - def to_dict(self, **kwargs): + def to_dict(self, include_outputs=False): + thrust_source = self.thrust_source + + if isinstance(thrust_source, str): + thrust_source = self.thrust.source + elif callable(thrust_source) and not isinstance(thrust_source, Function): + thrust_source = Function(thrust_source) + data = { - "thrust_source": self.thrust, + "thrust_source": thrust_source, "dry_I_11": self.dry_I_11, "dry_I_22": self.dry_I_22, "dry_I_33": self.dry_I_33, @@ -1258,94 +1282,31 @@ def to_dict(self, **kwargs): "reference_pressure": self.reference_pressure, } - if kwargs.get("include_outputs", False): - total_mass = self.total_mass - propellant_mass = self.propellant_mass - mass_flow_rate = self.total_mass_flow_rate - center_of_mass = self.center_of_mass - center_of_propellant_mass = self.center_of_propellant_mass - exhaust_velocity = self.exhaust_velocity - I_11 = self.I_11 - I_22 = self.I_22 - I_33 = self.I_33 - I_12 = self.I_12 - I_13 = self.I_13 - I_23 = self.I_23 - propellant_I_11 = self.propellant_I_11 - propellant_I_22 = self.propellant_I_22 - propellant_I_33 = self.propellant_I_33 - propellant_I_12 = self.propellant_I_12 - propellant_I_13 = self.propellant_I_13 - propellant_I_23 = self.propellant_I_23 - if kwargs.get("discretize", False): - total_mass = total_mass.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_mass = propellant_mass.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - mass_flow_rate = mass_flow_rate.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - center_of_mass = center_of_mass.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - center_of_propellant_mass = ( - center_of_propellant_mass.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - ) - exhaust_velocity = exhaust_velocity.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - I_11 = I_11.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_22 = I_22.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_33 = I_33.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_12 = I_12.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_13 = I_13.set_discrete_based_on_model(self.thrust, mutate_self=False) - I_23 = I_23.set_discrete_based_on_model(self.thrust, mutate_self=False) - propellant_I_11 = propellant_I_11.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_22 = propellant_I_22.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_33 = propellant_I_33.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_12 = propellant_I_12.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_13 = propellant_I_13.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) - propellant_I_23 = propellant_I_23.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) + if include_outputs: data.update( { "vacuum_thrust": self.vacuum_thrust, - "total_mass": total_mass, - "propellant_mass": propellant_mass, - "mass_flow_rate": mass_flow_rate, - "center_of_mass": center_of_mass, - "center_of_propellant_mass": center_of_propellant_mass, + "total_mass": self.total_mass, + "propellant_mass": self.propellant_mass, + "mass_flow_rate": self.mass_flow_rate, + "center_of_mass": self.center_of_mass, + "center_of_propellant_mass": self.center_of_propellant_mass, "total_impulse": self.total_impulse, - "exhaust_velocity": exhaust_velocity, + "exhaust_velocity": self.exhaust_velocity, "propellant_initial_mass": self.propellant_initial_mass, "structural_mass_ratio": self.structural_mass_ratio, - "I_11": I_11, - "I_22": I_22, - "I_33": I_33, - "I_12": I_12, - "I_13": I_13, - "I_23": I_23, - "propellant_I_11": propellant_I_11, - "propellant_I_22": propellant_I_22, - "propellant_I_33": propellant_I_33, - "propellant_I_12": propellant_I_12, - "propellant_I_13": propellant_I_13, - "propellant_I_23": propellant_I_23, + "I_11": self.I_11, + "I_22": self.I_22, + "I_33": self.I_33, + "I_12": self.I_12, + "I_13": self.I_13, + "I_23": self.I_23, + "propellant_I_11": self.propellant_I_11, + "propellant_I_22": self.propellant_I_22, + "propellant_I_33": self.propellant_I_33, + "propellant_I_12": self.propellant_I_12, + "propellant_I_13": self.propellant_I_13, + "propellant_I_23": self.propellant_I_23, } ) @@ -1573,7 +1534,7 @@ def center_of_propellant_mass(self): Function Function representing the center of mass of the motor. """ - return Function(self.chamber_position).set_discrete_based_on_model(self.thrust) + return self.chamber_position @funcify_method("Time (s)", "Inertia I_11 (kg m²)") def propellant_I_11(self): @@ -1595,11 +1556,11 @@ def propellant_I_11(self): ---------- https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - return Function( + return ( self.propellant_mass * (3 * self.chamber_radius**2 + self.chamber_height**2) / 12 - ).set_discrete_based_on_model(self.thrust) + ) @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def propellant_I_22(self): @@ -1643,21 +1604,19 @@ def propellant_I_33(self): ---------- https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - return Function( - self.propellant_mass * self.chamber_radius**2 / 2 - ).set_discrete_based_on_model(self.thrust) + return self.propellant_mass * self.chamber_radius**2 / 2 @funcify_method("Time (s)", "Inertia I_12 (kg m²)") def propellant_I_12(self): - return Function(0).set_discrete_based_on_model(self.thrust) + return Function(0) @funcify_method("Time (s)", "Inertia I_13 (kg m²)") def propellant_I_13(self): - return Function(0).set_discrete_based_on_model(self.thrust) + return Function(0) @funcify_method("Time (s)", "Inertia I_23 (kg m²)") def propellant_I_23(self): - return Function(0).set_discrete_based_on_model(self.thrust) + return Function(0) @staticmethod def load_from_eng_file( @@ -2076,8 +2035,8 @@ def all_info(self): self.prints.all() self.plots.all() - def to_dict(self, **kwargs): - data = super().to_dict(**kwargs) + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs) data.update( { "chamber_radius": self.chamber_radius, diff --git a/rocketpy/motors/solid_motor.py b/rocketpy/motors/solid_motor.py index f5e89c2f8..d13612841 100644 --- a/rocketpy/motors/solid_motor.py +++ b/rocketpy/motors/solid_motor.py @@ -4,6 +4,15 @@ from scipy import integrate from ..mathutils.function import Function, funcify_method, reset_funcified_methods +from ..mathutils.inertia import ( + parallel_axis_theorem_I11, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I22, + parallel_axis_theorem_I23, + parallel_axis_theorem_I33, +) +from ..mathutils.vector_matrix import Vector from ..plots.solid_motor_plots import _SolidMotorPlots from ..prints.solid_motor_prints import _SolidMotorPrints from .motor import Motor @@ -198,7 +207,6 @@ class SolidMotor(Motor): Grain length remains constant throughout the burn. Default is False. """ - # pylint: disable=too-many-arguments def __init__( self, thrust_source, @@ -331,6 +339,22 @@ class Function. Thrust units are Newtons. # Store before calling super().__init__() since it calls evaluate_geometry() self.only_radial_burn = only_radial_burn + self.rse_motor_data = None + self.description_eng_file = None + + if isinstance(thrust_source, str): + if thrust_source.endswith(".eng"): + comments, description, thrust_source_data = Motor.import_eng( + thrust_source + ) + self.description_eng_file = description + self.comments_eng_file = comments + thrust_source = thrust_source_data + elif thrust_source.endswith(".rse"): + rse_data, thrust_source_data = Motor.import_rse(thrust_source) + self.rse_motor_data = rse_data + thrust_source = thrust_source_data + super().__init__( thrust_source=thrust_source, dry_inertia=dry_inertia, @@ -344,11 +368,10 @@ class Function. Thrust units are Newtons. coordinate_system_orientation=coordinate_system_orientation, reference_pressure=reference_pressure, ) - # 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 @@ -357,7 +380,6 @@ class Function. Thrust units are Newtons. 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 @@ -367,9 +389,60 @@ class Function. Thrust units are Newtons. self.evaluate_geometry() - # Initialize plots and prints object self.prints = _SolidMotorPrints(self) self.plots = _SolidMotorPlots(self) + self.propellant_I_11_from_propellant_CM = self.propellant_I_11 + self.propellant_I_22_from_propellant_CM = self.propellant_I_22 + self.propellant_I_33_from_propellant_CM = self.propellant_I_33 + self.propellant_I_12_from_propellant_CM = self.propellant_I_12 + self.propellant_I_13_from_propellant_CM = self.propellant_I_13 + self.propellant_I_23_from_propellant_CM = self.propellant_I_23 + propellant_com_func = self.center_of_propellant_mass + + propellant_com_vector_func = Function( + lambda t: Vector([0, 0, propellant_com_func(t)]), + inputs="t", + outputs="Vector (m)", + ) + + # Utiliser les nouvelles fonctions PAT + self.propellant_I_11 = parallel_axis_theorem_I11( + self.propellant_I_11_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_22 = parallel_axis_theorem_I22( + self.propellant_I_22_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_33 = parallel_axis_theorem_I33( + self.propellant_I_33_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_12 = parallel_axis_theorem_I12( + self.propellant_I_12_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_13 = parallel_axis_theorem_I13( + self.propellant_I_13_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + self.propellant_I_23 = parallel_axis_theorem_I23( + self.propellant_I_23_from_propellant_CM, + self.propellant_mass, + propellant_com_vector_func, + ) + + self.I_11 = Function(lambda t: self.dry_I_11) + self.propellant_I_11 + self.I_22 = Function(lambda t: self.dry_I_22) + self.propellant_I_22 + self.I_33 = Function(lambda t: self.dry_I_33) + self.propellant_I_33 + self.I_12 = Function(lambda t: self.dry_I_12) + self.propellant_I_12 + self.I_13 = Function(lambda t: self.dry_I_13) + self.propellant_I_13 + self.I_23 = Function(lambda t: self.dry_I_23) + self.propellant_I_23 @funcify_method("Time (s)", "Mass (kg)") def propellant_mass(self): @@ -380,7 +453,8 @@ def propellant_mass(self): Function Mass of the motor, in kg. """ - return self.grain_volume * self.grain_density * self.grain_number + prop_mass = self.grain_volume * self.grain_density * self.grain_number + return prop_mass.set_discrete_based_on_model(self.grain_height) @funcify_method("Time (s)", "Grain volume (m³)") def grain_volume(self): @@ -736,7 +810,7 @@ def propellant_I_11(self): # Calculate inertia for all grains I_11 = grain_number * grain_inertia11 + grain_mass * np.sum(d**2) - return I_11 + return I_11.set_discrete_based_on_model(self.propellant_mass) @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def propellant_I_22(self): @@ -785,19 +859,22 @@ def propellant_I_33(self): * self.propellant_mass * (self.grain_outer_radius**2 + self.grain_inner_radius**2) ) - return I_33 + return I_33.set_discrete_based_on_model(self.propellant_mass) @funcify_method("Time (s)", "Inertia I_12 (kg m²)") def propellant_I_12(self): - return 0 + zero = Function(0) + return zero.set_discrete_based_on_model(self.propellant_mass) @funcify_method("Time (s)", "Inertia I_13 (kg m²)") def propellant_I_13(self): - return 0 + zero = Function(0) + return zero.set_discrete_based_on_model(self.propellant_mass) @funcify_method("Time (s)", "Inertia I_23 (kg m²)") def propellant_I_23(self): - return 0 + zero = Function(0) + return zero.set_discrete_based_on_model(self.propellant_mass) def draw(self, *, filename=None): """Draw a representation of the SolidMotor. @@ -816,8 +893,8 @@ def draw(self, *, filename=None): """ self.plots.draw(filename=filename) - def to_dict(self, **kwargs): - data = super().to_dict(**kwargs) + def to_dict(self, include_outputs=False): + data = super().to_dict(include_outputs) data.update( { "nozzle_radius": self.nozzle_radius, @@ -833,18 +910,13 @@ def to_dict(self, **kwargs): } ) - if kwargs.get("include_outputs", False): - burn_rate = self.burn_rate - if kwargs.get("discretize", False): - burn_rate = burn_rate.set_discrete_based_on_model( - self.thrust, mutate_self=False - ) + if include_outputs: data.update( { "grain_inner_radius": self.grain_inner_radius, "grain_height": self.grain_height, "burn_area": self.burn_area, - "burn_rate": burn_rate, + "burn_rate": self.burn_rate, "Kn": self.Kn, } ) diff --git a/rocketpy/plots/plot_helpers.py b/rocketpy/plots/plot_helpers.py index f795d46e9..0aeaeaa0a 100644 --- a/rocketpy/plots/plot_helpers.py +++ b/rocketpy/plots/plot_helpers.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt from matplotlib.figure import Figure -from ..tools import get_matplotlib_supported_file_endings +from rocketpy.tools import get_matplotlib_supported_file_endings SAVEFIG_DPI = 300 diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py index 8eaaded16..587848cd9 100644 --- a/rocketpy/plots/rocket_plots.py +++ b/rocketpy/plots/rocket_plots.py @@ -1,7 +1,13 @@ import matplotlib.pyplot as plt import numpy as np -from rocketpy.motors import EmptyMotor, HybridMotor, LiquidMotor, SolidMotor +from rocketpy.motors import ( + ClusterMotor, + EmptyMotor, + HybridMotor, + LiquidMotor, + SolidMotor, +) from rocketpy.rocket.aero_surface import Fins, NoseCone, Tail from rocketpy.rocket.aero_surface.generic_surface import GenericSurface @@ -201,7 +207,7 @@ def draw(self, vis_args=None, plane="xz", *, filename=None): drawn_surfaces = self._draw_aerodynamic_surfaces(ax, vis_args, plane) last_radius, last_x = self._draw_tubes(ax, drawn_surfaces, vis_args) - self._draw_motor(last_radius, last_x, ax, vis_args) + self._draw_motor(last_radius, last_x, ax, vis_args, plane) self._draw_rail_buttons(ax, vis_args) self._draw_center_of_mass_and_pressure(ax) self._draw_sensors(ax, self.rocket.sensors, plane) @@ -224,15 +230,7 @@ def __validate_aerodynamic_surfaces(self): def _draw_aerodynamic_surfaces(self, ax, vis_args, plane): """Draws the aerodynamic surfaces and saves the position of the points of interest for the tubes.""" - # List of drawn surfaces with the position of points of interest - # and the radius of the rocket at that point drawn_surfaces = [] - # Idea is to get the shape of each aerodynamic surface in their own - # coordinate system and then plot them in the rocket coordinate system - # using the position of each surface - # For the tubes, the surfaces need to be checked in order to check for - # diameter changes. The final point of the last surface is the final - # point of the last tube for surface, position in self.rocket.aerodynamic_surfaces: if isinstance(surface, NoseCone): @@ -264,14 +262,14 @@ def _draw_nose_cone(self, ax, surface, position, drawn_surfaces, vis_args): color=vis_args["nose"], linewidth=vis_args["line_width"], ) - # close the nosecone + ax.plot( [x_nosecone[-1], x_nosecone[-1]], [y_nosecone[-1], -y_nosecone[-1]], color=vis_args["nose"], linewidth=vis_args["line_width"], ) - # Add the nosecone to the list of drawn surfaces + drawn_surfaces.append( (surface, x_nosecone[-1], surface.rocket_radius, x_nosecone[-1]) ) @@ -287,7 +285,7 @@ def _draw_tail(self, ax, surface, position, drawn_surfaces, vis_args): ax.plot( x_tail, -y_tail, color=vis_args["tail"], linewidth=vis_args["line_width"] ) - # close above and below the tail + ax.plot( [x_tail[-1], x_tail[-1]], [y_tail[-1], -y_tail[-1]], @@ -300,7 +298,7 @@ def _draw_tail(self, ax, surface, position, drawn_surfaces, vis_args): color=vis_args["tail"], linewidth=vis_args["line_width"], ) - # Add the tail to the list of drawn surfaces + drawn_surfaces.append((surface, position, surface.bottom_radius, x_tail[-1])) def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args): @@ -312,16 +310,12 @@ def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args): rotation_angles = [2 * np.pi * i / num_fins for i in range(num_fins)] for angle in rotation_angles: - # Create a rotation matrix for the current angle around the x-axis rotation_matrix = np.array([[1, 0], [0, np.cos(angle)]]) - # Apply the rotation to the original fin points rotated_points_2d = np.dot(rotation_matrix, np.vstack((x_fin, y_fin))) - # Extract x and y coordinates of the rotated points x_rotated, y_rotated = rotated_points_2d - # Project points above the XY plane back into the XY plane (set z-coordinate to 0) x_rotated = np.where( rotated_points_2d[1] > 0, rotated_points_2d[0], x_rotated ) @@ -343,20 +337,18 @@ def _draw_generic_surface( surface, position, drawn_surfaces, - vis_args, # pylint: disable=unused-argument + vis_args, plane, ): """Draws the generic surface and saves the position of the points of interest for the tubes.""" if plane == "xz": - # z position of the sensor is the x position in the plot x_pos = position[2] - # x position of the surface is the y position in the plot + y_pos = position[0] elif plane == "yz": - # z position of the surface is the x position in the plot x_pos = position[2] - # y position of the surface is the y position in the plot + y_pos = position[1] else: # pragma: no cover raise ValueError("Plane must be 'xz' or 'yz'.") @@ -375,22 +367,16 @@ def _draw_tubes(self, ax, drawn_surfaces, vis_args): radius = 0 last_x = 0 for i, d_surface in enumerate(drawn_surfaces): - # Draw the tubes, from the end of the first surface to the beginning - # of the next surface, with the radius of the rocket at that point surface, position, radius, last_x = d_surface if i == len(drawn_surfaces) - 1: - # If the last surface is a tail, do nothing if isinstance(surface, Tail): continue - # Else goes to the end of the surface + x_tube = [position, last_x] y_tube = [radius, radius] y_tube_negated = [-radius, -radius] else: - # If it is not the last surface, the tube goes to the beginning - # of the next surface - # [next_surface, next_position, next_radius, next_last_x] next_position = drawn_surfaces[i + 1][1] x_tube = [last_x, next_position] y_tube = [radius, radius] @@ -410,45 +396,122 @@ def _draw_tubes(self, ax, drawn_surfaces, vis_args): ) return radius, last_x - def _draw_motor(self, last_radius, last_x, ax, vis_args): + def _draw_motor(self, last_radius, last_x, ax, vis_args, plane="xz"): """Draws the motor from motor patches""" - total_csys = self.rocket._csys * self.rocket.motor._csys - nozzle_position = ( - self.rocket.motor_position + self.rocket.motor.nozzle_position * total_csys - ) - # Get motor patches translated to the correct position - motor_patches = self._generate_motor_patches(total_csys, ax) + motor_patches = self._generate_motor_patches(ax, plane) - # Draw patches if not isinstance(self.rocket.motor, EmptyMotor): - # Add nozzle last so it is in front of the other patches - nozzle = self.rocket.motor.plots._generate_nozzle( - translate=(nozzle_position, 0), csys=self.rocket._csys - ) - motor_patches += [nozzle] + if isinstance(self.rocket.motor, ClusterMotor): + for patch in motor_patches: + ax.add_patch(patch) + + if self.rocket.motor.positions: + cluster_z_positions = [ + pos[2] for pos in self.rocket.motor.positions + ] + z_connector = ( + min(cluster_z_positions) + if self.rocket._csys == 1 + else max(cluster_z_positions) + ) + self._draw_nozzle_tube( + last_radius, last_x, z_connector, ax, vis_args + ) - outline = self.rocket.motor.plots._generate_motor_region( - list_of_patches=motor_patches - ) - # add outline first so it is behind the other patches - ax.add_patch(outline) - for patch in motor_patches: - ax.add_patch(patch) + else: + total_csys = self.rocket._csys * self.rocket.motor._csys + nozzle_position = ( + self.rocket.motor_position + + self.rocket.motor.nozzle_position * total_csys + ) + + nozzle = self.rocket.motor.plots._generate_nozzle( + translate=(nozzle_position, 0), csys=self.rocket._csys + ) + motor_patches += [nozzle] + + outline = self.rocket.motor.plots._generate_motor_region( + list_of_patches=motor_patches + ) - self._draw_nozzle_tube(last_radius, last_x, nozzle_position, ax, vis_args) + ax.add_patch(outline) + for patch in motor_patches: + ax.add_patch(patch) - def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-argument + self._draw_nozzle_tube( + last_radius, last_x, nozzle_position, ax, vis_args + ) + + def _generate_motor_patches(self, ax, plane="xz"): """Generates motor patches for drawing""" motor_patches = [] + total_csys = self.rocket._csys * self.rocket.motor._csys + + if isinstance(self.rocket.motor, ClusterMotor): + cluster = self.rocket.motor + all_sub_patches = [] + + for sub_motor, sub_pos in zip(cluster.motors, cluster.positions): + motor_longitudinal_pos = sub_pos[2] + + offset = 0 + if plane == "xz": + offset = sub_pos[0] + elif plane == "yz": + offset = sub_pos[1] + + sub_total_csys = self.rocket._csys * sub_motor._csys + + if isinstance(sub_motor, SolidMotor): + current_motor_patches = [] + + grains_cm_pos = ( + motor_longitudinal_pos + + sub_motor.grains_center_of_mass_position * sub_total_csys + ) + ax.scatter( + grains_cm_pos.real, + offset, + color="brown", + label=f"Grains CM ({sub_pos[0]:.2f}, {sub_pos[1]:.2f})", + s=8, + zorder=10, + ) + + chamber = sub_motor.plots._generate_combustion_chamber( + translate=(grains_cm_pos.real, offset), label=None + ) + grains = sub_motor.plots._generate_grains( + translate=(grains_cm_pos.real, offset) + ) + + nozzle_position = ( + motor_longitudinal_pos + + sub_motor.nozzle_position * sub_total_csys + ) + + nozzle = sub_motor.plots._generate_nozzle( + translate=(nozzle_position.real, offset) + ) + + current_motor_patches += [chamber, *grains, nozzle] + all_sub_patches.extend(current_motor_patches) + + if all_sub_patches: + outline = self.rocket.motor.motors[0].plots._generate_motor_region( + list_of_patches=all_sub_patches + ) + motor_patches.append(outline) # Mettre l'outline en premier + motor_patches.extend(all_sub_patches) - if isinstance(self.rocket.motor, SolidMotor): + elif isinstance(self.rocket.motor, SolidMotor): grains_cm_position = ( self.rocket.motor_position + self.rocket.motor.grains_center_of_mass_position * total_csys ) ax.scatter( - grains_cm_position, + grains_cm_position.real, 0, color="brown", label="Grains Center of Mass", @@ -457,10 +520,10 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg ) chamber = self.rocket.motor.plots._generate_combustion_chamber( - translate=(grains_cm_position, 0), label=None + translate=(grains_cm_position.real, 0), label=None ) grains = self.rocket.motor.plots._generate_grains( - translate=(grains_cm_position, 0) + translate=(grains_cm_position.real, 0) ) motor_patches += [chamber, *grains] @@ -471,7 +534,7 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg + self.rocket.motor.grains_center_of_mass_position * total_csys ) ax.scatter( - grains_cm_position, + grains_cm_position.real, 0, color="brown", label="Grains Center of Mass", @@ -483,16 +546,16 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg translate=(self.rocket.motor_position, 0), csys=total_csys ) chamber = self.rocket.motor.plots._generate_combustion_chamber( - translate=(grains_cm_position, 0), label=None + translate=(grains_cm_position.real, 0), label=None ) grains = self.rocket.motor.plots._generate_grains( - translate=(grains_cm_position, 0) + translate=(grains_cm_position.real, 0) ) motor_patches += [chamber, *grains] for tank, center in tanks_and_centers: ax.scatter( - center[0], - center[1], + center[0].real, + center[1].real, color="black", alpha=0.2, s=5, @@ -506,8 +569,8 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg ) for tank, center in tanks_and_centers: ax.scatter( - center[0], - center[1], + center[0].real, + center[1].real, color="black", alpha=0.2, s=4, @@ -519,8 +582,6 @@ def _generate_motor_patches(self, total_csys, ax): # pylint: disable=unused-arg def _draw_nozzle_tube(self, last_radius, last_x, nozzle_position, ax, vis_args): """Draws the tube from the last surface to the nozzle position.""" - # Check if nozzle is beyond the last surface, if so draw a tube - # to it, with the radius of the last surface if self.rocket._csys == 1: if nozzle_position < last_x: x_tube = [last_x, nozzle_position] @@ -539,7 +600,7 @@ def _draw_nozzle_tube(self, last_radius, last_x, nozzle_position, ax, vis_args): color=vis_args["body"], linewidth=vis_args["line_width"], ) - else: # if self.rocket._csys == -1: + else: if nozzle_position > last_x: x_tube = [last_x, nozzle_position] y_tube = [last_radius, last_radius] @@ -575,13 +636,16 @@ def _draw_rail_buttons(self, ax, vis_args): def _draw_center_of_mass_and_pressure(self, ax): """Draws the center of mass and center of pressure of the rocket.""" - # Draw center of mass and center of pressure - cm = self.rocket.center_of_mass(0) - ax.scatter(cm, 0, color="#1565c0", label="Center of Mass", s=10) - cp = self.rocket.cp_position(0) + cm_vector = self.rocket.center_of_mass(0) + + cm_z = float(cm_vector.z.real) + + ax.scatter(cm_z, 0, color="#1565c0", label="Center of Mass", s=10) + + cp_z = float(self.rocket.cp_position(0).real) ax.scatter( - cp, 0, label="Static Center of Pressure", color="red", s=10, zorder=10 + cp_z, 0, label="Static Center of Pressure", color="red", s=10, zorder=10 ) def _draw_sensors(self, ax, sensors, plane): @@ -593,23 +657,20 @@ def _draw_sensors(self, ax, sensors, plane): sensor = sensor_pos[0] pos = sensor_pos[1] if plane == "xz": - # z position of the sensor is the x position in the plot x_pos = pos[2] normal_x = sensor.normal_vector.z - # x position of the sensor is the y position in the plot + y_pos = pos[0] normal_y = sensor.normal_vector.x elif plane == "yz": - # z position of the sensor is the x position in the plot x_pos = pos[2] normal_x = sensor.normal_vector.z - # y position of the sensor is the y position in the plot + y_pos = pos[1] normal_y = sensor.normal_vector.y - else: # pragma: no cover + else: raise ValueError("Plane must be 'xz' or 'yz'.") - # line length is 2/5 of the rocket radius line_length = self.rocket.radius / 2.5 ax.plot( @@ -644,34 +705,30 @@ def all(self): None """ - # Rocket draw if len(self.rocket.aerodynamic_surfaces) > 0: print("\nRocket Draw") print("-" * 40) self.draw() - # Mass Plots print("\nMass Plots") print("-" * 40) self.total_mass() self.reduced_mass() - # Aerodynamics Plots print("\nAerodynamics Plots") print("-" * 40) # Drag Plots print("Drag Plots") - print("-" * 20) # Separator for Drag Plots + print("-" * 20) self.drag_curves() # Stability Plots print("\nStability Plots") - print("-" * 20) # Separator for Stability Plots + print("-" * 20) self.static_margin() self.stability_margin() - # Thrust-to-Weight Plot print("\nThrust-to-Weight Plot") print("-" * 40) self.thrust_to_weight() diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index c3b1f364d..ac3669d7d 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -1,12 +1,20 @@ import math import warnings -from typing import Iterable import numpy as np from rocketpy.control.controller import _Controller from rocketpy.mathutils.function import Function +from rocketpy.mathutils.inertia import ( + parallel_axis_theorem_I11, + parallel_axis_theorem_I12, + parallel_axis_theorem_I13, + parallel_axis_theorem_I22, + parallel_axis_theorem_I23, + parallel_axis_theorem_I33, +) from rocketpy.mathutils.vector_matrix import Matrix, Vector +from rocketpy.motors.cluster_motor import ClusterMotor from rocketpy.motors.empty_motor import EmptyMotor from rocketpy.plots.rocket_plots import _RocketPlots from rocketpy.prints.rocket_prints import _RocketPrints @@ -23,11 +31,6 @@ from rocketpy.rocket.aero_surface.generic_surface import GenericSurface from rocketpy.rocket.components import Components from rocketpy.rocket.parachute import Parachute -from rocketpy.tools import ( - deprecated, - find_obj_from_hash, - parallel_axis_theorem_from_com, -) # pylint: disable=too-many-instance-attributes, too-many-public-methods, too-many-instance-attributes @@ -365,24 +368,48 @@ def __init__( # pylint: disable=too-many-statements self.reduced_mass = None self.total_mass = None self.dry_mass = None + self.propellant_initial_mass = 0 # Ajout pour la cohérence + + # Initialize plots and prints object + self.prints = _RocketPrints(self) + self.plots = _RocketPlots(self) + self._full_rocket_update() + + def _full_rocket_update(self): + """ + Centralized method to update all rocket properties after adding or modifying a motor. + """ + if hasattr(self.motor, "propellant_initial_mass"): + self.propellant_initial_mass = self.motor.propellant_initial_mass + else: + self.propellant_initial_mass = ( + self.motor.propellant_mass.get_value_opt(0) + if self.motor.propellant_mass + else 0 + ) - # calculate dynamic inertial quantities self.evaluate_dry_mass() - self.evaluate_structural_mass_ratio() self.evaluate_total_mass() self.evaluate_center_of_dry_mass() + + if hasattr(self.motor, "nozzle_position"): + self.nozzle_position = self.motor.nozzle_position + else: + self.nozzle_position = self.motor_position + + self.evaluate_nozzle_to_cdm() self.evaluate_center_of_mass() + self.evaluate_dry_inertias() + self.evaluate_inertias() self.evaluate_reduced_mass() self.evaluate_thrust_to_weight() - - # Evaluate stability (even though no aerodynamic surfaces are present yet) self.evaluate_center_of_pressure() + self.evaluate_surfaces_cp_to_cdm() self.evaluate_stability_margin() self.evaluate_static_margin() - - # Initialize plots and prints object - self.prints = _RocketPrints(self) - self.plots = _RocketPlots(self) + self.evaluate_com_to_cdm_function() + self.evaluate_nozzle_gyration_tensor() + self.evaluate_structural_mass_ratio() def _check_missing_components(self): """Check if the rocket is missing any essential components and issue a warning. @@ -483,13 +510,23 @@ def evaluate_structural_mass_ratio(self): divided by total mass (Rocket + Motor + Propellant) (kg). """ try: - self.structural_mass_ratio = self.dry_mass / ( - self.dry_mass + self.motor.propellant_initial_mass + dry_mass = self.dry_mass if self.dry_mass is not None else self.mass + propellant_mass_initial = self.propellant_initial_mass + total_mass_initial = dry_mass + propellant_mass_initial + + if total_mass_initial <= 1e-9: + print( + "Warning: Total initial mass is near zero, setting structural mass ratio to 1.0" + ) + self.structural_mass_ratio = 1.0 + else: + self.structural_mass_ratio = dry_mass / total_mass_initial + except (AttributeError, TypeError) as e: + print( + f"Warning: Could not calculate structural mass ratio, setting to 1.0. Details: {e}" ) - except ZeroDivisionError as e: - raise ValueError( - "Total rocket mass (dry + propellant) cannot be zero" - ) from e + self.structural_mass_ratio = 1.0 + return self.structural_mass_ratio def evaluate_center_of_mass(self): @@ -504,10 +541,18 @@ def evaluate_center_of_mass(self): See :doc:`Positions and Coordinate Systems ` for more information. """ + + com_without_motor_vec = Function( + lambda t: Vector([0, 0, self.center_of_mass_without_motor]), + inputs="t", + outputs="Vector (m)", + ) + self.center_of_mass = ( - self.center_of_mass_without_motor * self.mass + com_without_motor_vec * self.mass + self.motor_center_of_mass_position * self.motor.total_mass ) / self.total_mass + self.center_of_mass.set_inputs("Time (s)") self.center_of_mass.set_outputs("Center of Mass Position (m)") self.center_of_mass.set_title( @@ -522,13 +567,19 @@ def evaluate_center_of_dry_mass(self): Returns ------- - self.center_of_dry_mass_position : int, float + self.center_of_dry_mass_position : Vector Rocket's center of dry mass position (with unloaded motor) """ + + com_without_motor_vec = Vector([0, 0, self.center_of_mass_without_motor]) + + # motor_center_of_dry_mass_position + motor_com_dry_vec = self.motor_center_of_dry_mass_position + self.center_of_dry_mass_position = ( - self.center_of_mass_without_motor * self.mass - + self.motor_center_of_dry_mass_position * self.motor.dry_mass + com_without_motor_vec * self.mass + motor_com_dry_vec * self.motor.dry_mass ) / self.dry_mass + return self.center_of_dry_mass_position def evaluate_reduced_mass(self): @@ -633,7 +684,7 @@ def __evaluate_single_surface_cp_to_cdm(self, surface, position): [ (position.x - self.cm_eccentricity_x) * self._csys - surface.cpx, (position.y - self.cm_eccentricity_y) - surface.cpy, - (position.z - self.center_of_dry_mass_position) * self._csys + (position.z - self.center_of_dry_mass_position.z) * self._csys - surface.cpz, ] ) @@ -642,21 +693,10 @@ def __evaluate_single_surface_cp_to_cdm(self, surface, position): def evaluate_stability_margin(self): """Calculates the stability margin of the rocket as a function of mach number and time. - - Returns - ------- - stability_margin : Function - Stability margin of the rocket, in calibers, as a function of mach - number and time. Stability margin is defined as the distance between - the center of pressure and the center of mass, divided by the - rocket's diameter. """ self.stability_margin.set_source( lambda mach, time: ( - ( - self.center_of_mass.get_value_opt(time) - - self.cp_position.get_value_opt(mach) - ) + (self.center_of_mass(time).z - self.cp_position.get_value_opt(mach)) / (2 * self.radius) ) * self._csys @@ -664,31 +704,21 @@ def evaluate_stability_margin(self): return self.stability_margin def evaluate_static_margin(self): - """Calculates the static margin of the rocket as a function of time. - - Returns - ------- - static_margin : Function - Static margin of the rocket, in calibers, as a function of time. - Static margin is defined as the distance between the center of - pressure and the center of mass, divided by the rocket's diameter. - """ - # Calculate static margin + """Calculates the static margin of the rocket as a function of time.""" self.static_margin.set_source( lambda time: ( - self.center_of_mass.get_value_opt(time) - - self.cp_position.get_value_opt(0) + self.center_of_mass(time).z - self.cp_position.get_value_opt(0) ) / (2 * self.radius) ) - # Change sign if coordinate system is upside down self.static_margin *= self._csys self.static_margin.set_inputs("Time (s)") self.static_margin.set_outputs("Static Margin (c)") self.static_margin.set_title("Static Margin") - self.static_margin.set_discrete( - lower=0, upper=self.motor.burn_out_time, samples=200 - ) + if hasattr(self.motor, "burn_out_time"): + self.static_margin.set_discrete( + lower=0, upper=self.motor.burn_out_time, samples=200 + ) return self.static_margin def evaluate_dry_inertias(self): @@ -696,76 +726,49 @@ def evaluate_dry_inertias(self): the rocket's center of dry mass. The inertias are saved and returned in units of kg*m². This does not consider propellant mass but does take into account the motor dry mass. - - Returns - ------- - self.dry_I_11 : float - Float value corresponding to rocket inertia tensor 11 - component, which corresponds to the inertia relative to the - e_1 axis, centered at the center of dry mass. - self.dry_I_22 : float - Float value corresponding to rocket inertia tensor 22 - component, which corresponds to the inertia relative to the - e_2 axis, centered at the center of dry mass. - self.dry_I_33 : float - Float value corresponding to rocket inertia tensor 33 - component, which corresponds to the inertia relative to the - e_3 axis, centered at the center of dry mass. - self.dry_I_12 : float - Float value corresponding to rocket inertia tensor 12 - component, which corresponds to the inertia relative to the - e_1 and e_2 axes, centered at the center of dry mass. - self.dry_I_13 : float - Float value corresponding to rocket inertia tensor 13 - component, which corresponds to the inertia relative to the - e_1 and e_3 axes, centered at the center of dry mass. - self.dry_I_23 : float - Float value corresponding to rocket inertia tensor 23 - component, which corresponds to the inertia relative to the - e_2 and e_3 axes, centered at the center of dry mass. - - Notes - ----- - #. The ``e_1`` and ``e_2`` directions are assumed to be the directions \ - perpendicular to the rocket axial direction. - #. The ``e_3`` direction is assumed to be the direction parallel to the \ - axis of symmetry of the rocket. - #. RocketPy follows the definition of the inertia tensor that includes \ - the minus sign for all products of inertia. - - See Also - -------- - `Inertia Tensor `_ """ # Get masses motor_dry_mass = self.motor.dry_mass mass = self.mass - # Compute axes distances (CDM: Center of Dry Mass) + center_of_mass_without_motor_vec = Vector( + [0, 0, self.center_of_mass_without_motor] + ) + + # Compute axes distances (CDM: Center of Dry Mass) using vector subtraction center_of_mass_without_motor_to_CDM = ( - self.center_of_mass_without_motor - self.center_of_dry_mass_position + center_of_mass_without_motor_vec - self.center_of_dry_mass_position ) motor_center_of_dry_mass_to_CDM = ( self.motor_center_of_dry_mass_position - self.center_of_dry_mass_position ) - # Compute dry inertias - self.dry_I_11 = parallel_axis_theorem_from_com( - self.I_11_without_motor, mass, center_of_mass_without_motor_to_CDM - ) + parallel_axis_theorem_from_com( - self.motor.dry_I_11, motor_dry_mass, motor_center_of_dry_mass_to_CDM - ) + d_rocket = center_of_mass_without_motor_to_CDM + d_motor = motor_center_of_dry_mass_to_CDM - self.dry_I_22 = parallel_axis_theorem_from_com( - self.I_22_without_motor, mass, center_of_mass_without_motor_to_CDM - ) + parallel_axis_theorem_from_com( - self.motor.dry_I_22, motor_dry_mass, motor_center_of_dry_mass_to_CDM - ) + self.dry_I_11 = parallel_axis_theorem_I11( + self.I_11_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I11(self.motor.dry_I_11, motor_dry_mass, d_motor) + + self.dry_I_22 = parallel_axis_theorem_I22( + self.I_22_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I22(self.motor.dry_I_22, motor_dry_mass, d_motor) + + self.dry_I_33 = parallel_axis_theorem_I33( + self.I_33_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I33(self.motor.dry_I_33, motor_dry_mass, d_motor) - self.dry_I_33 = self.I_33_without_motor + self.motor.dry_I_33 - self.dry_I_12 = self.I_12_without_motor + self.motor.dry_I_12 - self.dry_I_13 = self.I_13_without_motor + self.motor.dry_I_13 - self.dry_I_23 = self.I_23_without_motor + self.motor.dry_I_23 + self.dry_I_12 = parallel_axis_theorem_I12( + self.I_12_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I12(self.motor.dry_I_12, motor_dry_mass, d_motor) + + self.dry_I_13 = parallel_axis_theorem_I13( + self.I_13_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I13(self.motor.dry_I_13, motor_dry_mass, d_motor) + + self.dry_I_23 = parallel_axis_theorem_I23( + self.I_23_without_motor, mass, d_rocket + ) + parallel_axis_theorem_I23(self.motor.dry_I_23, motor_dry_mass, d_motor) return ( self.dry_I_11, @@ -780,58 +783,43 @@ def evaluate_inertias(self): """Calculates and returns the rocket's inertias relative to the rocket's center of dry mass. The inertias are saved and returned in units of kg*m². - - Returns - ------- - self.I_11 : float - Float value corresponding to rocket inertia tensor 11 - component, which corresponds to the inertia relative to the - e_1 axis, centered at the center of dry mass. - self.I_22 : float - Float value corresponding to rocket inertia tensor 22 - component, which corresponds to the inertia relative to the - e_2 axis, centered at the center of dry mass. - self.I_33 : float - Float value corresponding to rocket inertia tensor 33 - component, which corresponds to the inertia relative to the - e_3 axis, centered at the center of dry mass. - - Notes - ----- - #. The ``e_1`` and ``e_2`` directions are assumed to be the directions \ - perpendicular to the rocket axial direction. - #. The ``e_3`` direction is assumed to be the direction parallel to the \ - axis of symmetry of the rocket. - #. RocketPy follows the definition of the inertia tensor that includes \ - the minus sign for all products of inertia. - - See Also - -------- - `Inertia Tensor `_ """ - # Get masses - prop_mass = self.motor.propellant_mass # Propellant mass as a function of time - # Compute axes distances - CDM_to_CPM = ( - self.center_of_dry_mass_position - self.center_of_propellant_position + prop_mass = self.motor.propellant_mass + + cdm_position_func = Function( + lambda t: self.center_of_dry_mass_position, inputs="t", outputs="Vector (m)" ) - # Compute inertias - self.I_11 = self.dry_I_11 + parallel_axis_theorem_from_com( + CDM_to_CPM = cdm_position_func - self.center_of_propellant_position + + propellant_inertia_term_11 = parallel_axis_theorem_I11( self.motor.propellant_I_11, prop_mass, CDM_to_CPM ) - - self.I_22 = self.dry_I_22 + parallel_axis_theorem_from_com( + propellant_inertia_term_22 = parallel_axis_theorem_I22( self.motor.propellant_I_22, prop_mass, CDM_to_CPM ) + propellant_inertia_term_33 = parallel_axis_theorem_I33( + self.motor.propellant_I_33, prop_mass, CDM_to_CPM + ) + propellant_inertia_term_12 = parallel_axis_theorem_I12( + self.motor.propellant_I_12, prop_mass, CDM_to_CPM + ) + propellant_inertia_term_13 = parallel_axis_theorem_I13( + self.motor.propellant_I_13, prop_mass, CDM_to_CPM + ) + propellant_inertia_term_23 = parallel_axis_theorem_I23( + self.motor.propellant_I_23, prop_mass, CDM_to_CPM + ) + + self.I_11 = Function(lambda t: self.dry_I_11) + propellant_inertia_term_11 + self.I_22 = Function(lambda t: self.dry_I_22) + propellant_inertia_term_22 + self.I_33 = Function(lambda t: self.dry_I_33) + propellant_inertia_term_33 - self.I_33 = self.dry_I_33 + self.motor.propellant_I_33 - self.I_12 = self.dry_I_12 + self.motor.propellant_I_12 - self.I_13 = self.dry_I_13 + self.motor.propellant_I_13 - self.I_23 = self.dry_I_23 + self.motor.propellant_I_23 + self.I_12 = Function(lambda t: self.dry_I_12) + propellant_inertia_term_12 + self.I_13 = Function(lambda t: self.dry_I_13) + propellant_inertia_term_13 + self.I_23 = Function(lambda t: self.dry_I_23) + propellant_inertia_term_23 - # Return inertias return ( self.I_11, self.I_22, @@ -844,15 +832,9 @@ def evaluate_inertias(self): def evaluate_nozzle_to_cdm(self): """Evaluates the distance between the nozzle exit and the rocket's center of dry mass. - - Returns - ------- - self.nozzle_to_cdm : float - Distance between the nozzle exit and the rocket's center of dry - mass position, in meters. """ self.nozzle_to_cdm = ( - -(self.nozzle_position - self.center_of_dry_mass_position) * self._csys + -(self.nozzle_position - self.center_of_dry_mass_position.z) * self._csys ) return self.nozzle_to_cdm @@ -881,29 +863,21 @@ def evaluate_nozzle_gyration_tensor(self): def evaluate_com_to_cdm_function(self): """Evaluates the z-coordinate of the center of mass (COM) relative to the center of dry mass (CDM). + """ - Notes - ----- - 1. The `com_to_cdm_function` plus `center_of_mass` should be equal - to `center_of_dry_mass_position` at every time step. - 2. The `com_to_cdm_function` is a function of time and will usually - already be discretized. + cpm_z_func = Function( + lambda t: self.center_of_propellant_position(t).z, inputs="t" + ) + + cdm_z_scalar = self.center_of_dry_mass_position.z - Returns - ------- - self.com_to_cdm_function : Function - Function of time expressing the z-coordinate of the center of mass - relative to the center of dry mass. - """ self.com_to_cdm_function = ( -1 - * ( - (self.center_of_propellant_position - self.center_of_dry_mass_position) - * self._csys - ) + * ((cpm_z_func - cdm_z_scalar) * self._csys) * self.motor.propellant_mass / self.total_mass ) + self.com_to_cdm_function.set_inputs("Time (s)") self.com_to_cdm_function.set_outputs("Z Coordinate COM to CDM (m)") self.com_to_cdm_function.set_title("Z Coordinate COM to CDM") @@ -969,62 +943,77 @@ def get_inertia_tensor_derivative_at_time(self, t): ] ) - def add_motor(self, motor, position): # pylint: disable=too-many-statements - """Adds a motor to the rocket. + def add_motor(self, motor, position=0): + """ + Adds a motor or a motor cluster to the rocket. This method is universal + and handles standard motors, EmptyMotor, and ClusterMotor objects. + """ + if hasattr(self, "motor") and not isinstance(self.motor, EmptyMotor): + print("A motor is already attached. Overwriting previous motor.") - Parameters - ---------- - motor : Motor, SolidMotor, HybridMotor, LiquidMotor, GenericMotor - Motor to be added to the rocket. - position : int, float - Position, in meters, of the motor's coordinate system origin - relative to the user defined rocket coordinate system. + self.motor = motor + self.motor_position = position - See Also - -------- - :ref:`addsurface` + if isinstance(motor, EmptyMotor): + self.center_of_propellant_position = Function( + lambda t: Vector([0, 0, position]) + ) + self.motor_center_of_mass_position = Function( + lambda t: Vector([0, 0, position]) + ) + self.motor_center_of_dry_mass_position = Vector([0, 0, position]) + self.nozzle_position = position + self.total_mass_flow_rate = Function(0) + return - Returns - ------- - None + _ = self._csys * getattr(motor, "_csys", 1) + + is_cluster = hasattr(motor, "get_total_thrust_vector") + + if is_cluster: + self.center_of_propellant_position = motor.center_of_propellant_mass + self.motor_center_of_mass_position = motor.center_of_mass + self.motor_center_of_dry_mass_position = motor.center_of_dry_mass_position + else: + self.motor_center_of_mass_position = Function( + lambda t: Vector([0, 0, motor.center_of_mass(t) * _ + position]), + inputs="t", + outputs="Vector (m)", + ) + self.center_of_propellant_position = Function( + lambda t: Vector( + [0, 0, motor.center_of_propellant_mass(t) * _ + position] + ), + inputs="t", + outputs="Vector (m)", + ) + self.motor_center_of_dry_mass_position = Vector( + [0, 0, motor.center_of_dry_mass_position * _ + position] + ) + + self.nozzle_position = motor.nozzle_position * _ + position + self.total_mass_flow_rate = motor.total_mass_flow_rate + + self._full_rocket_update() + + def add_cluster_motor(self, motors, positions, orientations=None): """ - if hasattr(self, "motor"): - # pylint: disable=access-member-before-definition - if not isinstance(self.motor, EmptyMotor): - print( - "Only one motor per rocket is currently supported. " - + "Overwriting previous motor." - ) - self.motor = motor - self.motor_position = position - _ = self._csys * self.motor._csys - self.center_of_propellant_position = ( - self.motor.center_of_propellant_mass * _ + self.motor_position - ) - self.motor_center_of_mass_position = ( - self.motor.center_of_mass * _ + self.motor_position - ) - self.motor_center_of_dry_mass_position = ( - self.motor.center_of_dry_mass_position * _ + self.motor_position - ) - self.nozzle_position = self.motor.nozzle_position * _ + self.motor_position + Creates a ClusterMotor and adds it to the rocket. + """ + + cluster = ClusterMotor(motors, positions, orientations) + + self.motor = cluster + self.motor_position = 0 # Positions are handled within the cluster + + self.center_of_propellant_position = self.motor.center_of_propellant_mass + self.motor_center_of_mass_position = self.motor.center_of_mass + self.motor_center_of_dry_mass_position = self.motor.center_of_dry_mass_position self.total_mass_flow_rate = self.motor.total_mass_flow_rate - self.evaluate_dry_mass() - self.evaluate_structural_mass_ratio() - self.evaluate_total_mass() - self.evaluate_center_of_dry_mass() - self.evaluate_nozzle_to_cdm() - self.evaluate_center_of_mass() - self.evaluate_dry_inertias() - self.evaluate_inertias() - self.evaluate_reduced_mass() - self.evaluate_thrust_to_weight() - self.evaluate_center_of_pressure() - self.evaluate_surfaces_cp_to_cdm() - self.evaluate_stability_margin() - self.evaluate_static_margin() - self.evaluate_com_to_cdm_function() - self.evaluate_nozzle_gyration_tensor() + + # Call the centralized update method + self._full_rocket_update() + return cluster def __add_single_surface(self, surface, position): """Adds a single aerodynamic surface to the rocket. Makes checks for @@ -1207,16 +1196,16 @@ def add_nose( self.add_surfaces(nose, position) return nose - @deprecated( - reason="This method is set to be deprecated in version 1.0.0 and fully " - "removed by version 2.0.0", - alternative="Rocket.add_trapezoidal_fins", - ) def add_fins(self, *args, **kwargs): # pragma: no cover """See Rocket.add_trapezoidal_fins for documentation. This method is set to be deprecated in version 1.0.0 and fully removed by version 2.0.0. Use Rocket.add_trapezoidal_fins instead. It keeps the same arguments and signature.""" + warnings.warn( + "This method is set to be deprecated in version 1.0.0 and fully " + "removed by version 2.0.0. Use Rocket.add_trapezoidal_fins instead", + DeprecationWarning, + ) return self.add_trapezoidal_fins(*args, **kwargs) def add_trapezoidal_fins( @@ -1468,16 +1457,7 @@ def add_free_form_fins( return fin_set def add_parachute( - self, - name, - cd_s, - trigger, - sampling_rate=100, - lag=0, - noise=(0, 0, 0), - radius=1.5, - height=None, - porosity=0.0432, + self, name, cd_s, trigger, sampling_rate=100, lag=0, noise=(0, 0, 0) ): """Creates a new parachute, storing its parameters such as opening delay, drag coefficients and trigger function. @@ -1536,39 +1516,16 @@ def add_parachute( The values are used to add noise to the pressure signal which is passed to the trigger function. Default value is (0, 0, 0). Units are in pascal. - radius : float, optional - Length of the non-unique semi-axis (radius) of the inflated hemispheroid - parachute. Default value is 1.5. - Units are in meters. - height : float, optional - Length of the unique semi-axis (height) of the inflated hemispheroid - parachute. Default value is the radius of the parachute. - Units are in meters. - porosity : float, optional - Geometric porosity of the canopy (ratio of open area to total canopy area), - in [0, 1]. Affects only the added-mass scaling during descent; it does - not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass - of 1.0 (“neutral” behavior). Returns ------- parachute : Parachute - Parachute containing trigger, sampling_rate, lag, cd_s, noise, radius, - height, porosity and name. Furthermore, it stores clean_pressure_signal, + Parachute containing trigger, sampling_rate, lag, cd_s, noise + and name. Furthermore, it stores clean_pressure_signal, noise_signal and noisyPressureSignal which are filled in during Flight simulation. """ - parachute = Parachute( - name, - cd_s, - trigger, - sampling_rate, - lag, - noise, - radius, - height, - porosity, - ) + parachute = Parachute(name, cd_s, trigger, sampling_rate, lag, noise) self.parachutes.append(parachute) return self.parachutes[-1] @@ -1957,16 +1914,7 @@ def all_info(self): self.info() self.plots.all() - # pylint: disable=too-many-statements - def to_dict(self, **kwargs): - discretize = kwargs.get("discretize", False) - - power_off_drag = self.power_off_drag - power_on_drag = self.power_on_drag - if discretize: - power_off_drag = power_off_drag.set_discrete(0, 4, 50, mutate_self=False) - power_on_drag = power_on_drag.set_discrete(0, 4, 50, mutate_self=False) - + def to_dict(self, include_outputs=False): rocket_dict = { "radius": self.radius, "mass": self.mass, @@ -1976,8 +1924,8 @@ def to_dict(self, **kwargs): "I_12_without_motor": self.I_12_without_motor, "I_13_without_motor": self.I_13_without_motor, "I_23_without_motor": self.I_23_without_motor, - "power_off_drag": power_off_drag, - "power_on_drag": power_on_drag, + "power_off_drag": self.power_off_drag, + "power_on_drag": self.power_on_drag, "center_of_mass_without_motor": self.center_of_mass_without_motor, "coordinate_system_orientation": self.coordinate_system_orientation, "motor": self.motor, @@ -1990,51 +1938,7 @@ def to_dict(self, **kwargs): "sensors": self.sensors, } - if kwargs.get("include_outputs", False): - thrust_to_weight = self.thrust_to_weight - cp_position = self.cp_position - stability_margin = self.stability_margin - center_of_mass = self.center_of_mass - motor_center_of_mass_position = self.motor_center_of_mass_position - reduced_mass = self.reduced_mass - total_mass = self.total_mass - total_mass_flow_rate = self.total_mass_flow_rate - center_of_propellant_position = self.center_of_propellant_position - - if discretize: - thrust_to_weight = thrust_to_weight.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - cp_position = cp_position.set_discrete(0, 4, 25, mutate_self=False) - stability_margin = stability_margin.set_discrete( - (0, self.motor.burn_time[0]), - (2, self.motor.burn_time[1]), - (10, 10), - mutate_self=False, - ) - center_of_mass = center_of_mass.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - motor_center_of_mass_position = ( - motor_center_of_mass_position.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - ) - reduced_mass = reduced_mass.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - total_mass = total_mass.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - total_mass_flow_rate = total_mass_flow_rate.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - center_of_propellant_position = ( - center_of_propellant_position.set_discrete_based_on_model( - self.motor.thrust, mutate_self=False - ) - ) - + if include_outputs: rocket_dict["area"] = self.area rocket_dict["center_of_dry_mass_position"] = ( self.center_of_dry_mass_position @@ -2042,26 +1946,30 @@ def to_dict(self, **kwargs): rocket_dict["center_of_mass_without_motor"] = ( self.center_of_mass_without_motor ) - rocket_dict["motor_center_of_mass_position"] = motor_center_of_mass_position + rocket_dict["motor_center_of_mass_position"] = ( + self.motor_center_of_mass_position + ) rocket_dict["motor_center_of_dry_mass_position"] = ( self.motor_center_of_dry_mass_position ) - rocket_dict["center_of_mass"] = center_of_mass - rocket_dict["reduced_mass"] = reduced_mass - rocket_dict["total_mass"] = total_mass - rocket_dict["total_mass_flow_rate"] = total_mass_flow_rate - rocket_dict["thrust_to_weight"] = thrust_to_weight + rocket_dict["center_of_mass"] = self.center_of_mass + rocket_dict["reduced_mass"] = self.reduced_mass + rocket_dict["total_mass"] = self.total_mass + rocket_dict["total_mass_flow_rate"] = self.total_mass_flow_rate + rocket_dict["thrust_to_weight"] = self.thrust_to_weight rocket_dict["cp_eccentricity_x"] = self.cp_eccentricity_x rocket_dict["cp_eccentricity_y"] = self.cp_eccentricity_y rocket_dict["thrust_eccentricity_x"] = self.thrust_eccentricity_x rocket_dict["thrust_eccentricity_y"] = self.thrust_eccentricity_y - rocket_dict["cp_position"] = cp_position - rocket_dict["stability_margin"] = stability_margin + rocket_dict["cp_position"] = self.cp_position + rocket_dict["stability_margin"] = self.stability_margin rocket_dict["static_margin"] = self.static_margin rocket_dict["nozzle_position"] = self.nozzle_position rocket_dict["nozzle_to_cdm"] = self.nozzle_to_cdm rocket_dict["nozzle_gyration_tensor"] = self.nozzle_gyration_tensor - rocket_dict["center_of_propellant_position"] = center_of_propellant_position + rocket_dict["center_of_propellant_position"] = ( + self.center_of_propellant_position + ) return rocket_dict @@ -2104,29 +2012,17 @@ def from_dict(cls, data): for parachute in data["parachutes"]: rocket.parachutes.append(parachute) - for sensor, position in data["sensors"]: - rocket.add_sensor(sensor, position) - - for air_brake in data["air_brakes"]: - rocket.air_brakes.append(air_brake) - - for controller in data["_controllers"]: - interactive_objects_hash = getattr(controller, "_interactive_objects_hash") - if interactive_objects_hash is not None: - is_iterable = isinstance(interactive_objects_hash, Iterable) - if not is_iterable: - interactive_objects_hash = [interactive_objects_hash] - for hash_ in interactive_objects_hash: - if (hashed_obj := find_obj_from_hash(data, hash_)) is not None: - if not is_iterable: - controller.interactive_objects = hashed_obj - else: - controller.interactive_objects.append(hashed_obj) - else: - warnings.warn( - "Could not find controller interactive objects." - "Deserialization will proceed, results may not be accurate." - ) - rocket._add_controllers(controller) + for air_brakes in data["air_brakes"]: + rocket.add_air_brakes( + drag_coefficient_curve=air_brakes["drag_coefficient_curve"], + controller_function=air_brakes["controller_function"], + sampling_rate=air_brakes["sampling_rate"], + clamp=air_brakes["clamp"], + reference_area=air_brakes["reference_area"], + initial_observed_variables=air_brakes["initial_observed_variables"], + override_rocket_drag=air_brakes["override_rocket_drag"], + name=air_brakes["name"], + controller_name=air_brakes["controller_name"], + ) return rocket diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 30ea66466..b294c7ec0 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -7,8 +7,6 @@ import numpy as np from scipy.integrate import BDF, DOP853, LSODA, RK23, RK45, OdeSolver, Radau -from rocketpy.simulation.flight_data_exporter import FlightDataExporter - from ..mathutils.function import Function, funcify_method from ..mathutils.vector_matrix import Matrix, Vector from ..motors.point_mass_motor import PointMassMotor @@ -17,7 +15,6 @@ from ..rocket import PointMassRocket from ..tools import ( calculate_cubic_hermite_coefficients, - deprecated, euler313_to_quaternions, find_closest, find_root_linear_interpolation, @@ -26,6 +23,7 @@ quaternions_to_precession, quaternions_to_spin, ) +from .flight_data_exporter import FlightDataExporter ODE_SOLVER_MAP = { "RK23": RK23, @@ -712,8 +710,9 @@ def __simulate(self, verbose): # u_dot for all sensors u_dot = phase.derivative(self.t, self.y_sol) for sensor, position in node._component_sensors: + cdm_z = self.rocket.center_of_dry_mass_position.z relative_position = position - self.rocket._csys * Vector( - [0, 0, self.rocket.center_of_dry_mass_position] + [0, 0, cdm_z] ) sensor.measure( self.t, @@ -1908,90 +1907,64 @@ def u_dot_generalized_3dof(self, t, u, post_processing=False): return u_dot - def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too-many-locals,too-many-statements - """Calculates derivative of u state vector with respect to time when the - rocket is flying in 6 DOF motion in space and significant mass variation - effects exist. Typical flight phases include powered ascent after launch - rail. - - Parameters - ---------- - t : float - Time in seconds - u : list - State vector defined by u = [x, y, z, vx, vy, vz, q0, q1, - q2, q3, omega1, omega2, omega3]. - post_processing : bool, optional - If True, adds flight data information directly to self variables - such as self.angle_of_attack, by default False. - - Returns - ------- - u_dot : list - State vector defined by u_dot = [vx, vy, vz, ax, ay, az, - e0_dot, e1_dot, e2_dot, e3_dot, alpha1, alpha2, alpha3]. + def u_dot_generalized(self, t, u, post_processing=False): """ - # Retrieve integration data - _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u + Calculates the derivative of the state vector u (6-DOF) using a + generalized formalism suitable for motor clusters and including all + physical forces. + """ + # State variable extraction + _, _, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3 = u - # Create necessary vectors - # r = Vector([x, y, z]) # CDM position vector - v = Vector([vx, vy, vz]) # CDM velocity vector - e = [e0, e1, e2, e3] # Euler parameters/quaternions - w = Vector([omega1, omega2, omega3]) # Angular velocity vector + # Create vectors for calculations + velocity_vector = Vector([vx, vy, vz]) + euler_params = [e0, e1, e2, e3] + angular_velocity_vector = Vector([w1, w2, w3]) + w = angular_velocity_vector - # Retrieve necessary quantities - ## Rocket mass + # Rocket properties at time t total_mass = self.rocket.total_mass.get_value_opt(t) - total_mass_dot = self.rocket.total_mass_flow_rate.get_value_opt(t) - total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate_complex_step(t) - ## CM position vector and time derivatives relative to CDM in body frame - r_CM_z = self.rocket.com_to_cdm_function - r_CM_t = r_CM_z.get_value_opt(t) - r_CM = Vector([0, 0, r_CM_t]) - r_CM_dot = Vector([0, 0, r_CM_z.differentiate_complex_step(t)]) - r_CM_ddot = Vector([0, 0, r_CM_z.differentiate(t, order=2)]) - ## Nozzle position vector - r_NOZ = Vector([0, 0, self.rocket.nozzle_to_cdm]) - ## Nozzle gyration tensor - S_nozzle = self.rocket.nozzle_gyration_tensor - ## Inertia tensor - inertia_tensor = self.rocket.get_inertia_tensor_at_time(t) - ## Inertia tensor time derivative in the body frame - I_dot = self.rocket.get_inertia_tensor_derivative_at_time(t) - - # Calculate the Inertia tensor relative to CM - H = (r_CM.cross_matrix @ -r_CM.cross_matrix) * total_mass - I_CM = inertia_tensor - H - # Prepare transformation matrices - K = Matrix.transformation(e) + # Transformation matrix from body frame to inertial frame + K = Matrix.transformation(euler_params) Kt = K.transpose - # Compute aerodynamic forces and moments + # --- Forces and Moments Calculation --- + + # Gravity in body frame + gravity_force_inertial = Vector( + [0, 0, -total_mass * self.env.gravity.get_value_opt(z)] + ) + gravity_force_body = Kt @ gravity_force_inertial + + # --- START OF INTEGRATED AERODYNAMIC CALCULATION BLOCK --- R1, R2, R3, M1, M2, M3 = 0, 0, 0, 0, 0, 0 - ## Drag force - rho = self.env.density.get_value_opt(z) + # Get freestream speed wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z) wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z) - wind_velocity = Vector([wind_velocity_x, wind_velocity_y, 0]) - free_stream_speed = abs((wind_velocity - Vector(v))) speed_of_sound = self.env.speed_of_sound.get_value_opt(z) + rho = self.env.density.get_value_opt(z) + + # Relative wind speed in inertial frame + stream_velocity_inertial = Vector( + [wind_velocity_x - vx, wind_velocity_y - vy, -vz] + ) + free_stream_speed = abs(stream_velocity_inertial) free_stream_mach = free_stream_speed / speed_of_sound - if self.rocket.motor.burn_start_time < t < self.rocket.motor.burn_out_time: - pressure = self.env.pressure.get_value_opt(z) - net_thrust = max( - self.rocket.motor.thrust.get_value_opt(t) - + self.rocket.motor.pressure_thrust(pressure), - 0, - ) + # Determine the base drag of the rocket + if t < self.rocket.motor.burn_out_time: drag_coeff = self.rocket.power_on_drag.get_value_opt(free_stream_mach) else: - net_thrust = 0 drag_coeff = self.rocket.power_off_drag.get_value_opt(free_stream_mach) - R3 += -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff + + R3_base_drag = ( + -0.5 * rho * (free_stream_speed**2) * self.rocket.area * drag_coeff + ) + R3 += R3_base_drag + + # Add air brakes drag if active for air_brakes in self.rocket.air_brakes: if air_brakes.deployment_level > 0: air_brakes_cd = air_brakes.drag_coefficient.get_value_opt( @@ -2005,36 +1978,40 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too * air_brakes_cd ) if air_brakes.override_rocket_drag: - R3 = air_brakes_force # Substitutes rocket drag coefficient + R3 = air_brakes_force else: R3 += air_brakes_force - # Get rocket velocity in body frame - velocity_in_body_frame = Kt @ v - # Calculate lift and moment for each component of the rocket + + # Moment due to CoP eccentricity + M1 += self.rocket.cp_eccentricity_y * R3 + M2 -= self.rocket.cp_eccentricity_x * R3 + + # Rocket velocity in body frame + velocity_body_frame = Kt @ velocity_vector + + # Calculate lift and moment for each aerodynamic surface for aero_surface, _ in self.rocket.aerodynamic_surfaces: - # Component cp relative to CDM in body frame comp_cp = self.rocket.surfaces_cp_to_cdm[aero_surface] # Component absolute velocity in body frame - comp_vb = velocity_in_body_frame + (w ^ comp_cp) + comp_vb = velocity_body_frame + (w ^ comp_cp) # Wind velocity at component altitude comp_z = z + (K @ comp_cp).z comp_wind_vx = self.env.wind_velocity_x.get_value_opt(comp_z) comp_wind_vy = self.env.wind_velocity_y.get_value_opt(comp_z) - # Component freestream velocity in body frame + # Wind velocity in body frame comp_wind_vb = Kt @ Vector([comp_wind_vx, comp_wind_vy, 0]) comp_stream_velocity = comp_wind_vb - comp_vb comp_stream_speed = abs(comp_stream_velocity) comp_stream_mach = comp_stream_speed / speed_of_sound - # Reynolds at component altitude - # TODO: Reynolds is only used in generic surfaces. This calculation - # should be moved to the surface class for efficiency + comp_reynolds = ( self.env.density.get_value_opt(comp_z) * comp_stream_speed * aero_surface.reference_length / self.env.dynamic_viscosity.get_value_opt(comp_z) ) - # Forces and moments + + # Calculate forces and moments for the surface X, Y, Z, M, N, L = aero_surface.compute_forces_and_moments( comp_stream_velocity, comp_stream_speed, @@ -2051,70 +2028,87 @@ def u_dot_generalized(self, t, u, post_processing=False): # pylint: disable=too M2 += N M3 += L - # Off center moment - M1 += ( - self.rocket.cp_eccentricity_y * R3 - + self.rocket.thrust_eccentricity_y * net_thrust - ) - M2 -= ( - self.rocket.cp_eccentricity_x * R3 - + self.rocket.thrust_eccentricity_x * net_thrust - ) + # Moment due to CoP eccentricity M3 += self.rocket.cp_eccentricity_x * R2 - self.rocket.cp_eccentricity_y * R1 - weight_in_body_frame = Kt @ Vector( - [0, 0, -total_mass * self.env.gravity.get_value_opt(z)] - ) + # Create final vectors for aero forces and moments + R_aero_body = Vector([R1, R2, R3]) + M_aero_body = Vector([M1, M2, M3]) - T00 = total_mass * r_CM - T03 = 2 * total_mass_dot * (r_NOZ - r_CM) - 2 * total_mass * r_CM_dot - T04 = ( - Vector([0, 0, net_thrust]) - - total_mass * r_CM_ddot - - 2 * total_mass_dot * r_CM_dot - + total_mass_ddot * (r_NOZ - r_CM) - ) - T05 = total_mass_dot * S_nozzle - I_dot - - T20 = ( - ((w ^ T00) ^ w) - + (w ^ T03) - + T04 - + weight_in_body_frame - + Vector([R1, R2, R3]) - ) + # --- END OF AERODYNAMIC CALCULATION BLOCK --- - T21 = ( - ((inertia_tensor @ w) ^ w) - + T05 @ w - - (weight_in_body_frame ^ r_CM) - + Vector([M1, M2, M3]) - ) + # Thrust and moment from cluster/motor + if hasattr(self.rocket.motor, "get_total_thrust_vector"): # It's a ClusterMotor + thrust_vector_body = self.rocket.motor.get_total_thrust_vector(t) + rocket_cm = self.rocket.center_of_mass.get_value_opt(t) + moment_vector_body = self.rocket.motor.get_total_moment(t, rocket_cm) + net_thrust_scalar = abs(thrust_vector_body) + else: # It's a standard motor + pressure = self.env.pressure.get_value_opt(z) + net_thrust_scalar = max( + self.rocket.motor.thrust.get_value_opt(t) + + self.rocket.motor.pressure_thrust(pressure), + 0, + ) + thrust_vector_body = Vector([0, 0, net_thrust_scalar]) # Axial thrust + moment_vector_body = Vector( + [ # Moment due to eccentricity + self.rocket.thrust_eccentricity_y * net_thrust_scalar, + -self.rocket.thrust_eccentricity_x * net_thrust_scalar, + 0, + ] + ) - # Angular velocity derivative - w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM)) + # Total force and moment in body frame + total_force_body = gravity_force_body + thrust_vector_body + R_aero_body + total_moment_body = moment_vector_body + M_aero_body - # Euler parameters derivative - e_dot = [ - 0.5 * (-omega1 * e1 - omega2 * e2 - omega3 * e3), - 0.5 * (omega1 * e0 + omega3 * e2 - omega2 * e3), - 0.5 * (omega2 * e0 - omega3 * e1 + omega1 * e3), - 0.5 * (omega3 * e0 + omega2 * e1 - omega1 * e2), - ] + # --- Equations of Motion --- - # Velocity vector derivative + Coriolis acceleration - w_earth = Vector(self.env.earth_rotation_vector) - v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2 * (w_earth ^ v) + # Linear acceleration + linear_acceleration_body = total_force_body / total_mass + linear_acceleration_inertial = K @ linear_acceleration_body - # Position vector derivative - r_dot = [vx, vy, vz] + # Angular acceleration (Euler's equation for rigid body with variable mass) + inertia_tensor = self.rocket.get_inertia_tensor_at_time(t) + I_w = inertia_tensor @ angular_velocity_vector + w_cross_Iw = angular_velocity_vector.cross(I_w) + I_dot = self.rocket.get_inertia_tensor_derivative_at_time(t) - # Create u_dot - u_dot = [*r_dot, *v_dot, *e_dot, *w_dot] + angular_acceleration_body = inertia_tensor.inverse @ ( + total_moment_body - w_cross_Iw - (I_dot @ angular_velocity_vector) + ) + + # Euler parameters derivatives + e0dot = 0.5 * (-w1 * e1 - w2 * e2 - w3 * e3) + e1dot = 0.5 * (w1 * e0 + w3 * e2 - w2 * e3) + e2dot = 0.5 * (w2 * e0 - w3 * e1 + w1 * e3) + e3dot = 0.5 * (w3 * e0 + w2 * e1 - w1 * e2) + + # Assemble the derivative vector + u_dot = [ + vx, + vy, + vz, + *linear_acceleration_inertial, + e0dot, + e1dot, + e2dot, + e3dot, + *angular_acceleration_body, + ] + # Post-processing (if necessary) if post_processing: self.__post_processed_variables.append( - [t, *v_dot, *w_dot, R1, R2, R3, M1, M2, M3, net_thrust] + [ + t, + *linear_acceleration_inertial, + *angular_acceleration_body, + *R_aero_body, + *total_moment_body, + net_thrust_scalar, + ] ) return u_dot @@ -3276,13 +3270,10 @@ def __calculate_rail_button_forces(self): # TODO: complex method. ) lower_button_position = rail_buttons_tuple.position.z angular_position_rad = rail_buttons_tuple.component.angular_position_rad - D1 = ( - upper_button_position - self.rocket.center_of_dry_mass_position - ) * self.rocket._csys + cdm_z = self.rocket.center_of_dry_mass_position.z + D1 = (upper_button_position - cdm_z) * self.rocket._csys # Distance from Rail Button 2 (lower) to CM - D2 = ( - lower_button_position - self.rocket.center_of_dry_mass_position - ) * self.rocket._csys + D2 = (lower_button_position - cdm_z) * self.rocket._csys F11 = (self.R1 * D2 - self.M2) / (D1 + D2) F11.set_outputs("Upper button force direction 1 (m)") F12 = (self.R2 * D2 + self.M1) / (D1 + D2) diff --git a/rocketpy/tools.py b/rocketpy/tools.py index bc80b469c..00098b46c 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -1081,6 +1081,7 @@ def wrapper(*args, **kwargs): return decorator +# TODO: move it to inertia.py file when possible def parallel_axis_theorem_from_com(com_inertia_moment, mass, distance): """Calculates the moment of inertia of a object relative to a new axis using the parallel axis theorem. The new axis is parallel to and at a distance