AeroelasticityFluid-Structure InteractionTurbomachineryStructural DynamicsCascade Flutter35 min read

Aeroelasticity & Fluid-Structure Interaction (FSI) in Turbomachinery

An advanced engineering analysis of fluid-structure interaction (FSI) in axial-flow turbomachines. This textbook-quality guide covers Campbell diagrams, cascade flutter, vortex-induced vibrations, aerodynamic work loops, and a rigorous 2-DOF mathematical derivation with a worked numerical example.

Aeroelasticity & Fluid-Structure Interaction (FSI) in Turbomachinery

Aeroelasticity is the multidisciplinary field of engineering that studies the physical phenomena resulting from the interaction of aerodynamic, elastic, and inertial forces. In turbomachinery, this interaction takes place within high-performance, complex axial-flow configurations such as those in aircraft gas turbines, industrial steam turbines, and rocket engine turbopumps. Unlike the classical aeroelasticity of aircraft wings, which is characterized by isolated, high-aspect-ratio lifting surfaces in an unbounded medium, turbomachinery aeroelasticity operates in a highly confined, multi-body cascade environment. The proximity of adjacent blades, high structural frequencies, and extremely high fluid densities create unique fluid-structure interaction (FSI) dynamics that represent one of the most critical design boundaries for modern turbomachines.

Historically, aeroelastic instabilities in turbomachinery, such as blade flutter and forced response resonance, have led to catastrophic structural damage, performance degradation, and engine loss. Consequently, understanding and modeling the coupling between the unsteady aerodynamics of the cascade and the structural dynamics of the rotor blade assembly is paramount. In this article, we present a comprehensive textbook-level study of turbomachinery aeroelasticity. We analyze the mechanisms of dynamic flutter, vortex-induced vibrations (VIV), aerodynamic work loops, and Campbell diagrams. Furthermore, we derive the governing 2-degree-of-freedom (2-DOF) equations of motion for a pitching and heaving airfoil, incorporating Theodorsen's unsteady aerodynamic theory, and solve a representative engineering problem to calculate the critical flutter speed of a turbine blade cascade.

1. The Aeroelastic Framework: Collar's Triangle and Cascade Coupling

The fundamental structure of aeroelasticity is represented by Collar's Triangle, which details the coupling between three types of forces:

Inertial Forces: Governed by the mass distribution, density, and rotational accelerations of the blade structure.
Elastic Forces: Determined by the structural stiffness, material properties (Young's modulus, shear modulus), and blade geometry.
Aerodynamic Forces: Produced by the pressure and shear stress distributions acting on the blade surfaces due to the fluid flow.

In external aerodynamics (e.g., aircraft wings), the fluid-to-structure mass ratio $ \mu = \frac{m}{\pi \rho b^2} $ is high, typically between 50 and 200, which means that the structural inertia dominates the system, and added mass effects are negligible. In internal turbomachinery flow, however, the working fluid is highly compressed, leading to very low mass ratios, often in the range of 2 to 10. Consequently, the fluid inertial forces (added mass) and aerodynamic damping are of the same order of magnitude as the structural inertial and elastic forces. The fluid cannot be treated as a passive loading mechanism; it acts as an additional mass and damping element coupled to the structure.

Furthermore, turbomachinery blades are arranged in a circular cascade. The aerodynamic field of any individual blade is strongly influenced by the position and velocity of its neighboring blades. When one blade vibrates, it sheds vorticity into the flow that alters the angle of attack and pressure distribution of the adjacent blades. This cascade coupling requires that aeroelastic analysis be conducted on the entire blade assembly, accounting for the phase differences between blades, rather than modeling a single isolated blade.

2. Structural Dynamics: Natural Modes, Centrifugal Effects, and Campbell Diagrams

Turbomachinery blades are complex, three-dimensional twisted cantilevered structures that exhibit a variety of vibrational modes. The primary modes of interest are:

Flexural (Bending) Modes: Flapwise (out-of-plane, perpendicular to the chord) and edgewise (in-plane, parallel to the chord) bending.
Torsional Modes: Twist of the blade about its elastic axis.
Coupled Modes: Combinations of bending and torsion, which are particularly susceptible to flutter due to energy transfer between the modes.

The natural frequencies of these modes are heavily modified by the rotor's rotational speed. As the rotor spins, centrifugal forces pull the blades radially outward, exerting a large tensile stress field. This centrifugal tension acts as a geometric stiffness effect, which increases the bending natural frequencies. This centrifugal stiffening is mathematically modeled using the Southwell relation:

$$ f_n^2 = f_{n,0}^2 + K_s \left( \frac{\Omega}{60} \right)^2 $$

where $ f_{n,0} $ is the natural frequency of the $ n $-th mode at standstill, $ \Omega $ is the rotational speed in RPM, and $ K_s $ is the Southwell coefficient, a non-dimensional constant that depends on the mode shape and the blade's radial mass distribution. Typically, $ K_s $ is positive for bending modes, indicating a frequency increase with speed, but is near zero or slightly negative for pure torsional modes.

Conversely, high operating temperatures in the turbine sections cause thermal softening of the blade material. The elasticity modulus $ E $ decreases with temperature, reducing the structural stiffness and thus lowering the natural frequencies:

$$ E(T) = E_0 - \beta (T - T_0) $$

The Campbell Diagram is a critical design tool that plots these speed-dependent natural frequencies alongside Engine Order (EO) excitation lines. An Engine Order excitation is a periodic forcing frequency caused by flow non-uniformities (such as the wakes from upstream vanes or struts). The frequency of the $ k $-th engine order excitation is given by:

$$ f_{EO} = k \frac{\Omega}{60} $$

Intersections between the blade natural frequency curves and the EO lines indicate resonance conditions. If the operating range of the engine falls near these crossing points, the blade will experience forced response resonance, leading to rapid fatigue accumulation. To prevent this, engineers use Campbell diagrams to adjust the blade stiffness or vane count, moving these resonance points outside the operating speed envelope.

To account for the spatial variation of these modes across the entire bladed disc assembly, engineers utilize the SAFE (Singh's Advanced Frequency Evaluation) Diagram (also known as the Interference Diagram). The SAFE diagram plots frequency against the number of nodal diameters (ND), which represent the spatial phase variation of the mode shape around the rotor circumference. Resonance can only occur if two conditions are met simultaneously:

Frequency Match: The excitation frequency matches the natural frequency of a blade mode ($ f_{excitation} = f_{mode} $).
Spatial Shape Match (Nodal Diameter Match): The number of nodal diameters of the excitation pattern matches the nodal diameter of the structural mode ($ ND_{excitation} = ND_{mode} $).

This spatial match requirement is defined by the stator vane and rotor blade counts. An upstream stator row with $ N_v $ vanes produces an excitation with spatial harmonics. The excitation spatial harmonic order is:

$$ ND = m N_v \pm n N_b $$

where $ N_b $ is the rotor blade count, and $ m, n $ are integers. The SAFE diagram acts as a three-dimensional extension of the Campbell diagram, ensuring that spatial mode shapes and excitation shapes do not couple constructively to drive resonance.

3. Fluid-Structure Interaction (FSI) in Turbomachinery

Fluid-Structure Interaction (FSI) refers to the coupled physics where the fluid flow deforms the solid structure, and the structural deformation in turn modifies the fluid flow field. FSI numerical simulations are classified based on the level of coupling between the fluid and structural domains:

One-Way FSI

In one-way FSI, the flow field is solved independently on a rigid or statically deformed geometry using Computational Fluid Dynamics (CFD). The resulting pressure and shear stress distributions at the blade surface are extracted and mapped as boundary conditions onto a structural finite element method (FEM) solver. The structural solver then computes the stresses and deformations. One-way FSI assumes that the structural displacements are small enough that they do not feedback and significantly alter the flow field. While computationally efficient, this method is incapable of predicting self-excited instabilities like flutter, where the structural motion is the primary driver of the unsteady aerodynamic forces.

Two-Way Partitioned FSI

Two-way partitioned FSI maintains separate solvers for the fluid and structural domains, using a coupling interface to transfer data at each time step. The CFD solver computes the fluid forces and applies them to the structural boundary. The structural solver (CSD) computes the displacements, which are then passed back to the CFD solver, requiring the fluid mesh to deform to fit the new solid boundary. Partitioned solvers can be weakly coupled (explicit) or strongly coupled (implicit):

Weak Coupling (Explicit): The fluid and structural solvers are executed once per time step. The fluid solver passes forces, the structural solver solves for displacements, and the mesh is updated for the next time step. This approach is prone to numerical instability (known as the artificial mass effect) when the fluid density is comparable to the solid density, a common scenario in high-pressure steam and gas turbines.
Strong Coupling (Implicit): The fluid and structural equations are solved iteratively within each time step. Sub-iterations are performed, transferring forces and displacements back and forth until the interface variables converge. This method is numerically stable for low mass ratios but demands substantial computational resources.

Monolithic FSI

Monolithic FSI solves the fluid and structural equations simultaneously in a single mathematical framework, using a unified system of equations. This avoids the interface lag and coupling iterations of partitioned methods, providing maximum numerical stability. However, monolithic solvers are highly complex to formulate and program, as they require a single solver capable of handling both the hyperbolic-elliptic equations of fluid flow and the elliptic equations of structural mechanics.

A critical challenge in two-way FSI is mesh deformation. As the blade vibrates, the CFD grid must be updated dynamically without producing cell crossover or highly skewed elements. Advanced algorithms such as Radial Basis Functions (RBF) interpolation, transfinite interpolation (TFI), and torsional spring analogy are used to smooth the deformation field throughout the CFD domain, preserving mesh quality and convergence.

4. Dynamic Blade Flutter in Cascades

Flutter is a self-excited, self-sustained aeroelastic instability. Unlike forced response resonance, which occurs at specific rotational speeds matching engine order lines, flutter can occur at any operating condition where the aerodynamic damping becomes negative. The blade extracts energy from the steady flow, and the amplitude of vibration grows exponentially until a limit cycle oscillation (LCO) is reached or structural failure occurs.

The likelihood of flutter is characterized by the Reduced Frequency ($ k $), which is a non-dimensional measure of the flow unsteadiness:

$$ k = \frac{\omega c}{2 U} $$

where $ \omega $ is the vibrational frequency, $ c $ is the blade chord, and $ U $ is the flow velocity.

$ k < 0.05 $: The flow is quasi-steady; the aerodynamic forces adjust instantaneously to the blade motion. Flutter is highly unlikely.
$ 0.05 < k < 0.2 $: The flow is moderately unsteady. Phase lags between the blade motion and the unsteady aerodynamic forces begin to emerge.
$ k > 0.2 $: The flow is highly unsteady. The phase lag between blade motion and the resulting aerodynamic forces is large enough to drive self-excited flutter.

In turbomachinery cascades, flutter is a collective phenomenon. Because the blades are arranged in a circular rotor, their vibration is described by a travelling wave around the wheel, characterized by the Inter-Blade Phase Angle (IBPA, $ \sigma_n $). In a tuned rotor (where all blades are identical), the blades vibrate at the same frequency and amplitude, but with a constant phase difference $ \sigma_n $ between adjacent blades:

$$ \sigma_n = \frac{2\pi n}{N}, \quad n = 0, 1, \dots, N-1 $$

where $ N $ is the number of blades and $ n $ is the nodal diameter index. A positive $ n $ indicates a forward-travelling wave (moving in the direction of rotor rotation), while a negative $ n $ indicates a backward-travelling wave. For a given structural mode, the aerodynamic stability must be evaluated for all possible values of IBPA. Flutter occurs if the aerodynamic damping is negative for even a single IBPA value.

Aeroelastic stability in turbomachinery is often mapped across different flow regimes:

Subsonic/Transonic Stall Flutter: Occurs when the engine operates at off-design conditions with high incidence angles. The flow separates periodically from the blade suction surface (dynamic stall), creating a strong phase-lagged force that drives torsional oscillations.
Supersonic Unstalled Flutter: Occurs at design speeds with supersonic inlet velocities and low incidence. The shock waves within the blade passages move in response to the blade vibrations. The phase lag between the shock motion and the blade displacement can feed acoustic energy back into the blade, driving instabilities.
Choke Flutter: Occurs at negative incidence angles under choking conditions. The passage of a sonic throat through the cascade is highly sensitive to blade displacements, creating unstable pressure distributions.

To mitigate flutter, engineers use mistuning, which is the intentional introduction of small variations in the structural properties (mass, stiffness, or natural frequency) of individual blades around the rotor wheel. Mistuning disrupts the spatial symmetry of the rotor, preventing the formation of coherent travelling wave modes (IBPA). The vibrational energy becomes localized, which prevents the aerodynamic forces from doing constructive work on the entire cascade, thereby stabilizing the rotor.

5. Aerodynamic Work Loops and the Energy Method

In turbomachinery design, the stability of a cascade is commonly assessed using the Energy Method of Carta (1967). Because turbomachinery blades are structurally very stiff compared to the aerodynamic forces acting on them, the aerodynamic forces do not significantly alter the structural mode shape. Thus, we can assume that the blade vibrates in a single, predetermined structural mode shape (e.g., pure bending or pure torsion) at a given IBPA, and calculate the aerodynamic work done on the blade over one vibration cycle.

The aerodynamic work per cycle $ W_{aero} $ is defined as the integral of the product of the unsteady pressure $ p(t) $ and the blade surface normal velocity $ \mathbf{v}(t) $ over the blade area and one period of oscillation $ T $:

$$ W_{aero} = \int_0^T \oint_{Blade} p(t, \mathbf{s}) \left( \mathbf{v}(t, \mathbf{s}) \cdot \mathbf{n} \right) dA \, dt $$

The stability of the system is determined by the sign of $ W_{aero} $:

Stable ($ W_{aero} < 0 $): The unsteady aerodynamic forces oppose the motion, extracting energy from the structure and dissipating it into the fluid. The aerodynamic forces act as positive damping.
Unstable/Flutter ($ W_{aero} > 0 $): The unsteady aerodynamic forces act in the direction of the blade motion, doing positive work and feeding energy into the structure. The aerodynamic damping is negative. If this negative aerodynamic damping exceeds the positive structural damping, flutter occurs.

This energy transfer is governed by the phase angle $ \phi $ between the blade displacement and the unsteady aerodynamic force. If we assume a harmonic displacement $ x(t) = x_0 \cos(\omega t) $ and a harmonic aerodynamic force $ F(t) = F_0 \cos(\omega t + \phi) $, the work done per cycle is:

$$ W = \oint F \, dx = \int_0^{2\pi/\omega} F_0 \cos(\omega t + \phi) \left( -\omega x_0 \sin(\omega t) \right) dt = \pi x_0 F_0 \sin(\phi) $$

If the force leads the displacement ($ 0 < \phi < \pi $), then $ \sin(\phi) > 0 $, resulting in positive work ($ W > 0 $) and flutter. If the force lags the displacement ($ -\pi < \phi < 0 $), then $ \sin(\phi) < 0 $, and the aerodynamic forces stabilize the blade.

This behavior is visualized through aerodynamic work loops. By plotting the unsteady force or moment against the displacement or twist angle over one cycle, a closed loop is formed. The area of this loop represents the work per cycle:

A loop traversed in a direction where force and displacement are in phase (producing positive area) represents a net energy input to the blade structure.
A loop traversed in the opposite direction represents net energy dissipation, which stabilizes the structure.

In compressor or fan blade design, the local work distribution is evaluated along the span. Unsteady CFD simulations calculate the local work coefficient at multiple radial sections. This allows designers to identify which regions of the blade are stabilizing (dissipating energy) and which are destabilizing (generating positive work). For instance, if the tip region of a fan blade is producing positive work due to shock-motion hysteresis, the blade shape can be optimized locally (e.g., through sweep or lean adjustments) to force the local work back into a negative, stabilizing regime.

6. Vortex-Induced Vibrations (VIV) and Lock-in

Vortex-Induced Vibration (VIV) is a resonance-based fluid-structure interaction phenomenon. When fluid flows past a blade, alternating vortices are shed from the trailing edge, forming a von Kármán vortex street. The shedding of these vortices creates a periodic pressure asymmetry, which exerts a fluctuating lift force on the blade transverse to the flow direction.

The frequency of vortex shedding $ f_s $ is dictated by the Strouhal relation:

$$ f_s = \frac{St \cdot U}{d} $$

where $ St $ is the Strouhal number, $ U $ is the local flow velocity, and $ d $ is the characteristic trailing-edge thickness of the blade. For typical turbomachinery profiles, $ St $ is approximately 0.2.

As the flow velocity $ U $ increases, the vortex shedding frequency $ f_s $ increases linearly. When $ f_s $ approaches one of the natural frequencies of the blade $ f_n $, the structural vibration begins to influence the vortex shedding. Instead of continuing to rise linearly with flow speed, the vortex shedding frequency becomes "captured" or **locked-in** to the structural natural frequency:

$$ f_s \approx f_n $$

This lock-in (or synchronization) phenomenon persists over a discrete range of flow velocities. During lock-in, the vortex shedding becomes highly coherent and aligned along the span, which dramatically increases the amplitude of the transverse aerodynamic force. The resulting large-amplitude vibrations can cause high-cycle fatigue (HCF) failure in a very short time. Unlike flutter, which is self-excited and grows indefinitely, VIV is self-limiting because the structural displacement eventually disrupts the vortex organization at very large amplitudes, but the stresses experienced during lock-in are often well above the fatigue limit of turbine alloys.

In turbomachinery, VIV is particularly critical for the trailing edges of turbine blades, which are often relatively thick to accommodate internal cooling passages. The shed vortices can produce intense high-frequency excitation. To mitigate VIV, engineers employ trailing-edge modifications such as asymmetric bevels, cutbacks, or chamfers. These geometrical features disrupt the alternate shedding of vortices, weakening the coherent vortex structures and shifting the Strouhal shedding frequencies away from the structural modes.

7. Mathematical Formulations: 2-DOF Airfoil Aeroelastic Model

To analyze the stability boundary of an airfoil, we model the system with two degrees of freedom: heave displacement $ h $ (vertical translational displacement, positive downward) and pitch angle $ \alpha $ (rotational displacement about the elastic axis, positive nose-up).

Let $ b = c/2 $ be the semi-chord, $ m $ the mass per unit span, $ I_{\alpha} $ the mass moment of inertia per unit span about the elastic axis (EA), and $ S_{\alpha} $ the static imbalance (static coupling term) defined as:

$$ S_{\alpha} = m x_{\alpha} b $$

where $ x_{\alpha} $ is the non-dimensional distance from the elastic axis to the center of gravity (CG), positive aft. The governing structural equations of motion are:

$$ m \ddot{h} + S_{\alpha} \ddot{\alpha} + c_h \dot{h} + k_h h = -L_{aero} $$

$$ I_{\alpha} \ddot{\alpha} + S_{\alpha} \ddot{h} + c_{\alpha} \dot{\alpha} + k_{\alpha} \alpha = M_{aero} $$

where $ k_h $ and $ k_{\alpha} $ are the heave and pitch stiffnesses, and $ c_h $ and $ c_{\alpha} $ are the structural damping coefficients. $ L_{aero} $ is the aerodynamic lift (positive upward, hence the negative sign on the right-hand side of the heave equation), and $ M_{aero} $ is the aerodynamic moment about the elastic axis (positive nose-up).

8. Derivation of Lift and Moment: Theodorsen's Unsteady Aerodynamic Theory

Under the assumptions of thin airfoil theory in an incompressible, potential flow, the unsteady lift and moment can be derived by combining non-circulatory (added mass) and circulatory forces. The non-circulatory forces represent the inertial reaction of the fluid displaced by the accelerating airfoil, while the circulatory forces represent the aerodynamic lift due to the circulation generated by the boundary layer and trailing-edge vortex wake.

Step 1: Non-Circulatory Lift and Moment

The non-circulatory lift $ L_{NC} $ and moment $ M_{NC} $ are derived from the potential flow around a thin plate undergoing acceleration. The acceleration of the fluid creates pressure distributions that integrate to:

$$ L_{NC} = \pi \rho b^2 \left( \ddot{h} + U \dot{\alpha} - b a \ddot{\alpha} \right) $$

$$ M_{NC} = \pi \rho b^2 \left[ b a \ddot{h} - U b \left( \frac{1}{2} - a \right) \dot{\alpha} - b^2 \left( \frac{1}{8} + a^2 \right) \ddot{\alpha} \right] $$

where $ a $ is the non-dimensional location of the elastic axis measured from the mid-chord (positive aft, with the leading edge at $ -1 $ and the trailing edge at $ +1 $). The term $ \pi \rho b^2 $ represents the mass of a cylinder of air of diameter equal to the chord length.

Step 2: Circulatory Lift and Moment

The circulatory lift depends on the boundary layer circulation, which is governed by the Kutta condition at the trailing edge. Because the airfoil is oscillating, vortices are continuously shed into the wake. The wake vortices induce a velocity field (downwash) on the airfoil surface that reduces the effective angle of attack. The circulatory lift is driven by the downwash velocity evaluated at the $ 3/4 $-chord point:

$$ w_{3/4} = \dot{h} + U \alpha + b \left( \frac{1}{2} - a \right) \dot{\alpha} $$

Using Theodorsen's function $ C(k) = F(k) + i G(k) $ to account for the phase lag and attenuation caused by the wake vortex sheet, the circulatory lift $ L_C $ is:

$$ L_C = 2 \pi \rho U b C(k) w_{3/4} = 2 \pi \rho U b C(k) \left[ \dot{h} + U \alpha + b \left( \frac{1}{2} - a \right) \dot{\alpha} \right] $$

The circulatory lift acts at the aerodynamic center (AC), which for thin airfoil theory is located at the $ 1/4 $-chord point ($ x_{AC} = -b/2 $). The moment about the elastic axis due to this lift force is:

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

Step 3: Total Aerodynamic Lift and Moment

Summing the circulatory and non-circulatory contributions, we obtain the complete expressions for the aerodynamic forces acting on the oscillating airfoil:

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

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

These equations show that the lift and moment contain:

Inertial (Added Mass) Terms: Proportional to $ \ddot{h} $ and $ \ddot{\alpha} $.
Damping Terms: Proportional to $ \dot{h} $ and $ \dot{\alpha} $. Note that some terms depend on velocity $ U $ (quasi-steady damping).
Stiffness Terms: Proportional to $ \alpha $. These represent the aerodynamic lift and moment due to steady incidence, acting as an aerodynamic stiffness. There is no aerodynamic stiffness term proportional to heave $ h $, because a purely vertical translation does not change the angle of attack in a uniform flow.

9. Comparative Analysis of Aeroelastic Phenomena in Turbomachinery

Phenomenon	Key Driving Force	Primary Mode Shape	Reduced Frequency $ k $	Sensitivity Parameters	Mitigation Strategies
Classical Flutter	Coherent phase-lagged pressure waves; potential flow	Coupled Bending-Torsion	$ k \approx 0.1 - 0.35 $	Mass ratio, frequency ratio, EA position	Mistuning, increased structural stiffness, mass balance
Subsonic Stall Flutter	Dynamic stall, flow separation, dynamic hysteresis	Pure Torsion	$ k < 0.1 $	Stagger angle, mean flow incidence, solidity	Blade profiling, stagger adjustment, casing treatment
Supersonic Flutter	Shock wave oscillations, shock-boundary layer interaction	Bending or Torsion	$ k > 0.3 $	Inlet Mach number, shock position, pressure ratio	Leading-edge sweep, thin airfoils, friction damping
Vortex-Induced Vibration	Von Kármán vortex shedding resonance (lock-in)	Pure Bending	$ k \approx 0.15 - 0.25 $	Trailing-edge thickness, Strouhal number, damping	Asymmetric trailing-edge beveling, friction dampers
Forced Response Resonance	Rotor-stator wake intersection, potential field coupling	Bending, Torsion, Plate modes	Wide Range	Vane count, rotor speed, spatial harmonics	Mistuning, shroud coupling, under-platform dampers

10. Two-Degree-of-Freedom Airfoil Model Diagram

The schematic below shows a 2-DOF airfoil cross-section. It illustrates the coordinate systems and key variables used in the mathematical formulation: the elastic axis (EA), the center of gravity (CG), the aerodynamic center (AC), the heave displacement $ h $, the pitch angle $ \alpha $, and the aerodynamic forces (lift $ L $, drag $ D $, and moment $ M $).

11. Worked Numerical Example: Critical Flutter Speed of a Turbine Blade

In this worked example, we calculate the critical flutter speed $ U_{crit} $ and the flutter frequency $ f_F $ for a turbine blade model using the quasi-steady aeroelastic formulation. The blade is modeled as a 2-DOF system with bending (heave) and torsional (pitch) displacements.

Problem Specification

The turbine blade model possesses the following structural and geometric properties:

Chord Length: $ c = 0.15 \text{ m} $ (hence semi-chord $ b = c/2 = 0.075 \text{ m} $).
Mass per unit span: $ m = 4.5 \text{ kg/m} $.
Radius of Gyration: $ r_{\alpha} = 0.5 $ about the Elastic Axis (EA).
Mass Moment of Inertia: $ I_{\alpha} = m b^2 r_{\alpha}^2 = 4.5 \times 0.075^2 \times 0.5^2 = 0.006328125 \text{ kg}\cdot\text{m} $.
Static Imbalance Parameter: $ x_{\alpha} = 0.2 $ (center of gravity is 20% semi-chord aft of the elastic axis).
Static Imbalance: $ S_{\alpha} = m x_{\alpha} b = 4.5 \times 0.2 \times 0.075 = 0.0675 \text{ kg}\cdot\text{m/m} $.
Bending (Heave) Natural Frequency: $ \omega_h = 120 \text{ rad/s} $ ($ f_h \approx 19.1 \text{ Hz} $).
Torsional (Pitch) Natural Frequency: $ \omega_{\alpha} = 240 \text{ rad/s} $ ($ f_{\alpha} \approx 38.2 \text{ Hz} $).
Heave Stiffness: $ k_h = m \omega_h^2 = 4.5 \times 120^2 = 64,800 \text{ N/m} $.
Pitch Stiffness: $ k_{\alpha} = I_{\alpha} \omega_{\alpha}^2 = 0.006328125 \times 240^2 = 364.5 \text{ N}\cdot\text{m/rad} $.
Structural Damping Ratio: $ \zeta_h = \zeta_{\alpha} = 0.01 $.
Structural Damping Coefficients:
$$ c_h = 2 \zeta_h m \omega_h = 2 \times 0.01 \times 4.5 \times 120 = 10.8 \text{ N}\cdot\text{s/m} $$

$$ c_{\alpha} = 2 \zeta_{\alpha} I_{\alpha} \omega_{\alpha} = 2 \times 0.01 \times 0.006328125 \times 240 \approx 0.030375 \text{ N}\cdot\text{m}\cdot\text{s} $$
Elastic Axis Position: $ a = -0.4 $ (located at 30% chord, which is 10% chord forward of the mid-chord).
Aerodynamic Center: $ a_{ac} = -0.5 $ (located at 25% chord).
Aerodynamic Center to EA Distance: $ e = a - a_{ac} = -0.4 - (-0.5) = 0.1 $.
Air Density: $ \rho = 1.225 \text{ kg/m}^3 $.

Step 1: Governing Aeroelastic Equations

Under the quasi-steady aerodynamic model, the lift $ L $ and moment $ M $ are modeled as:

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

$$ M = L \cdot e \cdot b = 2 \pi \rho U^2 b^2 e \alpha + 2 \pi \rho U b^2 e \dot{h} + 2 \pi \rho U b^3 e \left( \frac{1}{2} - a \right) \dot{\alpha} $$

Let $ C_{coeff} = 2 \pi \rho b = 2 \pi \times 1.225 \times 0.075 \approx 0.577268 \text{ kg/m}^2 $. Substituting these definitions into the structural equations of motion yields the coupled system:

$$ \begin{bmatrix} m & S_{\alpha} \\ S_{\alpha} & I_{\alpha} \end{bmatrix} \begin{bmatrix} \ddot{h} \\ \ddot{\alpha} \end{bmatrix} + \begin{bmatrix} c_{11} & c_{12} \\ c_{21} & c_{22} \end{bmatrix} \begin{bmatrix} \dot{h} \\ \dot{\alpha} \end{bmatrix} + \begin{bmatrix} k_{11} & k_{12} \\ k_{21} & k_{22} \end{bmatrix} \begin{bmatrix} h \\ \alpha \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} $$

where the coefficients of the damping and stiffness matrices are speed-dependent:

$ c_{11} = C_{coeff} U + c_h = 0.577268 U + 10.8 $
$ c_{12} = C_{coeff} U b \left( \frac{1}{2} - a \right) = 0.577268 U \times 0.075 \times 0.9 = 0.038966 U $
$ c_{21} = -C_{coeff} U b e = -0.577268 U \times 0.075 \times 0.1 = -0.004330 U $
$ c_{22} = -C_{coeff} U b^2 e \left( \frac{1}{2} - a \right) + c_{\alpha} = -0.577268 U \times 0.075^2 \times 0.1 \times 0.9 + 0.030375 = -0.000292 U + 0.030375 $
$ k_{11} = k_h = 64800 $
$ k_{12} = C_{coeff} U^2 = 0.577268 U^2 $
$ k_{21} = 0 $
$ k_{22} = k_{\alpha} - C_{coeff} U^2 b e = 364.5 - 0.577268 U^2 \times 0.075 \times 0.1 = 364.5 - 0.004330 U^2 $

Step 2: Characteristic Equation

We assume solutions of the form $ h(t) = \bar{h} e^{s t} $ and $ \alpha(t) = \bar{\alpha} e^{s t} $. The system has non-trivial solutions if the determinant of the system matrix vanishes:

$$ \det \left( \begin{bmatrix} m & S_{\alpha} \\ S_{\alpha} & I_{\alpha} \end{bmatrix} s^2 + \begin{bmatrix} c_{11} & c_{12} \\ c_{21} & c_{22} \end{bmatrix} s + \begin{bmatrix} k_{11} & k_{12} \\ k_{21} & k_{22} \end{bmatrix} \right) = 0 $$

Expanding the determinant yields the fourth-order characteristic polynomial:

$$ A_4 s^4 + A_3 s^3 + A_2 s^2 + A_1 s + A_0 = 0 $$

where the coefficients are:

$ A_4 = m I_{\alpha} - S_{\alpha}^2 = 4.5 \times 0.006328125 - 0.0675^2 = 0.02392031 \text{ kg}^2\cdot\text{m} $
$ A_3 = m c_{22} + I_{\alpha} c_{11} - S_{\alpha} (c_{12} + c_{21}) $
$ A_2 = m k_{22} + I_{\alpha} k_{11} - S_{\alpha} (k_{12} + k_{21}) + (c_{11} c_{22} - c_{12} c_{21}) $
$ A_1 = c_{11} k_{22} + c_{22} k_{11} - c_{12} k_{21} - c_{21} k_{12} $
$ A_0 = k_{11} k_{22} - k_{12} k_{21} = k_h k_{22} $

Step 3: Routh-Hurwitz Stability Criterion

For a fourth-order system, the boundary of stability occurs when the third Routh determinant $ \Delta_3 $ becomes zero:

$$ \Delta_3 = A_3 A_2 A_1 - A_4 A_1^2 - A_3^2 A_0 = 0 $$

At this boundary, the characteristic roots are purely imaginary, $ s = \pm i \omega_F $, where $ \omega_F $ is the flutter frequency. The frequency is related to the coefficients by:

$$ \omega_F = \sqrt{\frac{A_1}{A_3}} $$

By performing a numerical bisection on $ \Delta_3(U) $ between the stable speed $ U = 40 \text{ m/s} $ and the unstable speed $ U = 45 \text{ m/s} $, we locate the root:

$$ U_{crit} = 41.48 \text{ m/s} $$

Step 4: Coefficient Verification and Frequency Calculation

Evaluating the matrices and coefficients at the exact critical velocity $ U = 41.4779 \text{ m/s} $:

$ c_{11} = 0.577268 \times 41.4779 + 10.8 = 34.7438 \text{ N}\cdot\text{s/m} $
$ c_{12} = 0.038966 \times 41.4779 = 1.6162 \text{ N}\cdot\text{s} $
$ c_{21} = -0.004330 \times 41.4779 = -0.1796 \text{ N}\cdot\text{m}\cdot\text{s/m} $
$ c_{22} = -0.000292 \times 41.4779 + 0.030375 = 0.01826 \text{ N}\cdot\text{m}\cdot\text{s} $
$ k_{11} = 64800 \text{ N/m} $
$ k_{12} = 0.577268 \times 41.4779^2 = 993.125 \text{ N/rad} $
$ k_{21} = 0 $
$ k_{22} = 364.5 - 0.004330 \times 41.4779^2 = 357.049 \text{ N}\cdot\text{m/rad} $

The polynomial coefficients evaluate to:

$ A_4 = 0.02392031 $
$ A_3 = 4.5 \times 0.01826 + 0.006328 \times 34.7438 - 0.0675 \times (1.6162 - 0.1796) = 0.20503125 $
$ A_2 = 4.5 \times 357.049 + 0.006328 \times 64800 - 0.0675 \times 993.125 + (34.7438 \times 0.01826 - 1.6162 \times (-0.1796)) = 1950.6817 $
$ A_1 = 34.7438 \times 357.049 + 0.01826 \times 64800 - 0 + 0.1796 \times 993.125 = 13766.5027 $
$ A_0 = 64800 \times 357.049 = 23,136,934.7 $

Checking the Routh determinant $ \Delta_3 $:

$$ \Delta_3 = (0.20503125 \times 1950.6817 \times 13766.5027) - (0.02392031 \times 13766.5027^2) - (0.20503125^2 \times 23136934.7) $$ $$ \Delta_3 = 5,505,963.09 - 4,533,300.99 - 972,626.54 \approx 0 $$

This confirms that the system is exactly at the stability boundary. The flutter frequency $ \omega_F $ is:

$$ \omega_F = \sqrt{\frac{A_1}{A_3}} = \sqrt{\frac{13766.5027}{0.20503125}} = 259.12 \text{ rad/s} $$

$$ f_F = \frac{\omega_F}{2 \pi} = \frac{259.12}{2 \pi} = 41.24 \text{ Hz} $$

At speeds below $ 41.48 \text{ m/s} $, $ \Delta_3 > 0 $, and the system is stable with damped oscillations. At speeds above $ 41.48 \text{ m/s} $, $ \Delta_3 < 0 $, and the system is unstable, leading to exponentially growing bending-torsion coupled vibrations.

References

Sisto, F., "Aeroelasticity in Turbomachinery: An Overview", Journal of Aircraft, Vol. 14, No. 11, 1977, pp. 1025-1031.
Bölcs, A., and Fransson, T. H., "Aeroelasticity in Turbomachines", Communication du Laboratoire de Thermique Appliquee et de Turbomachines, EPFL, Lausanne, 1986.
Dowell, E. H., Clark, R. L., Cox, D. E., Curtiss, H. C., Edwards, J. W., Hall, K. C., Peters, D. A., Scanlan, R. H., Simiu, E., Sisto, F., and Strganac, T. W., A Modern Course in Aeroelasticity, 5th Edition, Springer, 2015.
Platzer, M. F., and Carta, F. O. (Eds.), "AGARD Manual on Aeroelasticity in Axial-Flow Turbomachines", AGARD-AG-298, 1988.
Marshall, J. G., and Imregun, M., "A Review of Aeroelastic Calculations with Special Reference to Double-Row Turbomachinery Cascades", Journal of Fluids and Structures, Vol. 10, No. 3, 1996, pp. 237-256.