From 87c4e6bbb3ef23b2435d73de5c9602da20a3b76e Mon Sep 17 00:00:00 2001 From: David Rundus <114828622+rundvd@users.noreply.github.com> Date: Sat, 26 Sep 2026 02:22:28 +0000 Subject: [PATCH] [PATCH] BUG: remove extra factor of 2 from the fin flutter equation fin_flutter_analysis used the form of Martin's flutter boundary (NACA TN 4197) given in Apogee Peak of Flight issue 291. Its constant 1.337 is Martin's 39.3 psi divided by the sea-level pressure and by 2, and the form multiplies the numerator by 2 again, so the flutter velocity came out sqrt(2) too high for every fin. Bennett identified the error in Peak of Flight issue 615 (2023). Add tests against Bennett's reference calculator for fins with the centroid at mid-root-chord, where it reduces to Martin's equation, run through fin_flutter_analysis with TrapezoidalFins, and update the recorded Calisto values. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_017KkQuxkQhvUHzSJbU9Va8a --- rocketpy/utilities.py | 9 +++-- tests/unit/test_utilities.py | 67 ++++++++++++++++++++++++++++++++---- 2 files changed, 66 insertions(+), 10 deletions(-) diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py index 97ebe3f45..3ab67927b 100644 --- a/rocketpy/utilities.py +++ b/rocketpy/utilities.py @@ -304,8 +304,8 @@ def fin_flutter_analysis( filename=None, ): """Calculate and plot the Fin Flutter velocity using the pressure profile - provided by the selected atmospheric model. It considers the Flutter - Boundary Equation that published in NACA Technical Paper 4197. + provided by the selected atmospheric model. It uses the flutter boundary + equation of Martin, published in NACA Technical Note 4197 (1958). These results are only estimates of a real problem and may not be useful for fins made from non-isotropic materials. Currently, this function works if only a single set of fins is added, @@ -390,8 +390,11 @@ def fin_flutter_analysis( def _flutter_mach_number( fin_thickness, shear_modulus, flight, root_chord, aspect_ratio, lambda_ ): + # Martin: (Vf / a)^2 = G / (39.3 A^3 / ((t/c)^3 (A + 2)) * (lambda + 1) / 2 + # * p / p0), with 39.3 in psi and p0 = 14.696 psi. 1.337 is 39.3 / p0 / 2, + # so Martin's factor of 1/2 is already included in it. flutter_mach = ( - (shear_modulus * 2 * (aspect_ratio + 2) * (fin_thickness / root_chord) ** 3) + (shear_modulus * (aspect_ratio + 2) * (fin_thickness / root_chord) ** 3) / (1.337 * (aspect_ratio**3) * (lambda_ + 1) * flight.pressure) ) ** 0.5 flutter_mach.set_title("Fin Flutter Mach Number") diff --git a/tests/unit/test_utilities.py b/tests/unit/test_utilities.py index a6ed1f3eb..1e29ea240 100644 --- a/tests/unit/test_utilities.py +++ b/tests/unit/test_utilities.py @@ -1,11 +1,12 @@ import logging import os +from types import SimpleNamespace from unittest.mock import patch import numpy as np import pytest -from rocketpy import Function, utilities +from rocketpy import Function, TrapezoidalFins, utilities @pytest.mark.parametrize( @@ -110,12 +111,64 @@ def test_fin_flutter_analysis(flight_calisto_custom_wind): see_prints=False, see_graphs=False, ) - assert np.isclose(flutter_mach(0), 1.00482, atol=5e-3) - assert np.isclose(flutter_mach(10), 1.1413572089696549, atol=5e-3) - assert np.isclose(flutter_mach(np.inf), 1.0048188594647927, atol=5e-3) - assert np.isclose(safety_factor(0), 64.78797, atol=5e-3) - assert np.isclose(safety_factor(10), 2.1948620401502072, atol=5e-3) - assert np.isclose(safety_factor(np.inf), 61.669562809629035, atol=5e-3) + assert np.isclose(flutter_mach(0), 0.7105140845699177, atol=5e-3) + assert np.isclose(flutter_mach(10), 0.8070387292674097, atol=5e-3) + assert np.isclose(flutter_mach(np.inf), 0.7105140859542808, atol=5e-3) + assert np.isclose(safety_factor(0), 45.81200340574091, atol=5e-3) + assert np.isclose(safety_factor(10), 1.551202495333104, atol=5e-3) + assert np.isclose(safety_factor(np.inf), 43.610793287739654, atol=5e-3) + + +# Fin geometry in inches, shear modulus and pressure in psi, and the flutter +# Mach number Vf / a from John K. Bennett's reference calculator (Fin Flutter +# Boundary Calculator v1.3, github.com/jkb-git/Fin-Flutter-Velocity-Calculator). +# Each fin has its centroid at mid-root-chord, where the calculator reduces to +# Martin's equation (NACA TN 4197) with its constant of 39.3 psi. +MARTIN_REFERENCE_CASES = [ + # thickness, sweep, tip chord, root chord, span, G, p, flutter Mach + (1 / 8, 0.0, 4.0, 4.0, 4.0, 380000, 14.173548105504, 0.9581278421552433), + (1 / 8, 2.5, 2.5, 7.5, 3.0, 600000, 11.474749839986, 1.2785504155701972), + (1 / 8, 3.0, 0.0, 6.0, 2.5, 600000, 12.6957004220342, 1.2509974714673138), + (3 / 16, 6.0, 4.0, 16.0, 5.0, 700000, 14.0274369504818, 0.9803746521458), + (1 / 16, 1.0, 2.0, 4.0, 3.0, 90000, 8.81559972329207, 0.24137409108946878), +] + + +@pytest.mark.parametrize( + "thickness, sweep, tip_chord, root_chord, span, shear_modulus, pressure, expected", + MARTIN_REFERENCE_CASES, +) +def test_fin_flutter_analysis_matches_martin( + thickness, sweep, tip_chord, root_chord, span, shear_modulus, pressure, expected +): + """The flutter Mach number must match Martin's equation. The form from + Apogee Peak of Flight issue 291 used previously applied Martin's factor of + 1/2 twice and returned sqrt(2) times these values.""" + inch, psi = 0.0254, 6894.757293168361 + fins = TrapezoidalFins( + n=4, + root_chord=root_chord * inch, + tip_chord=tip_chord * inch, + span=span * inch, + rocket_radius=0.05, + sweep_length=sweep * inch, + ) + flight = SimpleNamespace( + rocket=SimpleNamespace(fins=[fins]), + pressure=Function(pressure * psi), + mach_number=Function(0.8), + ) + + flutter_mach, safety_factor = utilities.fin_flutter_analysis( + fin_thickness=thickness * inch, + shear_modulus=shear_modulus * psi, + flight=flight, + see_prints=False, + see_graphs=False, + ) + + assert flutter_mach(0) == pytest.approx(expected, rel=1e-4) + assert safety_factor(0) == pytest.approx(expected / 0.8, rel=1e-4) def test_calculate_stall_wind_velocity_returns_value(flight_calisto_custom_wind):