All projects
P14 / Course assignment · Aeroelasticity2025

Aeroelastic modelling & flutter prediction.

An individual model of a three-degree-of-freedom airfoil section with a control surface, built from first principles in Python to predict divergence, control reversal and flutter. The predicted flutter speed agrees with a NASA reference to within 1%.

TU DELFT · COURSE AE4ASM506Fundamentals of Aeroelasticity
PythonNumPy / SciPyTheodorsen theoryWagner functionState-space modellingEigenvalue analysis
Airfoil section with a trailing-edge flap, showing heave, pitch and flap degrees of freedom restrained by springs.
The typical section, with heave (h), pitch (θ) and flap (β) degrees of freedom restrained by the springs K_h, K_θ and K_β.

Predicted flutter speed

299 m/s

Within 1% of the NASA reference value of 302 m/s for the same typical section.

Assignment grade

9.0 / 10

Individual take-home assignment, TU Delft.

The problem

When a wing deforms under aerodynamic load, the deformation changes the load in turn. At high enough speed this feedback becomes unstable: the wing can twist until it fails (divergence), a control surface can start acting in the wrong direction (control reversal), or the structure can enter a self-excited oscillation that grows with every cycle (flutter).

This individual take-home assignment for TU Delft’s Fundamentals of Aeroelasticity course (AE4ASM506) asked for all three to be predicted for a classic test case: a two-dimensional airfoil section that can heave, pitch and deflect a trailing-edge flap, with structural parameters from a NASA reference study. I derived the model and implemented it in Python.

Structural model

Each degree of freedom is restrained by a spring, and the derivation gives a coupled mass matrix and a diagonal stiffness matrix. Solving the resulting eigenvalue problem shows how inertial coupling changes the natural frequencies:

ModeUncoupledCoupled
1 · Heave7.96 Hz (50 rad/s)7.75 Hz (48.7 rad/s)
2 · Pitch15.92 Hz (100 rad/s)17.64 Hz (110.8 rad/s)
3 · Flap47.75 Hz (300 rad/s)55.51 Hz (348.8 rad/s)

Natural frequencies of the typical section with and without inertial coupling.

Coupling softens the heave mode slightly and stiffens the pitch and flap modes, and the mode shapes no longer belong to a single degree of freedom.

From steady to unsteady aerodynamics

For the static problems, Theodorsen’s theory reduced to steady flow gives an aerodynamic stiffness matrix and an external force vector, and from them closed-form divergence and control-reversal speeds.

Flutter needs the unsteady aerodynamics. Using the indicial-function method, the lag in lift build-up after a change in motion is described by a two-term exponential approximation of the Wagner function. This adds six aerodynamic lag states, giving a linear 12-state model: three velocities, three displacements and six lag states. The same structure can take gust loads by adding lag states based on the Küssner function, which the assignment discussed but did not implement.

The derivation

Every matrix in the model was derived by hand before being implemented. First, substituting the effective angle of attack and the Wagner approximation into the convolution integral gives the circulatory lift in terms of the lag states:

Derivation of the unsteady circulatory lift expressed with six aerodynamic lag states, Equation 52 of the report.

Unsteady lift with lag states (report, Eq. 52). Click any equation block to enlarge it.

The lift and moments are then assembled into the aeroelastic equations of motion and rewritten as a first-order state-space system:

Assembly of the unsteady aeroelastic equations of motion and the 12-state state-space model, Equations 56 to 62 of the report.

Assembly into the state-space model ż = Az (report, Eqs. 56–62).

Finally, the aerodynamic mass, damping and stiffness matrices, the lag-state coupling matrices and the lag dynamics that make up the state matrix A:

All matrices of the state-space model: aerodynamic mass, circulatory and non-circulatory damping and stiffness, lag-state coupling, lag dynamics and input matrices, Equations 63 to 73 of the report.

The matrices of the state-space model (report, Eqs. 63–73).

Finding flutter

The eigenvalues of the state matrix were computed across a range of speeds. Flutter begins where a complex pair crosses into the right half-plane, so that an oscillation grows instead of decaying. Here, the pitch-dominated mode crosses at 299 m/s, against 302 m/s in the NASA reference, a deviation of less than 1%. The flutter mode is a coupled heave–pitch (bending–torsion) motion, with only a small flap contribution.

InstabilityCritical speed
Flutter299.0 m/s
Control reversal367.3 m/s
Divergence707.1 m/s

Critical aeroelastic speeds of the typical section.

Flutter occurs first and is therefore the design driver for this section. The assignment closed by identifying which single parameter would make each of the other instabilities critical instead: the position of the elastic axis relative to the aerodynamic centre for divergence, and the aerodynamic effectiveness of the flap for control reversal.