ENH: model unbonded solid-motor grain CM shift (#340)

rocketpy/motors/solid_motor.py

@@ -196,6 +196,10 @@ class SolidMotor(Motor):
     SolidMotor.only_radial_burn : bool
         If True, grain regression is restricted to radial burn only (inner radius growth).
         Grain length remains constant throughout the burn. Default is False.
+    SolidMotor.grains_bonded : bool
+        If True (default), grain axial positions stay fixed at the assembled layout
+        (BATES bonded/glued grains). If False, a first-order packing model shifts the
+        propellant center of mass toward the nozzle as grain height regresses.
     """
 
     # pylint: disable=too-many-arguments
@@ -221,6 +225,7 @@ def __init__(
         coordinate_system_orientation="nozzle_to_combustion_chamber",
         reference_pressure=None,
         only_radial_burn=False,
+        grains_bonded=True,
     ):
         """Initialize Motor class, process thrust curve and geometrical
         parameters and store results.
@@ -323,6 +328,24 @@ class Function. Thrust units are Newtons.
             radial burn. If False, allows the grain to also burn
             axially. May be useful for axially inhibited grains or hybrid motors.
             Default is False.
+        grains_bonded : bool, optional
+            If True (default), grains keep fixed axial positions about
+            ``grains_center_of_mass_position`` (bonded / BATES-style assembly).
+            If False, grains are treated as freestanding and packed against the
+            nozzle-side (aft) face of the initial grain stack once acceleration
+            settles them. The propellant CM then moves toward the nozzle as
+            ``grain_height`` regresses:
+
+            ``CM(t) = grains_center_of_mass_position
+            - _csys * (grain_number / 2) * (grain_initial_height - grain_height(t))``.
+
+            This is a first-order inertial packing model: it does not integrate
+            grain rigid-body dynamics, friction, DEM contacts, or discontinuous
+            rattling. Inter-grain ``grain_separation`` (e.g. spacers) is kept
+            while the stack shortens from grain-height loss only. With
+            ``only_radial_burn=True``, height is constant so the CM does not
+            shift. Follow-ups may add acceleration-dependent settling or a
+            full multi-body grain dynamics model.
 
         Returns
         -------
@@ -356,6 +379,7 @@ class Function. Thrust units are Newtons.
         self.grain_outer_radius = grain_outer_radius
         self.grain_initial_inner_radius = grain_initial_inner_radius
         self.grain_initial_height = grain_initial_height
+        self.grains_bonded = grains_bonded
 
         # Grains initial geometrical parameters
         self.grain_initial_volume = (
@@ -478,11 +502,47 @@ def center_of_propellant_mass(self):
         -------
         Function
             Position of the propellant center of mass as a function of time.
+
+        Notes
+        -----
+        When ``grains_bonded`` is True, the CM stays at
+        ``grains_center_of_mass_position`` (fixed grain layout).
+
+        When ``grains_bonded`` is False, grains are packed against the aft
+        (nozzle-side) face of the initial grain stack. As grain height
+        regresses, the packed stack shortens and the CM shifts toward the
+        nozzle by ``(grain_number / 2) * (grain_initial_height -
+        grain_height(t))`` along the motor axis (signed by ``_csys``).
+        This is a first-order packing model, not a discrete-element
+        simulation of grain motion.
         """
-        time_source = self.grain_inner_radius.x_array
-        center_of_mass = np.full_like(time_source, self.grains_center_of_mass_position)
-        return np.column_stack((time_source, center_of_mass))
+        if self.grains_bonded:
+            time_source = self.grain_inner_radius.x_array
+            center_of_mass = np.full_like(
+                time_source, self.grains_center_of_mass_position
+            )
+            return np.column_stack((time_source, center_of_mass))
+
+        # First-order packing: fixed aft face, stack shortens with grain height.
+        return self.grains_center_of_mass_position - self._csys * (
+            self.grain_number / 2.0
+        ) * (self.grain_initial_height - self.grain_height)
 
+    def _grain_pitch_squared_sum(self):
+        """Return ``pitch**2 * sum(index_offsets**2)`` for parallel-axis inertia.
+
+        Bonded grains use fixed initial pitch; unbonded grains use the
+        instantaneous packed pitch ``grain_height + grain_separation``.
+        """
+        grain_number = self.grain_number
+        initial_value = (grain_number - 1) / 2.0
+        index_offsets = np.linspace(-initial_value, initial_value, grain_number)
+        sum_sq_index = float(np.sum(index_offsets**2))
+        if self.grains_bonded:
+            pitch = self.grain_initial_height + self.grain_separation
+            return (pitch**2) * sum_sq_index
+        pitch = self.grain_height + self.grain_separation
+        return (pitch**2) * sum_sq_index
     # pylint: disable=too-many-statements
     def evaluate_geometry(self):
         """Calculates grain inner radius and grain height as a function of time
@@ -728,14 +788,9 @@ def propellant_I_11(self):
             + (1 / 12) * self.grain_height**2
         )
 
-        # Calculate each grain's distance d to propellant center of mass
-        # Assuming each grain's COM are evenly spaced
-        initial_value = (grain_number - 1) / 2
-        d = np.linspace(-initial_value, initial_value, grain_number)
-        d = d * (self.grain_initial_height + self.grain_separation)
-
-        # Calculate inertia for all grains
-        I_11 = grain_number * grain_inertia11 + grain_mass * np.sum(d**2)
+        # Parallel-axis term from grain COM offsets about the propellant COM.
+        # Bonded: fixed initial pitch. Unbonded: packed pitch tracks grain_height.
+        I_11 = grain_number * grain_inertia11 + grain_mass * self._grain_pitch_squared_sum()
 
         return I_11
 
@@ -831,6 +886,7 @@ def to_dict(self, **kwargs):
                 "grain_separation": self.grain_separation,
                 "grains_center_of_mass_position": self.grains_center_of_mass_position,
                 "only_radial_burn": self.only_radial_burn,
+                "grains_bonded": self.grains_bonded,
             }
         )
 
@@ -881,4 +937,5 @@ def from_dict(cls, data):
             coordinate_system_orientation=data["coordinate_system_orientation"],
             reference_pressure=data.get("reference_pressure"),
             only_radial_burn=data.get("only_radial_burn", False),
+            grains_bonded=data.get("grains_bonded", True),
         )

tests/unit/motors/test_solidmotor.py

@@ -4,7 +4,7 @@
 import numpy as np
 import pytest
 
-from rocketpy import Function
+from rocketpy import Function, SolidMotor
 
 BURN_TIME = 3.9
 GRAIN_NUMBER = 5
@@ -266,6 +266,86 @@ def test_burn_area_asserts_extreme_values(cesaroni_m1670):
     )
 
 
+def _cesaroni_like_kwargs():
+    """Shared Cesaroni M1670-like SolidMotor constructor kwargs."""
+    return {
+        "thrust_source": "data/motors/cesaroni/Cesaroni_M1670.eng",
+        "burn_time": BURN_TIME,
+        "dry_mass": 1.815,
+        "dry_inertia": (0.125, 0.125, 0.002),
+        "center_of_dry_mass_position": 0.317,
+        "nozzle_position": 0,
+        "grain_number": GRAIN_NUMBER,
+        "grain_density": GRAIN_DENSITY,
+        "nozzle_radius": NOZZLE_RADIUS,
+        "throat_radius": THROAT_RADIUS,
+        "grain_separation": GRAIN_SEPARATION,
+        "grain_outer_radius": GRAIN_OUTER_RADIUS,
+        "grain_initial_height": GRAIN_INITIAL_HEIGHT,
+        "grains_center_of_mass_position": 0.397,
+        "grain_initial_inner_radius": GRAIN_INITIAL_INNER_RADIUS,
+        "interpolation_method": "linear",
+        "coordinate_system_orientation": "nozzle_to_combustion_chamber",
+    }
+
+
+def test_grains_bonded_default_matches_prior_cm(cesaroni_m1670):
+    """Default grains_bonded=True keeps a fixed propellant CM (prior behavior)."""
+    assert cesaroni_m1670.grains_bonded is True
+    assert np.allclose(
+        cesaroni_m1670.center_of_propellant_mass(0),
+        cesaroni_m1670.grains_center_of_mass_position,
+    )
+    assert np.allclose(
+        cesaroni_m1670.center_of_propellant_mass(2.0),
+        cesaroni_m1670.grains_center_of_mass_position,
+    )
+
+    bonded = SolidMotor(**_cesaroni_like_kwargs(), grains_bonded=True)
+    assert np.allclose(
+        bonded.center_of_propellant_mass.get_source()[:, 1],
+        cesaroni_m1670.center_of_propellant_mass.get_source()[:, 1],
+    )
+
+
+def test_grains_unbonded_shifts_cm_aft_as_height_regresses():
+    """Unbonded multi-grain motors pack aft; CM moves toward the nozzle."""
+    kwargs = _cesaroni_like_kwargs()
+    bonded = SolidMotor(**kwargs, grains_bonded=True)
+    unbonded = SolidMotor(**kwargs, grains_bonded=False)
+
+    assert unbonded.grains_bonded is False
+    # At ignition the packed and bonded layouts share the same CM.
+    assert np.allclose(
+        unbonded.center_of_propellant_mass(0),
+        bonded.center_of_propellant_mass(0),
+    )
+
+    t = 2.0
+    height = unbonded.grain_height(t)
+    # nozzle_to_combustion_chamber: _csys = +1, aft (toward nozzle) is smaller z.
+    expected_cm = kwargs["grains_center_of_mass_position"] - (
+        GRAIN_NUMBER / 2.0
+    ) * (GRAIN_INITIAL_HEIGHT - height)
+
+    assert np.allclose(unbonded.center_of_propellant_mass(t), expected_cm)
+    assert unbonded.center_of_propellant_mass(t) < bonded.center_of_propellant_mass(t)
+    # Packing also shrinks grain pitch, so transverse propellant inertia drops.
+    assert unbonded.propellant_I_11(t) < bonded.propellant_I_11(t)
+
+
+def test_grains_unbonded_roundtrip_serialization():
+    """grains_bonded persists through to_dict / from_dict."""
+    unbonded = SolidMotor(**_cesaroni_like_kwargs(), grains_bonded=False)
+    restored = SolidMotor.from_dict(unbonded.to_dict())
+    assert restored.grains_bonded is False
+    assert np.allclose(
+        restored.center_of_propellant_mass(2.0),
+        unbonded.center_of_propellant_mass(2.0),
+        atol=1e-6,
+    )
+
+
 @pytest.mark.parametrize("tuple_parametric", [(5, 3000)])
 def test_reshape_thrust_curve_asserts_resultant_thrust_curve_correct(
     cesaroni_m1670_shifted, tuple_parametric, linear_func