Feature/cluster motor structure

rocketpy/motors/cluster_motor.py

@@ -0,0 +1,228 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from rocketpy import Function
+from rocketpy.motors import Motor
+
+
+class ClusterMotor(Motor):
+    """
+    A class representing a cluster of N identical motors arranged symmetrically.
+
+    This class aggregates the physical properties (thrust, mass, inertia) of
+    multiple motors using the Parallel Axis Theorem (Huygens-Steiner theorem).
+
+    Attributes
+    ----------
+    motor : SolidMotor
+        The single motor instance used in the cluster.
+    number : int
+        The number of motors in the cluster.
+    radius : float
+        The radial distance from the rocket's central axis to the center of each motor.
+    """
+
+    def __init__(self, motor, number, radius):
+        """
+        Initialize the ClusterMotor.
+
+        Parameters
+        ----------
+        motor : SolidMotor
+            The base motor to be clustered.
+        number : int
+            Number of motors. Must be >= 2.
+        radius : float
+            Distance from center of rocket to center of motor (m).
+        """
+        self.motor = motor
+        self.number = number
+        self.radius = radius
+        dry_inertia_cluster = self._calculate_dry_inertia()
+
+        super().__init__(
+            thrust_source=motor.thrust_source,
+            nozzle_radius=motor.nozzle_radius,
+            burn_time=motor.burn_time,
+            dry_mass=motor.dry_mass * number,
+            dry_inertia=dry_inertia_cluster,
+            center_of_dry_mass_position=motor.center_of_dry_mass_position,
+            coordinate_system_orientation=motor.coordinate_system_orientation,
+            interpolation_method="linear",
+        )
+
+        self.throat_radius = motor.throat_radius
+        self.grain_number = motor.grain_number
+        self.grain_density = motor.grain_density
+        self.grain_outer_radius = motor.grain_outer_radius
+        self.grain_initial_inner_radius = motor.grain_initial_inner_radius
+        self.grain_initial_height = motor.grain_initial_height
+        self.grains_center_of_mass_position = motor.grains_center_of_mass_position
+        self._thrust = self.motor.thrust * self.number
+        self._propellant_mass = self.motor.propellant_mass * self.number
+        self._propellant_initial_mass = self.number * self.motor.propellant_initial_mass
+        self._center_of_propellant_mass = self.motor.center_of_propellant_mass
+        Ixx_term1 = self.motor.propellant_I_11 * self.number
+        Ixx_term2 = self.motor.propellant_mass * (0.5 * self.number * self.radius**2)
+        self._propellant_I_11 = Ixx_term1 + Ixx_term2
+        self._propellant_I_22 = self._propellant_I_11
+
+        Izz_term1 = self.motor.propellant_I_33 * self.number
+        Izz_term2 = self.motor.propellant_mass * (self.number * self.radius**2)
+        self._propellant_I_33 = Izz_term1 + Izz_term2
+
+        zero_func = Function(0)
+        self._propellant_I_12 = zero_func
+        self._propellant_I_13 = zero_func
+        self._propellant_I_23 = zero_func
+
+    @property
+    def thrust(self):
+        return self._thrust
+
+    @thrust.setter
+    def thrust(self, value):
+        self._thrust = value
+
+    @property
+    def propellant_mass(self):
+        return self._propellant_mass
+
+    @propellant_mass.setter
+    def propellant_mass(self, value):
+        self._propellant_mass = value
+
+    @property
+    def propellant_initial_mass(self):
+        return self._propellant_initial_mass
+
+    @propellant_initial_mass.setter
+    def propellant_initial_mass(self, value):
+        self._propellant_initial_mass = value
+
+    @property
+    def center_of_propellant_mass(self):
+        return self._center_of_propellant_mass
+
+    @center_of_propellant_mass.setter
+    def center_of_propellant_mass(self, value):
+        self._center_of_propellant_mass = value
+
+    @property
+    def propellant_I_11(self):
+        return self._propellant_I_11
+
+    @propellant_I_11.setter
+    def propellant_I_11(self, value):
+        self._propellant_I_11 = value
+
+    @property
+    def propellant_I_22(self):
+        return self._propellant_I_22
+
+    @propellant_I_22.setter
+    def propellant_I_22(self, value):
+        self._propellant_I_22 = value
+
+    @property
+    def propellant_I_33(self):
+        return self._propellant_I_33
+
+    @propellant_I_33.setter
+    def propellant_I_33(self, value):
+        self._propellant_I_33 = value
+
+    @property
+    def propellant_I_12(self):
+        return self._propellant_I_12
+
+    @propellant_I_12.setter
+    def propellant_I_12(self, value):
+        self._propellant_I_12 = value
+
+    @property
+    def propellant_I_13(self):
+        return self._propellant_I_13
+
+    @propellant_I_13.setter
+    def propellant_I_13(self, value):
+        self._propellant_I_13 = value
+
+    @property
+    def propellant_I_23(self):
+        return self._propellant_I_23
+
+    @propellant_I_23.setter
+    def propellant_I_23(self, value):
+        self._propellant_I_23 = value
+
+    def exhaust_velocity(self, t):
+        return self.motor.exhaust_velocity(t)
+
+    def _calculate_dry_inertia(self):
+        Ixx_loc = self.motor.dry_I_11
+        Iyy_loc = self.motor.dry_I_22
+        Izz_loc = self.motor.dry_I_33
+        m_dry = self.motor.dry_mass
+
+        Izz_cluster = self.number * Izz_loc + self.number * m_dry * (self.radius**2)
+        I_transverse = self.number * Ixx_loc + (self.number / 2) * m_dry * (
+            self.radius**2
+        )
+
+        return (I_transverse, I_transverse, Izz_cluster)
+
+    def info(self):
+        print(f"Cluster Configuration:")
+        print(f" - Motors: {self.number} x {type(self.motor).__name__}")
+        print(f" - Radial Distance: {self.radius} m")
+        return self.motor.info()
+
+    def draw_cluster_layout(self, rocket_radius=None):
+        fig, ax = plt.subplots(figsize=(6, 6))
+        ax.plot(0, 0, "k+", markersize=10, label="Central axis")
+        if rocket_radius:
+            rocket_tube = plt.Circle(
+                (0, 0),
+                rocket_radius,
+                color="black",
+                fill=False,
+                linestyle="--",
+                linewidth=2,
+                label="Rocket",
+            )
+            ax.add_patch(rocket_tube)
+            limit = rocket_radius * 1.2
+        else:
+            limit = self.radius * 2
+        motor_outer_radius = self.grain_outer_radius
+        angles = np.linspace(0, 2 * np.pi, self.number, endpoint=False)
+
+        for i, angle in enumerate(angles):
+            x = self.radius * np.cos(angle)
+            y = self.radius * np.sin(angle)
+            motor_circle = plt.Circle(
+                (x, y),
+                motor_outer_radius,
+                color="red",
+                alpha=0.5,
+                label="Engine" if i == 0 else "",
+            )
+            ax.add_patch(motor_circle)
+            ax.text(
+                x,
+                y,
+                str(i + 1),
+                color="white",
+                ha="center",
+                va="center",
+                fontweight="bold",
+            )
+        ax.set_aspect("equal", "box")
+        ax.set_xlim(-limit, limit)
+        ax.set_ylim(-limit, limit)
+        ax.set_xlabel("Position X (m)")
+        ax.set_ylabel("Position Y (m)")
+        ax.set_title(f"Cluster Configuration : {self.number} engines")
+        ax.grid(True, linestyle=":", alpha=0.6)
+        ax.legend(loc="upper right")
+        plt.show()

tests/integration/motors/test_cluster_motor.py

@@ -0,0 +1,115 @@
+import pytest
+import numpy as np
+from rocketpy import SolidMotor, Function
+from rocketpy.motors.cluster_motor import ClusterMotor
+
+@pytest.fixture
+def base_motor():
+    """
+    Creates a simplified SolidMotor for testing purposes.
+    Properties:
+    - Constant Thrust: 1000 N
+    - Burn time: 5 s
+    - Dry mass: 10 kg
+    - Dry Inertia: (1.0, 1.0, 0.1)
+    """
+    thrust_curve = Function(lambda t: 1000 if t < 5 else 0, "Time (s)", "Thrust (N)")
+
+    return SolidMotor(
+        thrust_source=thrust_curve,
+        burn_time=5,
+        dry_mass=10.0,
+        dry_inertia=(1.0, 1.0, 0.1), # Ixx, Iyy, Izz
+        grain_number=1,
+        grain_density=1000,
+        grain_outer_radius=0.05,
+        grain_initial_inner_radius=0.02,
+        grain_initial_height=0.5,
+        coordinate_system_orientation="nozzle_to_combustion_chamber",
+        nozzle_radius=0.02,
+        grain_separation=0.001,
+        grains_center_of_mass_position=0.25,
+        center_of_dry_mass_position=0.25
+    )
+
+def test_cluster_initialization(base_motor):
+    """
+    Tests if the ClusterMotor initializes basic attributes correctly.
+    """
+    N = 3
+    R = 0.5
+    cluster = ClusterMotor(motor=base_motor, number=N, radius=R)
+
+    assert cluster.number == N
+    assert cluster.radius == R
+    assert cluster.grain_outer_radius == base_motor.grain_outer_radius
+
+def test_cluster_mass_and_thrust_scaling(base_motor):
+    """
+    Tests if scalar properties (Thrust, Mass) are correctly multiplied by N.
+    """
+    N = 4
+    R = 0.2
+    cluster = ClusterMotor(motor=base_motor, number=N, radius=R)
+
+    # 1. Check Thrust Scaling
+    # Thrust at t=1 should be N * single_motor_thrust
+    assert np.isclose(cluster.thrust(1), base_motor.thrust(1) * N)
+
+    # 2. Check Dry Mass Scaling
+    assert np.isclose(cluster.dry_mass, base_motor.dry_mass * N)
+
+    # 3. Check Propellant Mass Scaling
+    assert np.isclose(cluster.propellant_mass(0), base_motor.propellant_mass(0) * N)
+
+def test_cluster_dry_inertia_steiner_theorem(base_motor):
+    """
+    Tests the implementation of the Parallel Axis Theorem (Huygens-Steiner)
+    for the static (dry) mass of the cluster.
+
+    Theoretical Formulas:
+    I_zz_cluster = N * I_zz_local + N * m * R^2
+    I_xx_cluster = N * I_xx_local + (N/2) * m * R^2  (Radial symmetry approximation)
+    """
+    N = 3
+    R = 1.0 # 1 meter radius for simpler checking
+    cluster = ClusterMotor(motor=base_motor, number=N, radius=R)
+
+    m_dry = base_motor.dry_mass
+    Ixx_loc = base_motor.dry_I_11
+    Izz_loc = base_motor.dry_I_33
+
+    # Expected Izz (Longitudinal / Roll)
+    expected_Izz = N * Izz_loc + N * m_dry * (R**2)
+
+    # Expected Ixx/Iyy (Transverse / Pitch / Yaw)
+    expected_I_trans = N * Ixx_loc + (N / 2) * m_dry * (R**2)
+
+    assert np.isclose(cluster.dry_I_33, expected_Izz)
+    assert np.isclose(cluster.dry_I_11, expected_I_trans)
+    assert np.isclose(cluster.dry_I_22, expected_I_trans)
+
+def test_cluster_propellant_inertia_dynamic(base_motor):
+    """
+    Tests if the Steiner theorem is correctly applied dynamically 
+    to the changing propellant mass over time.
+    """
+    N = 2
+    R = 0.5
+    cluster = ClusterMotor(motor=base_motor, number=N, radius=R)
+
+    t = 0 # Check at t=0
+
+    m_prop = base_motor.propellant_mass(t)
+    Ixx_prop_loc = base_motor.propellant_I_11(t)
+    Izz_prop_loc = base_motor.propellant_I_33(t)
+
+    # Expected Dynamic Ixx
+    # Ixx_term1 (Local rotation) + Ixx_term2 (Parallel axis offset)
+    expected_Ixx = (Ixx_prop_loc * N) + (m_prop * 0.5 * N * R**2)
+
+    # Expected Dynamic Izz
+    expected_Izz = (Izz_prop_loc * N) + (m_prop * N * R**2)
+
+    assert np.isclose(cluster.propellant_I_11(t), expected_Ixx)
+    assert np.isclose(cluster.propellant_I_33(t), expected_Izz)