Remove the extra factor of 2 from the fin flutter equation

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")

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):