BUG: use the angular rates in the solid propulsion equations of motion

rocketpy/simulation/flight.py

@@ -2191,7 +2191,6 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
         # Retrieve integration data
         _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u
         # Determine lift force and moment
-        omega1, omega2, omega3 = 0, 0, 0
         R1, R2, M1, M2, M3 = 0, 0, 0, 0, 0
         # Thrust correction parameters
         pressure = self.env.pressure.get_value_opt(z)

tests/integration/simulation/test_flight.py

@@ -264,6 +264,43 @@ def test_simpler_parachute_triggers(mock_show, example_plain_env, calisto_robust
     assert test_flight.all_info() is None
 
 
+def test_solid_propulsion_equations_damp_rotation(example_plain_env, calisto_robust):
+    """The solid-propulsion equations must see the angular rates.
+
+    They once zeroed the rates before the aerodynamics, which removed all
+    aerodynamic damping: a pitch rate produced no restoring moment and a
+    canted fin spun the rocket up without bound. A rolling flight must reach
+    the same roll rate with either set of equations.
+    """
+    calisto_robust.aerodynamic_surfaces.clear()  # rebuild with canted fins
+    calisto_robust.add_nose(length=0.55829, kind="vonkarman", position=1.160)
+    calisto_robust.add_tail(
+        top_radius=0.0635, bottom_radius=0.0435, length=0.060, position=-1.313
+    )
+    calisto_robust.add_trapezoidal_fins(
+        4,
+        span=0.100,
+        root_chord=0.120,
+        tip_chord=0.040,
+        position=-1.168,
+        cant_angle=1.5,
+    )
+    roll = {}
+    for equations in ("solid_propulsion", "standard"):
+        rolling = Flight(
+            rocket=calisto_robust,
+            environment=example_plain_env,
+            rail_length=5.2,
+            inclination=85,
+            heading=0,
+            equations_of_motion=equations,
+            terminate_on_apogee=True,
+        )
+        roll[equations] = np.max(np.abs(rolling.w3[:, 1]))
+    assert roll["solid_propulsion"] == pytest.approx(roll["standard"], rel=0.02)
+    assert roll["standard"] < 100  # rad/s; unbounded spin-up would give thousands
+
+
 @patch("matplotlib.pyplot.show")
 def test_rolling_flight(  # pylint: disable=unused-argument
     mock_show,
@@ -302,6 +339,18 @@ def test_rolling_flight(  # pylint: disable=unused-argument
 
     assert test_flight.all_info() is None
 
+    # The canted fins spin the rocket up and their damping balances it, so the
+    # roll rate tracks the airspeed during the coast. Without damping it grows
+    # past a thousand rad/s.
+    times = test_flight.w3[:, 0]
+    coast = (times > test_flight.rocket.motor.burn_out_time + 1) & (
+        times < test_flight.apogee_time - 3
+    )
+    roll_rate = np.abs(test_flight.w3[coast, 1])
+    assert 1.0 < roll_rate.max() < 100.0
+    roll_per_speed = roll_rate / test_flight.speed(times[coast])
+    assert roll_per_speed.std() < 0.15 * roll_per_speed.mean()
+
 
 @patch("matplotlib.pyplot.show")
 def test_eccentricity_on_flight(  # pylint: disable=unused-argument

tests/unit/simulation/test_flight.py

@@ -831,3 +831,58 @@ def test_max_acceleration_power_off_time_with_controllers(
     assert test.max_acceleration_power_off > 0, (
         "max_acceleration_power_off should be greater than zero"
     )
+
+
+# ---------------------------------------------------------------------------
+# The equations of motion themselves, evaluated on a fixed state
+# ---------------------------------------------------------------------------
+
+
+def _coasting_state(flight, seconds_after_burnout=3.0):
+    """Return ``(t, state)`` a few seconds into the coast of a flown flight."""
+    t = flight.rocket.motor.burn_out_time + seconds_after_burnout
+    return t, list(flight.get_solution_at_time(t)[1:])
+
+
+def _with_rates(state, w1=0.0, w2=0.0, w3=0.0):
+    return [*state[:10], w1, w2, w3]
+
+
+@pytest.mark.parametrize("derivative", [Flight.u_dot, Flight.u_dot_generalized])
+def test_an_angular_rate_is_aerodynamically_damped(flight_calisto_robust, derivative):
+    """A rate about any body axis must produce an angular acceleration against it.
+
+    The solid propulsion equations once zeroed the rates before the
+    aerodynamics, so no rate was ever damped: a canted fin spun the rocket up
+    without bound. This checks each axis on a coasting state, for both sets of
+    equations.
+    """
+    flight = flight_calisto_robust
+    t, state = _coasting_state(flight)
+    at_rest = derivative(flight, t, state)
+    for axis, rate in ((0, 2.0), (1, 2.0), (2, 20.0)):
+        slot = 10 + axis
+        positive = derivative(flight, t, _with_rates(state, **{f"w{axis + 1}": rate}))
+        negative = derivative(flight, t, _with_rates(state, **{f"w{axis + 1}": -rate}))
+        assert positive[slot] < at_rest[slot] < negative[slot], f"axis {axis + 1}"
+        # damping is close to linear in the rate at these magnitudes
+        assert positive[slot] - at_rest[slot] == pytest.approx(
+            -(negative[slot] - at_rest[slot]), rel=0.2
+        )
+
+
+def test_solid_propulsion_and_generalized_equations_agree_when_rotating(
+    flight_calisto_robust,
+):
+    """With a pitch rate applied, both models give the same damping moment.
+
+    They differ in how variable mass is handled, not in the aerodynamics, so
+    the angular acceleration a rate produces must agree closely.
+    """
+    flight = flight_calisto_robust
+    t, state = _coasting_state(flight)
+    pitching = _with_rates(state, w1=2.0)
+    solid = Flight.u_dot(flight, t, pitching)
+    generalized = Flight.u_dot_generalized(flight, t, pitching)
+    assert solid[10] == pytest.approx(generalized[10], rel=0.1)
+    assert solid[3:6] == pytest.approx(generalized[3:6], rel=0.05)