diff --git a/CHANGELOG.md b/CHANGELOG.md index 1cdd8479c..0d38e07b4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,8 @@ 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) @@ -39,7 +41,7 @@ Attention: The newest changes should be on top --> ### 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 diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py index 42024fb9d..9c1275294 100644 --- a/rocketpy/__init__.py +++ b/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 ( diff --git a/rocketpy/rocket/__init__.py b/rocketpy/rocket/__init__.py index 94ad5e8ee..7a5081e26 100644 --- a/rocketpy/rocket/__init__.py +++ b/rocketpy/rocket/__init__.py @@ -2,6 +2,7 @@ from rocketpy.rocket.aero_surface import ( AeroSurface, AirBrakes, + BodyTube, ControllableGenericSurface, EllipticalFin, EllipticalFins, diff --git a/rocketpy/rocket/aero_surface/__init__.py b/rocketpy/rocket/aero_surface/__init__.py index 7a6e7ac2d..eb3602b64 100644 --- a/rocketpy/rocket/aero_surface/__init__.py +++ b/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, ) diff --git a/rocketpy/rocket/aero_surface/_barrowman_surface.py b/rocketpy/rocket/aero_surface/_barrowman_surface.py index dc161b9e6..7c14f689b 100644 --- a/rocketpy/rocket/aero_surface/_barrowman_surface.py +++ b/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,6 +35,8 @@ 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 @@ -41,6 +44,13 @@ class _BarrowmanSurface(LinearGenericSurface): # 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) diff --git a/rocketpy/rocket/aero_surface/body_tube.py b/rocketpy/rocket/aero_surface/body_tube.py new file mode 100644 index 000000000..3d48a31cf --- /dev/null +++ b/rocketpy/rocket/aero_surface/body_tube.py @@ -0,0 +1,240 @@ +import numpy as np + +from rocketpy.mathutils.function import Function +from rocketpy.plots.aero_surface_plots import _NoseConePlots +from rocketpy.prints.aero_surface_prints import _NoseConePrints + +from ._barrowman_surface import _BarrowmanSurface + + +class BodyTube(_BarrowmanSurface): + """Keeps body tube information. + + A constant-radius cylindrical section of the airframe. By slender-body + theory it produces no linear normal force (``clalpha = 0``); its entire + in-flight normal force comes from the nonlinear Galejs body-lift term + ``K · (A_plan / A_ref) · sin²α`` applied at the planform centroid, exactly + like OpenRocket's ``SymmetricComponentCalc`` for a straight tube. + + Note + ---- + Local coordinate system: + - the origin at the top of the tube (the portion closest to the rocket's nose) and + - the Z axis along the longitudinal axis of symmetry, positive downwards (top -> bottom). + + Attributes + ---------- + BodyTube.length : float + Body tube length. Has units of length and must be given in meters. + BodyTube.radius : float + Body tube outer radius. Has units of length and must be given in meters. + BodyTube.rocket_radius : float + The reference rocket radius used for lift coefficient normalization, + in meters. Defaults to the tube radius. + BodyTube.name : string + Body tube name. Has no impact in simulation, as it is only used to + display data in a more organized matter. + BodyTube.cp : tuple + Tuple with the x, y and z local coordinates of the body tube center of + pressure. Has units of length and is given in meters. + BodyTube.clalpha : float + Normal-force coefficient slope. Identically zero for a constant-radius + tube. + BodyTube.plots : plots.aero_surface_plots._NoseConePlots + This contains all the plots methods. Use help(BodyTube.plots) to know + more about it. + BodyTube.prints : prints.aero_surface_prints._NoseConePrints + This contains all the prints methods. Use help(BodyTube.prints) to know + more about it. + """ + + def __init__( + self, + length, + radius, + rocket_radius=None, + name="Body Tube", + ): + """Initializes the body tube object by computing and storing the most + important values. + + Parameters + ---------- + length : float + Body tube length. Has units of length and must be given in meters. + radius : float + Body tube outer radius. Has units of length and must be given in + meters. + rocket_radius : int, float, optional + The reference rocket radius used for lift coefficient normalization. + Defaults to the tube radius when not given. + name : str, optional + Body tube name. Has no impact in simulation, as it is only used to + display data in a more organized matter. + + Returns + ------- + None + """ + rocket_radius = rocket_radius or radius + self.name = name + self.reference_area = np.pi * rocket_radius**2 + self.reference_length = 2 * rocket_radius + + self._length = length + self._radius = radius + self._rocket_radius = rocket_radius + + self.evaluate_geometrical_parameters() + self.evaluate_lift_coefficient() + self.evaluate_center_of_pressure() + self.evaluate_body_lift_geometry() + + # Translate the Barrowman geometry into the linear generic-surface + # coefficient model and build the shared compute path. + super().__init__( + reference_area=self.reference_area, + reference_length=self.reference_length, + coefficients={}, + center_of_pressure=(self.cpx, self.cpy, self.cpz), + name=name, + ) + + self.plots = _NoseConePlots(self) + self.prints = _NoseConePrints(self) + + @property + def radius(self): + return self._radius + + @radius.setter + def radius(self, value): + self._radius = value + self.evaluate_geometrical_parameters() + self.evaluate_body_lift_geometry() + + @property + def length(self): + return self._length + + @length.setter + def length(self, value): + self._length = value + self.evaluate_center_of_pressure() + self.evaluate_body_lift_geometry() + + @property + def rocket_radius(self): + return self._rocket_radius + + @rocket_radius.setter + def rocket_radius(self, value): + self._rocket_radius = value + self.reference_area = np.pi * value**2 + self.reference_length = 2 * value + + def evaluate_geometrical_parameters(self): + """Calculates and saves the body tube's surface area. + + Returns + ------- + None + """ + self.surface_area = 2 * np.pi * self.radius * self.length + self.fineness_ratio = self.length / (2 * self.radius) + + def evaluate_lift_coefficient(self): + """A constant-radius tube produces no slender-body normal force, so + its lift-curve slope is identically zero; all normal force comes from + the Galejs body-lift term. + + Returns + ------- + None + """ + self.clalpha = Function( + lambda mach: 0.0, + "Mach", + f"Normal-force coefficient derivative for {self.name}", + ) + + def evaluate_center_of_pressure(self): + """The geometric center of pressure sits at the tube midpoint in local + coordinates. With zero slender-body lift this point carries no force; + the Galejs term is applied at the planform centroid instead. + + Returns + ------- + self.cp : tuple + Tuple containing cpx, cpy, cpz. + """ + self.cpx = 0 + self.cpy = 0 + self.cpz = self.length / 2 + self.cp = (self.cpx, self.cpy, self.cpz) + return self.cp + + def evaluate_body_lift_geometry(self): + """Compute the planform (side-projection) geometry used by the Galejs + body-lift term: a rectangle of width 2·R and height L, with centroid + at the midpoint. The slender-body CP equals the same midpoint. + + Returns + ------- + None + """ + self._planform_area = 2 * self.radius * self.length + self._planform_centroid = self.length / 2 + self._cp_slender = self.cpz + + def info(self): + """Prints and plots summarized information of the body tube. + + Return + ------ + None + """ + self.prints.geometry() + self.prints.lift() + + def all_info(self): + """Prints and plots all the available information of the body tube. + + Returns + ------- + None + """ + self.prints.all() + self.plots.all() + + def to_dict(self, **kwargs): + data = { + "length": self._length, + "radius": self._radius, + "rocket_radius": self._rocket_radius, + "name": self.name, + } + + if kwargs.get("include_outputs", False): + clalpha = self.clalpha + if kwargs.get("discretize", False): + clalpha = Function(clalpha).set_discrete(0, 4, 50) + + data.update( + { + "cp": self.cp, + "clalpha": clalpha, + "surface_area": self.surface_area, + } + ) + + return data + + @classmethod + def from_dict(cls, data): + return cls( + length=data["length"], + radius=data["radius"], + rocket_radius=data["rocket_radius"], + name=data["name"], + ) diff --git a/rocketpy/rocket/aero_surface/nose_cone.py b/rocketpy/rocket/aero_surface/nose_cone.py index 7f17c159b..643a72e55 100644 --- a/rocketpy/rocket/aero_surface/nose_cone.py +++ b/rocketpy/rocket/aero_surface/nose_cone.py @@ -165,6 +165,7 @@ def __init__( # pylint: disable=too-many-statements self.evaluate_lift_coefficient() self.evaluate_center_of_pressure() + self.evaluate_body_lift_geometry() # Translate the Barrowman geometry (clalpha, cpz) into the linear # generic-surface coefficient model and build the shared compute path. @@ -200,6 +201,7 @@ def base_radius(self, value): self.evaluate_geometrical_parameters() self.evaluate_lift_coefficient() self.evaluate_nose_shape() + self.evaluate_body_lift_geometry() @property def length(self): @@ -210,6 +212,7 @@ def length(self, value): self._length = value self.evaluate_center_of_pressure() self.evaluate_nose_shape() + self.evaluate_body_lift_geometry() @property def power(self): @@ -464,6 +467,30 @@ def final_shape(x): ) self.fineness_ratio = self.length / (2 * self.base_radius) + def evaluate_body_lift_geometry(self): + """Compute the planform (side-projection) geometry used by the Galejs + body-lift term of ``_BarrowmanSurface.compute_forces_and_moments``. + + The planform area is the lateral projection of the nose contour + ``y_nosecone`` and its centroid is measured from the nose tip, in the + same convention as ``cpz``. The slender-body CP required by the CP + blend is stored as well. + + Returns + ------- + None + """ + # Numerical integration of the contour handles every nose kind, + # including the ones without a simple closed-form planform (ogive, + # tangent, lvhaack, vonkarman). + x_samples = np.linspace(0.0, self._length, 401) + y_samples = np.array([self.y_nosecone.get_value_opt(x) for x in x_samples]) + self._planform_area = float(np.trapezoid(y_samples, x_samples)) + self._planform_centroid = float( + np.trapezoid(x_samples * y_samples, x_samples) / self._planform_area + ) + self._cp_slender = self.cpz + def evaluate_lift_coefficient(self): """Calculates and returns nose cone's lift coefficient. The lift coefficient is saved and returned. This function diff --git a/rocketpy/rocket/aero_surface/tail.py b/rocketpy/rocket/aero_surface/tail.py index 7ccd2e1e3..53ba4804a 100644 --- a/rocketpy/rocket/aero_surface/tail.py +++ b/rocketpy/rocket/aero_surface/tail.py @@ -90,6 +90,7 @@ def __init__(self, top_radius, bottom_radius, length, rocket_radius, name="Tail" self.evaluate_geometrical_parameters() self.evaluate_lift_coefficient() self.evaluate_center_of_pressure() + self.evaluate_body_lift_geometry() # Translate the Barrowman geometry into the linear generic-surface # coefficient model and build the shared compute path. @@ -114,6 +115,7 @@ def top_radius(self, value): self.evaluate_geometrical_parameters() self.evaluate_lift_coefficient() self.evaluate_center_of_pressure() + self.evaluate_body_lift_geometry() @property def bottom_radius(self): @@ -125,6 +127,7 @@ def bottom_radius(self, value): self.evaluate_geometrical_parameters() self.evaluate_lift_coefficient() self.evaluate_center_of_pressure() + self.evaluate_body_lift_geometry() @property def length(self): @@ -135,6 +138,7 @@ def length(self, value): self._length = value self.evaluate_geometrical_parameters() self.evaluate_center_of_pressure() + self.evaluate_body_lift_geometry() @property def rocket_radius(self): @@ -167,6 +171,27 @@ def evaluate_shape(self): np.array([self.top_radius, self.bottom_radius]), ] + def evaluate_body_lift_geometry(self): + """Compute the planform (side-projection) geometry used by the Galejs + body-lift term of ``_BarrowmanSurface.compute_forces_and_moments``. + + The tail is a conical frustum, so its planform is a trapezoid with + parallel sides ``2 * top_radius`` and ``2 * bottom_radius``. The + centroid is measured from the top of the tail, in the same convention + as ``cpz``. The slender-body CP required by the CP blend is stored as + well. + + Returns + ------- + None + """ + self._planform_area = (self.top_radius + self.bottom_radius) * self.length + self._planform_centroid = ( + self.length / 3 * (self.top_radius + 2 * self.bottom_radius) + / (self.top_radius + self.bottom_radius) + ) + self._cp_slender = self.cpz + def evaluate_lift_coefficient(self): """Calculates and returns tail's lift coefficient. The lift coefficient is saved and returned. This function diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 6367009d8..8a8b97845 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -13,6 +13,7 @@ from rocketpy.prints.rocket_prints import _RocketPrints from rocketpy.rocket.aero_surface import ( AirBrakes, + BodyTube, EllipticalFins, Fin, Fins, @@ -1911,6 +1912,47 @@ def add_tail( self.add_surfaces(tail, position) return tail + def add_body_tube( + self, length, position, radius=None, rocket_radius=None, name="Body Tube" + ): + """Create a new constant-radius body tube, storing it as part of the + aerodynamic_surfaces list. The tube produces no slender-body normal + force; its in-flight lift comes entirely from the nonlinear Galejs + body-lift term applied at the planform centroid. + + Parameters + ---------- + length : int, float + Body tube length in meters. Must be a positive value. + position : int, float + Tube position relative to the rocket's coordinate system. By tube + position, understand the point belonging to the tube which is + highest in the rocket coordinate system (i.e. the point closest to + the nose cone). + radius : int, float, optional + Body tube outer radius in meters. If None, which is default, the + rocket radius will be used. + rocket_radius : int, float, optional + Reference radius used for lift coefficient normalization. If None, + which is default, the rocket radius will be used. + name : string + Body tube name. Default is "Body Tube". + + See Also + -------- + :ref:`addsurface` + + Returns + ------- + body_tube : BodyTube + BodyTube object created. + """ + radius = self.radius if radius is None else radius + rocket_radius = self.radius if rocket_radius is None else rocket_radius + body_tube = BodyTube(length, radius, rocket_radius, name) + self.add_surfaces(body_tube, position) + return body_tube + def add_nose( self, length, diff --git a/tests/unit/rocket/aero_surface/test_barrowman_body_lift.py b/tests/unit/rocket/aero_surface/test_barrowman_body_lift.py new file mode 100644 index 000000000..05880da9c --- /dev/null +++ b/tests/unit/rocket/aero_surface/test_barrowman_body_lift.py @@ -0,0 +1,325 @@ +"""Unit tests for the Galejs body-lift extension of ``_BarrowmanSurface``. + +The body-lift hook adds a nonlinear ``sin²α`` term (Galejs, K = 1.1) to the +classic Barrowman normal force, applied at a blended center of pressure +between the slender-body CP and the planform centroid. The tests pin down: + +- the closed-form magnitude of the body-lift contribution, +- the low-speed / high-α (apogee) damping factor ``(M / 0.05)²``, +- the force-weighted CP blend between slender-body CP and planform centroid, +- the pure-tube limit where the slender-body term vanishes, +- backward compatibility: surfaces that do not opt in are unchanged. + +Reference: Galejs, R., "Body Lift Extension to Barrowman's CP Calculation", +Apogee Rockets Newsletter, 2003; OpenRocket ``SymmetricComponentCalc``. +""" + +import numpy as np +import pytest + +from rocketpy import BodyTube, Function, NoseCone, Tail +from rocketpy.mathutils.vector_matrix import Vector +from rocketpy.rocket.aero_surface._barrowman_surface import _BarrowmanSurface + +RHO = 1.225 +SPEED = 100.0 +MACH = 0.3 + + +class _BodyLiftStub(_BarrowmanSurface): + """Minimal Barrowman surface with constant clalpha and opt-in body lift. + + Defined only for exercising ``compute_forces_and_moments`` against + closed-form expectations; mirrors how geometry-defined subclasses + (nose cones, tails, future BodyTube) populate the planform attributes. + """ + + # pylint: disable=super-init-not-called + def __init__( + self, + reference_area, + clalpha_value=2.0, + cpz_slender=0.0, + planform_area=0.0, + planform_centroid=0.0, + cp_slender=0.0, + ): + self.name = "body lift stub" + self.reference_area = reference_area + self.reference_length = 2 * (reference_area / np.pi) ** 0.5 + self.clalpha = Function(lambda mach: clalpha_value) + self.cpz = cpz_slender + self._planform_area = planform_area + self._planform_centroid = planform_centroid + self._cp_slender = cp_slender + + # Attributes expected by GenericSurface.__init__ machinery. + self._unsteady_aero = False + self.control_variables = {} + self.evaluate_coefficients() + super().__init__( + reference_area=reference_area, + reference_length=self.reference_length, + coefficients={}, + center_of_pressure=(0.0, 0.0, cpz_slender), + name=self.name, + ) + + +def _velocity(alpha): + """Stream velocity at a total angle of attack alpha, transverse along y.""" + return [0.0, SPEED * np.sin(alpha), -SPEED * np.cos(alpha)] + + +def _forces(surface, velocity, mach=MACH, cp=(0.0, 0.0, 0.0)): + """Call compute_forces_and_moments with standard test conditions.""" + speed = float(np.linalg.norm(velocity)) + return surface.compute_forces_and_moments( + Vector(velocity), speed, mach, RHO, Vector(cp), Vector([0, 0, 0]) + ) + + +def test_default_surface_has_no_body_lift(): + """Surfaces that do not set a planform area must produce exactly the + linear Barrowman force (backward compatibility).""" + surface = _BodyLiftStub(reference_area=np.pi * 0.0635**2) + + _, r2, *_ = _forces(surface, _velocity(0.15)) + + expected = 0.5 * RHO * SPEED**2 * surface.reference_area * 2.0 * 0.15 + assert r2 == pytest.approx(expected, rel=1e-12) + + +def test_body_lift_magnitude_matches_galejs_formula(): + """The total lift coefficient must equal the linear term plus + K · (A_plan/A_ref) · sin²α.""" + ref_area = np.pi * 0.0635**2 + plan_area = 2 * 0.0635 * 0.5 # tube-like planform: diameter × length + k = 1.1 + alpha = 0.35 + surface = _BodyLiftStub(reference_area=ref_area, planform_area=plan_area) + + _, r2, *_ = _forces(surface, _velocity(alpha)) + + c_linear = 2.0 * alpha + c_body = k * plan_area / ref_area * np.sin(alpha) ** 2 + expected = 0.5 * RHO * SPEED**2 * ref_area * (c_linear + c_body) + + assert r2 == pytest.approx(expected, rel=1e-12) + + +def test_low_speed_high_alpha_damping(): + """Below M = 0.05 with α > 45°, the body-lift term is damped by + (M / 0.05)².""" + ref_area = np.pi * 0.0635**2 + plan_area = 2 * 0.0635 * 0.5 + k = 1.1 + alpha = np.deg2rad(60) + mach = 0.03 + slow_speed = mach * 340.0 + velocity = [0.0, slow_speed * np.sin(alpha), -slow_speed * np.cos(alpha)] + surface = _BodyLiftStub(reference_area=ref_area, planform_area=plan_area) + + _, r2, *_ = _forces(surface, velocity, mach=mach) + + damping = (mach / 0.05) ** 2 + c_linear = 2.0 * alpha + c_body = k * plan_area / ref_area * np.sin(alpha) ** 2 * damping + expected = 0.5 * RHO * slow_speed**2 * ref_area * (c_linear + c_body) + + assert r2 == pytest.approx(expected, rel=1e-12) + + +def test_no_damping_below_45_degrees(): + """The damping factor applies only when α > 45°, even at very low Mach.""" + ref_area = np.pi * 0.0635**2 + plan_area = 2 * 0.0635 * 0.5 + k = 1.1 + alpha = np.deg2rad(30) + mach = 0.03 + slow_speed = mach * 340.0 + velocity = [0.0, slow_speed * np.sin(alpha), -slow_speed * np.cos(alpha)] + surface = _BodyLiftStub(reference_area=ref_area, planform_area=plan_area) + + _, r2, *_ = _forces(surface, velocity, mach=mach) + + c_linear = 2.0 * alpha + c_body = k * plan_area / ref_area * np.sin(alpha) ** 2 # no damping + expected = 0.5 * RHO * slow_speed**2 * ref_area * (c_linear + c_body) + + assert r2 == pytest.approx(expected, rel=1e-12) + + +def test_blended_cp_is_force_weighted_average(): + """The moment must match the force-weighted blend of the slender-body CP + and the planform centroid.""" + ref_area = np.pi * 0.0635**2 + plan_area = 2 * 0.0635 * 0.5 + k = 1.1 + cp_slender = 0.10 + centroid = 0.25 # aft of the slender CP in local (nose→tail) z + alpha = 0.30 + surface = _BodyLiftStub( + reference_area=ref_area, + planform_area=plan_area, + planform_centroid=centroid, + cp_slender=cp_slender, + ) + + # Pass the slender-body CP as the geometric application point, in the + # body frame (the surface-local z runs nose→tail, hence the flip). + forces = _forces(surface, _velocity(alpha), cp=(0.0, 0.0, -cp_slender)) + _, r2, _, m1, _, _ = forces + + c_linear = 2.0 * alpha + c_body = k * plan_area / ref_area * np.sin(alpha) ** 2 + cp_blend = (c_linear * cp_slender + c_body * centroid) / (c_linear + c_body) + # M1 = -(cp_z_body) × F2 transport about the origin. + expected_m1 = cp_blend * r2 + + assert m1 == pytest.approx(expected_m1, rel=1e-12) + + +def test_blended_cp_stays_between_the_two_contributions(): + """For any α the effective CP implied by the moment must lie between the + slender-body CP and the planform centroid.""" + ref_area = np.pi * 0.0635**2 + plan_area = 2 * 0.0635 * 0.5 + cp_slender, centroid = 0.05, 0.40 + surface = _BodyLiftStub( + reference_area=ref_area, + planform_area=plan_area, + planform_centroid=centroid, + cp_slender=cp_slender, + ) + + lo, hi = sorted([cp_slender, centroid]) + for alpha_deg in (5, 15, 30, 60, 85): + alpha = np.deg2rad(alpha_deg) + forces = _forces(surface, _velocity(alpha), cp=(0.0, 0.0, -cp_slender)) + _, r2, _, m1, _, _ = forces + + c_linear = 2.0 * alpha + c_body = 1.1 * plan_area / ref_area * np.sin(alpha) ** 2 + cp_blend = ( + c_linear * cp_slender + c_body * centroid + ) / (c_linear + c_body) + + assert lo <= cp_blend <= hi + assert m1 == pytest.approx(cp_blend * r2, rel=1e-12) + + +def test_pure_tube_no_nan_when_linear_term_vanishes(): + """A constant-radius tube has zero slender-body lift; the body-lift-only + limit must not produce NaN/inf and must apply the force at the planform + centroid.""" + cp_slender, centroid = 0.10, 0.25 + surface = _BodyLiftStub( + reference_area=np.pi * 0.0635**2, + clalpha_value=0.0, # pure tube: no linear term at all + planform_area=2 * 0.0635 * 0.5, + planform_centroid=centroid, + cp_slender=cp_slender, + ) + alpha = np.deg2rad(80) + + forces = _forces(surface, _velocity(alpha), cp=(0.0, 0.0, -cp_slender)) + assert all(np.isfinite(f) for f in forces) + + # Force is purely from the Galejs term. + _, r2, _, m1, _, _ = forces + q_s = 0.5 * RHO * SPEED**2 * surface.reference_area + c_body = 1.1 * (2 * 0.0635 * 0.5) / surface.reference_area * np.sin(alpha) ** 2 + assert r2 == pytest.approx(q_s * c_body, rel=1e-12) + # Applied entirely at the planform centroid: M1 = -(cp_z_body) × R2 + # with the blended CP landing exactly on the centroid in the body frame. + assert m1 == pytest.approx(centroid * r2, rel=1e-12) + + +def test_nose_cone_planform_matches_closed_form(): + """A real geometry subclass populates the planform attributes from its + contour: for a conical nose, A_plan = R·L/2 with centroid at 2L/3 from + the tip, and zeroing the planform recovers the legacy linear-only force.""" + nose = NoseCone(length=0.5, kind="conical", base_radius=0.05, rocket_radius=0.05) + # The planform comes from trapezoidal integration of the contour (401 + # samples), so allow its discretization error (~1e-6 relative). + assert nose._planform_area == pytest.approx( # noqa: SLF001 + 0.05 * 0.5 / 2, rel=1e-4 + ) + assert nose._planform_centroid == pytest.approx( # noqa: SLF001 + 2 * 0.5 / 3, rel=1e-4 + ) + assert nose._cp_slender == pytest.approx(nose.cpz) # noqa: SLF001 + + alpha = 0.20 + with_body = _forces(nose, _velocity(alpha)) + nose._planform_area = 0.0 # noqa: SLF001 + without_body = _forces(nose, _velocity(alpha)) + + # Body lift adds force beyond the linear term at this alpha. + assert with_body[1] > without_body[1] + # And zeroing the planform reproduces the legacy linear Barrowman force. + q_s = 0.5 * RHO * SPEED**2 * nose.reference_area + assert without_body[1] == pytest.approx( + q_s * 2 * (nose.radius_ratio**2) * alpha, rel=1e-12 + ) + + +def test_tail_planform_matches_closed_form(): + """The tail's trapezoid planform is (r_top + r_bot)·L with centroid at + L/3·(r_top + 2·r_bot)/(r_top + r_bot) from the top.""" + tail = Tail( + top_radius=0.0635, bottom_radius=0.0435, length=0.06, rocket_radius=0.0635 + ) + assert tail._planform_area == pytest.approx( # noqa: SLF001 + (0.0635 + 0.0435) * 0.06 + ) + expected_centroid = ( + 0.06 / 3 * (0.0635 + 2 * 0.0435) / (0.0635 + 0.0435) + ) + assert tail._planform_centroid == pytest.approx(expected_centroid) # noqa: SLF001 + assert tail._cp_slender == pytest.approx(tail.cpz) # noqa: SLF001 + + +def test_body_tube_geometry(): + """A constant-radius tube has a rectangular planform of width 2·R and + height L with the centroid at the midpoint, zero clalpha and midpoint CP.""" + tube = BodyTube(length=0.5, radius=0.0635) + assert tube._planform_area == pytest.approx(2 * 0.0635 * 0.5) # noqa: SLF001 + assert tube._planform_centroid == pytest.approx(0.25) # noqa: SLF001 + assert tube._cp_slender == pytest.approx(0.25) # noqa: SLF001 + assert tube.cpz == pytest.approx(0.25) + assert tube.clalpha.get_value_opt(0.3) == pytest.approx(0.0) + assert tube.surface_area == pytest.approx(2 * np.pi * 0.0635 * 0.5) + + +def test_body_tube_force_is_pure_galejs(): + """The BodyTube's in-flight normal force must equal exactly the Galejs + term K·(A_plan/A_ref)·sin²α applied at the planform centroid.""" + tube = BodyTube(length=0.5, radius=0.0635, rocket_radius=0.0635) + alpha = np.deg2rad(30) + + forces = _forces(tube, _velocity(alpha), cp=(0.0, 0.0, -tube.cpz)) + _, r2, _, m1, _, _ = forces + + q_s = 0.5 * RHO * SPEED**2 * tube.reference_area + c_body = 1.1 * (2 * 0.0635 * 0.5) / tube.reference_area * np.sin(alpha) ** 2 + assert r2 == pytest.approx(q_s * c_body, rel=1e-12) + assert m1 == pytest.approx(tube._planform_centroid * r2, rel=1e-12) # noqa: SLF001 + + +def test_rocket_add_body_tube(calisto_motorless): + """``Rocket.add_body_tube`` registers the tube as an aerodynamic surface + and it shifts the rocket's aerodynamic center aft.""" + rocket = calisto_motorless + rocket.add_nose(length=0.55829, kind="vonkarman", position=1.278) + + cp_before = rocket.aerodynamic_center.get_value_opt(0.2) + tube = rocket.add_body_tube(length=0.8, position=-0.4) + cp_after = rocket.aerodynamic_center.get_value_opt(0.2) + + assert isinstance(tube, BodyTube) + assert any(s.component is tube for s in rocket.aerodynamic_surfaces) + # The tube has zero slender-body slope, so the linear aerodynamic-center + # diagnostic is unchanged; its Galejs lift acts only in-flight. + assert cp_after == pytest.approx(cp_before)