enh:Airfoil fins normal force coefficient derivative calculation.

rocketpy/Flight.py

@@ -1324,7 +1324,7 @@ def uDot(self, t, u, postProcessing=False):
                 compStreamVzBn = compStreamVzB / compStreamSpeed
                 if -1 * compStreamVzBn < 1:
                     compAttackAngle = np.arccos(-compStreamVzBn)
-                    cLift = abs(aerodynamicSurface[1](compAttackAngle))
+                    cLift = abs(aerodynamicSurface[1](compAttackAngle, freestreamMach))
                     # Component lift force magnitude
                     compLift = (
                         0.5 * rho * (compStreamSpeed ** 2) * self.rocket.area * cLift

rocketpy/Rocket.py

@@ -370,11 +370,13 @@ def evaluateStaticMargin(self):
         # Calculate total lift coefficient derivative and center of pressure
         if len(self.aerodynamicSurfaces) > 0:
             for aerodynamicSurface in self.aerodynamicSurfaces:
-                self.totalLiftCoeffDer += aerodynamicSurface[1].differentiate(
-                    x=1e-2, dx=1e-3
-                )
+                self.totalLiftCoeffDer += Function(
+                    lambda alpha: aerodynamicSurface[1](alpha, 0)
+                ).differentiate(x=1e-2, dx=1e-3)
                 self.cpPosition += (
-                    aerodynamicSurface[1].differentiate(x=1e-2, dx=1e-3)
+                    Function(
+                        lambda alpha: aerodynamicSurface[1](alpha, 0)
+                    ).differentiate(x=1e-2, dx=1e-3)
                     * aerodynamicSurface[0][2]
                 )
             self.cpPosition /= self.totalLiftCoeffDer
@@ -436,11 +438,9 @@ def addTail(self, topRadius, bottomRadius, length, distanceToCM):
         # Calculate clalpha
         clalpha = -2 * (1 - r ** (-2)) * (topRadius / rref) ** 2
         cldata = Function(
-            lambda x: clalpha * x,
-            "Alpha (rad)",
+            lambda alpha, mach: clalpha * alpha,
+            ["Alpha (rad)", "Mach"],
             "Cl",
-            interpolation="linear",
-            extrapolation="natural",
         )
 
         # Store values as new aerodynamic surface
@@ -500,11 +500,9 @@ def addNose(self, length, kind, distanceToCM):
         # Calculate clalpha
         clalpha = 2
         cldata = Function(
-            lambda x: clalpha * x,
-            "Alpha (rad)",
+            lambda alpha, mach: clalpha * alpha,
+            ["Alpha (rad)", "Mach"],
             "Cl",
-            interpolation="linear",
-            extrapolation="natural",
         )
 
         # Store values
@@ -526,7 +524,7 @@ def addFins(
         distanceToCM,
         radius=0,
         cantAngle=0,
-        airfoil=None,
+        airfoil=False,
     ):
         """Create a fin set, storing its parameters as part of the
         aerodynamicSurfaces list. Its parameters are the axial position
@@ -556,11 +554,10 @@ def addFins(
         cantAngle : int, float, optional
             Fins cant angle with respect to the rocket centerline. Must
             be given in degrees.
-        airfoil : string
-            Fin's lift curve. It must be a .csv file. The .csv file shall
-            contain no headers and the first column must specify time in
-            seconds, while the second column specifies lift coefficient. Lift
-            coefficient is dimensionaless.
+        airfoil : bool, optional
+            Fin's airfoil shape. If True, generic airfoil lift
+            calculations will be performed. If False, calculations for
+            the trapezoildal shape will be perfomed
 
         Returns
         -------
@@ -576,6 +573,7 @@ def addFins(
         Yr = rootChord + tipChord
         s = span
         Af = Yr * s / 2  # fin area
+        gamac = np.arctan((Cr - Ct) / (2 * span))  # mid chord angle
         Ymac = (
             (s / 3) * (Cr + 2 * Ct) / Yr
         )  # span wise position of fin's mean aerodynamic chord
@@ -593,6 +591,25 @@ def addFins(
         self.span = s
         self.distanceRocketFins = distanceToCM
 
+        # Auxiliary functions
+
+        # Defines beta parameter
+        def beta(mach):
+            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 correction
+        def finNumCorrection(n):
+            correctorFactor = [2.37, 2.74, 2.99, 3.24]
+            if n >= 5 and n <= 8:
+                return correctorFactor[n - 5]
+            else:
+                return n / 2
+
         # Calculate cp position relative to cm
         if distanceToCM < 0:
             cpz = distanceToCM - (
@@ -605,102 +622,66 @@ def addFins(
                 + (1 / 6) * (Cr + Ct - Cr * Ct / (Cr + Ct))
             )
 
-        # Calculate lift parameters for planar fins
-        if not airfoil:
-            # Calculate clalpha
-            clalpha = (4 * n * (s / d) ** 2) / (1 + np.sqrt(1 + (2 * Lf / Yr) ** 2))
-            clalpha *= 1 + radius / (s + radius)
+        if airfoil:  # Calculate lift parameters for generic airfoil.
+            # Documented at https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Fins_Lift_Coefficient.pdf
 
-            # # Create a function of lift values by attack angle
-            cldata = Function(
-                lambda x: clalpha * x, "Alpha (rad)", "Cl", interpolation="linear"
-            )
-            # Parameters for Roll Moment. Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf
-            clfDelta = n * (Ymac + radius) * clalpha / d
-            cldOmega = (
-                n * clalpha * np.cos(cantAngleRad) * trapezoidalConstant / (Af * d)
-            )
-            rollParameters = (
-                [clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0]
-            )
+            # Fin–body interference correction
+            const = 1 + radius / (s + radius)
 
-            # Store values
-            fin = [(0, 0, cpz), cldata, rollParameters, "Fins"]
-            self.aerodynamicSurfaces.append(fin)
+            # Aplies number of fins correction to lift coefficient
+            const *= finNumCorrection(n)
 
-            # Refresh static margin calculation
-            self.evaluateStaticMargin()
+            # Calculates clalpha * alpha
+            cldata = Function(
+                lambda alpha, mach: alpha
+                * const
+                * 2
+                * np.pi
+                * (s ** 2 / Af)
+                / (
+                    1 + np.sqrt(1 + ((beta(mach) * s ** 2) / (Af * np.cos(gamac))) ** 2)
+                ),
+                ["Alpha (rad)", "Mach"],
+                "Cl",
+            )
 
-            # Return self
-            return self.aerodynamicSurfaces[-1]
+            # Calculates clalpha
+            clalpha = Function(
+                lambda alpha: cldata(alpha, 0),
+            ).differentiate(x=1e-2, dx=1e-3)
 
-        else:
+        else:  # Calculate lift parameters for trapezoildal planar fins
 
-            def cnalfa1(cn):
-                """Calculates the normal force coefficient derivative of a 3D
-                airfoil for a given Cnalfa0
-
-                Parameters
-                ----------
-                cn : int
-                    Normal force coefficient derivative of a 2D airfoil.
-
-                Returns
-                -------
-                Cnalfa1 : int
-                    Normal force coefficient derivative of a 3D airfoil.
-                """
-
-                # Retrieve parameters for calculations
-                Af = (Cr + Ct) * span / 2
-                # fin area
-                AR = 2 * (span ** 2) / Af  # Aspect ratio
-                gamac = np.arctan((Cr - Ct) / (2 * span))
-                # mid chord angle
-                FD = 2 * np.pi * AR / (cn * np.cos(gamac))
-                Cnalfa1 = (
-                    cn
-                    * FD
-                    * (Af / self.area)
-                    * np.cos(gamac)
-                    / (2 + FD * (1 + (4 / FD ** 2)) ** 0.5)
-                )
-                return Cnalfa1
+            # Calculates clalpha
+            clalpha = (4 * n * (s / d) ** 2) / (1 + np.sqrt(1 + (2 * Lf / Yr) ** 2))
 
-            # Import the lift curve as a function of lift values by attack angle
-            read = genfromtxt(airfoil, delimiter=",")
+            # Fin–body interference correction
+            clalpha *= 1 + radius / (s + radius)
 
-            # Applies number of fins to lift coefficient data
-            data = [[cl[0], (n / 2) * cnalfa1(cl[1])] for cl in read]
+            # Create a function of lift values by attack angle
             cldata = Function(
-                data,
-                "Alpha (rad)",
+                lambda alpha, mach: clalpha * alpha,
+                ["Alpha (rad)", "Mach"],
                 "Cl",
-                interpolation="linear",
-                extrapolation="natural",
             )
 
-            # Takes an approximation to an angular coefficient
-            clalpha = cldata.differentiate(x=0, dx=1e-2)
-
-            # Parameters for Roll Moment. Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf
-            clfDelta = n * (Ymac + radius) * clalpha / d
-            cldOmega = (
-                n * clalpha * np.cos(cantAngleRad) * trapezoidalConstant / (Af * d)
-            )
-            rollParameters = (
-                [clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0]
-            )
+        # Parameters for Roll Moment.
+        # Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf
+        clfDelta = n * (Ymac + radius) * clalpha / d
+        cldOmega = n * clalpha * np.cos(cantAngleRad) * trapezoidalConstant / (Af * d)
+        rollParameters = (
+            [clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0]
+        )
 
-            # Store values
-            fin = [(0, 0, cpz), cldata, rollParameters, "Fins"]
-            self.aerodynamicSurfaces.append(fin)
+        # Store values
+        fin = [(0, 0, cpz), cldata, rollParameters, "Fins"]
+        self.aerodynamicSurfaces.append(fin)
 
-            # Refresh static margin calculation
-            self.evaluateStaticMargin()
+        # Refresh static margin calculation
+        self.evaluateStaticMargin()
 
-            # Return self
-            return self.aerodynamicSurfaces[-1]
+        # Return self
+        return self.aerodynamicSurfaces[-1]
 
     def addParachute(
         self, name, CdS, trigger, samplingRate=100, lag=0, noise=(0, 0, 0)
@@ -987,7 +968,9 @@ def allInfo(self):
         print("\nAerodynamics Lift Coefficient Derivatives")
         for aerodynamicSurface in self.aerodynamicSurfaces:
             name = aerodynamicSurface[-1]
-            clalpha = aerodynamicSurface[1].differentiate(x=1e-2, dx=1e-3)
+            clalpha = Function(
+                lambda alpha: aerodynamicSurface[1](alpha, 0),
+            ).differentiate(x=1e-2, dx=1e-3)
             print(
                 name + " Lift Coefficient Derivative: {:.3f}".format(clalpha) + "/rad"
             )