ENH: General Axisymmetric Aerodynamic Surfaces

rocketpy/rocket/aero_surface.py

@@ -83,6 +83,182 @@ def all_info(self):
         pass
 
 
+class AxisymmetricAeroSurface(AeroSurface):
+    """Class that holds common methods for the axisymmetric aero surfaces.
+
+    Attributes
+    ----------
+    AxisymmetricAeroSurface.c_A : Function
+        Axial force coefficient as a function of angle of attack (rad) and
+        Mach number as independent variables.
+    AxisymmetricAeroSurface.c_N : Function
+        Normal force coefficient as a function of angle of attack (rad) and
+        Mach number as independent variables.
+    AxisymmetricAeroSurface.c_m : Function
+        Pitch moment coefficient as a function of angle of attack (rad) and
+        Mach number as independent variables.
+    AxisymmetricAeroSurface.c_l : Function
+        Roll moment coefficient as a function of angle of attack (rad) and
+        Mach number as independent variables.
+    AxisymmetricAeroSurface.c_m_d : Function
+        Pitch moment damping coefficient as a function of angle of attack
+        (rad) and Mach number as independent variables.
+    AxisymmetricAeroSurface.c_l_d : Function
+        Roll moment damping coefficient as a function of angle of attack
+        (rad) and Mach number as independent variables.
+    AxisymmetricAeroSurface.origin : float
+        Origin of the axisymmetric aero surface. Has units of length and
+        must be given in meters.
+    AxisymmetricAeroSurface.name : string
+        Name of the axisymmetric aero surface. Has no impact in simulation,
+        as it is only used to display data in a more organized matter.
+    """
+
+    def __init__(
+        self,
+        c_A,
+        c_N,
+        c_m,
+        c_l,
+        c_m_d,
+        c_l_d,
+        origin=0,
+        reference_area=None,
+        reference_length=None,
+        name="Axisymmetric Aero Surface",
+    ):
+        """Initializes the axisymmetric aero surface.
+
+        Parameters
+        ----------
+        c_A : Function
+            Axial force coefficient as a function of angle of attack (rad) and
+            Mach number as independent variables.
+        c_N : Function
+            Normal force coefficient as a function of angle of attack (rad) and
+            Mach number as independent variables.
+        c_m : Function
+            Pitch moment coefficient as a function of angle of attack (rad) and
+            Mach number as independent variables.
+        c_l : Function
+            Roll moment coefficient as a function of angle of attack (rad) and
+            Mach number as independent variables.
+        c_m_d : Function
+            Pitch moment damping coefficient as a function of angle of attack
+            (rad) and Mach number as independent variables.
+        c_l_d : Function
+            Roll moment damping coefficient as a function of angle of attack
+            (rad) and Mach number as independent variables.
+        origin : Function, tuple
+            Origin of the axisymmetric aero surface as a function of angle of
+            attack (rad) and Mach number as independent variables. Has units of
+            length and must be given in meters. Default is constant 0. If a
+            tuple is given, it must be a tuple of three Functions, each one
+            corresponding to the x, y and z coordinates of the origin,
+            respectively (z aligned with the rocket axis).
+        reference_area : float, optional
+            Reference area of the surface. Has units of length squared and must
+            be given in meters squared. If not given, the reference area will
+            be assumed as 1.
+        reference_length : float, optional
+            Reference length of the surface. Has units of length and must be
+            given in meters. If not given, the reference length will be assumed
+            as 1.
+        name : string, optional
+            Name of the axisymmetric aero surface. Has no impact in simulation,
+            as it is only used to display data in a more organized matter.
+
+        Returns
+        -------
+        None
+        """
+        super().__init__(name)
+        self.c_A = c_A
+        self.c_N = c_N
+        self.c_m = c_m
+        self.c_l = c_l
+        self.c_m_d = c_m_d
+        self.c_l_d = c_l_d
+        self.reference_area = reference_area
+        self.reference_length = reference_length
+
+        if isinstance(origin, (tuple, list)):
+            self.cpz = origin
+            self.cpx = Function(lambda a, M: 0)
+            self.cpy = Function(lambda a, M: 0)
+        else:
+            self.cpx, self.cpy, self.cpz = origin
+
+        return None
+
+    @abstractmethod
+    def evaluate_center_of_pressure(self):
+        """Evaluates the center of pressure of the aerodynamic surface in local
+        coordinates. This corresponds to the point where the pitching moment
+        coefficient is null. Assumes zero pitching rate.
+
+        Returns
+        -------
+        None
+        """
+        pass
+
+    @abstractmethod
+    def evaluate_lift_coefficient(self):
+        """Evaluates the lift coefficient curve of the aerodynamic surface.
+
+        Returns
+        -------
+        None
+        """
+        pass
+
+    @abstractmethod
+    def evaluate_geometrical_parameters(self):
+        """Evaluates the geometrical parameters of the aerodynamic surface.
+
+        Returns
+        -------
+        None
+        """
+        pass
+
+    def clalpha(self, mach):
+        """Calculates the lift coefficient slope.
+
+        Parameters
+        ----------
+        mach : float
+            Mach number.
+
+        Returns
+        -------
+        clalpha : float
+            Lift coefficient slope. Has units of 1/rad.
+        """
+        cl = Function(lambda alpha: self.c_N(alpha, mach))
+        return self.c_N.differentiate(0.01, dx=0.01)
+
+    @abstractmethod
+    def info(self):
+        """Prints and plots summarized information of the aerodynamic surface.
+
+        Returns
+        -------
+        None
+        """
+
+    @abstractmethod
+    def all_info(self):
+        """Prints and plots all the available information of the aero surface.
+
+        Returns
+        -------
+        None
+        """
+        pass
+
+
 class NoseCone(AeroSurface):
     """Keeps nose cone information.
 

rocketpy/simulation/flight.py

@@ -1712,8 +1712,12 @@ def u_dot_generalized(self, t, u, post_processing=False):
                 position - self.rocket.center_of_dry_mass_position
             ) * self.rocket._csys - aero_surface.cpz
             comp_cp = Vector([0, 0, comp_cpz])
-            surface_radius = aero_surface.rocket_radius
-            reference_area = np.pi * surface_radius**2
+            if hasattr(aero_surface, "reference_area"):
+                reference_area = aero_surface.reference_area
+                reference_length = aero_surface.reference_length
+            else:
+                reference_length = 2 * aero_surface.rocket_radius
+                reference_area = np.pi * aero_surface.rocket_radius**2
             # Component absolute velocity in body frame
             comp_vb = vB + (w ^ comp_cp)
             # Wind velocity at component altitude
@@ -1728,15 +1732,19 @@ def u_dot_generalized(self, t, u, post_processing=False):
             comp_stream_mach = (
                 comp_stream_speed / self.env.speed_of_sound.get_value_opt(z)
             )
-            # Component attack angle and lift force
+            comp_dynamic_pressure = 0.5 * rho * (comp_stream_speed**2)
+            # Component attack angle, lift force and drag force
             comp_attack_angle = 0
             comp_lift, comp_lift_xb, comp_lift_yb = 0, 0, 0
             if comp_stream_vx_b**2 + comp_stream_vy_b**2 != 0:
                 # Normalize component stream velocity in body frame
                 comp_stream_vz_bn = comp_stream_vz_b / comp_stream_speed
                 if -1 * comp_stream_vz_bn < 1:
                     comp_attack_angle = np.arccos(-comp_stream_vz_bn)
-                    c_lift = aero_surface.cl(comp_attack_angle, comp_stream_mach)
+                    if hasattr(aero_surface, "cl"):
+                        c_lift = aero_surface.cl(comp_attack_angle, comp_stream_mach)
+                    else:
+                        c_lift = aero_surface.c_N(comp_attack_angle, comp_stream_mach)
                     # Component lift force magnitude
                     comp_lift = (
                         0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift
@@ -1753,28 +1761,75 @@ def u_dot_generalized(self, t, u, post_processing=False):
                     # Add to total moment
                     M1 -= (comp_cpz + r_CM_z.get_value_opt(t)) * comp_lift_yb
                     M2 += (comp_cpz + r_CM_z.get_value_opt(t)) * comp_lift_xb
+                    # Add to total drag force
+                    R3 += (
+                        -0.5
+                        * rho
+                        * (comp_stream_speed**2)
+                        * reference_area
+                        * (aero_surface.c_A(comp_attack_angle, comp_stream_mach))
+                    )
             # Calculates Roll Moment
-            try:
+            if hasattr(aero_surface, "roll_parameters"):
                 clf_delta, cld_omega, cant_angle_rad = aero_surface.roll_parameters
                 M3f = (
-                    (1 / 2 * rho * comp_stream_speed**2)
+                    comp_dynamic_pressure
                     * reference_area
-                    * 2
-                    * surface_radius
+                    * reference_length
                     * clf_delta(comp_stream_mach)
                     * cant_angle_rad
                 )
                 M3d = (
-                    (1 / 2 * rho * comp_stream_speed)
+                    (comp_dynamic_pressure / comp_stream_speed)
                     * reference_area
-                    * (2 * surface_radius) ** 2
+                    * (reference_length) ** 2
                     * cld_omega(comp_stream_mach)
                     * omega3
                     / 2
                 )
                 M3 += M3f - M3d
-            except AttributeError:
-                pass
+            elif hasattr(aero_surface, "c_m"):
+                # Pitching moment direction
+                pitching_moment_dir = comp_stream_velocity ^ Vector([0, 0, 1])
+                pitching_moment_dir = pitching_moment_dir.unit_vector
+                # Pitching moment magnitude
+                c_m = aero_surface.c_m(comp_attack_angle, comp_stream_mach)
+                comp_pitching_moment = (
+                    comp_dynamic_pressure * reference_area * reference_length * c_m
+                )
+                # Pitch damping moment magnitude
+                c_m_d = aero_surface.c_m_d(comp_attack_angle, comp_stream_mach)
+                pitching_rate = Vector([omega1, omega2, omega3]) @ pitching_moment_dir
+                comp_pitching_moment -= (
+                    (comp_dynamic_pressure / comp_stream_speed)
+                    * reference_area
+                    * reference_length**2
+                    * c_m_d
+                    * pitching_rate
+                    / 2
+                )
+                # Add to total pitching moment
+                comp_pitching_moment *= pitching_moment_dir
+                M1 += comp_pitching_moment.x
+                M2 += comp_pitching_moment.y
+                # Roll moment magnitude
+                c_l = aero_surface.c_l(comp_attack_angle, comp_stream_mach)
+                comp_roll_moment = (
+                    comp_dynamic_pressure * reference_area * reference_length * c_l
+                )
+                # Roll damping moment magnitude
+                c_l_d = aero_surface.c_l_d(comp_attack_angle, comp_stream_mach)
+                comp_roll_moment = -(
+                    (comp_dynamic_pressure / comp_stream_speed)
+                    * reference_area
+                    * reference_length**2
+                    * c_l_d
+                    * omega3
+                    / 2
+                )
+                # Add to total roll moment
+                M3 += comp_roll_moment
+
         weightB = Kt @ Vector([0, 0, -total_mass * self.env.gravity(z)])
         T00 = total_mass * r_CM
         T03 = (