Merge pull request #674 from RocketPy-Team/bug/rotational-eoms-not-relative-to-cdm

BUG: Rotational EOMs Not Relative To CDM

Author: Gui-FernandesBR

.pylintrc

@@ -225,8 +225,7 @@ good-names=FlightPhases,
            motor_I_11_derivative_at_t,
            M3_forcing,
            M3_damping,
-           CM_to_CDM,
-           CM_to_CPM,
+           CDM_to_CPM,
            center_of_mass_without_motor_to_CDM,
            motor_center_of_dry_mass_to_CDM,
            generic_motor_cesaroni_M1520,

CHANGELOG.md

@@ -43,6 +43,7 @@ Attention: The newest changes should be on top -->
 
 ### Fixed
 
+- BUG: Rotational EOMs Not Relative To CDM [#674](https://github.com/RocketPy-Team/RocketPy/pull/674)
 - BUG: Pressure ISA Extrapolation as "linear" [#675](https://github.com/RocketPy-Team/RocketPy/pull/675)
 - BUG: fix the Frequency Response plot of Flight class [#653](https://github.com/RocketPy-Team/RocketPy/pull/653)
 

docs/user/rocket/rocket_usage.rst

@@ -41,6 +41,8 @@ and radius:
     Pay special attention to the following:
 
     - ``mass`` is the rocket's mass, **without the motor**, in kg.
+    - All ``inertia`` values are given in relation to the rocket's center of 
+      mass without motor.
     - ``inertia`` is defined as a tuple of the form ``(I11, I22, I33)``.
       Where ``I11`` and ``I22`` are the inertia of the mass around the
       perpendicular axes to the rocket, and ``I33`` is the inertia around the
@@ -432,10 +434,10 @@ The lets check all the information available about the rocket:
 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 relative to the center of dry mass 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:
 

rocketpy/motors/motor.py

@@ -1352,6 +1352,12 @@ def __init__(self):
         self.dry_I_12 = 0
         self.dry_I_13 = 0
         self.dry_I_23 = 0
+        self.propellant_I_11 = Function(0, "Time (s)", "Propellant I_11 (kg m²)")
+        self.propellant_I_22 = Function(0, "Time (s)", "Propellant I_22 (kg m²)")
+        self.propellant_I_33 = Function(0, "Time (s)", "Propellant I_33 (kg m²)")
+        self.propellant_I_12 = Function(0, "Time (s)", "Propellant I_12 (kg m²)")
+        self.propellant_I_13 = Function(0, "Time (s)", "Propellant I_13 (kg m²)")
+        self.propellant_I_23 = Function(0, "Time (s)", "Propellant I_23 (kg m²)")
         self.I_11 = Function(0)
         self.I_22 = Function(0)
         self.I_33 = Function(0)

rocketpy/rocket/rocket.py

@@ -596,7 +596,7 @@ def evaluate_static_margin(self):
 
     def evaluate_dry_inertias(self):
         """Calculates and returns the rocket's dry inertias relative to
-        the rocket's center of mass. The inertias are saved and returned
+        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.
 
@@ -605,27 +605,27 @@ def evaluate_dry_inertias(self):
         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 instantaneous center of mass.
+            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 instantaneous center of mass.
+            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 instantaneous center of mass.
+            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 instantaneous center of mass.
+            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 instantaneous center of mass.
+            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 instantaneous center of mass.
+            e_2 and e_3 axes, centered at the center of dry mass.
 
         Notes
         -----
@@ -681,23 +681,23 @@ def evaluate_dry_inertias(self):
 
     def evaluate_inertias(self):
         """Calculates and returns the rocket's inertias relative to
-        the rocket's center of mass. The inertias are saved and returned
+        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 instantaneous center of mass.
+            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 instantaneous center of mass.
+            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 instantaneous center of mass.
+            e_3 axis, centered at the center of dry mass.
 
         Notes
         -----
@@ -714,25 +714,25 @@ def evaluate_inertias(self):
         """
         # Get masses
         prop_mass = self.motor.propellant_mass  # Propellant mass as a function of time
-        dry_mass = self.dry_mass  # Constant rocket mass with motor, without propellant
 
         # Compute axes distances
-        CM_to_CDM = self.center_of_mass - self.center_of_dry_mass_position
-        CM_to_CPM = self.center_of_mass - self.center_of_propellant_position
+        CDM_to_CPM = (
+            self.center_of_dry_mass_position - self.center_of_propellant_position
+        )
 
         # Compute inertias
-        self.I_11 = parallel_axis_theorem_from_com(
-            self.dry_I_11, dry_mass, CM_to_CDM
-        ) + parallel_axis_theorem_from_com(self.motor.I_11, prop_mass, CM_to_CPM)
+        self.I_11 = self.dry_I_11 + parallel_axis_theorem_from_com(
+            self.motor.propellant_I_11, prop_mass, CDM_to_CPM
+        )
 
-        self.I_22 = parallel_axis_theorem_from_com(
-            self.dry_I_22, dry_mass, CM_to_CDM
-        ) + parallel_axis_theorem_from_com(self.motor.I_22, prop_mass, CM_to_CPM)
+        self.I_22 = self.dry_I_22 + parallel_axis_theorem_from_com(
+            self.motor.propellant_I_22, prop_mass, CDM_to_CPM
+        )
 
-        self.I_33 = self.dry_I_33 + self.motor.I_33
-        self.I_12 = self.dry_I_12 + self.motor.I_12
-        self.I_13 = self.dry_I_13 + self.motor.I_13
-        self.I_23 = self.dry_I_23 + self.motor.I_23
+        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
 
         # Return inertias
         return (
@@ -814,7 +814,7 @@ def evaluate_com_to_cdm_function(self):
 
     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
+        respect to the rocket's center of dry mass at a given time. It evaluates
         each inertia tensor component at the given time and returns a Matrix
         with the computed values.
 
@@ -844,8 +844,8 @@ def get_inertia_tensor_at_time(self, t):
 
     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
+        tensor of the rocket with respect to the rocket's center of dry 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

rocketpy/simulation/flight.py

@@ -1400,7 +1400,6 @@ def u_dot(
             * self.rocket._csys
         )
         c = self.rocket.nozzle_to_cdm
-        a = self.rocket.com_to_cdm_function.get_value_opt(t)
         nozzle_radius = self.rocket.motor.nozzle_radius
         # Prepare transformation matrix
         a11 = 1 - 2 * (e2**2 + e3**2)
@@ -1503,8 +1502,8 @@ def u_dot(
                     R1 += comp_lift_xb
                     R2 += comp_lift_yb
                     # Add to total moment
-                    M1 -= (comp_cp + a) * comp_lift_yb
-                    M2 += (comp_cp + a) * comp_lift_xb
+                    M1 -= (comp_cp) * comp_lift_yb
+                    M2 += (comp_cp) * comp_lift_xb
             # Calculates Roll Moment
             try:
                 clf_delta, cld_omega, cant_angle_rad = aero_surface.roll_parameters
@@ -1779,8 +1778,8 @@ def u_dot_generalized(
                     R1 += comp_lift_xb
                     R2 += comp_lift_yb
                     # Add to total moment
-                    M1 -= (comp_cpz + r_CM_t) * comp_lift_yb
-                    M2 += (comp_cpz + r_CM_t) * comp_lift_xb
+                    M1 -= (comp_cpz) * comp_lift_yb
+                    M2 += (comp_cpz) * comp_lift_xb
             # Calculates Roll Moment
             try:
                 clf_delta, cld_omega, cant_angle_rad = aero_surface.roll_parameters

tests/unit/test_flight.py

@@ -167,9 +167,9 @@ def test_out_of_rail_stability_margin(flight_calisto_custom_wind):
 @pytest.mark.parametrize(
     "flight_time, expected_values",
     [
-        ("t_initial", (0.17179073815516033, -0.431117, 0)),
-        ("out_of_rail_time", (0.543760, -1.364593, 0)),
-        ("apogee_time", (-0.5874848151271623, -0.7563596, 0)),
+        ("t_initial", (0.258818, -0.649515, 0)),
+        ("out_of_rail_time", (0.788918, -1.979828, 0)),
+        ("apogee_time", (-0.522394, -0.744154, 0)),
         ("t_final", (0, 0, 0)),
     ],
 )
@@ -208,7 +208,7 @@ def test_aerodynamic_moments(flight_calisto_custom_wind, flight_time, expected_v
     [
         ("t_initial", (1.6542528, 0.65918, -0.067107)),
         ("out_of_rail_time", (5.05334, 2.01364, -1.7541)),
-        ("apogee_time", (2.366258, -1.830744, -0.875342)),
+        ("apogee_time", (2.354663, -1.652953, -0.936126)),
         ("t_final", (0, 0, 159.2212)),
     ],
 )
@@ -247,7 +247,10 @@ def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_va
     [
         ("t_initial", (0, 0, 0)),
         ("out_of_rail_time", (0, 2.248727, 25.703072)),
-        ("apogee_time", (-13.204789, 15.990903, -0.000138)),
+        (
+            "apogee_time",
+            (-14.485655, 15.580647, -0.000240),
+        ),
         ("t_final", (5, 2, -5.65998)),
     ],
 )
@@ -334,10 +337,10 @@ def test_rail_buttons_forces(flight_calisto_custom_wind):
     """
     test = flight_calisto_custom_wind
     atol = 5e-3
-    assert pytest.approx(3.876749, abs=atol) == test.max_rail_button1_normal_force
-    assert pytest.approx(1.544799, abs=atol) == test.max_rail_button1_shear_force
-    assert pytest.approx(1.178420, abs=atol) == test.max_rail_button2_normal_force
-    assert pytest.approx(0.469574, abs=atol) == test.max_rail_button2_shear_force
+    assert pytest.approx(1.825283, abs=atol) == test.max_rail_button1_normal_force
+    assert pytest.approx(0.727335, abs=atol) == test.max_rail_button1_shear_force
+    assert pytest.approx(3.229578, abs=atol) == test.max_rail_button2_normal_force
+    assert pytest.approx(1.286915, abs=atol) == test.max_rail_button2_shear_force
 
 
 def test_max_values(flight_calisto_robust):

tests/unit/test_rocket.py

@@ -458,9 +458,9 @@ def test_evaluate_com_to_cdm_function(calisto):
 def test_get_inertia_tensor_at_time(calisto):
     # Expected values (for t = 0)
     # TODO: compute these values by hand or using CAD.
-    I_11 = 10.31379
-    I_22 = 10.31379
-    I_33 = 0.039942
+    I_11 = 10.516647727227216
+    I_22 = 10.516647727227216
+    I_33 = 0.0379420341586346
 
     # Set tolerance threshold
     atol = 1e-5
@@ -484,9 +484,9 @@ def test_get_inertia_tensor_at_time(calisto):
 def test_get_inertia_tensor_derivative_at_time(calisto):
     # Expected values (for t = 2s)
     # TODO: compute these values by hand or using CAD.
-    I_11_dot = -0.634805230901143
-    I_22_dot = -0.634805230901143
-    I_33_dot = -0.000671493662305
+    I_11_dot = -0.7164327431607691
+    I_22_dot = -0.7164327431607691
+    I_33_dot = -0.0006714936623050
 
     # Set tolerance threshold
     atol = 1e-3