Aeroelastic Flutter Analysis

Q: Divergent vibration... does that mean it won't stop until the structure fails?

Precisely. The 1940 collapse of the Tacoma Narrows Bridge is a famous example. In aircraft design, sufficient margin against the flutter speed must be ensured across the entire flight envelope. FAA regulations (FAR 25.629) require a flutter margin of at least 1.15 times the flight speed.

Q: How is flutter described mathematically?

Let's start with a typical 2-degree-of-freedom wing section flutter model. Taking the bending displacement $h$ and torsion angle $\alpha$ as the degrees of freedom, the equations of motion are as follows: $$m\ddot{h} + S_\alpha\ddot{\alpha} + K_h h = -L$$ $$S_\alpha\ddot{h} + I_\alpha\ddot{\alpha} + K_\alpha \alpha = M_{EA}$$ Here, $m$ is the mass per unit span, $S_\alpha$ is the static unbalance moment, $I_\alpha$ is the moment of inertia about the elastic axis, and $K_h$ and $K_\alpha$ are the bending and torsional spring constants, respectively.

Q: How are the aerodynamic terms $L$ and $M_{EA}$ on the right-hand side expressed?

The unsteady aerodynamic forces are expressed using the Theodorsen function $C(k)$. For the reduced frequency $k = \omega b / U$, $$L = \pi\rho b^2[\ddot{h} + U\dot{\alpha} - ba\ddot{\alpha}] + 2\pi\rho U b C(k)\left[\dot{h} + U\alpha + b\left(\frac{1}{2}-a\right)\dot{\alpha}\right]$$ The Theodorsen function $C(k) = F(k) + iG(k)$, defined as a ratio of Hankel functions, represents the amplitude and phase lag of the circulatory aerodynamic force. It reduces to the steady aerodynamic force as $k \to 0$ and to the non-circulatory force as $k \to \infty$.

Category: 解析 | Integrated 2026-04-06

Theory and Physics

🎓

Flutter is a self-excited vibration caused by the dynamic coupling of fluid forces and structural elastic forces. When a certain critical speed (flutter speed $U_F$) is exceeded, the aerodynamic damping becomes negative, and the vibration grows divergently. It is a problem in aircraft wings, turbine blades, bridge cables, etc.

🧑‍🎓

Divergent vibration... does that mean it doesn't stop until it breaks?

🎓

That's right. The collapse of the Tacoma Narrows Bridge in 1940 is a famous example. In aircraft design, sufficient margin against the flutter speed must be ensured across the entire flight envelope. FAA regulations (FAR 25.629) require a flutter margin of 1.15 times the flight speed or more.

Governing Equations

🧑‍🎓

How is flutter described mathematically?

🎓

Let's start with a typical 2-degree-of-freedom airfoil flutter model. Taking the bending displacement $h$ and the torsion angle $\alpha$ as degrees of freedom, the equations of motion are as follows.

$$m\ddot{h} + S_\alpha\ddot{\alpha} + K_h h = -L$$

$$S_\alpha\ddot{h} + I_\alpha\ddot{\alpha} + K_\alpha \alpha = M_{EA}$$

Here, $m$ is the mass per unit span, $S_\alpha$ is the static unbalance moment, $I_\alpha$ is the moment of inertia about the elastic axis, and $K_h$ and $K_\alpha$ are the spring constants for bending and torsion, respectively.

🧑‍🎓

How are the aerodynamic terms $L$ and $M_{EA}$ on the right-hand side expressed?

🎓

Unsteady aerodynamic forces are expressed using the Theodorsen function $C(k)$. For the reduced frequency $k = \omega b / U$,

$$L = \pi\rho b^2[\ddot{h} + U\dot{\alpha} - ba\ddot{\alpha}] + 2\pi\rho U b C(k)\left[\dot{h} + U\alpha + b\left(\frac{1}{2}-a\right)\dot{\alpha}\right]$$

The Theodorsen function $C(k) = F(k) + iG(k)$, defined as the ratio of Hankel functions, represents the amplitude and phase lag of circulatory aerodynamic forces. It reduces to steady aerodynamic forces as $k \to 0$ and to non-circulatory forces as $k \to \infty$.

🧑‍🎓

How is it extended to general structures?

🎓

In the generalized form using modal coordinates $\mathbf{q}$, it becomes:

$$\mathbf{M}\ddot{\mathbf{q}} + \mathbf{C}\dot{\mathbf{q}} + \mathbf{K}\mathbf{q} = q_\infty \mathbf{Q}(k)\mathbf{q}$$

For dynamic pressure $q_\infty = \frac{1}{2}\rho U^2$, the generalized aerodynamic matrix $\mathbf{Q}(k)$ is calculated using the Doublet Lattice Method (DLM) or CFD.

Flutter Solution Methods

🧑‍🎓

What methods are there to determine the flutter speed?

🎓

Let's organize the representative solution methods.

Method	Characteristics	Application Scenario
V-g Method (American Method)	Introduces artificial structural damping $g$ and performs eigenvalue analysis	Preliminary design stage
p-k Method	Iteratively updates the reduced frequency; more physical	Detailed design
p-Method (State-Space Method)	Time-domain representation of aerodynamics via rational function approximation	Coupled with control systems
CFD-CSD Coupling	Direct calculation of nonlinear aerodynamic forces	Transonic flutter

In the p-k method, the following eigenvalue problem is solved.

$$\left[p^2\mathbf{M} + p\left(\mathbf{C} - \frac{q_\infty}{k}\text{Im}(\mathbf{Q})\right) + \mathbf{K} - q_\infty\text{Re}(\mathbf{Q})\right]\mathbf{q} = 0$$

The speed at which the real part $\gamma$ of $p = \gamma \pm i\omega$ crosses zero is the flutter speed.

🧑‍🎓

So linear theory can't be used in the transonic regime?

🎓

Exactly. In the transonic regime, the nonlinearity of aerodynamics becomes strong due to shock wave and boundary layer interaction. In this regime, a CFD-CSD approach that directly couples CFD (Euler or RANS equations) with structural FEM is necessary.

Implementation in Commercial Tools

🧑‍🎓

What software can be used for flutter analysis?

Tool	Method	Characteristics
MSC Nastran (SOL 145/146)	DLM + p-k/V-g Method	Aerospace industry standard
Ansys Mechanical + Fluent	CFD-CSD Coupling	Transonic capability
ZAERO (ZONA Technology)	ZONA6/7 Unsteady Panel Method	High-speed flutter
STAR-CCM+ + Abaqus	Co-simulation	Nonlinear FSI

🧑‍🎓

Nastran's SOL 145 is the industry standard. I'll remember that!

Coffee Break Yomoyama Talk

"V-n Diagram" and "Flutter Margin" — The Aeroelastic Wall That Keeps Designers Awake at Night

One of the most nerve-wracking tests in aircraft design is the "flutter flight test." International standards (FAR Part 25) mandate a margin of 15% or more above the design flutter speed (VF). In actual flight tests, speed is gradually increased up to the design dive speed (VD) to detect the onset of flutter. The problem is that flutter is a "suddenly occurring instability phenomenon" — even if signs are detected, it might be too late to turn back, making the test pilot's risk extremely high. In the 1990s, during the certification tests for a next-generation European airliner, theoretical calculations overlooked the effect of center-of-gravity shift due to external fuel tank attachment, and the test aircraft showed signs of flutter at a speed 20% lower than the design VF. Theoretical "oversights" manifesting in flight tests — that's why flutter theory must always be questioned with humility.

Physical Meaning of Each Term

Structural-Thermal Coupling Term: Thermal expansion due to temperature change induces structural deformation, and deformation affects the temperature field. $\sigma = D(\varepsilon - \alpha \Delta T)$. 【Everyday Example】Railroad tracks expanding in summer, narrowing the gaps — temperature rise → thermal expansion → stress generation is a typical example. Warping of circuit boards after soldering is also due to differences in thermal expansion coefficients of different materials. Engine cylinder blocks develop thermal stress due to temperature differences between hot and cold sections, potentially leading to cracks.
Fluid-Structure Interaction (FSI) Term: Bidirectional interaction where fluid pressure/shear forces deform the structure, and structural deformation changes the fluid domain. 【Everyday Example】Suspension bridge cables vibrating in strong winds (vortex-induced vibration) — wind forces shake the structure, the shaken structure alters the airflow, further amplifying vibration. Blood flow in the heart and elastic deformation of blood vessel walls, aircraft wing flutter (aeroelastic instability) are also typical FSI problems. One-way coupling may suffice in some cases, but bidirectional coupling is essential for large deformations.
Electromagnetic-Thermal Coupling Term: Feedback loop where Joule heating $Q = J^2/\sigma$ causes temperature rise, and temperature change alters electrical resistance. 【Everyday Example】Nichrome wire in an electric heater heats up (Joule heat) and glows red when current flows — temperature rise changes resistance, altering current distribution. Eddy current heating in IH cooking heaters, increased sag in power lines due to temperature rise are also examples of this coupling.
Data Transfer Term: Interpolation resolves mesh mismatch between different physical fields. 【Everyday Example】When calculating "feels-like temperature" by combining "air temperature data" and "wind data" in weather forecasting, interpolation is needed if observation points differ — similarly in CAE coupled analysis, structural mesh and CFD mesh generally don't match, so data transfer (interpolation) accuracy at the interface directly affects result reliability.

Assumptions and Applicability Limits

Weak Coupling Assumption (One-way coupling): Valid when one physical field affects the other but the reverse is negligible.
Cases Requiring Strong Coupling: Large deformations in FSI, strong temperature dependence in electromagnetic-thermal coupling.
Time Scale Separation: If characteristic times of each physical field differ significantly, efficiency can be improved via subcycling.
Interface Condition Consistency: Ensure energy/momentum conservation at the coupling interface is satisfied numerically.
Non-applicable Cases: When three or more physical fields are strongly coupled simultaneously, monolithic methods may be required.

Dimensional Analysis and Unit Systems

Variable	SI Unit	Notes / Conversion Memo
Thermal Expansion Coefficient $\alpha$	1/K	Steel: ~12×10⁻⁶, Aluminum: ~23×10⁻⁶
Coupled Interface Force	N/m² (pressure) or N (concentrated force)	Verify force balance between fluid and structure sides.
Data Transfer Error	Dimensionless (%)	Interpolation accuracy depends on mesh density ratio. Below 5% is a guideline.

Numerical Methods and Implementation

🎓

The wing surface is divided into trapezoidal panels, and a doublet of acceleration potential is placed on each panel. A system of equations is solved to satisfy the downwash boundary condition at the 3/4 chord point of each panel to obtain the pressure distribution.

$$w_j = \sum_{i=1}^{N} A_{ij}(k, M) \Delta c_{p,i}$$

$A_{ij}$ is the aerodynamic influence coefficient containing the kernel function, dependent on Mach number $M$ and reduced frequency $k$.

🧑‍🎓

Where are the accuracy limits of DLM?

🎓

It is based on linear potential theory, so accuracy degrades in the following cases:

Transonic regime ($M \approx 0.8$〜1.2): Shock wave effects.
High angle of attack: Flow involving separation.
When wing thickness effects are significant.

For these cases, CFD-based methods are necessary.

Mode Extraction and Aerodynamic Spline

🧑‍🎓

Structural modes and aerodynamic panels have completely different meshes, right? How are they connected?

🎓

Spline interpolation is used. Representative methods are as follows.

Spline	Method	Application
Infinite Plate Spline (IPS)	Fundamental solution for thin plate bending	Displacement interpolation within wing surface
Thin Plate Spline (TPS)	Improved version of IPS	General wing structures
RBF (Radial Basis Function)	Radial basis function	3D shapes
Beam Spline	1D interpolation	High aspect ratio wings

In Nastran, it is defined using SPLINE1/SPLINE2 cards. It interpolates structural node displacements to aerodynamic panel control points and conversely maps aerodynamic loads to structural nodes.

🎓

Using the transformation matrix $\mathbf{G}$,

$$\mathbf{u}_{aero} = \mathbf{G} \mathbf{u}_{struct}$$

$$\mathbf{F}_{struct} = \mathbf{G}^T \mathbf{F}_{aero}$$

This force transformation ensures conservation of virtual work.

p-k Method Implementation Procedure

🧑‍🎓

Please explain the specific calculation steps of the p-k method.

🎓

1. Modal analysis from structural model (SOL 103) → Eigenmodes $\phi_i$, natural frequencies $\omega_i$

2. Calculate the generalized aerodynamic matrix $\mathbf{Q}(M, k)$ for each Mach number and reduced frequency using DLM.

3. Iterative solution for each speed $U$:

Interpolate $\mathbf{Q}$ using initial estimate of $k$.
Solve eigenvalue problem to obtain $p = \gamma + i\omega$.
Update $\mathbf{Q}$ with new $k = \omega b / U$.
Repeat until $k$ converges.

4. From the plot of $\gamma(U)$, the speed where $\gamma = 0$ is crossed is the flutter speed.

🧑‍🎓

About how many iterations does it take to converge?

🎓

It usually converges in 3 to 5 iterations. However, in cases where modes are close, tracking becomes difficult, so mode tracking algorithms (e.g., tracking via MAC values) are important.