BUG: Airfoil fins normal force coefficient derivative calculation.

rocketpy/Rocket.py

@@ -564,11 +564,20 @@ 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 : 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
+            tuple's first item specefies the aifoil's lift coefficient
+            by angle of attack and must be either a .csv, .txt, ndarray
+            or callable. The .csv and .txt files must contain no headers
+            and the first column must specify the angle of attack, while
+            the second column must specify the lift coefficient. The
+            ndarray should be as [(x0, y0), (x1, y1), (x2, y2), ...]
+            where x0 is the angle of attack and y0 is the lift coefficient.
+            If callable, it should take an angle of attack as input and
+            return the lift coefficient at that angle of attack.
+            The tuple's second item is the unit of the angle of attack,
+            accepting either "radians" or "degrees".
 
         Returns
         -------
@@ -590,8 +599,8 @@ def addFins(
         Yma = (
             (s / 3) * (Cr + 2 * Ct) / Yr
         )  # span wise position of fin's mean aerodynamic chord
-        gamac = np.arctan((Cr - Ct) / (2 * span))
-        Lf = np.sqrt((rootChord / 2 - tipChord / 2) ** 2 + span**2)
+        gamac = np.arctan((Cr - Ct) / (2 * s))
+        Lf = np.sqrt((Cr / 2 - Ct / 2) ** 2 + s**2)
         radius = self.radius if radius == 0 else radius
         d = 2 * radius
         Aref = np.pi * radius**2
@@ -689,74 +698,47 @@ def finNumCorrection(n):
                 + (1 / 6) * (Cr + Ct - Cr * Ct / (Cr + Ct))
             )
 
-        # Calculate lift parameters for planar fins
         if not airfoil:
-
-            clalphaSingleFin = Function(
-                lambda mach: 2
-                * np.pi
-                * AR
-                * (Af / Aref)
-                / (2 + np.sqrt(4 + ((beta(mach) * AR) / (np.cos(gamac))) ** 2)),
+            # Defines clalpha2D as 2*pi for planar fins
+            clalpha2D = Function(lambda mach: 2 * np.pi / beta(mach))
+        else:
+            # Defines clalpha2D as the derivative of the
+            # lift coefficient curve for a specific airfoil
+            airfoilCl = Function(
+                airfoil[0],
+                interpolation="linear",
             )
 
-            clalphaMultipleFins = (
-                liftInterferenceFactor * finNumCorrection(n) * clalphaSingleFin
-            )  # Function of mach number
+            # Differentiating at x = 0 to get cl_alpha
+            clalpha2D_Mach0 = airfoilCl.differentiate(x=1e-3, dx=1e-3)
 
-            # Calculates clalpha * alpha
-            cl = Function(
-                lambda alpha, mach: alpha * clalphaMultipleFins(mach),
-                ["Alpha (rad)", "Mach"],
-                "Cl",
-            )
+            # Convert to radians if needed
+            if airfoil[1] == "degrees":
+                clalpha2D_Mach0 *= 180 / np.pi
 
-        else:
+            # Correcting for compressible flow
+            clalpha2D = Function(lambda mach: clalpha2D_Mach0 / beta(mach))
 
-            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
+        # Diederich's Planform Correlation Parameter
+        FD = 2 * np.pi * AR / (clalpha2D * np.cos(gamac))
 
-            # Import the lift curve as a function of lift values by attack angle
-            read = genfromtxt(airfoil, delimiter=",")
+        # Lift coefficient derivative for a single fin
+        clalphaSingleFin = Function(
+            lambda mach: (clalpha2D(mach) * FD(mach) * (Af / Aref) * np.cos(gamac))
+            / (2 + FD(mach) * np.sqrt(1 + (2 / FD(mach)) ** 2))
+        )
 
-            # Applies number of fins to lift coefficient data
-            data = [[cl[0], (n / 2) * cnalfa1(cl[1])] for cl in read]
-            cl = Function(
-                data,
-                "Alpha (rad)",
-                "Cl",
-                interpolation="linear",
-                extrapolation="natural",
-            )
+        # Lift coefficient derivative for a number of n fins corrected for Fin-Body interference
+        clalphaMultipleFins = (
+            liftInterferenceFactor * finNumCorrection(n) * clalphaSingleFin
+        )  # Function of mach number
 
-            # Takes an approximation to an angular coefficient
-            clalpha = cl.differentiate(x=0, dx=1e-2)
+        # Calculates clalpha * alpha
+        cl = Function(
+            lambda alpha, mach: alpha * clalphaMultipleFins(mach),
+            ["Alpha (rad)", "Mach"],
+            "Cl",
+        )
 
         # Parameters for Roll Moment.
         # Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf

tests/fixtures/airfoils/NACA0012-radians.txt

@@ -0,0 +1,59 @@
+0.0,0.0
+0.017453293,0.11
+0.034906585,0.22
+0.052359878,0.33
+0.06981317,0.44
+0.087266463,0.55
+0.104719755,0.66
+0.122173048,0.746
+0.13962634,0.8274
+0.157079633,0.8527
+0.174532925,0.1325
+0.191986218,0.1095
+0.20943951,0.1533
+0.226892803,0.203
+0.244346095,0.2546
+0.261799388,0.3082
+0.27925268,0.362
+0.296705973,0.42
+0.314159265,0.47689
+0.331612558,0.5322
+0.34906585,0.587
+0.366519143,0.6414
+0.383972435,0.6956
+0.401425728,0.7497
+0.41887902,0.8034
+0.436332313,0.8572
+0.453785606,0.9109
+0.471238898,0.9646
+0.523598776,0.915
+0.610865238,1.02
+0.698131701,1.075
+0.785398163,1.085
+0.872664626,1.04
+0.959931089,0.965
+1.047197551,0.875
+1.134464014,0.765
+1.221730476,0.65
+1.308996939,0.515
+1.396263402,0.37
+1.483529864,0.22
+1.570796327,0.07
+1.658062789,-0.07
+1.745329252,-0.22
+1.832595715,-0.37
+1.919862177,-0.51
+2.00712864,-0.625
+2.094395102,-0.735
+2.181661565,-0.84
+2.268928028,-0.91
+2.35619449,-0.945
+2.443460953,-0.945
+2.530727415,-0.91
+2.617993878,-0.85
+2.705260341,-0.74
+2.792526803,-0.66
+2.879793266,-0.675
+2.967059728,-0.085
+3.054326191,-0.069
+3.141592654,0.0

tests/fixtures/airfoils/e473-10e6-degrees.csv

@@ -0,0 +1,142 @@
+-18.500,-1.2376
+-18.250,-1.2660
+-18.000,-1.2964
+-17.750,-1.3297
+-17.500,-1.3599
+-17.250,-1.3903
+-17.000,-1.4125
+-16.750,-1.4275
+-16.500,-1.4384
+-16.250,-1.4489
+-16.000,-1.4526
+-15.750,-1.4533
+-15.500,-1.4563
+-15.250,-1.4537
+-15.000,-1.4517
+-14.750,-1.4485
+-14.500,-1.4420
+-14.250,-1.4396
+-14.000,-1.4323
+-13.750,-1.4302
+-13.500,-1.4243
+-13.250,-1.4136
+-13.000,-1.4014
+-12.750,-1.3881
+-12.500,-1.3734
+-12.250,-1.3574
+-12.000,-1.3413
+-11.750,-1.3234
+-11.500,-1.3062
+-11.250,-1.2865
+-11.000,-1.2684
+-10.750,-1.2483
+-10.500,-1.2284
+-10.250,-1.2083
+-9.750,-1.1669
+-9.500,-1.1449
+-9.250,-1.1240
+-9.000,-1.1024
+-8.750,-1.0798
+-8.500,-1.0583
+-8.250,-1.0218
+-8.000,-0.9837
+-7.750,-0.9463
+-7.500,-0.9069
+-7.250,-0.8715
+-7.000,-0.8349
+-6.750,-0.7967
+-6.500,-0.7556
+-6.250,-0.7247
+-5.500,-0.6225
+-5.250,-0.5452
+-5.000,-0.5173
+-4.750,-0.4969
+-4.500,-0.4742
+-4.250,-0.4500
+-4.000,-0.4249
+-3.750,-0.3991
+-3.500,-0.3742
+-3.250,-0.3487
+-3.000,-0.3225
+-2.750,-0.2960
+-2.500,-0.2692
+-2.000,-0.2162
+-1.750,-0.1894
+-1.500,-0.1624
+-1.250,-0.1351
+-1.000,-0.1083
+-0.750,-0.0814
+-0.500,-0.0544
+-0.250,-0.0268
+0.000,0.0000
+0.250,0.0268
+0.500,0.0544
+0.750,0.0814
+1.000,0.1082
+1.250,0.1351
+1.500,0.1624
+1.750,0.1894
+2.000,0.2162
+2.500,0.2692
+2.750,0.2960
+3.000,0.3225
+3.250,0.3487
+3.500,0.3742
+3.750,0.3991
+4.000,0.4249
+4.250,0.4500
+4.500,0.4742
+4.750,0.4968
+5.000,0.5174
+5.250,0.5452
+5.500,0.6225
+6.250,0.7247
+6.500,0.7556
+6.750,0.7968
+7.000,0.8350
+7.250,0.8718
+7.500,0.9069
+7.750,0.9464
+8.000,0.9838
+8.250,1.0225
+8.500,1.0582
+8.750,1.0798
+9.000,1.1023
+9.250,1.1239
+9.500,1.1449
+9.750,1.1669
+10.000,1.1867
+10.250,1.2083
+10.500,1.2284
+10.750,1.2483
+11.000,1.2684
+11.250,1.2865
+11.500,1.3062
+11.750,1.3235
+12.000,1.3414
+12.250,1.3574
+12.500,1.3735
+12.750,1.3879
+13.000,1.4017
+13.250,1.4137
+13.500,1.4247
+13.750,1.4305
+14.000,1.4327
+14.250,1.4391
+14.500,1.4433
+14.750,1.4492
+15.000,1.4531
+15.250,1.4547
+15.500,1.4569
+15.750,1.4546
+16.000,1.4541
+16.250,1.4502
+16.500,1.4407
+16.750,1.4295
+17.000,1.4155
+17.250,1.3928
+17.500,1.3645
+17.750,1.3317
+18.000,1.2977
+18.250,1.2672
+18.500,1.2406