Merge pull request #609 from RocketPy-Team/enh/add-prandtl-nosecone-tail

ENH: Add Prandtl-Glauert transformation to NoseCone and Tail

CHANGELOG.md

@@ -44,6 +44,7 @@ You can install this version by running `pip install rocketpy==1.3.0`
 
 ### Added
 
+- ENH: Add Prandtl-Gauss transformation to NoseCone and Tail [#609](https://github.com/RocketPy-Team/RocketPy/pull/609)
 - DOC: Adds prometheus data, Spaceport America 2022 [#601](https://github.com/RocketPy-Team/RocketPy/pull/601)
 - ENH: Pre-calculate attributes in Rocket class [#595](https://github.com/RocketPy-Team/RocketPy/pull/595)
 - ENH: Complex step differentiation [#594](https://github.com/RocketPy-Team/RocketPy/pull/594)

rocketpy/rocket/aero_surface.py

@@ -34,6 +34,37 @@ def __init__(self, name):
         self.name = name
         return None
 
+    # Defines beta parameter
+    @staticmethod
+    def _beta(mach):
+        """Defines a parameter that is often used in aerodynamic
+        equations. It is commonly used in the Prandtl factor which
+        corrects subsonic force coefficients for compressible flow.
+        This is applied to the lift coefficient of the nose cone,
+        fins and tails/transitions as in [1].
+
+        Parameters
+        ----------
+        mach : int, float
+            Number of mach.
+
+        Returns
+        -------
+        beta : int, float
+            Value that characterizes flow speed based on the mach number.
+
+        References
+        ----------
+        [1] Barrowman, James S. https://arc.aiaa.org/doi/10.2514/6.1979-504
+        """
+
+        if mach < 0.8:
+            return np.sqrt(1 - mach**2)
+        elif mach < 1.1:
+            return np.sqrt(1 - 0.8**2)
+        else:
+            return np.sqrt(mach**2 - 1)
+
     @abstractmethod
     def evaluate_center_of_pressure(self):
         """Evaluates the center of pressure of the aerodynamic surface in local
@@ -465,7 +496,7 @@ def evaluate_lift_coefficient(self):
         # It must be set as a Function because it will be called and treated
         # as a function of mach in the simulation.
         self.clalpha = Function(
-            lambda mach: 2 * self.radius_ratio**2,
+            lambda mach: 2 / self._beta(mach) * self.radius_ratio**2,
             "Mach",
             f"Lift coefficient derivative for {self.name}",
         )
@@ -781,7 +812,7 @@ def evaluate_lift_coefficient(self):
                 clalpha2D_incompressible *= 180 / np.pi
 
         # Correcting for compressible flow (apply Prandtl-Glauert correction)
-        clalpha2D = Function(lambda mach: clalpha2D_incompressible / self.__beta(mach))
+        clalpha2D = Function(lambda mach: clalpha2D_incompressible / self._beta(mach))
 
         # Diederich's Planform Correlation Parameter
         FD = 2 * np.pi * self.AR / (clalpha2D * np.cos(self.gamma_c))
@@ -857,30 +888,6 @@ def evaluate_roll_parameters(self):
         self.roll_parameters = [clf_delta, cld_omega, self.cant_angle_rad]
         return self.roll_parameters
 
-    # Defines beta parameter
-    def __beta(_, mach):
-        """Defines a parameter that is often used in aerodynamic
-        equations. It is commonly used in the Prandtl factor which
-        corrects subsonic force coefficients for compressible flow.
-
-        Parameters
-        ----------
-        mach : int, float
-            Number of mach.
-
-        Returns
-        -------
-        beta : int, float
-            Value that characterizes flow speed based on the mach number.
-        """
-
-        if mach < 0.8:
-            return np.sqrt(1 - mach**2)
-        elif mach < 1.1:
-            return np.sqrt(1 - 0.8**2)
-        else:
-            return np.sqrt(mach**2 - 1)
-
     # Defines number of fins  factor
     def __fin_num_correction(_, n):
         """Calculates a correction factor for the lift coefficient of multiple
@@ -1748,6 +1755,7 @@ def evaluate_lift_coefficient(self):
         # as a function of mach in the simulation.
         self.clalpha = Function(
             lambda mach: 2
+            / self._beta(mach)
             * (
                 (self.bottom_radius / self.rocket_radius) ** 2
                 - (self.top_radius / self.rocket_radius) ** 2

tests/test_flight.py

@@ -640,10 +640,10 @@ def test_rail_buttons_forces(flight_calisto_custom_wind):
     """
     test = flight_calisto_custom_wind
     atol = 5e-3
-    assert pytest.approx(3.80358, abs=atol) == test.max_rail_button1_normal_force
-    assert pytest.approx(1.63602, abs=atol) == test.max_rail_button1_shear_force
-    assert pytest.approx(1.19353, abs=atol) == test.max_rail_button2_normal_force
-    assert pytest.approx(0.51337, abs=atol) == test.max_rail_button2_shear_force
+    assert pytest.approx(3.833613, abs=atol) == test.max_rail_button1_normal_force
+    assert pytest.approx(1.648938, abs=atol) == test.max_rail_button1_shear_force
+    assert pytest.approx(1.165307, abs=atol) == test.max_rail_button2_normal_force
+    assert pytest.approx(0.501229, abs=atol) == test.max_rail_button2_shear_force
 
 
 @pytest.mark.parametrize(
@@ -690,7 +690,7 @@ def test_accelerations(flight_calisto_custom_wind, flight_time, expected_values)
     [
         ("t_initial", (0, 0, 0)),
         ("out_of_rail_time", (0, 2.248727, 25.703072)),
-        ("apogee_time", (-13.209436, 16.05115, -0.000257)),
+        ("apogee_time", (-13.204789, 15.990903, -0.000138)),
         ("t_final", (5, 2, -5.65998)),
     ],
 )
@@ -728,7 +728,7 @@ def test_velocities(flight_calisto_custom_wind, flight_time, expected_values):
     [
         ("t_initial", (1.6542528, 0.65918, -0.067107)),
         ("out_of_rail_time", (5.05334, 2.01364, -1.7541)),
-        ("apogee_time", (2.35291, -1.8275, -0.87851)),
+        ("apogee_time", (2.366258, -1.830744, -0.875342)),
         ("t_final", (0, 0, 159.2212)),
     ],
 )
@@ -766,7 +766,7 @@ def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_va
     "flight_time, expected_values",
     [
         ("t_initial", (0.17179073815516033, -0.431117, 0)),
-        ("out_of_rail_time", (0.547026, -1.3727895, 0)),
+        ("out_of_rail_time", (0.543760, -1.364593, 0)),
         ("apogee_time", (-0.5874848151271623, -0.7563596, 0)),
         ("t_final", (0, 0, 0)),
     ],