ENH: Galejs body lift with planform geometry (nose cone, tail) and new BodyTube surface

CHANGELOG.md

@@ -32,14 +32,16 @@ Attention: The newest changes should be on top -->
 
 ### Added
 
+- ENH: Add nonlinear Galejs body lift (sin²α, K=1.1) to Barrowman surfaces, with planform geometry for nose cones and tails, matching OpenRocket's SymmetricComponentCalc
+- ENH: Add BodyTube aerodynamic surface and `Rocket.add_body_tube`, producing normal force purely from the Galejs body-lift term
 - ENH: MNT: introduce pressure unit conversion when using forecast/reanalysis/ensemble data [#955](https://github.com/RocketPy-Team/RocketPy/pull/955)
 - ENH: Auto Populate Changelog [#919](https://github.com/RocketPy-Team/RocketPy/pull/919)
 - ENH: Adaptive Monte Carlo via Convergence Criteria [#922](https://github.com/RocketPy-Team/RocketPy/pull/922)
 - TST: Add acceptance tests for 3DOF flight simulation based on Bella Lui rocket [#914](https://github.com/RocketPy-Team/RocketPy/pull/914)
 
 ### Changed
 
-- 
+- ENH: Nose cones and tails now include the Galejs body-lift term, migrating the center of pressure aft at high angle of attack; results for near-apogee / high-α flight conditions differ from previous versions
 
 ### Fixed
 

rocketpy/__init__.py

@@ -29,6 +29,7 @@
 from .rocket import (
     AeroSurface,
     AirBrakes,
+    BodyTube,
     Components,
     ControllableGenericSurface,
     EllipticalFin,
@@ -48,6 +49,7 @@
     TrapezoidalFin,
     TrapezoidalFins,
 )
+
 from .sensitivity import SensitivityModel
 from .sensors import Accelerometer, Barometer, GnssReceiver, Gyroscope
 from .simulation import (

rocketpy/rocket/__init__.py

@@ -2,6 +2,7 @@
 from rocketpy.rocket.aero_surface import (
     AeroSurface,
     AirBrakes,
+    BodyTube,
     ControllableGenericSurface,
     EllipticalFin,
     EllipticalFins,

rocketpy/rocket/aero_surface/__init__.py

@@ -1,5 +1,6 @@
 from rocketpy.rocket.aero_surface.aero_surface import AeroSurface
 from rocketpy.rocket.aero_surface.air_brakes import AirBrakes
+from rocketpy.rocket.aero_surface.body_tube import BodyTube
 from rocketpy.rocket.aero_surface.controllable_generic_surface import (
     ControllableGenericSurface,
 )

rocketpy/rocket/aero_surface/_barrowman_surface.py

@@ -6,22 +6,23 @@
 
 
 class _BarrowmanSurface(LinearGenericSurface):
-    """Intermediate base for Barrowman-defined aerodynamic surfaces
+    r"""Intermediate base for Barrowman-defined aerodynamic surfaces
     such as nose cones, tails/transitions and fin sets.
 
     These surfaces expose a lift-curve slope ``clalpha`` (a ``Function`` of
     Mach), a geometric center of pressure ``cpz`` and, for fins, a pair of roll
     forcing/damping coefficients.
 
-    The in-flight normal force and its moment are computed with the classic
-    Barrowman method (see :meth:`compute_forces_and_moments`): the normal force
-    uses the true total angle of attack and acts at the geometric center of
-    pressure, and its moment about the center of dry mass is the geometric
-    transport (``cp ^ force``). This reproduces the formulation used in
-    RocketPy's flight-test validation. The resultant force is therefore reported
-    at the geometric center of pressure (:attr:`force_application_point`), which
-    the surface-local frame maps to the body frame through
-    :meth:`_default_surface_rotation`.
+    The in-flight normal force and its moment are computed with the Barrowman
+    method (see :meth:`compute_forces_and_moments`): the normal force uses the
+    true total angle of attack, acting at the geometric center of pressure, and
+    its moment about the center of dry mass is the geometric transport
+    (``cp ^ force``).  When the subclass provides planform geometry (see
+    :attr:`_planform_area`, :attr:`_planform_centroid`, :attr:`_cp_slender`), a
+    non-linear Galejs body-lift term :math:`K \cdot (A_\text{plan} /
+    A_\text{ref}) \cdot \sin^2\alpha` is added and the CP is blended
+    accordingly.  The resultant force is reported at the blended center of
+    pressure in the body frame through :meth:`_default_surface_rotation`.
 
     The class also derives the linear normal-force slopes ``cN_alpha`` (pitch
     plane) and ``cY_beta`` (yaw plane), which feed the stability and
@@ -34,13 +35,22 @@ class _BarrowmanSurface(LinearGenericSurface):
     center of pressure before calling ``super().__init__`` (which passes the
     geometric cp through ``center_of_pressure``), and, for fins, set
     ``self.roll_parameters = [clf_delta, cld_omega, cant_angle_rad]``.
+    Subclasses that wish to enable body lift must also set
+    :attr:`_planform_area`, :attr:`_planform_centroid` and :attr:`_cp_slender`.
     """
 
     # Geometry-defined Barrowman surfaces are axisymmetric by construction
     # (``cY_beta = -cN_alpha``, etc.), so they contribute identically to the
     # pitch and yaw planes. The individual ``Fin`` overrides this back to False.
     is_axisymmetric = True
 
+    # Galejs body-lift parameters.  Subclasses may override these to enable
+    # the nonlinear sin²α body-lift term (see :meth:`compute_forces_and_moments`).
+    _body_lift_k = 1.1          # Galejs constant K
+    _planform_area = 0.0        # projected (planform) area, m²
+    _planform_centroid = 0.0    # planform centroid local z, m
+    _cp_slender = 0.0           # slender-body CP local z, m
+
     @staticmethod
     def _beta(mach):
         """Prandtl-Glauert compressibility factor used to correct subsonic
@@ -126,16 +136,24 @@ def compute_forces_and_moments(
         omega,
         *args,  # pylint: disable=unused-argument
     ):
-        """Compute the surface's forces and moments with the classic Barrowman
-        method. Called at each simulation step.
+        """Compute the surface's forces and moments with the Barrowman method
+        plus the optional Galejs body-lift extension. Called at each
+        simulation step.
+
+        The normal force has two contributions:
+
+        1. **Slender-body linear term**:
+           ``0.5 ρ V² A_ref · clalpha(Mach) · α``
+
+        2. **Galejs body-lift term** (nonlinear, when :attr:`_planform_area`
+           > 0):
+           ``0.5 ρ V² A_ref · K · (A_plan / A_ref) · sin²α``,
+           with ``K = 1.1``.  At very low speed and high α (apogee) the
+           term is damped by a factor ``(M / 0.05)²``.
 
-        The normal force uses the true total angle of attack between the flow
-        and the body axis, ``attack_angle = arccos(-v_z / |v|)``, giving
-        ``0.5 * rho * V**2 * A_ref * clalpha(Mach) * attack_angle``. It is
-        applied perpendicular to the body axis (along the transverse flow) at the
-        geometric center of pressure, and its moment about the rocket's center of
-        dry mass is the geometric transport ``cp ^ force``. Fin sets add their
-        roll moment on top.
+        The total force is applied at the blended center of pressure of the
+        two contributions.  Fin sets (including canards) add their roll moment
+        on top.
 
         Parameters
         ----------
@@ -150,7 +168,8 @@ def compute_forces_and_moments(
         cp : Vector
             Surface center of pressure relative to the center of dry mass, in
             the body frame (the force-application point; see
-            :attr:`force_application_point`).
+            :attr:`force_application_point`).  When body lift is active this
+            is the *slender-body* CP; the blended CP is computed internally.
         omega : tuple of float
             Body angular velocity about the x, y, z axes. Only the roll
             component (``omega[2]``) is used, by fin sets.
@@ -172,18 +191,53 @@ def compute_forces_and_moments(
             stream_vzn = stream_vz / stream_speed
             if -stream_vzn < 1:
                 attack_angle = np.arccos(-stream_vzn)
-                c_lift = self.clalpha.get_value_opt(stream_mach) * attack_angle
-                lift = 0.5 * rho * stream_speed**2 * self.reference_area * c_lift
+
+                # --- Slender-body linear term ---
+                c_lift_linear = (
+                    self.clalpha.get_value_opt(stream_mach) * attack_angle
+                )
+
+                # --- Galejs body-lift term (nonlinear sin²α) ---
+                c_lift_body = 0.0
+                if self._planform_area > 0:
+                    sin2_alpha = np.sin(attack_angle) ** 2
+                    c_lift_body = (
+                        self._body_lift_k
+                        * self._planform_area
+                        / self.reference_area
+                        * sin2_alpha
+                    )
+                    # Low-speed / high-α damping (avoids apogee CP anomaly)
+                    if stream_mach < 0.05 and attack_angle > np.pi / 4:
+                        c_lift_body *= (stream_mach / 0.05) ** 2
+
+                c_lift = c_lift_linear + c_lift_body
+                lift = (
+                    0.5
+                    * rho
+                    * stream_speed**2
+                    * self.reference_area
+                    * c_lift
+                )
                 # Normal force, perpendicular to the body axis, directed along
                 # the transverse component of the flow.
                 transverse_norm = (stream_vx**2 + stream_vy**2) ** 0.5
                 R1 = lift * stream_vx / transverse_norm
                 R2 = lift * stream_vy / transverse_norm
-                # The normal force acts at the geometric center of pressure,
-                # which ``cp`` already locates relative to the center of dry
-                # mass; transport its moment from there.
+                # The total force acts at the blended center of pressure:
+                # slender-body CP + Galejs offset.
                 force = Vector([R1, R2, R3])
-                M1, M2, M3 = cp ^ force
+                if c_lift_body > 0 and c_lift != 0:
+                    # Body-frame offset between the two CPs (the
+                    # _default_surface_rotation flips the local z axis, hence
+                    # the minus sign).
+                    dz_body = -(self._planform_centroid - self._cp_slender)
+                    cp_effective = cp + Vector(
+                        [0.0, 0.0, dz_body * c_lift_body / c_lift]
+                    )
+                else:
+                    cp_effective = cp
+                M1, M2, M3 = cp_effective ^ force
 
         # Fin roll (cant forcing + rate damping); zero for non-fin surfaces.
         M3 += self._roll_moment(stream_speed, stream_mach, rho, omega)