Fixed lift coefficient, roll coefficients and roll moment calculations

rocketpy/Flight.py

@@ -1256,6 +1256,8 @@ def uDot(self, t, u, postProcessing=False):
         c = -self.rocket.distanceRocketNozzle
         a = b * Mt / M
         rN = self.rocket.motor.nozzleRadius
+        Aref = self.rocket.area
+        d = self.rocket.radius * 2
         # Prepare transformation matrix
         a11 = 1 - 2 * (e2 ** 2 + e3 ** 2)
         a12 = 2 * (e1 * e2 - e0 * e3)
@@ -1324,7 +1326,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
@@ -1342,13 +1344,25 @@ def uDot(self, t, u, postProcessing=False):
             # Calculates Roll Moment
             if aerodynamicSurface[-1] == "Fins":
                 Clfdelta, Cldomega, cantAngleRad = aerodynamicSurface[2]
-                Clf = Clfdelta * cantAngleRad
-                Cld = (
-                        Cldomega * omega3 * min(1, 1 / freestreamSpeed)
-                        if freestreamSpeed != 0
-                        else Cldomega * omega3
+                M3f = (
+                    (1 / 2 * rho * freestreamSpeed ** 2)
+                    * Aref
+                    * d
+                    * Clfdelta(freestreamMach)
+                    * cantAngleRad
+                )
+                M3d = (
+                    (
+                        (1 / 2 * rho * freestreamSpeed)
+                        * Aref
+                        * d
+                        * Cldomega(freestreamMach)
+                        * omega3
+                        * d
+                        /   2
                     )
-                M3 += Clf - Cld
+                )
+                M3 += M3f - M3d
         # Calculate derivatives
         # Angular acceleration
         alpha1 = (

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
@@ -416,7 +418,7 @@ def addTail(self, topRadius, bottomRadius, length, distanceToCM):
 
         Returns
         -------
-        cldata : Function
+        cl : Function
             Object of the Function class. Contains tail's lift data.
         self : Rocket
             Object of the Rocket class.
@@ -435,16 +437,14 @@ 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)",
+        cl = Function(
+            lambda alpha, mach: clalpha * alpha,
+            ["Alpha (rad)", "Mach"],
             "Cl",
-            interpolation="linear",
-            extrapolation="natural",
         )
 
         # Store values as new aerodynamic surface
-        tail = [(0, 0, cpz), cldata, "Tail"]
+        tail = [(0, 0, cpz), cl, "Tail"]
         self.aerodynamicSurfaces.append(tail)
 
         # Refresh static margin calculation
@@ -476,7 +476,7 @@ def addNose(self, length, kind, distanceToCM):
 
         Returns
         -------
-        cldata : Function
+        cl : Function
             Object of the Function class. Contains nose's lift data.
         self : Rocket
             Object of the Rocket class.
@@ -499,16 +499,14 @@ def addNose(self, length, kind, distanceToCM):
 
         # Calculate clalpha
         clalpha = 2
-        cldata = Function(
-            lambda x: clalpha * x,
-            "Alpha (rad)",
+        cl = Function(
+            lambda alpha, mach: clalpha * alpha,
+            ["Alpha (rad)", "Mach"],
             "Cl",
-            interpolation="linear",
-            extrapolation="natural",
         )
 
         # Store values
-        nose = [(0, 0, cpz), cldata, "Nose Cone"]
+        nose = [(0, 0, cpz), cl, "Nose Cone"]
         self.aerodynamicSurfaces.append(nose)
 
         # Refresh static margin calculation
@@ -564,7 +562,7 @@ def addFins(
 
         Returns
         -------
-        cldata : Function
+        cl : Function
             Object of the Function class. Contains fin's lift data.
         self : Rocket
             Object of the Rocket class.
@@ -579,20 +577,67 @@ def addFins(
         Ymac = (
             (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)
         radius = self.radius if radius == 0 else radius
         d = 2 * radius
+        Aref = np.pi * radius ** 2
+        AR = 2 * s ** 2 / Af
         cantAngleRad = np.radians(cantAngle)
-        trapezoidalConstant = ((Yr) / 2) * (radius ** 2) * s
-        trapezoidalConstant += ((Cr + 2 * Ct) / 3) * radius * (s ** 2)
-        trapezoidalConstant += ((Cr + 3 * Ct) / 12) * (s ** 3)
+        trapezoidalConstant = (
+            (Cr + 3 * Ct) * s ** 3
+            + 4 * (Cr + 2 * Ct) * radius * s ** 2
+            + 6 * (Cr + Ct) * s * radius ** 2
+        ) / 12
+
+        # Fin–body interference correction parameters
+        tau = (s + radius) / radius
+        λ = Ct / Cr
+        liftInterferenceFactor = 1 + 1 / tau
+        rollForcingInterferenceFactor = (1 / np.pi ** 2) * (
+            (np.pi ** 2 / 4) * ((tau + 1) ** 2 / tau ** 2)
+            + ((np.pi * (tau ** 2 + 1) ** 2) / (tau ** 2 * (tau - 1) ** 2))
+            * np.arcsin((tau ** 2 - 1) / (tau ** 2 + 1))
+            - (2 * np.pi * (tau + 1)) / (tau * (tau - 1))
+            + ((tau ** 2 + 1) ** 2)
+            / (tau ** 2 * (tau - 1) ** 2)
+            * (np.arcsin((tau ** 2 - 1) / (tau ** 2 + 1))) ** 2
+            - (4 * (tau + 1))
+            / (tau * (tau - 1))
+            * np.arcsin((tau ** 2 - 1) / (tau ** 2 + 1))
+            + (8 / (tau - 1) ** 2) * np.log((tau ** 2 + 1) / (2 * tau))
+        )
+        rollDampingInterferenceFactor = 1 + (
+            ((tau - λ) / (tau)) - ((1 - λ) / (tau - 1)) * np.log(tau)
+        ) / (
+            ((tau + 1) * (tau - λ)) / (2) - ((1 - λ) * (tau ** 3 - 1)) / (3 * (tau - 1))
+        )
 
         # Save geometric parameters for later Fin Flutter Analysis and Roll Moment Calculation
         self.rootChord = Cr
         self.tipChord = Ct
         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 - (
@@ -607,44 +652,35 @@ def addFins(
 
         # 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)
 
-            # # 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)
+            clalphaSingleFin = Function(
+                lambda mach: 2
+                * np.pi
+                * AR
+                * (Af / Aref)
+                / (2 + np.sqrt(4 + ((beta(mach) * AR) / (np.cos(gamac))) ** 2)),
             )
-            rollParameters = (
-                [clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0]
-            )
-
-            # Store values
-            fin = [(0, 0, cpz), cldata, rollParameters, "Fins"]
-            self.aerodynamicSurfaces.append(fin)
 
-            # Refresh static margin calculation
-            self.evaluateStaticMargin()
+            clalphaMultipleFins = (
+                liftInterferenceFactor * finNumCorrection(n) * clalphaSingleFin
+            )  # Function of mach number
 
-            # Return self
-            return self.aerodynamicSurfaces[-1]
+            # Calculates clalpha * alpha
+            cl = Function(
+                lambda alpha, mach: alpha * clalphaMultipleFins(mach),
+                ["Alpha (rad)", "Mach"],
+                "Cl",
+            )
 
         else:
 
             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
@@ -672,7 +708,7 @@ def cnalfa1(cn):
 
             # Applies number of fins to lift coefficient data
             data = [[cl[0], (n / 2) * cnalfa1(cl[1])] for cl in read]
-            cldata = Function(
+            cl = Function(
                 data,
                 "Alpha (rad)",
                 "Cl",
@@ -681,26 +717,35 @@ def cnalfa1(cn):
             )
 
             # 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]
-            )
+            clalpha = cl.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 = (
+            rollForcingInterferenceFactor * n * (Ymac + radius) * clalphaSingleFin / d
+        )  # Function of mach number
+        cldOmega = Function(
+            2
+            * rollDampingInterferenceFactor
+            * n
+            * clalphaSingleFin
+            * np.cos(cantAngleRad)
+            * trapezoidalConstant
+            / (Aref * d ** 2)
+        )  # Function of mach number
+        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), cl, 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 +1032,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"
             )