diff --git a/CHANGELOG.md b/CHANGELOG.md index 426a4fa6f..d2c854529 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/). ### Added +- ENH: Pre-calculate attributes in Rocket class [#595](https://github.com/RocketPy-Team/RocketPy/pull/595) - ENH: Complex step differentiation [#594](https://github.com/RocketPy-Team/RocketPy/pull/594) - ENH: Exponential backoff decorator (fix #449) [#588](https://github.com/RocketPy-Team/RocketPy/pull/588) - ENH: Function Validation Rework & Swap `np.searchsorted` to `bisect_left` [#582](https://github.com/RocketPy-Team/RocketPy/pull/582) diff --git a/docs/user/rocket.rst b/docs/user/rocket.rst index 28e091e8d..fa09f058e 100644 --- a/docs/user/rocket.rst +++ b/docs/user/rocket.rst @@ -12,6 +12,7 @@ Defining a Rocket in RocketPy is simple and requires a few steps: 4. Add, if desired, parachutes; 5. Set, if desired, rail guides; 6. See results. +7. Inertia Tensors. Lets go through each of these steps in detail. @@ -428,3 +429,54 @@ The lets check all the information available about the rocket: .. jupyter-execute:: calisto.all_info() + +7. Inertia Tensors +------------------ + +The inertia tensor of the rocket at a given time can be obtained using the +``get_inertia_tensor_at_time`` method. This method evaluates each component of +the inertia tensor at the specified time and returns a +:class:`rocketpy.mathutils.Matrix` object. + +The inertia tensor is a matrix that looks like this: + +.. math:: + :label: inertia_tensor + + \mathbf{I} = \begin{bmatrix} + I_{11} & I_{12} & I_{13} \\ + I_{21} & I_{22} & I_{23} \\ + I_{31} & I_{32} & I_{33} + \end{bmatrix} + +For example, to get the inertia tensor of the rocket at time 0.5 seconds, you +can use the following code: + +.. jupyter-execute:: + + calisto.get_inertia_tensor_at_time(0.5) + +Derivative of the Inertia Tensor +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +You can also get the derivative of the inertia tensor at a given time using the +``get_inertia_tensor_derivative_at_time`` method. Here's an example: + +.. jupyter-execute:: + + calisto.get_inertia_tensor_derivative_at_time(0.5) + +Implications from these results +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The inertia tensor reveals important information about the rocket's symmetry +and ease of rotation: + +1. **Axis Symmetry**: If I\ :sub:`11` and I\ :sub:`22` are equal, the rocket is symmetric around the axes perpendicular to the rocket's center axis. In our defined rocket, I\ :sub:`11` and I\ :sub:`22` are indeed equal, indicating that our rocket is axisymmetric. + +2. **Zero Products of Inertia**: The off-diagonal elements of the inertia tensor are zero, which means the products of inertia are zero. This indicates that the rocket is symmetric around its center axis. + +3. **Ease of Rotation**: The I\ :sub:`33` value is significantly lower than the other two. This suggests that the rocket is easier to rotate around its center axis than around the axes perpendicular to the rocket. This is an important factor when considering the rocket's stability and control. + +However, these conclusions are based on the assumption that the inertia tensor is calculated with respect to the rocket's center of mass and aligned with the principal axes of the rocket. If the inertia tensor is calculated with respect to a different point or not aligned with the principal axes, the conclusions may not hold. + diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py index d41240ac9..8f0d37f29 100644 --- a/rocketpy/rocket/aero_surface.py +++ b/rocketpy/rocket/aero_surface.py @@ -367,7 +367,7 @@ def evaluate_geometrical_parameters(self): return None def evaluate_nose_shape(self): - """Calculates and saves nose cone's shape as lists and reavulates the + """Calculates and saves nose cone's shape as lists and re-evaluates the nose cone's length for a given bluffness ratio. The shape is saved as two vectors, one for the x coordinates and one for the y coordinates. @@ -381,7 +381,9 @@ def evaluate_nose_shape(self): # Calculate a function to find the tangential intersection point between the circle and nosecone curve. def find_x_intercept(x): - return x + self.y_nosecone(x) * self.y_nosecone.differentiate(x) + return x + self.y_nosecone(x) * self.y_nosecone.differentiate_complex_step( + x + ) # Calculate a function to find the radius of the nosecone curve def find_radius(x): @@ -770,7 +772,9 @@ def evaluate_lift_coefficient(self): ) # Differentiating at alpha = 0 to get cl_alpha - clalpha2D_incompressible = self.airfoil_cl.differentiate(x=1e-3, dx=1e-3) + clalpha2D_incompressible = self.airfoil_cl.differentiate_complex_step( + x=1e-3, dx=1e-3 + ) # Convert to radians if needed if self.airfoil[1] == "degrees": diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index d09bd2acb..a0db26eff 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -4,6 +4,7 @@ from rocketpy.control.controller import _Controller from rocketpy.mathutils.function import Function +from rocketpy.mathutils.vector_matrix import Matrix from rocketpy.motors.motor import EmptyMotor from rocketpy.plots.rocket_plots import _RocketPlots from rocketpy.prints.rocket_prints import _RocketPrints @@ -72,6 +73,9 @@ class Rocket: for more information regarding the coordinate system. Expressed in meters as a function of time. + Rocket.com_to_cdm_function : Function + Function of time expressing the z-coordinate of the center of mass + relative to the center of dry mass. Rocket.reduced_mass : Function Function of time expressing the reduced mass of the rocket, defined as the product of the propellant mass and the mass @@ -146,6 +150,11 @@ class Rocket: defined rocket coordinate system. See :doc:`Positions and Coordinate Systems ` for more information. + Rocket.nozzle_to_cdm : float + Distance between the nozzle exit and the rocket's center of dry mass + position, in meters. + Rocket.nozzle_gyration_tensor: Matrix + Matrix representing the nozzle gyration tensor. Rocket.center_of_propellant_position : Function Position of the propellant's center of mass relative to the user defined rocket reference system. See @@ -734,8 +743,133 @@ def evaluate_inertias(self): self.I_23, ) + 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 + ) + return self.nozzle_to_cdm + def evaluate_nozzle_gyration_tensor(self): - pass + """Calculates and returns the nozzle gyration tensor relative to the + rocket's center of dry mass. The gyration tensor is saved and returned + in units of kg*m². + + Returns + ------- + self.nozzle_gyration_tensor : Matrix + Matrix containing the nozzle gyration tensor. + """ + S_noz_33 = 0.5 * self.motor.nozzle_radius**2 + S_noz_11 = S_noz_22 = 0.5 * S_noz_33 + 0.25 * self.nozzle_to_cdm**2 + S_noz_12, S_noz_13, S_noz_23 = 0, 0, 0 # Due to axis symmetry + self.nozzle_gyration_tensor = Matrix( + [ + [S_noz_11, S_noz_12, S_noz_13], + [S_noz_12, S_noz_22, S_noz_23], + [S_noz_13, S_noz_23, S_noz_33], + ] + ) + return self.nozzle_gyration_tensor + + 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. + + 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 + ) + * 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") + return self.com_to_cdm_function + + def get_inertia_tensor_at_time(self, t): + """Returns a Matrix representing the inertia tensor of the rocket with + respect to the rocket's center of mass at a given time. It evaluates + each inertia tensor component at the given time and returns a Matrix + with the computed values. + + Parameters + ---------- + t : float + Time at which the inertia tensor is to be evaluated. + + Returns + ------- + Matrix + Inertia tensor of the rocket at time t. + """ + I_11 = self.I_11.get_value_opt(t) + I_12 = self.I_12.get_value_opt(t) + I_13 = self.I_13.get_value_opt(t) + I_22 = self.I_22.get_value_opt(t) + I_23 = self.I_23.get_value_opt(t) + I_33 = self.I_33.get_value_opt(t) + return Matrix( + [ + [I_11, I_12, I_13], + [I_12, I_22, I_23], + [I_13, I_23, I_33], + ] + ) + + def get_inertia_tensor_derivative_at_time(self, t): + """Returns a Matrix representing the time derivative of the inertia + tensor of the rocket with respect to the rocket's center of mass at a + given time. It evaluates each inertia tensor component's derivative at + the given time and returns a Matrix with the computed values. + + Parameters + ---------- + t : float + Time at which the inertia tensor derivative is to be evaluated. + + Returns + ------- + Matrix + Inertia tensor time derivative of the rocket at time t. + """ + I_11_dot = self.I_11.differentiate_complex_step(t) + I_12_dot = self.I_12.differentiate_complex_step(t) + I_13_dot = self.I_13.differentiate_complex_step(t) + I_22_dot = self.I_22.differentiate_complex_step(t) + I_23_dot = self.I_23.differentiate_complex_step(t) + I_33_dot = self.I_33.differentiate_complex_step(t) + return Matrix( + [ + [I_11_dot, I_12_dot, I_13_dot], + [I_12_dot, I_22_dot, I_23_dot], + [I_13_dot, I_23_dot, I_33_dot], + ] + ) def add_motor(self, motor, position): """Adds a motor to the rocket. @@ -778,6 +912,7 @@ def add_motor(self, motor, position): self.evaluate_dry_mass() 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() @@ -786,6 +921,8 @@ def add_motor(self, motor, position): self.evaluate_center_of_pressure() self.evaluate_stability_margin() self.evaluate_static_margin() + self.evaluate_com_to_cdm_function() + self.evaluate_nozzle_gyration_tensor() def add_surfaces(self, surfaces, positions): """Adds one or more aerodynamic surfaces to the rocket. The aerodynamic diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index b4bc1a3be..4cbc752fb 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1341,12 +1341,8 @@ def u_dot(self, t, u, post_processing=False): ) * self.rocket._csys ) - # c = -self.rocket.distance_rocket_nozzle - c = ( - -(self.rocket.nozzle_position - self.rocket.center_of_dry_mass_position) - * self.rocket._csys - ) - a = b * Mt / M + c = self.rocket.nozzle_to_cdm + a = self.rocket.com_to_cdm_function.get_value_opt(t) rN = self.rocket.motor.nozzle_radius # Prepare transformation matrix a11 = 1 - 2 * (e2**2 + e3**2) @@ -1570,70 +1566,26 @@ def u_dot_generalized(self, t, u, post_processing=False): w = Vector([omega1, omega2, omega3]) # Angular velocity vector # Retrieve necessary quantities - rho = self.env.density.get_value_opt(z) + ## Rocket mass 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(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 = ( - -1 - * ( - ( - self.rocket.center_of_propellant_position - - self.rocket.center_of_dry_mass_position - ) - * self.rocket._csys - ) - * self.rocket.motor.propellant_mass - / total_mass - ) - r_CM = Vector([0, 0, r_CM_z.get_value_opt(t)]) - r_CM_dot = Vector([0, 0, r_CM_z.differentiate(t)]) + 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 - r_NOZ = ( - -(self.rocket.nozzle_position - self.rocket.center_of_dry_mass_position) - * self.rocket._csys - ) - S_noz_33 = 0.5 * self.rocket.motor.nozzle_radius**2 - S_noz_11 = S_noz_22 = 0.5 * S_noz_33 + 0.25 * r_NOZ**2 - S_noz_12, S_noz_13, S_noz_23 = 0, 0, 0 - S_nozzle = Matrix( - [ - [S_noz_11, S_noz_12, S_noz_13], - [S_noz_12, S_noz_22, S_noz_23], - [S_noz_13, S_noz_23, S_noz_33], - ] - ) + S_nozzle = self.rocket.nozzle_gyration_tensor ## Inertia tensor - I_11 = self.rocket.I_11.get_value_opt(t) - I_12 = self.rocket.I_12.get_value_opt(t) - I_13 = self.rocket.I_13.get_value_opt(t) - I_22 = self.rocket.I_22.get_value_opt(t) - I_23 = self.rocket.I_23.get_value_opt(t) - I_33 = self.rocket.I_33.get_value_opt(t) - I = Matrix( - [ - [I_11, I_12, I_13], - [I_12, I_22, I_23], - [I_13, I_23, I_33], - ] - ) + I = self.rocket.get_inertia_tensor_at_time(t) ## Inertia tensor time derivative in the body frame - I_11_dot = self.rocket.I_11.differentiate(t) - I_12_dot = self.rocket.I_12.differentiate(t) - I_13_dot = self.rocket.I_13.differentiate(t) - I_22_dot = self.rocket.I_22.differentiate(t) - I_23_dot = self.rocket.I_23.differentiate(t) - I_33_dot = self.rocket.I_33.differentiate(t) - I_dot = Matrix( - [ - [I_11_dot, I_12_dot, I_13_dot], - [I_12_dot, I_22_dot, I_23_dot], - [I_13_dot, I_23_dot, I_33_dot], - ] - ) - ## Inertia tensor relative to CM + 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 = I - H @@ -1723,8 +1675,8 @@ def u_dot_generalized(self, t, u, post_processing=False): R1 += comp_lift_xb R2 += comp_lift_yb # Add to total moment - M1 -= (comp_cpz + r_CM_z.get_value_opt(t)) * comp_lift_yb - M2 += (comp_cpz + r_CM_z.get_value_opt(t)) * comp_lift_xb + M1 -= (comp_cpz + r_CM_t) * comp_lift_yb + M2 += (comp_cpz + r_CM_t) * comp_lift_xb # Calculates Roll Moment try: clf_delta, cld_omega, cant_angle_rad = aero_surface.roll_parameters @@ -1749,13 +1701,12 @@ def u_dot_generalized(self, t, u, post_processing=False): pass weightB = Kt @ Vector([0, 0, -total_mass * self.env.gravity.get_value_opt(z)]) T00 = total_mass * r_CM - r_NOZ_vector = Vector([0, 0, r_NOZ]) - T03 = 2 * total_mass_dot * (r_NOZ_vector - r_CM) - 2 * total_mass * r_CM_dot + T03 = 2 * total_mass_dot * (r_NOZ - r_CM) - 2 * total_mass * r_CM_dot T04 = ( - self.rocket.motor.thrust.get_value_opt(t) * Vector([0, 0, 1]) + Vector([0, 0, self.rocket.motor.thrust.get_value_opt(t)]) - total_mass * r_CM_ddot - 2 * total_mass_dot * r_CM_dot - + total_mass_ddot * (r_NOZ_vector - r_CM) + + total_mass_ddot * (r_NOZ - r_CM) ) T05 = total_mass_dot * S_nozzle - I_dot diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py index 7d942ee60..4d934efef 100644 --- a/tests/unit/test_rocket.py +++ b/tests/unit/test_rocket.py @@ -413,3 +413,86 @@ def test_evaluate_center_of_mass(calisto): A predefined instance of the calisto Rocket with a motor, used as a base for testing. """ assert isinstance(calisto.evaluate_center_of_mass(), Function) + + +def test_evaluate_nozzle_to_cdm(calisto): + expected_distance = 1.255 + atol = 1e-3 # Equivalent to 1mm + assert pytest.approx(expected_distance, atol) == calisto.nozzle_to_cdm + # Test if calling the function returns the same result + res = calisto.evaluate_nozzle_to_cdm() + assert pytest.approx(expected_distance, atol) == res + + +def test_evaluate_nozzle_gyration_tensor(calisto): + expected_gyration_tensor = np.array( + [[0.3940207, 0, 0], [0, 0.3940207, 0], [0, 0, 0.0005445]] + ) + atol = 1e-3 * 1e-2 * 1e-2 # Equivalent to 1g * 1cm^2 + assert np.allclose( + expected_gyration_tensor, np.array(calisto.nozzle_gyration_tensor), atol=atol + ) + # Test if calling the function returns the same result + res = calisto.evaluate_nozzle_gyration_tensor() + assert np.allclose(expected_gyration_tensor, np.array(res), atol=atol) + + +def test_evaluate_com_to_cdm_function(calisto): + atol = 1e-3 # Equivalent to 1mm + assert np.allclose( + (calisto.center_of_dry_mass_position - calisto.center_of_mass).source, + calisto.com_to_cdm_function.source, + atol=atol, + ) + + +def test_get_inertia_tensor_at_time(calisto): + # Expected values (for t = 0) + # TODO: compute these values by hand or using CAD. + Ixx = 10.31379 + Iyy = 10.31379 + Izz = 0.039942 + + # Set tolerance threshold + atol = 1e-5 + + # Get inertia tensor at t = 0 + inertia_tensor = calisto.get_inertia_tensor_at_time(0) + + # Check if the values are close to the expected ones + assert pytest.approx(Ixx, atol) == inertia_tensor.x[0] + assert pytest.approx(Iyy, atol) == inertia_tensor.y[1] + assert pytest.approx(Izz, atol) == inertia_tensor.z[2] + # Check if products of inertia are zero + assert pytest.approx(0, atol) == inertia_tensor.x[1] + assert pytest.approx(0, atol) == inertia_tensor.x[2] + assert pytest.approx(0, atol) == inertia_tensor.y[0] + assert pytest.approx(0, atol) == inertia_tensor.y[2] + assert pytest.approx(0, atol) == inertia_tensor.z[0] + assert pytest.approx(0, atol) == inertia_tensor.z[1] + + +def test_get_inertia_tensor_derivative_at_time(calisto): + # Expected values (for t = 2s) + # TODO: compute these values by hand or using CAD. + Ixx_dot = -0.634805230901143 + Iyy_dot = -0.634805230901143 + Izz_dot = -0.000671493662305 + + # Set tolerance threshold + atol = 1e-3 + + # Get inertia tensor at t = 2s + inertia_tensor = calisto.get_inertia_tensor_derivative_at_time(2) + + # Check if the values are close to the expected ones + assert pytest.approx(Ixx_dot, atol) == inertia_tensor.x[0] + assert pytest.approx(Iyy_dot, atol) == inertia_tensor.y[1] + assert pytest.approx(Izz_dot, atol) == inertia_tensor.z[2] + # Check if products of inertia are zero + assert pytest.approx(0, atol) == inertia_tensor.x[1] + assert pytest.approx(0, atol) == inertia_tensor.x[2] + assert pytest.approx(0, atol) == inertia_tensor.y[0] + assert pytest.approx(0, atol) == inertia_tensor.y[2] + assert pytest.approx(0, atol) == inertia_tensor.z[0] + assert pytest.approx(0, atol) == inertia_tensor.z[1]