BUG: fin_flutter_analysis overestimates the flutter velocity by sqrt(2) #1200

Describe the bug

fin_flutter_analysis uses the form of Martin's flutter boundary (NACA TN 4197, 1958) given by Howard in Apogee Peak of Flight issue 291 (2011):

flutter_mach = (
    (shear_modulus * 2 * (aspect_ratio + 2) * (fin_thickness / root_chord) ** 3)
    / (1.337 * (aspect_ratio**3) * (lambda_ + 1) * flight.pressure)
) ** 0.5

Martin's equation is

(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. The constant 1.337 is 39.3 / 14.696 / 2, so it already contains Martin's factor of 1/2, and the extra 2 * in the numerator applies it a second time. The flutter Mach number comes out sqrt(2) too high for every fin, so every reported safety factor is 41% too optimistic. John K. Bennett identified this error in Peak of Flight issue 615 (December 2023).

To Reproduce

A 4 x 4 in rectangular fin, 1/8 in thick, G = 380,000 psi, at p = 14.1735 psi:

from types import SimpleNamespace
from rocketpy import Function
from rocketpy.utilities import _flutter_mach_number

inch, psi = 0.0254, 6894.757293168361
chord = span = 4 * inch
flight = SimpleNamespace(pressure=Function(14.173548 * psi))
mach = _flutter_mach_number(0.125 * inch, 380000 * psi, flight, chord, 1.0, 1.0)
print(mach(0))  # 1.3549

Martin's equation gives 0.9581, and so does Bennett's reference calculator (v1.3). A rocket at Mach 1.2 would get a safety factor of 1.13 from RocketPy against 0.80 from Martin's equation.

Expected behavior

The same expression without the extra 2 *. The ratio to Bennett's calculator is 1.41416 for every fin I checked.

Additional context

Peak of Flight issue 615 also proposes a second change that is a modeling choice rather than a bug, so I have left it out of the fix. Martin defines epsilon as the distance of a section's center of gravity behind its quarter-chord and assumes 0.25, which holds for a flat plate of uniform thickness. Bennett replaces it with the planform centroid behind the root quarter-chord, which lowers the flutter velocity of fins whose centroid is behind mid-root-chord (including RocketPy's default sweep_length = root_chord - tip_chord) and raises it for fins whose centroid is ahead. For the issue 615 worked example, epsilon is 0.349, and with it the formula reproduces the printed 1425, 1991.7 and 2618.1 ft/s. It would also need a guard for fins swept forward far enough that epsilon reaches zero. An open OpenRocket feature request proposes following Bennett's method (openrocket/openrocket#2709). Would you want this as a follow-up, and if so as the default or as an option?

• Peak of Flight issue 615: https://www.apogeerockets.com/Peak-of-Flight/Newsletter615
• Bennett's reference calculator (spreadsheet and Python, BSD-2): https://github.com/jkb-git/Fin-Flutter-Velocity-Calculator

I have a fix with tests ready and will open a PR referencing this issue.