BUG: Add reference area factor correction to aero surfaces (solves #557)

CHANGELOG.md

@@ -39,6 +39,14 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 ### Fixed
 
 
+## [v1.2.1] - 2024-02-22
+
+You can install this version by running `pip install rocketpy==1.2.1`
+
+### Fixed
+
+- BUG: Add reference area factor correction to aero surfaces (solves #557) [#558](https://github.com/RocketPy-Team/RocketPy/pull/558)
+
 ## [v1.2.0] - 2024-02-12
 
 You can install this version by running `pip install rocketpy==1.2.0`

rocketpy/prints/rocket_prints.py

@@ -114,9 +114,13 @@ def rocket_aerodynamics_quantities(self):
         print("\nAerodynamics Lift Coefficient Derivatives\n")
         for surface, position in self.rocket.aerodynamic_surfaces:
             name = surface.name
+            # ref_factor corrects lift for different reference areas
+            ref_factor = (surface.rocket_radius / self.rocket.radius) ** 2
             print(
                 name
-                + " Lift Coefficient Derivative: {:.3f}".format(surface.clalpha(0))
+                + " Lift Coefficient Derivative: {:.3f}".format(
+                    ref_factor * surface.clalpha(0)
+                )
                 + "/rad"
             )
 
@@ -135,6 +139,9 @@ def rocket_aerodynamics_quantities(self):
         print(
             f"Center of Mass position (time=0): {self.rocket.center_of_mass(0):.3f} m"
         )
+        print(
+            f"Center of Pressure position (time=0): {self.rocket.cp_position(0):.3f} m"
+        )
         print(
             "Initial Static Margin (mach=0, time=0): "
             + "{:.3f}".format(self.rocket.static_margin(0))

rocketpy/rocket/aero_surface.py

@@ -349,12 +349,20 @@
         # If base radius is not given, the ratio between base radius and
         # rocket radius is assumed as 1, meaning that the nose cone has the
         # same radius as the rocket
-        if self.base_radius is None or self.rocket_radius is None:
+        if self.base_radius is None and self.rocket_radius is not None:
             self.radius_ratio = 1
+            self.base_radius = self.rocket_radius
+        elif self.base_radius is not None and self.rocket_radius is None:
+            self.radius_ratio = 1
+            self.rocket_radius = self.base_radius
         # If base radius is given, the ratio between base radius and rocket
         # radius is calculated
-        else:
+        elif self.base_radius is not None and self.rocket_radius is not None:
             self.radius_ratio = self.base_radius / self.rocket_radius
+        else:
+            raise ValueError(
+                "Either base radius or rocket radius must be given to calculate the nose cone radius ratio."
+            )
 
         self.fineness_ratio = self.length / (2 * self.base_radius)
         return None

rocketpy/rocket/rocket.py

@@ -517,9 +517,13 @@ def evaluate_center_of_pressure(self):
         # Calculate total lift coefficient derivative and center of pressure
         if len(self.aerodynamic_surfaces) > 0:
             for aero_surface, position in self.aerodynamic_surfaces:
-                self.total_lift_coeff_der += aero_surface.clalpha
-                self.cp_position += aero_surface.clalpha * (
-                    position - self._csys * aero_surface.cpz
+                # ref_factor corrects lift for different reference areas
+                ref_factor = (aero_surface.rocket_radius / self.radius) ** 2
+                self.total_lift_coeff_der += ref_factor * aero_surface.clalpha
+                self.cp_position += (
+                    ref_factor
+                    * aero_surface.clalpha
+                    * (position - self._csys * aero_surface.cpz)
                 )
             self.cp_position /= self.total_lift_coeff_der
 
@@ -871,7 +875,9 @@ def add_tail(
         self.add_surfaces(tail, position)
         return tail
 
-    def add_nose(self, length, kind, position, bluffness=0, name="Nose Cone"):
+    def add_nose(
+        self, length, kind, position, bluffness=0, name="Nose Cone", base_radius=None
+    ):
         """Creates a nose cone, storing its parameters as part of the
         aerodynamic_surfaces list. Its parameters are the axial position
         along the rocket and its derivative of the coefficient of lift
@@ -894,6 +900,9 @@ def add_nose(self, length, kind, position, bluffness=0, name="Nose Cone"):
             the radius of the base of the ogive.
         name : string
             Nose cone name. Default is "Nose Cone".
+        base_radius : int, float, optional
+            Nose cone base radius in meters. If not given, the rocket radius
+            will be used.
 
         See Also
         --------
@@ -907,8 +916,8 @@ def add_nose(self, length, kind, position, bluffness=0, name="Nose Cone"):
         nose = NoseCone(
             length=length,
             kind=kind,
-            base_radius=self.radius,
-            rocket_radius=self.radius,
+            base_radius=base_radius or self.radius,
+            rocket_radius=base_radius or self.radius,
             bluffness=bluffness,
             name=name,
         )
@@ -983,8 +992,9 @@ def add_trapezoidal_fins(
             with its base perpendicular to the rocket's axis. Cannot be used in
             conjunction with sweep_length.
         radius : int, float, optional
-            Reference radius to calculate lift coefficient. If None, which is
-            default, use rocket radius.
+            Reference fuselage radius where the fins are located. This is used
+            to calculate lift coefficient and to draw the rocket. If None,
+            which is default, the rocket radius will be used.
         airfoil : tuple, optional
             Default is null, in which case fins will be treated as flat plates.
             Otherwise, if tuple, fins will be considered as airfoils. The
@@ -1064,8 +1074,9 @@ def add_elliptical_fins(
             Fins cant angle with respect to the rocket centerline. Must be given
             in degrees.
         radius : int, float, optional
-            Reference radius to calculate lift coefficient. If None, which
-            is default, use rocket radius.
+            Reference fuselage radius where the fins are located. This is used
+            to calculate lift coefficient and to draw the rocket. If None,
+            which is default, the rocket radius will be used.
         airfoil : tuple, optional
             Default is null, in which case fins will be treated as flat plates.
             Otherwise, if tuple, fins will be considered as airfoils. The

tests/test_rocket.py

@@ -1,8 +1,10 @@
 from unittest.mock import patch
 
 import numpy as np
+import pytest
 
 from rocketpy import Rocket, SolidMotor
+from rocketpy.rocket import NoseCone
 
 
 @patch("matplotlib.pyplot.show")
@@ -205,3 +207,76 @@ def test_air_brakes_clamp_off(mock_show, calisto_air_brakes_clamp_off):
     assert air_brakes_clamp_off.deployment_level == 0
 
     assert air_brakes_clamp_off.all_info() == None
+
+
+def test_add_surfaces_different_noses(calisto):
+    """Test the add_surfaces method with different nose cone configurations.
+    More specifically, this will check the static margin of the rocket with
+    different nose cone configurations.
+
+    Parameters
+    ----------
+    calisto : Rocket
+        Pytest fixture for the calisto rocket.
+    """
+    length = 0.55829
+    kind = "vonkarman"
+    position = 1.16
+    bluffness = 0
+    base_radius = 0.0635
+    rocket_radius = 0.0635
+
+    # Case 1: base_radius == rocket_radius
+    nose1 = NoseCone(
+        length,
+        kind,
+        base_radius=base_radius,
+        bluffness=bluffness,
+        rocket_radius=rocket_radius,
+        name="Nose Cone 1",
+    )
+    calisto.add_surfaces(nose1, position)
+    assert nose1.radius_ratio == pytest.approx(1, 1e-8)
+    assert calisto.static_margin(0) == pytest.approx(-8.9053, 0.01)
+
+    # Case 2: base_radius == rocket_radius / 2
+    calisto.aerodynamic_surfaces.remove(nose1)
+    nose2 = NoseCone(
+        length,
+        kind,
+        base_radius=base_radius / 2,
+        bluffness=bluffness,
+        rocket_radius=rocket_radius,
+        name="Nose Cone 2",
+    )
+    calisto.add_surfaces(nose2, position)
+    assert nose2.radius_ratio == pytest.approx(0.5, 1e-8)
+    assert calisto.static_margin(0) == pytest.approx(-8.9053, 0.01)
+
+    # Case 3: base_radius == None
+    calisto.aerodynamic_surfaces.remove(nose2)
+    nose3 = NoseCone(
+        length,
+        kind,
+        base_radius=None,
+        bluffness=bluffness,
+        rocket_radius=rocket_radius * 2,
+        name="Nose Cone 3",
+    )
+    calisto.add_surfaces(nose3, position)
+    assert nose3.radius_ratio == pytest.approx(1, 1e-8)
+    assert calisto.static_margin(0) == pytest.approx(-8.9053, 0.01)
+
+    # Case 4: rocket_radius == None
+    calisto.aerodynamic_surfaces.remove(nose3)
+    nose4 = NoseCone(
+        length,
+        kind,
+        base_radius=base_radius,
+        bluffness=bluffness,
+        rocket_radius=None,
+        name="Nose Cone 4",
+    )
+    calisto.add_surfaces(nose4, position)
+    assert nose4.radius_ratio == pytest.approx(1, 1e-8)
+    assert calisto.static_margin(0) == pytest.approx(-8.9053, 0.01)