Introduction

Aeroelasticity is a specialized field in aerospace engineering that examines the interactions between aerodynamic forces, elastic deformations, and the structural dynamics of aerospace structures. The essence of the field can be described by a Collar Diagram^[1], which illustrates the interaction between the structural characteristics and the aerodynamic forces of a flight vehicle or another aeroelastic system. As shown in Figure 1, the Collar Diagram helps assess a structure’s stability and dynamic behavior under aerodynamic loads. This video provides an overview of this fundamental behavior as it affects the wings and tails of different airplanes.

The origins of aeroelasticity date back to the early days of aviation, when unexpected wing failures occurred on the first monoplanes. Many of these failures were traced to a behavior now known as aeroelastic divergence. As a critical airspeed was approached, the aerodynamic loads produced an increasing torsional moment about the wing’s elastic axis, as illustrated in Figure 2. If the resulting twist increased the local angle of attack, the lift would increase further, producing still more twist and an additional buildup of aerodynamic loads. This positive feedback could lead to excessively large torsional deformations and catastrophic structural failure unless the airspeed is quickly reduced.

Flutter, another behavior that can lead to wing or tail surface failure, is a self-excited oscillation arising from the interplay between unsteady aerodynamic forces and the structure’s dynamic response. As also shown in Figure 2, flutter is an oscillatory response that can grow rapidly in amplitude as a critical airspeed is approached. At or near this critical airspeed, a minor disturbance, such as a gust or control input, can trigger the onset of flutter. The flutter can reach a limiting amplitude, but it can also become a divergent oscillatory response in which aerodynamic work feeds energy into the structural motion. In these conditions, structural failure is also likely.

These observations, along with a history of structural failures in early airplanes, revealed that structural elasticity and flutter avoidance could not be ignored during the design process. Understanding aeroelastic phenomena became even more critical with the development of higher-speed and supersonic aircraft, in which aerodynamic forces and structural stresses are significantly more pronounced. Today, all aircraft must be designed for flutter-free operation over their entire flight envelopes, and thorough aeroelastic analyses and flight testing are integral to their design and certification process. This video from Airbus shows how engineers and pilots explore the flutter envelope of a new airplane.

History

The origins of aeroelastic behavior can be traced back to the earliest aircraft; see the review paper “Historical Development of Airplane Flutter” by I. E. Garrick and Wilmer H. Reed III. The Wright brothers and other early aviation pioneers sometimes encountered wing oscillations and buzzing. They also observed buzzing at the tips of their propellers, which they addressed by using wider-chord blades. They did not understand these effects at the time, but they were likely aeroelastic in origin, resulting from the interaction between aerodynamic forces and structural dynamics. In 1903, Samuel Langley attempted to fly an aircraft, but his monoplane’s wooden and fabric wings were structurally weak and failed under the aerodynamic loading, causing the aircraft to crash into the Potomac River. The episode is often cited as an early warning that aerodynamic loading, structural stiffness, and flight-vehicle stability could not be treated separately. The Wright brothers experimented with various structures and developed a more rigid biplane, which did not suffer as readily from aeroelastic effects. The catastrophic failure of the horizontal tail on a Handley Page O/400 bomber in 1916, caused by an oscillatory aeroelastic phenomenon later known as flutter, suggested that much remained to be learned about aircraft design.

Wing and control surface flutter became an increasingly critical concern as aircraft designs evolved, particularly with the transition to monoplanes. Monoplanes, with their single-wing configuration and lower structural stiffness,^[2] highlighted the importance of understanding and mitigating flutter during the design of aircraft. “The Blue Max” is a 1966 film set in WWI that tells the story of a German fighter pilot, Bruno Stachel, who aims to earn the prestigious Pour le Mérite, also known as the “Blue Max,” for outstanding acts of bravery in aerial combat. While wing divergence and flutter are not central themes in the film, it depicts the heightened risks of structural failure that early pilots faced as they transitioned from slow biplanes to faster, more agile monoplanes.

The rapid advancement in aircraft design and increased flight speeds after this period led to more frequent and severe aeroelastic problems. Engineers then systematically studied these issues, developing aeroelastic theories to predict and prevent flutter, often necessitating more rigorous structural design. The work of Theodore Theodorsen, among others, laid the groundwork for modern aeroelastic analyses. The severity of flutter issues in the decade following WWII is underscored by a 1956 state-of-the-art survey conducted by the NACA Subcommittee on Vibration and Flutter, which documented 54 instances of flutter-related difficulties with U.S. military aircraft. The introduction of thin swept wings for supersonic flight in the 1950s increased the challenges of predicting aeroelastic effects at higher Mach numbers. This period saw significant advances in theoretical, wind-tunnel, and flight testing, including a deeper understanding of unsteady aerodynamic effects.

The Lockheed L-188 Electra, a four-engine turboprop airliner introduced in the late 1950s, was found to exhibit a more complex but catastrophic form of flutter known as whirl-mode flutter. Under certain conditions, propeller oscillations were transmitted to the engine mountings, which, in turn, coupled to the wings, creating hazardous vibrations. After two midair breakups in which the wings detached from the airframe, wind tunnel tests (Figure 3) led to design changes, including reinforced engine mounts and thicker wing skins, making the aircraft flutter-free. While these failures damaged the reputation of the airplane and its manufacturer, the modified Electra became a reliable aircraft; the P-3 Orion, in particular, was used extensively by the U.S. Navy until quite recently. While flutter prediction methods continued to advance, wind-tunnel testing of aeroelastic Froude-scaled models remained crucial for validating these methods. Steps were also taken to establish formal criteria for predicting and mitigating the onset of flutter before the aircraft’s first flight.

The development of supersonic aircraft, including military (e.g., SR-71) and civil (e.g., Concorde), as well as hypersonic vehicles that experience thermal heating (e.g., the X-15), has further advanced aeroelastic research. Comprehensive computational methods emerged in the 1980s, enabling more confident predictions of aeroelastic effects and flutter. Today, computational tools such as Computational Fluid Dynamics (CFD) and the Finite Element Method (FEM), which can be coupled to include structural dynamics, are integral to most comprehensive aeroelastic analyses. These methods enable accurate simulations of aeroelastic interactions among airframe components, particularly for complex geometries and high-speed flight conditions. Unlike analytical solutions, which are more limited, CFD and FEM integration enables engineers to analyze the full spectrum of aeroelastic behavior for an entire flight vehicle, including wings, control surfaces, engines, undercarriage, and other components. Aeroelasticity is not solely a concern for airplanes; it also plays a crucial role in the design of helicopter rotors, wind turbines, and even bridges.

As flight speeds continue to increase and airplane construction evolves toward even lighter and more flexible aerostructures, understanding and controlling aeroelastic effects remain an ongoing engineering challenge. Using active control, either structural or aerodynamic, or both, can help mitigate undesirable aeroelastic effects with lightweight structures. For example, the Boeing 787 has a Flaps-Up Vertical Modal Suppression (F0VMS) system^[3] to enhance the damping of a stable but lightly damped low-frequency aeroelastic mode between the engines, wings, and fuselage. This system uses flight-control surface inputs, including elevators and flaperons, to augment damping rather than suppress a flutter mode. In the FAA special conditions for the 787-10, the agency stated that it was not prepared to accept an active flutter-suppression system that suppresses flutter within the airplane’s operational or design envelope.

Unsteady Airfoil Behavior

Fundamental to predicting flutter is the modeling of unsteady aerodynamics. Unsteady airfoil theory addresses the aerodynamic behavior of airfoils and wings under time-varying boundary conditions, such as oscillations in angle of attack, vertical and horizontal motions, or other situations in which the angle of attack varies with time.

In an unsteady flow, the wake behind an airfoil is time-varying in structure from the shedding of vorticity and circulation, as illustrated in Figure 4, thereby influencing the aerodynamic forces experienced by the airfoil. These wake effects alter the magnitude and phase of the lift response relative to the time-varying boundary conditions. Furthermore, at high angles of attack, the flow may separate from the airfoil surface, causing the leading-edge vortex to shed and resulting in aerodynamic hysteresis, commonly known as dynamic stall. Under these conditions, reduced aerodynamic damping can lead to another form of flutter, known as stall flutter.

Several methods exist for quantitatively calculating these unsteady effects. One instructive approach for those first learning about unsteady airfoil behavior is the classic linearized unsteady aerodynamics theories, including Theodorsen’s theory and the indicial response method. Unsteady panel methods discretize the airfoil surface into panels that extend into the wake, thereby capturing the same physical phenomena. The most computationally intensive CFD methods solve the Navier-Stokes equations for unsteady flow, capturing detailed phenomena such as wakes, vortex shedding, and dynamic stall, but at a high computational cost and time.

Reduced Frequency

A nondimensional parameter used to categorize whether a flow is steady or unsteady is the reduced frequency, denoted by the symbol . The reduced frequency is defined, in general, as

(1)  

where is a characteristic physical frequency of the unsteady flow, is a characteristic length scale, and is a reference flow velocity. The units must be used consistently so that is dimensionless. In practice, is normally expressed in radians per second, with and in compatible units.

For a wing or airfoil, such as one oscillating in angle of attack or an oscillatory vertical (heaving) motion as shown in Figure 5, the reduced frequency is often defined in terms of its semi-chord, i.e., , and the freestream velocity, , i.e., in this case, it is defined as

(2)  

The aerodynamics of non-steady and oscillating airfoil sections and wings must be addressed in the field of aeroelasticity and flutter. If not, flutter predictions will be inadequate. For example, a torsional wing motion introduces a pitching behavior at any wing section. Likewise, a bending motion of the wing introduces a plunging or heaving behavior at the wing section. Both types of motion alter the effective angle of attack and contribute to the aerodynamic loads. Increasing reduced frequency generally produces a greater departure from quasi-steady behavior, including changes in amplitude and phase. It does not necessarily increase the total aerodynamic load, because wake effects may attenuate the circulatory response while velocity- and acceleration-dependent terms may become more important.

For = 0, the flow is steady, and all the usual steady results will apply regarding the relationships between the aerodynamic quantities and the angle of attack. For 0 0.05, the flow can be considered quasi-steady; that is, unsteady effects are generally minor, and for some problems, they may be neglected completely. Under these conditions, the aerodynamic response is directly related to the instantaneous forcing; therefore, no significant time-history effect is observed. Nevertheless, additional noncirculatory apparent-mass terms associated with the airfoil section’s velocities and accelerations may also need to be included. These are unsteady local-in-time terms rather than quasi-steady aerodynamic terms, although they may be added to a quasi-steady approximation for the circulatory loading.

Flows with characteristic reduced frequencies of 0.05 and above are typically considered as unsteady, so all unsteady terms must be retained in the governing flow equations, including the effects of the shed wake. The manifestation of this is more significant changes in the frequency-dependent amplitude and phase responses of the aerodynamic loads with respect to motion. Such aerodynamic and aeroelastic problems are more challenging to model for two reasons. First, the local flow properties also depend on the previous time, i.e., on the history of the lift and other aerodynamic forces. Second, the effects of compressibility may need to be accounted for even when the freestream flow velocity is low, introducing additional challenges in modeling the unsteady aerodynamic response.

Quasi-Steady Aerodynamics

In quasi-steady aerodynamics, the behavior can be evaluated by applying the principles of steady flow under the instantaneous boundary conditions of flow tangency to the airfoil surface. Different forcing conditions, such as pure plunge, pitch, and a combination of pitch and plunge, may be considered, as shown in Figure 6.

For example, for a thin airfoil,^[4] the changes in the angle of attack of a harmonically oscillating airfoil can be expressed as

(3)  

where is the mean angle of attack and is the angle of attack amplitude of the oscillation. If the flow were quasi-steady, in the sense that flow adjustments were to take place instantaneously, the circulatory part of the lift coefficient would be given by

(4)  

recognizing that the pitch rate, i.e., , must affect the lift coefficient because it changes the angle of attack, i.e., to satisfy flow tangency at the rear-neutral point (3/4-chord). This equation then leads to

(5)  

where is a frequency-dependent phase angle given by

(6)  

Therefore, even based on quasi-steady arguments, the harmonic part of the circulatory lift response will no longer be in phase with the angle of attack; it leads by the angle in this case, and its amplitude is larger than the quasi-steady harmonic value by the factor . Notice that because the aerodynamic center is at the 1/4 chord for a thin airfoil, the incremental circulatory lift produced by the angle-of-attack and pitch-rate terms produces no incremental pitching moment about the 1/4-chord point. For a symmetric thin airfoil, this gives . For a cambered thin airfoil, however, a steady quarter-chord pitching moment may still be present and must be added separately.

Apparent Mass Terms

The apparent mass or “noncirculatory” terms^[5] account for the pressure forces required to accelerate the fluid near the airfoil. For example, consider a thin airfoil of chord moving normal to its surface at velocity ; the noncirculatory fluid force, , acting on the surface is

(7)  

The term is often known as the apparent mass or added mass, and in this case is given by

(8)  

Therefore, with the sign convention used here, the noncirculatory lift per unit span for a motion where the airfoil moves normal to its surface with velocity , as shown in Figure 7, is written as

(9)  

In terms of lift coefficient, then

(10)  

As also shown in Figure 7, an additional contribution arises when the airfoil undergoes pitching motion. The corresponding noncirculatory lift contribution depends on the location of the pitch axis as well as on the sign convention used for positive pitch, positive plunge, and positive lift. Therefore, it is better to write the pitch contribution first in its more general form. For a two-dimensional airfoil section, with and with locating the pitch axis relative to the mid-chord in semi-chord units, the noncirculatory lift per unit span may be written as

(11)  

for the sign convention used here. In coefficient form, this becomes

(12)  

For pitching about the quarter-chord, , so the pitching contribution becomes

(13)  

with the same sign convention. The corresponding moment about the leading edge from the pitching motion is

(14)  

In terms of the moment coefficient, then

(15)  

The moment about the 1/4-chord is related to the leading-edge moment by the application of statics using

(16)  

The noncirculatory aerodynamic terms introduce local-in-time forces associated with the airfoil section’s motion. These terms may be proportional to plunge acceleration, pitch rate, and pitch acceleration, depending on the adopted pitch-axis location and sign convention. Unlike the circulatory air loads, which depend on the effective angle of attack and the history of vorticity shed into the wake, the noncirculatory forces are instantaneous terms associated with the local velocity and acceleration of the airfoil section and the resulting acceleration of the surrounding fluid. These effects become increasingly important in high-frequency oscillations and in rapid transient motions, such as airplanes undergoing rapid pitch maneuvers.

Quasi-Steady Aerodynamics, Including Apparent Mass

Using the previous oscillatory pitching example and identifying the angle-of-attack variation with a small pitching motion about the quarter-chord, i.e., with when the pitch axis is measured from the mid-chord in semi-chord units, the noncirculatory contribution to the lift contains both pitch-rate and pitch-acceleration terms. Therefore, it is proportional to the harmonic part of the motion and its time derivatives, not to the mean angle of attack. For

(17)  

then

(18)  

and

(19)  

so the noncirculatory contribution may be written as

(20)  

or, using ,

(21)  

Adding this term to the previous quasi-steady circulatory lift expression gives

(22)  

Substituting leads to and . Therefore, using , the total lift coefficient is

(23)  

which, after some rearrangement, becomes

(24)  

where

(25)  

These results show that the apparent-mass terms modify both the amplitude and phase of the lift response, with the changes depending on reduced frequency. The pitch-rate apparent-mass term contributes a component proportional to , while the pitch-acceleration term contributes a component proportional to . Therefore, even in this quasi-steady-plus-apparent-mass approximation, the aerodynamic force need not be exactly in phase with the prescribed angle of attack. In a complete unsteady aerodynamic treatment, the shed wake must also be included; Theodorsen’s function represents this wake effect and generally attenuates and phase-shifts the circulatory part of the lift response.

Frequency Domain Theories

Although an extensive literature exists on unsteady airfoil aerodynamics, a physically correct theory must account for the vorticity shed from the trailing edge and convected downstream, which feeds back into the airfoil loads and introduces a history-dependent response. When the motion or excitation is harmonic in time, all unsteady quantities vary as , and the governing equations reduce to algebraic relations in the frequency domain. The effects of the shed wake then appear through complex, frequency-dependent transfer functions. Within thin-airfoil theory, this framework leads to two classical results: Theodorsen’s function for harmonic airfoil motion and Sears’s function for a convected sinusoidal gust.

Theodorsen’s Theory

Theodorsen’s theory^[6] is the classical frequency-domain theory for the unsteady aerodynamic loading on a thin airfoil undergoing small-amplitude harmonic motion in incompressible flow. It is one of the foundational results in aeroelasticity because it gives the amplitude and phase of the circulatory lift and moment when an airfoil pitches, plunges, or undergoes a combination of these motions.

The theory is based on linearized thin-airfoil assumptions, small disturbances, harmonic motion, and an incompressible freestream. Therefore, it should be regarded as a low-Mach-number baseline theory. At appreciable Mach numbers, compressibility modifies both the circulatory wake response and the noncirculatory apparent-mass response. Within its assumptions, however, Theodorsen’s theory gives an exact analytical result for the circulatory part of the unsteady thin-airfoil problem.

The central physical point is that an airfoil cannot change its bound circulation without shedding vorticity from the trailing edge. This shed vorticity is convected downstream in the wake, but it continues to induce velocities back on the airfoil. Therefore, the circulatory airload is not determined only by the instantaneous airfoil motion. It also depends on the phase-lagged influence of the wake. Theodorsen’s function is the frequency-domain transfer function that represents this wake effect.

Wake Vorticity & Bound Circulation

In thin-airfoil theory, the bound vorticity distribution on the chord is denoted by , and the shed vorticity in the wake is denoted by . The downwash on the airfoil chord is produced by both the bound vorticity and the wake vorticity, so that the governing integral equation may be written as

(26)  

where is the downwash required to satisfy the flow-tangency condition on the airfoil chord. The Kutta condition requires the bound vorticity to vanish at the trailing edge, i.e.,

(27)  

The total bound circulation is

(28)  

and conservation of circulation requires a change in bound circulation to be accompanied by vorticity shed into the wake. If the shed wake is convected downstream at the freestream speed , then

(29)  

This relationship is the origin of the wake-hereditary effect. Harmonic airfoil motion produces a harmonic change in bound circulation, which produces a harmonic shed wake, which then feeds back into the airfoil loading.

Theodorsen’s Function

For harmonic motion, the unsteady quantities vary as . The wake contribution can then be represented by a complex function of reduced frequency, called Theodorsen’s function. The reduced frequency is

(30)  

where is the airfoil chord, is the freestream speed, and is the angular frequency of the motion.

Theodorsen’s function is

(31)  

where and are Hankel functions of the second kind. With

(32)  

where and are Bessel functions of the first and second kind, respectively, the real and imaginary parts may be written as

(33)  

and

(34)  

where each Bessel function has the argument .

The magnitude and phase of Theodorsen’s function are

(35)  

The magnitude gives the attenuation of the circulatory response relative to the corresponding quasi-steady circulatory response. The phase angle gives the phase shift introduced by the shed wake.

At , Theodorsen’s function approaches unity, and the steady thin-airfoil circulatory result is recovered. As increases, the circulatory response is attenuated and phase-shifted because of the influence of the shed wake. In the limit as , approaches , so the circulatory part of the lift approaches one-half of the corresponding quasi-steady circulatory value.

Circulatory Lift from Pitch & Plunge

Theodorsen’s function applies to the circulatory component of the unsteady airloads. This is the part associated with bound circulation on the airfoil and vorticity shed into the wake. It does not include the noncirculatory or apparent-mass terms, which arise from the pressure field required to accelerate the surrounding fluid.

For a thin airfoil undergoing small-amplitude plunge and pitch , the lift coefficient is written as the sum of circulatory and noncirculatory parts, i.e.,

(36)  

where is the circulatory contribution and is the noncirculatory or apparent-mass contribution.

For the sign convention used here, and for pitching about the quarter-chord, the circulatory lift coefficient may be written in the compact form

(37)  

The bracketed quantity is the effective angle of attack that drives the circulatory loading. It contains the pitch angle, the plunge-velocity contribution, and the pitch-rate contribution associated with the motion-induced boundary condition. Theodorsen’s function then modifies this circulatory response by adding the wake-induced amplitude attenuation and phase shift.

The pitch-rate term inside the brackets is still part of the circulatory response. It should not be confused with the apparent-mass terms. The apparent-mass terms are local inertial terms and are added separately.

Representative results from Theodorsen’s theory are shown in Figure 10. For , the steady thin-airfoil circulatory result is recovered. As increases, the circulatory lift no longer follows the motion instantaneously, and a hysteresis loop appears in the -versus- plot. The orientation and shape of the loop depend on the sign convention used for the harmonic motion and on whether the plotted quantity is only the circulatory lift or the total lift including apparent mass.

Adding the Apparent-Mass Terms

The noncirculatory contribution is associated with the pressure field required to accelerate the fluid near the airfoil. It is instantaneous in the sense that it depends on the local velocity and acceleration of the airfoil section, not on the history of vorticity shed into the wake. Therefore, the noncirculatory terms are not multiplied by .

For a thin airfoil with , and with locating the pitch axis relative to the mid-chord in semi-chord units, the noncirculatory lift contribution may be written in coefficient form as

(38)  

for the sign convention used here. For pitching about the quarter-chord, , so the pitching contribution becomes

(39)  

For harmonic pitching motion,

(40)  

the derivatives are

(41)  

Using , the quarter-chord pitching apparent-mass contribution becomes

(42)  

This contribution is added to the circulatory lift. It is not multiplied by .

Total Lift for Quarter-Chord Pitching

For harmonic pitching about the quarter-chord with no plunge, the circulatory contribution is

(43)  

and the noncirculatory contribution is

(44)  

Therefore, the total lift coefficient is

(45)  

for the sign convention used here.

This equation is the practical form needed for many introductory aeroelastic calculations involving quarter-chord pitching. The first term is the circulatory lift modified by Theodorsen’s function. The last two terms are the apparent-mass lift terms and are added directly. If plunge is also present, then the plunge-velocity contribution is included in the circulatory effective angle of attack, and the plunge-acceleration contribution is included in the noncirculatory part.

The preceding example computed only the circulatory part of the unsteady lift response, i.e., the part modified by Theodorsen’s function. The next example adds the noncirculatory apparent-mass contribution to obtain the total unsteady lift coefficient.

Sinusoidal Gust: Sears’s Problem

Von Kármán & Sears^[7] analyzed the unsteady lift of a thin airfoil encountering a sinusoidal vertical gust convected by the free stream, as shown in Figure 11. The gust is represented as a prescribed upwash velocity field of the form

(46)  

where is the gust amplitude, is the gust (temporal) frequency, and is the free-stream speed. The corresponding gust wavelength is

(47)  

Using the identity then Eq. 46 may be written as

(48)  

Two reference conventions are common. If the gust is referenced to the leading edge, then and Eq. 46 reduces to . If the gust is referenced to the mid-chord, then and the forcing becomes , i.e., a frequency-dependent phase shift, where the (dimensionless) gust reduced frequency is

(49)  

The mid-chord convention was used in the original work of von Kármán & Sears, and their use of the mid-chord reference point has caused much confusion in the literature. For a thin airfoil, the harmonic lift response to the convected gust may be written in transfer-function form as

(50)  

where is called Sears’s function. An exact representation is

(51)  

where and are Bessel functions of the first kind and is Theodorsen’s function. Writing , and , Eq. 51 gives

(52)  

as plotted in Figure 12. The results can be calculated numerically, as in the case of the Theodorsen function. Sears’s function describes the complete linear unsteady lift response of a fixed thin airfoil to a prescribed convected sinusoidal gust under the assumptions of incompressible thin-airfoil theory. In the underlying general theory, the gust-induced lift may be decomposed into apparent-mass, quasi-steady, and wake-dependent contributions. The standard Sears function combines the contributions required for the complete sinusoidal-gust solution, so it should not be described as purely circulatory.

The apparent-mass terms derived separately for prescribed pitching or plunging motion should not be added to the Sears-function result, because those terms correspond to a different forcing mechanism and would double-count parts of the unsteady response.

Sears’s original convention references the sinusoidal gust to the mid-chord. With the convention used above,

the gust at the mid-chord lags the gust at the leading edge by , where . This mid-chord convention is the one associated with the familiar spiral locus of Sears’s function in the complex plane. If the same physical lift response is instead written relative to the leading-edge gust signal, then the transfer function differs by the corresponding convective phase factor. If is defined relative to the mid-chord gust signal, then

(53)  

Equivalently, if and , then

(54)  

The important point is that Sears’s function must always be used with the same chordwise reference point used to define the gust input. Mixing mid-chord and leading-edge conventions introduces a spurious phase error.

At low reduced frequencies, and the gust response becomes quasi-steady. At high reduced frequencies, Sears’s function decays in magnitude, with the asymptotic behavior given by

(55)  

with a phase that depends on the chosen reference point. Notice again that Sears’s original mid-chord convention is the one associated with the characteristic spiral locus in the complex plane. A leading-edge-referenced transfer function is obtained only by applying the appropriate convective phase factor, as in Eq. 53.

The Theodorsen and Sears functions describe different physical problems and are not interchangeable. Theodorsen’s function modifies the circulatory lift produced by prescribed airfoil motion, while the corresponding motion-induced noncirculatory apparent-mass terms must be added separately. Sears’s function gives the complete linear lift response caused by a sinusoidal vertical gust convected over a fixed airfoil. Therefore, no additional apparent-mass term should be formed from the gust velocity and added to the Sears-function response. If the airfoil is also pitching, plunging, or otherwise accelerating while encountering the gust, the apparent-mass terms caused by that independent airfoil motion must still be included. Under the assumptions of linear thin-airfoil theory, the motion-induced and gust-induced responses may then be superposed.

Time-Domain: The Indicial Response

The indicial aerodynamic response method provides a time-domain framework for calculating unsteady aerodynamic forces and moments. The indicial response refers to the system’s response to a unit step input applied at and held constant thereafter. Like all classic thin-airfoil theories, the response relies on the linearization of aerodynamic forces. For a step change in angle of attack, , from an initial angle , the final steady lift increment would be , but the transient circulatory lift does not reach this value immediately. Instead, the incremental circulatory indicial response may be written as

(56)  

or, including the initial steady circulatory lift,

(57)  

where is the circulatory indicial-response function. The noncirculatory apparent-mass contribution, if present, must be added separately.

For arbitrary variations in , Duhamel’s convolution integral applies to the circulatory part of the lift, i.e.,

(58)  

where is a dummy time variable of integration. In this dimensional-time form, the indicial-response function is understood to depend on the elapsed convective time after the disturbance. Therefore, the solution to the Duhamel integral accounts for the time-history effects of arbitrary forcing on the circulatory lift. Any noncirculatory apparent-mass contribution must be added separately.

Non-Dimensional Time

The nondimensional or reduced time is always used for transient problems where the reduced frequency becomes ambiguous. The reduced time is defined as

(59)  

where is the freestream velocity, so that represents the relative distance traveled by the airfoil through the flow in units of semi-chords. This parameter facilitates the analysis of unsteady flow problems that discrete frequencies cannot characterize. It will be apparent that the reduced time variable is the same as the distance traveled by the airfoil through the flow expressed in semi-chords.

Wagner’s Problem

Wagner determined the transient lift response to a step change in angle of attack.^[8] In Wagner’s problem, the airfoil experiences a sudden change in angle of attack caused by its own motion. The response accounts for the effects of the shed wake through Wagner’s function, , where

(60)  

is the reduced time, i.e., the distance traveled by the airfoil through the flow in units of semi-chords.

For a step change in angle of attack, , the lift coefficient increment may be written formally as the sum of a noncirculatory impulsive term and a circulatory term, i.e.,

(61)  

where is the Dirac delta function in nondimensional time. The first term is a distributional impulse associated with the apparent-mass response at the instant of the step; it is not a finite sustained lift coefficient and contributes only when integrated through the impulse. The second term represents the increment in circulatory lift that develops as vorticity is shed into the wake.

The Wagner function can be calculated exactly in terms of special functions and is shown in Figure 13. Because builds asymptotically from to as , the circulatory component of the lift increment increases from to The impulsive contribution at is separate and represents the noncirculatory apparent-mass response. Otherwise, Wagner’s function represents the gradual growth of circulatory lift as the starting vortex is shed into the downstream wake and circulation develops around the airfoil.

For arbitrary variations in angle of attack, the circulatory lift can be computed using Duhamel’s integral. In terms of nondimensional time , the circulatory lift coefficient is

(62)  

where is a dummy reduced-time variable. This expression gives only the circulatory part of the lift, i.e., the part associated with the gradual development of circulation and the shed wake. The noncirculatory apparent-mass terms are added separately as local velocity- or acceleration-dependent terms. For practical calculations, Wagner’s function is often approximated using a two-term exponential form, i.e.,

(63)  

where, for the common two-term approximation, , , , and . Therefore, Wagner’s function is the time-domain indicial response for airfoil-motion forcing. Recall that the corresponding frequency-domain response for harmonic airfoil motion is Theodorsen’s function.

Sharp-Edged Gust: Küssner’s Problem

The gust-response counterpart to Wagner’s problem is Küssner’s problem.^[9] In Wagner’s problem, the unsteadiness is generated by airfoil motion. In Küssner’s problem, the airfoil is fixed and encounters a sharp-edged vertical gust convected by the freestream, i.e., the unsteadiness is imposed by the flow field rather than generated by the motion of the airfoil.

For a sharp-edged gust of vertical velocity , the effective gust angle is

(64)  

for small angles. The corresponding gust-induced unsteady lift response may be written as

(65)  

where is Küssner’s function and is the nondimensional time. Unlike Wagner’s function, which begins at , Küssner’s function begins at and approaches unity as , i.e., . Therefore, the lift response to a sharp-edged gust builds from zero to the quasi-steady value, i.e.,

(66)  

The Küssner function can be obtained from the classical incompressible thin-airfoil gust problem and is shown in Figure 14.

For engineering calculations, Küssner’s function is often represented by a two-term exponential approximation of the form

(67)  

with constants chosen to match the exact or tabulated response over the range of reduced times of interest. A common approximation is , and , which satisfies the correct limiting values at and , but more accurate sets of coefficients may be used when better agreement with the classical Küssner function is required.

For an arbitrary gust history, the gust-induced lift can be computed by applying Duhamel’s integral to the gust angle, i.e.,

(68)  

where

(69)  

and is again a dummy reduced-time variable.

Therefore, it will be appreciated that Küssner’s function is the time-domain indicial response for a gust input, just as Wagner’s function is the time-domain indicial response for a sudden change in airfoil angle of attack. The corresponding frequency-domain response for a sinusoidal gust is Sears’s function. The two viewpoints are complementary: Küssner’s function is used for arbitrary gust histories in the time domain, whereas Sears’s function is used for harmonic gusts in the frequency domain.

Connection Between Apparent Mass & Indicial Response

The distinction between circulatory and noncirculatory lift can also be understood from the indicial-response viewpoint. Previously, the noncirculatory or apparent-mass terms were introduced as local velocity- and acceleration-dependent contributions that are added to the circulatory lift. These terms are not multiplied by Theodorsen’s function because they are not wake-hereditary effects. In the time domain, the same physical contribution appears as the impulsive part of the indicial response.

In Wagner’s problem, the circulatory lift following a step change in angle of attack develops gradually as vorticity is shed into the wake. This gradual development is represented by Wagner’s function, , where is the usual nondimensional time. The noncirculatory response, by contrast, is instantaneous because it is associated with the pressure field that accelerates the fluid around the airfoil. Therefore, it appears mathematically as an impulse.

The point can be made explicitly by considering a sudden step change in angle of attack, i.e.,

(70)  

where is the Heaviside step function. This function is zero before the step and one after the step, i.e.,

(71)  

Therefore, represents an angle of attack that jumps suddenly from zero to at . The derivative of the Heaviside function is the Dirac delta function, so

(72)  

The delta function is not an ordinary finite function at . It is defined by its integral, namely

(73)  

for any . Therefore,

(74)  

which recovers the finite jump in angle of attack. This is the mathematical meaning of an impulsive response: the instantaneous value is singular, but the integrated effect through the impulse is finite.

It is important to distinguish between an indicial-response kernel and the lift itself. Let the noncirculatory part of the indicial-response kernel be written as

(75)  

where is the appropriate apparent-mass coefficient for the motion being considered. The symbol denotes the kernel, not the lift coefficient. The corresponding noncirculatory lift coefficient is obtained from Duhamel’s integral, i.e.,

(76)  

Substituting gives

(77)  

Using the sifting property of the delta function,

(78)  

so that

(79)  

This result shows how the impulsive part of the indicial response becomes a local derivative term for a general smooth motion.

The angle-of-attack example illustrates how a delta-function kernel reduces to a local derivative term. For an actual airfoil section, the same mechanism acts on the motion variables that create normal velocity on the airfoil. Therefore, for a two-dimensional airfoil with semi-chord , freestream speed , plunge displacement , and pitch angle , the effective angle of attack contains contributions from both pitch and vertical motion. In a small-disturbance approximation, the plunge contribution gives an effective angle of attack that contains a term of the form

(80)  

with the sign depending on the convention used for positive plunge. Therefore, when the noncirculatory indicial response acts on the rate of change of the effective angle of attack, it naturally produces a term proportional to

(81)  

This is the time-domain mechanism by which the plunge-acceleration apparent-mass term appears. Physically, the airfoil must accelerate a surrounding mass of fluid as it plunges, so the resulting pressure field produces an instantaneous lift contribution proportional to .

Pitch produces two related noncirculatory effects. First, pitch rate changes the local normal velocity distribution over the chord, giving an apparent-mass contribution proportional to Second, pitch acceleration requires angular acceleration of the surrounding fluid, giving a contribution proportional to with a coefficient that depends on the location of the pitch axis. Therefore, the classical noncirculatory lift per unit span may be written as

(82)  

where locates the pitch axis relative to the mid-chord in semi-chord units. This form is written for a general pitch-axis location; simpler special cases follow by choosing the pitch-axis convention and sign definitions used in a particular example. The precise signs of the terms depend on the adopted sign convention for positive plunge, positive pitch, positive lift, and the definition of . However, the structure of the result is the important point: the noncirculatory lift is local in time and contains apparent-mass terms proportional to plunge acceleration, pitch rate, and pitch acceleration.

These terms are the physical counterparts of the delta-function contribution in the indicial response. The term is the apparent-mass lift associated with accelerating the surrounding fluid in plunge. The term is associated with the rotational velocity field produced by pitch rate in a forward stream. The term is associated with angular acceleration about an axis displaced from the mid-chord. These contributions arise from the instantaneous pressure field needed to accelerate the surrounding fluid.

The connection with nondimensional indicial time follows by introducing

(83)  

so that

(84)  

and

(85)  

Therefore,

(86)  

and

(87)  

If the nondimensional plunge displacement is defined as then

(88)  

Therefore, the dimensional apparent-mass terms , , and correspond in nondimensional indicial form to terms involving , , and For harmonic forcing, for example, , where is the nondimensional frequency based on , then

(89)  

and

(90)  

First-derivative apparent-mass terms therefore appear in the frequency domain as algebraic factors proportional to , while second-derivative apparent-mass terms appear as factors proportional to . The same reasoning applies to plunge. A harmonic plunge displacement produces in nondimensional form, and so the plunge contribution appears as an algebraic acceleration term rather than as a wake-hereditary term.

Consequently, the delta-function part of the indicial response is the time-domain representation of the same apparent-mass physics that appears in the classical unsteady airfoil equations as local terms proportional to , , and . The circulatory part contains the wake effects and is governed by Wagner’s function in the time domain or Theodorsen’s function in the frequency domain. The noncirculatory part has no hereditary effects and appears instead as instantaneous velocity- and acceleration-dependent terms.

Time- Vs. Frequency-Domain Functions

By now, it will be apparent that Theodorsen’s, Sears’s, Wagner’s, and Küssner’s functions are closely related, but they describe different forms of unsteady aerodynamic forcing. Theodorsen’s and Sears’s functions are frequency-domain transfer functions for harmonic forcing, whereas Wagner’s and Küssner’s functions are time-domain indicial response functions for step-type forcing.

For airfoil motions, Theodorsen’s function gives the circulatory lift response to harmonic pitching or plunging motion. In the time domain, the corresponding indicial response is Wagner’s function , which describes the growth of circulatory lift after a sudden change in angle of attack. Therefore, Theodorsen’s and Wagner’s functions represent the same circulatory wake physics viewed from two different domains, i.e.,

(91)  

where is the nondimensional time and is the reduced frequency.

For gust excitations, Sears’s function gives the circulatory lift response of a fixed airfoil to a sinusoidal vertical gust convected by the freestream. In the time domain, the corresponding indicial response is Küssner’s function , which describes the growth of lift after the airfoil encounters a sharp-edged vertical gust. Therefore,

(92)  

The distinction between the two pairs is physical. Wagner’s and Theodorsen’s functions apply when the unsteadiness is generated by airfoil motion. Küssner’s and Sears’s functions apply when the flow field imposes the unsteadiness as a gust. In both cases, the time-domain and frequency-domain descriptions are complementary. The indicial response gives the lift history following a step input, while the frequency-domain transfer function gives the amplitude and phase response to harmonic forcing.

In practical calculations, arbitrary motions or gust histories may be treated in the time domain by convolution with the appropriate indicial response function. Harmonic motions or gusts may be treated in the frequency domain using the corresponding complex transfer function. For prescribed airfoil motion, local noncirculatory apparent-mass terms are added separately from the circulatory Theodorsen or Wagner response. For the standard gust problems, Sears’s and Küssner’s functions represent the corresponding complete linear gust responses under their stated assumptions, so the motion-induced apparent-mass terms should not be added separately.

Numerical Superposition Using the Indicial Response

The indicial response gives the aerodynamic response to a unit step input. Duhamel’s integral is the corresponding superposition integral for an arbitrary time-varying input. It says that the response at the current time is obtained by adding the delayed effects of all previous small changes in the input. In other words, each increment in angle of attack, pitch, plunge-induced angle, or gust angle produces its own indicial response, and the total circulatory response is the sum of those contributions.

For a generic input angle , where is reduced time, this idea may be written schematically as

where is the appropriate indicial response function. For airfoil-motion inputs, is usually Wagner’s function. For a sharp-edged gust input, is usually Küssner’s function. If the input is formulated as an increment from a trimmed initial condition, then and only the convolution integral remains. This expression applies only to the circulatory lift response; the noncirculatory apparent-mass terms are added separately.

In practice, directly evaluating this integral at each time step can be inefficient because the previous input history must be retained and repeatedly integrated. A more efficient approach is to approximate the indicial response function by a sum of exponentials. This converts Duhamel’s integral into a small number of aerodynamic lag states that can be advanced in time using recurrence equations.

Consider a general two-term approximation to an indicial response function, i.e.,

(93)  

where may represent Wagner’s function , Küssner’s function , or another approximate indicial response function. Recall that for Wagner’s function, typical constants are , , , and , whereas for the simple Küssner approximation used above, the constants are , , , and . For a generic input angle , the effective angle associated with the circulatory response may be written as

(94)  

where and are aerodynamic lag states associated with the two exponential terms in the chosen indicial response function. The corresponding circulatory lift coefficient is then

(95)  

for a thin airfoil in incompressible flow.

The lag states are obtained from the Duhamel superposition. For example, the first lag state has the form

(96)  

with an analogous expression for . This expression shows that the lag states respond to changes in the input, not merely to the instantaneous value of the input.

For a piecewise-linear change in the input angle over the interval , a practical recurrence form is

(97)  

and

(98)  

where

(99)  

so that the effective angle at the new time step is

(100)  

The aerodynamic wake history is contained in the lag states and . For sufficiently small , the bracketed factors approach unity, giving the simpler limiting update in which each input increment contributes approximately or to the corresponding lag state.

For a true step input, the jump in must be treated as an input discontinuity. Immediately after a step change , the lag states receive the finite jumps

(101)  

so that the circulatory response begins at

(102)  

and then approaches the steady value as the lag states decay for the held final value of the input. This result recovers the usual indicial response behavior.

For airfoil-motion forcing, the input angle may be denoted by , and the Wagner-function constants are used. The corresponding circulatory lift contribution is

(103)  

where

(104)  

For gust forcing, the input angle may be written as

(105)  

and the Küssner-function constants are used. The corresponding gust-induced circulatory lift contribution is

(106)  

where

(107)  

For a general unsteady problem, both airfoil motion and gust excitation may be present simultaneously. Under the assumptions of linear thin-airfoil theory, their circulatory contributions may be superposed, so that

(108)  

or

(109)  

The noncirculatory terms are not included in the lag states because they have no wake-hereditary effect. They are evaluated directly from the instantaneous velocity or acceleration of the airfoil section at the current time step. Therefore, the total lift coefficient is assembled as

(110)  

or, equivalently,

(111)  

For example, for small-amplitude plunge and pitch motions of a thin airfoil, the noncirculatory lift contribution may be written in coefficient form as

(112)  

where locates the pitch axis relative to the mid-chord in semi-chord units. The signs of these terms must be interpreted according to the chosen definitions of positive , positive , positive lift, and the sign convention used for .

These local noncirculatory terms are evaluated directly from the imposed motion or from the structural equations of motion. If the motion history is prescribed or has already been computed, then finite-difference approximations may be used. For example, with known values on both sides of , a centered approximation is

(113)  

and

(114)  

In a forward time-marching aeroelastic solution, however, and may not yet be known when the loads are assembled at . In that case, the noncirculatory terms must be evaluated consistently with the chosen structural time-integration scheme, for example, by using the acceleration states solved at the current step or by using an appropriate backward or implicit approximation. Therefore, the numerical procedure is to compute the circulatory wake effects from the appropriate Wagner or Küssner lag-state recurrence and then add the local noncirculatory contribution without passing it through the wake lag states.

State-Space Form

The recurrence lag-state formulation developed above is well-suited to direct time-marching. The same Duhamel superposition can also be written in continuous state-space form, which is often more convenient for aeroelastic analysis because the aerodynamic lag states can be appended directly to the structural equations of motion.

Using the same lag-state definitions as in the preceding recurrence formulation, the states associated with the two exponential terms satisfy

(115)  

and

(116)  

where denotes the input angle for the particular forcing mechanism, such as airfoil motion or gust excitation. The corresponding effective angle remains

(117)  

so that, for a thin airfoil in incompressible flow,

(118)  

The initial values of the lag states must be set consistently with the initial input. For example, if the input has a finite initial value , then the corresponding lag-state initial conditions are

(119)  

unless the problem is formulated in terms of increments from a trimmed initial condition.

In matrix form, the aerodynamic lag-state equations may be written as

(120)  

with the output equation

(121)  

If dimensional time is used instead of nondimensional time , then for constant ,

(122)  

and the dimensional-time state equations become

(123)  

and

(124)  

The same state-space structure may be used separately for motion-induced and gust-induced inputs. Motion-induced states use the Wagner-function constants, while gust-induced states use the Küssner-function constants. Under the assumptions of linear thin-airfoil theory, the corresponding circulatory lift contributions are then superposed, i.e.,

(125)  

The total unsteady lift coefficient is obtained only after the noncirculatory terms are added, i.e.,

(126)  

The aerodynamic lag states model the circulatory wake-history effects, whereas the noncirculatory terms represent local inertial effects associated with the instantaneous motion of the airfoil section.

The important practical point is that the indicial-response superposition can be represented either by the recurrence equations developed above or by equivalent aerodynamic lags expressed as state-space equations. The recurrence form is convenient for direct time marching, while the state-space form is convenient for coupling the aerodynamic states to structural equations of motion. Flutter can then be studied by computing the eigenvalues of the complete coupled system or by integrating the coupled equations forward in time.

Compressibility Effects on Unsteady Airfoil Theory

The classical Theodorsen, Sears, Wagner, and Küssner functions discussed above are incompressible-flow results. They are useful reference functions because they reveal the basic physics of unsteady aerodynamic response, including wake history, amplitude attenuation, phase lag, and the distinction between circulatory and noncirculatory loading. There are no simple, exact compressible-flow equivalents of these functions. In particular, compressibility should not be introduced by multiplying Theodorsen’s, Sears’s, Wagner’s, or Küssner’s incompressible functions by a Prandtl-Glauert factor. The incompressible apparent-mass terms also do not apply in compressible flow, because the noncirculatory pressure field has a finite adjustment time due to the propagation of pressure waves.

For subsonic compressible flow, practical methods instead use compressible aerodynamic response models, such as transfer functions, aerodynamic influence coefficients, or approximate indicial functions. These models are obtained from linearized compressible-flow theory, unsteady airfoil calculations, wind-tunnel measurements, or CFD-based system identification. Their Mach-number dependence is built into the frequency-domain response or into the fitted time-domain coefficients; they are not exact compressible versions of the classical incompressible functions.

Approximate compressible indicial functions for two-dimensional subsonic flow have been developed and validated.^[10] In this approach, the unsteady aerodynamic response is obtained from oscillatory lift data over a range of reduced frequencies and free-stream Mach numbers. The frequency-domain response is then related to an equivalent time-domain indicial response using Fourier-transform relationships. The resulting indicial functions are approximate engineering representations, but they are not arbitrary curve fits. They are constructed to reproduce the principal Mach-number effects on both the circulatory and noncirculatory parts of the unsteady aerodynamic response.

The circulatory part represents the lift associated with the buildup of bound circulation and the shed wake. In this formulation, the circulatory indicial response is represented by a two-term exponential approximation of the form

(127)  

where is the reduced time, , and is the free-stream Mach number. The coefficients , , , and are not obtained by rescaling the incompressible Wagner function but from unsteady airfoil data in the frequency domain and then transformed into an equivalent time-domain indicial representation. A least-squares procedure is used to determine a single set of coefficients that best represents the measured or calculated unsteady response. For the examples shown here, the fitted circulatory coefficients are taken as , , , and .

The Mach-number dependence in this compact form enters the exponential arguments through the factor . Physically, compressibility changes the timescale over which the wake-induced pressure field influences the airfoil. In incompressible potential flow, the pressure field is determined instantaneously by the boundary conditions and the shed vorticity. In subsonic compressible flow, pressure disturbances propagate at a finite speed, so the circulatory response develops over a longer effective reduced-time scale. Therefore, as increases and decreases, the effective lag rates and decrease, so the response approaches its final value more slowly in reduced time. Equivalently,

(128)  

For a step change in the effective angle of attack, the corresponding circulatory lift response may be written as

(129)  

This form shows that compressibility affects the circulatory response through both the lift-curve slope and the Mach-number scaling of the indicial response in time.

The same distinction applies to the incompressible apparent-mass terms. In incompressible flow, the apparent-mass contribution is instantaneous because the pressure field adjusts everywhere without any delay. In subsonic compressible flow, these apparent-mass terms do not apply in the same form. The corresponding contribution is still noncirculatory in origin, but it is governed by the finite-time pressure-wave adjustment of the flow field and must be represented by a compressible unsteady aerodynamic model.

For a unit step change in angle of attack, the compressible indicial lift response can be written as the sum of a decaying noncirculatory part and a growing circulatory part, i.e.,

(130)  

or

(131)  

The initial value of the noncirculatory lift response is given by piston theory, i.e.,

(132)  

for a unit step change in angle of attack.^[11] The noncirculatory part is then approximated as

(133)  

where is the nondimensional noncirculatory time constant in reduced time . The prime distinguishes this reduced-time constant from the corresponding dimensional time constant .

The time constant is not obtained by fitting the noncirculatory term in isolation. Instead, it is chosen so that the approximate total indicial response has the correct initial value and initial slope given by exact linearized subsonic theory. This matching gives

(134)  

The corresponding dimensional time constant is

(135)  

or, equivalently,

(136)  

where is the speed of sound.

Therefore, for a step change in effective angle of attack, the total compressible indicial lift response may be written as

(137)  

This expression has three exponential terms: one decaying noncirculatory term and two growing circulatory wake-history terms.

Figure 15 shows the three-exponential compressible indicial lift response for a step change in angle of attack. At , the response is dominated by the piston-theory noncirculatory value . As , the noncirculatory term decays to zero and the lift approaches the steady compressible circulatory value .

This behavior is fundamentally different from applying a Prandtl-Glauert scaling to the incompressible apparent-mass term. The finite-Mach-number noncirculatory response has its own initial value and decay time; in the singular limit , it collapses to the impulsive apparent-mass contribution of incompressible indicial theory.

For the pitching motion about the quarter-chord, it is useful to define the nondimensional pitch-rate input as

(138)  

The lift response to a step change in pitch rate can be written as the sum of a noncirculatory part and a circulatory part, i.e.,

(139)  

For pitching about the quarter-chord, the effective angle of attack associated with the pitch-rate contribution at the three-quarter-chord point is . Therefore, the final circulatory lift per unit is

(140)  

A useful approximate circulatory indicial response is then

(141)  

where the same fitted circulatory coefficients introduced above are used.

The noncirculatory part has a different initial value and a different time constant. For pitch rate about the quarter-chord, the piston-theory initial value per unit is

(142)  

The noncirculatory contribution may be approximated as

(143)  

where is the nondimensional pitch-rate noncirculatory time constant in reduced time . The prime distinguishes this reduced-time constant from the corresponding dimensional time constant .

The pitch-rate noncirculatory time constant differs from the corresponding angle-of-attack value because it is obtained by matching the initial slope of the approximate pitch-rate indicial response to the corresponding exact linearized subsonic result. This matching gives

(144)  

The corresponding dimensional time constant is

(145)  

or, equivalently,

(146)  

Therefore, the lift response to pitch rate about the quarter-chord may be written as

(147)  

This expression has the same three-exponential structure as the angle-of-attack response: one decaying noncirculatory term and two growing circulatory wake-history terms.

Figure 16 shows the three-exponential compressible indicial lift response for a step change in pitch rate about the quarter-chord. At , the response is dominated by the piston-theory noncirculatory value . As , the noncirculatory term decays to zero and the lift approaches the steady compressible circulatory value .

Similar indicial-response forms can also be written for the pitching moment. The moment response has its own circulatory and noncirculatory coefficients and depends on the moment reference point, so the lift-response coefficients above should not be used directly for the moment response.

The main practical point is that incompressible unsteady aerodynamic functions and apparent-mass terms are reference results for incompressible flow. For subsonic compressible calculations, the circulatory response must include the compressible lift-curve-slope scaling and the Mach-number scaling of the indicial response in time. In contrast, the noncirculatory response must be represented by the corresponding compressible indicial model. For higher-fidelity aeroelastic analysis, the aerodynamic model may instead be based on experimentally identified aerodynamic influence coefficients, linearized compressible-flow methods, Euler calculations, or Reynolds-averaged Navier-Stokes CFD coupled to the structural model. The required model depends on the Mach number, reduced frequency, geometry, and whether the flow remains attached and linear or becomes nonlinear because of shocks, flow separation, or other effects.

Low-Mach-Number Limit of the Noncirculatory Response

In incompressible unsteady airfoil theory, a sudden change in effective angle of attack produces an impulsive noncirculatory lift associated with the apparent mass of the surrounding fluid. At finite Mach number, however, the pressure field cannot adjust instantaneously because pressure disturbances propagate at the speed of sound. The corresponding noncirculatory response is then a short acoustic transient rather than a mathematical impulse.

Using the reduced time , the acoustic time scale may be written as . The finite-Mach-number noncirculatory response to a step change in effective angle of attack has the asymptotic form

(148)  

where describes the shape of the short-time acoustic response. As the Mach number decreases, the amplitude of this transient increases in proportion to , while its duration in reduced time decreases in proportion to . At every fixed value of , the transient has already passed as , so its pointwise limit is zero. Nevertheless, its integrated effect remains finite. If

(149)  

then the family of increasingly tall and narrow acoustic transients approaches

(150)  

where is the causal Dirac-delta distribution. This impulse is the familiar incompressible apparent-mass contribution produced by a unit step in effective angle of attack.

Therefore, the incompressible apparent-mass impulse is not the pointwise low-Mach-number limit of the compressible noncirculatory response. It is the distributional limit of a finite acoustic transient whose amplitude increases while its duration decreases. This distinction is also important in approximate indicial and state-space models because matching the finite-Mach-number initial value alone is insufficient. The model must also preserve the integrated area of the short-time noncirculatory response.

Three-Dimensional Unsteady Aerodynamics

The indicial-response methods described above are based on two-dimensional thin-airfoil theory. They are useful because they provide exact solutions for the unsteady lift response that develops over time in response to changes in angle of attack, pitch rate, plunge velocity, or gust upwash. Numerical approximations to the indicial functions give sufficient accuracy for engineering purposes. However, an actual wing is finite, so its unsteady aerodynamic response is not purely two-dimensional. The spanwise lift distribution, induced downwash, tip vortices, sweep, taper, aspect ratio, and wing-body interactions will all affect the unsteady loading.

A practical engineering synthesis is to combine the two-dimensional indicial-response model with a spanwise discretization of the wing. The simplest version is a strip analysis. In this method, the wing is divided into spanwise strips, and each strip is treated as a local two-dimensional airfoil section. The local effective angle of attack is computed from the rigid-body motion, elastic deformation, pitch rate, plunge velocity, and any gust velocity at that spanwise station. A two-dimensional unsteady aerodynamic model, such as a Wagner-function lag-state model for airfoil motion or a Küssner-function lag-state model for gust penetration, is then applied to each strip.

For the th strip, the local motion-induced angle may be written as

(151)  

for the sign convention used here. In this expression, is the local plunge displacement, is the local pitch or torsional rotation, and is the local chord. If a vertical gust is present, then the gust-induced angle may be approximated as

(152)  

where is the gust velocity at the strip. The motion-induced and gust-induced inputs may then be passed through their corresponding indicial-response models and added by linear superposition.

For example, using the lag-state form of the indicial response, the local circulatory lift coefficient on the th strip may be written schematically as

(153)  

where is the local effective angle at each section after the appropriate wake-hereditary lag states have been applied. For example, may include Wagner-function lag states for section motion and Küssner-function lag states for gusts. The corresponding sectional circulatory lift per unit span is

(154)  

The total circulatory lift on the wing is then approximated by summing the strip contributions, i.e.,

(155)  

where is the spanwise width of the strip.

The same idea may be extended to moments and generalized aerodynamic forces. The sectional force or moment contribution is computed on each strip and then summed over the span. For example, a generalized aerodynamic force associated with a structural mode may be approximated as

(156)  

where represents the appropriate modal displacement or rotation at the strip location. This form is often used in elementary aeroelastic calculations because it connects local unsteady airfoil theory to the global bending and torsional response of a finite wing.

The sectional lift used in the strip summation must include both circulatory and noncirculatory contributions, i.e.,

(157)  

where is evaluated locally using the chord, pitch-axis location, and motion variables for the th strip. The wing-level forces, moments, and generalized forces are then obtained by summing the resulting sectional loads over the span.

This strip-based formulation naturally leads to a state-space representation. In this case, for a single strip, or for a collection of strips, the aerodynamic lag states can be assembled into a first-order system of equations of the form

(158)  

and

(159)  

where contains the aerodynamic lag states, contains the motion or gust inputs, and contains the aerodynamic outputs, such as sectional lift, pitching moment, or generalized aerodynamic forces. The matrices , , , and are determined by the chosen indicial-response approximation and by the spatial discretization or modal representation.

This state-space form is especially useful for time marching, stability analysis, flutter prediction, and control-system design. The structural equations of motion can be coupled to the aerodynamic lag-state equations, giving a combined aeroelastic system in which the hereditary effects of the aerodynamics are represented by a finite set of first-order differential equations rather than by an explicit convolution over the entire previous history. In this form, the two-dimensional section theory, the strip summation, and the structural dynamics are all brought into one computational framework.

The strip-analysis approach is useful and relatively inexpensive numerically, but it is still an approximation. It does not fully capture three-dimensional wake deformation, spanwise flow, tip-vortex dynamics, or the mutual aerodynamic influence between neighboring strips. More complete three-dimensional unsteady aerodynamic models use unsteady lifting-line theory, lifting-surface theory, doublet-lattice methods, vortex-lattice methods, panel methods, or CFD coupled with structural dynamics. Nevertheless, strip theory and its state-space form provide a valuable first engineering synthesis: they use well-understood two-dimensional unsteady airfoil theory locally, retain the essential wake-hereditary effects through aerodynamic lag states, and then integrate the resulting sectional loads over the finite wing.

Aeroelastic Analysis

The field of aeroelasticity is not easy to master, but understanding the basics, especially the underlying physics in a qualitative context, is a good place to start. Flutter is a self-excited oscillation caused by the interaction of aerodynamic, elastic, and inertial forces. In a fairly general way, the structural dynamics of an aeroelastic system can be described by

(160)  

where is the mass matrix, is the damping matrix, is the stiffness matrix, is the displacement vector, and is the aerodynamic force vector. Aerodynamic forces can be computed using computational fluid dynamics (CFD), but panel methods or unsteady airfoil theory may also be used, with their limitations acknowledged. The aerodynamic force on the structure using an inviscid panel method is usually obtained from the pressure distribution, , that acts over the surface, i.e.,

(161)  

where is the differential vector area and is the outward unit normal. In a viscous CFD calculation, surface shear stresses may also contribute to the total aerodynamic load, although pressure loads usually dominate most aircraft flutter calculations. Other methods, however, may be used to approximate aerodynamic behavior, including unsteady airfoil theory (using the Wagner function with the Duhamel integral), panel methods, and others. It is not unusual for flutter calculations to be performed using various aerodynamic methods (with different assumptions) and for the results to be compared.

For a quasi-steady aerodynamic model, or for an unsteady aerodynamic model that has been reduced to frequency-independent aerodynamic stiffness and damping matrices, the coupled equations may be written as

(162)  

where and represent aerodynamic stiffness and damping contributions under that approximation. The corresponding quadratic eigenvalue problem is

(163)  

In a general frequency-domain unsteady aerodynamic formulation, the aerodynamic influence matrices depend on airspeed and reduced frequency. The flutter problem is then nonlinear in and must be solved iteratively, or the aerodynamics must first be represented by an augmented state-space or rational-function model.

When the eigenvalue is written as a complex circular frequency, the real part of , i.e., , represents the oscillation frequency of the system. The imaginary part, , determines the temporal growth or decay of the motion. For the assumed convention , then

(164)  

so that if , the motion is damped and stable, whereas if , the motion grows with time, indicating an instability. Equivalently, the modal growth rate may be defined as

(165)  

so that is stable, is neutrally stable, and is unstable. Flutter occurs when , which corresponds to the transition from a stable damped response to an unstable oscillatory response.

The balance between stiffness (including aerodynamic stiffness) and inertial forces is contained within the term , which governs the oscillatory behavior of the system. The damping effects, including aerodynamic damping, are contained within the term . If the dynamic terms are neglected, the governing equation reduces to

(166)  

which corresponds to a static aeroelastic problem. A non-trivial solution exists when

(167)  

which defines the condition for static divergence. At this point, the structure loses stiffness and can undergo unbounded deformation under aerodynamic loading.

A typical flutter boundary is shown in Figure 17, which plots the modal growth or decay rate as a function of airspeed. Flutter speed is the point at which damping equals zero, indicating the onset of unstable oscillation. Notice that negative values indicate a stable system with decaying motion, the zero crossing marks the critical flutter speed, and positive values indicate an unstable system with growing oscillations.

The field of aeroelasticity also considers limit-cycle oscillations (LCOs), in which periodic, self-sustaining oscillations arise from nonlinearities in the aerodynamic or structural response. Engineers use sophisticated computational methods and experimental testing to predict and mitigate these aeroelastic phenomena, ensuring that aircraft and spacecraft designs meet stringent safety, performance, and certification standards in diverse flight conditions.

Torsional Divergence and Flutter of a Cantilevered Wing

Understanding aeroelastic instability is often obscured by complex mathematics or code that is treated as a “black box.” However, a classic example for those being introduced to the field of aeroelasticity is a cantilever wing encastré at its root with a single degree of freedom in pitch (torsional motion), as shown in Figure 18. This simple model can illustrate static torsional divergence, in which the effective torsional stiffness becomes zero. Classical flutter, however, is a dynamic instability involving oscillatory motion and aerodynamic damping, and it generally requires at least two coupled degrees of freedom or an equivalent unsteady aerodynamic coupling. Therefore, the maximum airspeed of many aircraft may be limited by aeroelastic effects, which involve compromises among wing strength, torsional rigidity, damping, weight, and cost.

In this simplified sectional idealization, the generalized coordinate is equivalent to , the pitch angle of the section about its elastic axis, so the equation of motion per unit span is

(168)  

where is the mass moment of inertia per unit span about the elastic axis, is the structural damping coefficient per unit span, is the torsional stiffness per unit span, and is the aerodynamic moment per unit span about the elastic axis. A simple linearized aerodynamic model for the moment about the elastic axis may be written as

(169)  

where is the aerodynamic moment slope about the elastic axis with respect to torsional displacement, is the aerodynamic moment slope with respect to pitch rate, and is the mean chord of the section. The first aerodynamic term represents an aerodynamic stiffness effect. The second aerodynamic term represents an aerodynamic damping or phase-lag effect. If the change in aerodynamic angle of attack, , is proportional to the change in , then may also be written as under this convention.

Therefore, the equation of motion becomes

(170)  

This equation shows that the aerodynamic terms modify both the effective damping and the effective torsional stiffness of the system.

Analytical Approach

Assume a trial solution of the form

(171)  

Differentiating once and twice gives

(172)  

Substituting these results into the equation of motion gives the characteristic equation

(173)  

The real and imaginary parts of this equation correspond to two distinct aeroelastic instability mechanisms.

Static Torsional Divergence

Static torsional divergence occurs when the effective torsional stiffness becomes zero. In this case, the instability is non-oscillatory, so

(174)  

The characteristic equation then reduces to

(175)  

or

(176)  

where is the divergence speed. Rearranging gives

(177)  

This result applies when the aerodynamic moment slope is destabilizing, i.e., when it reduces the effective torsional stiffness. If the aerodynamic moment is stabilizing, then this simple divergence condition is not obtained.

Single-Degree-of-Freedom Dynamic Instability

For this one-degree-of-freedom torsional model, an oscillatory dynamic instability can occur only if the aerodynamic damping becomes destabilizing. Therefore, this example is best interpreted as a simple negative-damping instability rather than a complete classical flutter model.

For a nonzero oscillatory response,

(178)  

the imaginary part of the characteristic equation gives

(179)  

Because , the term in parentheses must vanish. Therefore,

(180)  

where is the critical airspeed at which the net torsional damping becomes zero. Rearranging gives

(181)  

with the understanding that must have the sign that makes the aerodynamic damping destabilizing.

The corresponding flutter frequency is obtained from the real part of the characteristic equation, i.e.,

(182)  

so

(183)  

This exemplar illustrates the difference between two aeroelastic instabilities. Divergence is a static instability associated with loss of effective stiffness, whereas flutter is a dynamic instability associated with loss of effective damping at a nonzero oscillation frequency. Notice from Eq. 177 that increasing torsional stiffness raises the divergence speed. Notice also from Eq. 181 that, in this simplified model, the flutter speed depends on the balance between structural damping and aerodynamic damping. In more complete wing models, bending, torsion, structural mode shapes, and unsteady aerodynamic phase lags are coupled, so the flutter speed is normally obtained from a full eigenvalue analysis.

Time-Marching Solution

Another approach to determining flutter onset is a numerical integration method that marches the equations of motion forward in time. First, a time step must be chosen that is small enough to capture the response accurately. Second, initial values for and must be specified; these values represent disturbances used to excite the aeroelastic system. Physically, such a disturbance could be a gust or a sudden flight control input. During flight testing, mechanical shakers may also be used to excite the structure at selected frequencies to examine potential aeroelastic problems, such as flutter.

For a single torsional degree of freedom, the governing equation per unit span can be written as

(184)  

where and include both structural and aerodynamic contributions per unit span. Rearranging gives the angular acceleration at the current time step, i.e.,

(185)  

Then, a numerical integration method, such as the Euler method, which is only first-order accurate, or a Runge-Kutta method, which is much more accurate but more computationally expensive, can be used to update the solution at each step, e.g.,

(186)  

The aerodynamic moment at each time step may be computed from the current values of , , and . Unsteady aerodynamic effects can also be incorporated using the indicial method with superposition, as described previously for the numerical solution of the Duhamel integral. The process continues by analyzing the aerodynamic and structural responses until the desired end time is reached or a specified condition, such as flutter onset, is met.

When the aerodynamic lag states are included, the time-marching problem can also be written in augmented state-space form. For example, for a single torsional structural degree of freedom with indicial aerodynamic states, define the structural state vector as

(187)  

and let the aerodynamic lag states be collected into

(188)  

where and represent the hereditary states of the shed wake associated with the exponential approximation to the indicial response. The complete aeroelastic state vector may then be written as

(189)  

The coupled equations can then be written as

(190)  

where is the aeroelastic state matrix. This matrix contains the structural inertia, stiffness, and damping terms, as well as the aerodynamic lag-state equations and their coupling to the structural motion. For a forced response problem, gusts or control inputs may be added through an input matrix, giving

(191)  

where represents prescribed inputs such as gust velocity or control-surface motion.

This state-space form is useful because the same equations can be used for both time integration and stability analysis. In a time-marching calculation, the augmented state vector is advanced using a numerical integration method. For stability analysis, the eigenvalues of are computed at each airspeed. Flutter occurs when a complex-conjugate pair of eigenvalues crosses into the right-half plane, i.e., when its real part changes from negative to positive.

Figure 19 shows typical results obtained from an aeroelastic stability calculation. If the response decays sufficiently long after the disturbance, the system is stable at that airspeed. The process can then be repeated for other airspeeds or for variations in other parameters, such as the center-of-gravity location, which affects the torsional inertia about the elastic axis. Depending on the wing design, its aeroelastic response may be well damped, lightly damped, exhibit limit-cycle oscillations, or become dynamically unstable. The process would be repeated at different airspeeds and air-density values to estimate the stability boundaries for the structure or the entire aircraft.

Flutter of a Two-Dimensional Airfoil Section

The flutter of a two-dimensional (2D) airfoil section is a more tangible case to study. Consider a two-degree-of-freedom model for a typical 2D airfoil section, as shown in Figure 20. The section can undergo plunge motion, , representing the vertical displacement of the elastic axis, and pitch motion, , representing rotation about the elastic axis.

Let be the semi-chord. In this section, the parameter locates the elastic axis relative to the mid-chord in semi-chord units, i.e.,

(192)  

where is the mid-chord location. With this convention, the leading edge is at , the mid-chord is at , and the quarter-chord point is at . Therefore, corresponds approximately to the aerodynamic center of a thin airfoil in incompressible subsonic flow. The value of affects the aerodynamic moment arm and the coupling between plunge and pitch motions in the flutter analysis.

The governing equations for the plunge and pitch degrees of freedom may be written as

(193)  

where is the mass per unit span, is the static moment about the elastic axis, and is the mass moment of inertia about the elastic axis, where is the radius of gyration. The quantities and denote the plunge stiffness and torsional stiffness, respectively. The aerodynamic lift force and pitching moment per unit span are denoted by and .

Using a quasi-steady thin-airfoil approximation, the aerodynamic lift and pitching moment may be written as

(194)  

where, in a simplified quasi-steady form, the lift coefficient may be written as

(195)  

The pitching-moment coefficient depends on the selected reference point and sign convention. Therefore, it is better at this introductory level to write it in the linear form

(196)  

where the coefficients depend on the elastic-axis location and the adopted moment convention. These expressions include both displacement-dependent and velocity-dependent aerodynamic terms. The displacement-dependent terms modify the effective stiffness of the airfoil section, while the velocity-dependent terms modify the effective damping.

It is useful to write the equations in matrix form. Define

(197)  

The structural mass and stiffness matrices may be written as

(198)  

where the sign of depends on the adopted positive direction for , the positive sense of , and the location of the center of mass relative to the elastic axis. The aerodynamic forces can then be written as

(199)  

where , and

(200)  

with

(201)  

Therefore, the aeroelastic equations of motion become

(202)  

This form explicitly shows how the aerodynamic terms alter both the damping and stiffness of the system.

For flutter, assume a harmonic motion of the form

(203)  

so that

(204)  

Substituting into the equations of motion gives

(205)  

A non-trivial solution requires that

(206)  

This determinant is the flutter equation for the simplified two-degree-of-freedom airfoil section. It must generally be solved numerically for each airspeed value. As is increased, the roots of this equation change. Flutter occurs when one of the roots reaches zero damping, so that the corresponding oscillatory motion no longer decays with time. The associated airspeed is the flutter speed, , and the corresponding oscillation frequency is the flutter frequency, .

This formulation is still quasi-steady, so it is best interpreted as an introductory model. More complete flutter calculations replace the quasi-steady aerodynamic matrices with unsteady aerodynamic influence matrices, often involving Theodorsen’s function, indicial-response approximations, panel methods, or CFD-based aerodynamic models. Nevertheless, this two-degree-of-freedom airfoil section illustrates the essential mechanism of flutter: the coupling of plunge, pitch, inertia, stiffness, and aerodynamic phase effects.

Flutter Speeds Using Theodorsen’s Function

The quasi-steady two-degree-of-freedom airfoil-section model developed above is useful for illustrating the basic flutter mechanism, but it does not fully represent the wake history and frequency-dependent phase lag of the aerodynamic loads. A more complete two-dimensional incompressible analysis uses Theodorsen’s function to represent the circulatory part of the unsteady aerodynamic response.

The key difference from the quasi-steady model is that the circulatory aerodynamic terms are no longer real constants. They depend on the complex Theodorsen function , where

(207)  

and . Therefore, the aerodynamic stiffness and damping contributions depend on the oscillation frequency and the airspeed. The noncirculatory apparent-mass terms must also be retained separately because they are not multiplied by .

The resulting flutter problem is a coupled eigenvalue problem. This may be written as

(208)  

where is the structural mass matrix, while and contain both structural and unsteady aerodynamic contributions. The aerodynamic contributions depend on , and hence on the reduced frequency.

A non-trivial solution requires

(209)  

with

(210)  

No general closed-form flutter speed follows from this formulation because , , and are coupled through

(211)  

In practice, the airspeed is varied, the corresponding complex eigenvalues are computed, and the flutter speed is identified as the speed at which one oscillatory mode first reaches zero damping.

For low-speed, two-dimensional incompressible problems, Theodorsen’s function provides an instructive analytical aerodynamic model. For finite wings, compressible flow, transonic flow, separated flow, or complex aircraft geometry, higher-fidelity aerodynamic models such as doublet-lattice methods, experimentally identified aerodynamic influence coefficients, Euler methods, or CFD-based methods are normally required.

Flutter Using FEM/CFD Coupling

Aeroelastic analyses may be conducted using time-marching methods, frequency-domain eigenvalue methods, or iterative coupling between aerodynamic and structural models. An iterative coupling process will typically proceed along these lines:

1. Define the geometry and grids for FEM and CFD. They will differ because the grids required to address aerodynamics differ from those for the structural shape.
2. Apply the material properties and boundary conditions in the FEM. Using isotropic materials, such as aluminum, will yield different properties from those of composites.
3. Set up the CFD with the specified flow conditions, or use an alternative aerodynamic model, e.g., a panel method.
4. Compute the aerodynamic forces, , from the selected aerodynamic model.
5. Update the FEM model with the aerodynamic forces .
6. Iterate between CFD and FEM to capture aeroelastic interactions. In a linearized aeroelastic formulation, this coupling may be represented as

   (212)  

   where and include both structural and aerodynamic contributions. In a time-marching CFD/FEM calculation, however, the aerodynamic loads are usually transferred directly between the flow solver and the structural solver rather than being represented by fixed aerodynamic damping and stiffness matrices.

7. Repeat the FEM analysis to obtain the structural response of the deformed structure.
8. Conduct CFD simulations of the deformed structure to obtain new aerodynamic loads.
9. Exchange data between FEM and CFD to formalize the coupled analysis.
10. Analyze displacements, stresses, and aerodynamic forces.
11. Evaluate the results against the stability criteria for flutter or divergence.

By integrating aerodynamic and structural models, engineers can predict and analyze aeroelastic phenomena, thereby supporting the safety and performance of aerospace structures. For example, this process has been applied to large aircraft such as the Boeing 747, as shown in Figure 21. These results provide valuable information on mode shapes, natural frequencies, damping, and aeroelastic deformation patterns, helping engineers identify and prevent aeroelastic and flutter issues before the first flight. However, flight testing is essential to validate these calculations and ensure the aircraft is flutter-free across its operational flight envelope.

Flutter Occurrences in Practice

In practice, aircraft structures can exhibit coupled responses between different “modes,” e.g., between wing torsion and bending, between the engine mount and the wing, or between the wing and a control surface. Two fatal crashes of the Lockheed Electra were attributed to a phenomenon known as “pylon whirl flutter,” in which aerodynamic forces on the propellers, gyroscopic moments, nacelle motion, engine-mount flexibility, and wing structural dynamics coupled to produce an unstable oscillation. Extensive structural modifications were made to the Electra’s design, including changes to the wings and the engine mounts to increase their stiffness.

During its first flight tests, the Boeing Dash 80 (predecessor to the Boeing 707) experienced aileron flutter, after which the control system was modified. This incident underscored the importance of flutter testing for jet airliners and led to improvements in control-surface design and testing protocols. The Boeing 747 experienced unexpected flutter in its horizontal stabilizers during flight testing, and the Boeing 747-8 experienced wingtip flutter before it was resolved.^[12]

The interactions among the relevant structural modes must be accounted for in an aircraft to capture the structure’s overall coupled dynamic behavior. This requires a structural dynamics model coupled to an appropriate aerodynamic model, which may range from linear unsteady aerodynamic methods, such as panel or doublet-lattice methods, to higher-fidelity CFD when nonlinear, transonic, or separated-flow effects are important. The process involves solving the coupled structural and aerodynamic equations and determining the eigenvalues that represent the system’s natural frequencies and damping ratios, as previously outlined. Figure 22 illustrates the information that can be extracted from this type of aeroelastic analysis: frequency and damping as functions of airspeed or dynamic pressure.

The frequency and damping of two structural modes as a function of airspeed can indicate that the likelihood of flutter increases as their modal frequencies approach each other. The undamped zero airspeed frequencies are distinct at low speeds, and the system remains stable. As the airspeed increases, the aerodynamic forces acting on the structure also increase in magnitude. These forces alter the structure’s frequencies by introducing effective aerodynamic stiffness and damping, which, in some cases, may cause the modal frequencies to converge, a phenomenon known as frequency coalescence. Frequency coalescence is often associated with classical flutter, although it is not a universal flutter criterion.

Flutter arises when the net structural-aerodynamic damping of an aeroelastic mode becomes negative. In this condition, the unsteady aerodynamic forces transfer energy to the coupled structural motion at a rate that exceeds the energy dissipated by structural and aerodynamic damping. The resulting flutter frequency is an eigenfrequency of the coupled aeroelastic system, not necessarily one of the natural frequencies at zero airspeed. If the net energy input over an oscillation cycle is positive, the oscillations grow exponentially, leading to flutter.

Flight Control Reversal

Control reversal occurs when the deflection of a control surface, such as an aileron, rudder, or elevator, produces an aircraft response opposite to the pilot’s intended input as a result of adverse aerodynamic or structural interactions. This behavior typically occurs at higher airspeeds, when aerodynamic forces on the control surface may induce excessive structural deformation, particularly in the wing or tail section. Suppose a pilot applies an aileron input to roll the aircraft. On the wing where the aileron deflection is downward and intended to increase lift, the aerodynamic loads from the deflected aileron may twist the wing in the opposite sense nose-down, decreasing the local effective angle of attack instead. If this aeroelastic twist is large enough, it can offset or exceed the intended aileron effect, producing a rolling moment opposite to the pilot’s command, as shown in Figure 23.

This undesirable effect is particularly pronounced in aircraft with long, slender wings, such as sailplanes, or those with wings that have insufficient torsional stiffness. In an attempt to minimize structural weight, wing skins are sometimes too thin to provide the necessary torsional stiffness. Historically, roll or aileron reversal was a significant issue in high-performance fighter airplanes during WWII and in early jet airplanes. Some aircraft, such as the Supermarine Spitfire, exhibited aileron reversal at high speeds, thereby limiting maneuverability in combat. Another contributing factor is the change in aerodynamic loading and pressure distribution as speed and Mach number increase, especially near transonic conditions. If the wing is not adequately designed to resist the resulting torsional moments, control effectiveness can degrade as dynamic pressure increases, eventually leading to reversal. This issue was common in early high-speed aircraft before engineers developed techniques to counteract it.

To mitigate control reversal, engineers may employ several strategies. Increasing the torsional stiffness of the aircraft’s structure, particularly in the wing and empennage, is usually the most direct way to resist excessive structural twisting. Other approaches include changing the size or spanwise placement of the control surface, using aerodynamic balance or geared tabs to reduce hinge moments, limiting control deflection at high dynamic pressure, or using spoilers for roll control. Mass balancing is important for suppressing control-surface flutter and buzzing, but it does not by itself prevent static control reversal. On modern high-speed aircraft, powered control systems, such as hydraulic or fly-by-wire actuators, can help manage aeroelastic effects by scheduling control-surface deflections, limiting commands at high dynamic pressure, or using alternative control surfaces. However, powered actuation does not by itself eliminate structural twisting or prevent static control reversal; adequate torsional stiffness and proper aeroelastic design remain essential.

Dynamic Stall & Stall Flutter

Stall flutter is an aeroelastic phenomenon that occurs when flow separation induces unsteady aerodynamic forces that couple with structural dynamics, leading to oscillatory instability. It is relevant in various aerospace and engineering applications, where predicting and mitigating its effects are essential for ensuring structural integrity and aerodynamic performance. Understanding stall flutter requires nonlinear aerodynamic modeling and experimental validation to develop practical mitigation steps and design improvements.

Physics of Dynamic Stall

Under unsteady conditions, airfoils and wings behave differently from those in steady flow. In particular, their stall characteristics differ significantly, a phenomenon known as dynamic stall. This phenomenon has attracted considerable research interest due to its unsteady aerodynamic effects, including the shedding of a leading-edge vortex, as illustrated by the simulation in Figure 24. The onset of a dynamic stall can lead to another type of aeroelastic behavior called stall flutter, which can occur on helicopter blades.

Vortex shedding increases lift and drag, resulting in powerful nose-down (negative) pitching moments induced by the aft-moving center of pressure, as illustrated in Figure 25, based on measurements and flow visualization. Dynamic stall is characterized by higher values of maximum lift, drag, and pitching moment, as well as hysteresis effects that can lead to flutter. It is particularly relevant in the design and analysis of helicopter rotors, wind turbine blades, and specific aircraft configurations.

A defining feature of dynamic stall is the presence of strong hysteresis and phase lag between the airfoil motion and the resulting aerodynamic loads. Unlike steady aerodynamics, in which lift and moment depend only on the instantaneous angle of attack, unsteady forces during dynamic stall depend on the prior history of flow separation and reattachment. Consequently, the lift and pitching moment do not remain in phase with the airfoil motion and generally lag behind it because of delayed separation, vortex formation, and reattachment.

This phase lag is especially pronounced in the pitching moment because of the formation and downstream convection of a leading-edge vortex and the associated aft movement of the center of pressure. During portions of an oscillation cycle, the aerodynamic moment can act in the same sense as the airfoil’s angular velocity rather than opposing it. In energy terms, unsteady aerodynamics can perform positive net work on the airfoil over part of a cycle rather than dissipating energy.

This history-dependent behavior provides the essential aerodynamic mechanism by which dynamic stall can reduce or even reverse aerodynamic damping. When such unsteady separated-flow forces act on a flexible lifting surface, they can couple with the structural dynamics, triggering an aeroelastic instability known as stall flutter.

Engineers and researchers study dynamic stall to understand its effects better and develop strategies to mitigate its adverse impacts on the performance of aerospace systems. In particular, helicopter rotor blades often experience dynamic stall during higher-speed forward flight or during maneuvers, and its onset effectively limits the helicopter’s operational flight envelope. Wind turbine blades may also experience dynamic stall during sudden changes in wind conditions, leading to fluctuations in power output and high blade loads that can result in alarmingly high vibration levels. Aircraft wings can also experience dynamic stall during aggressive maneuvers, a phenomenon that must be accounted for when establishing the maneuvering flight envelope for military combat airplanes.

Leishman-Beddoes Dynamic-Stall Model

A widely used engineering model for dynamic stall is the Leishman-Beddoes model.^[13] It is a semi-empirical model that represents the main unsteady aerodynamic processes observed during dynamic stall, including attached-flow unsteady response, trailing-edge separation, leading-edge vortex formation, vortex convection, and reattachment. The model was originally developed for practical rotorcraft and airfoil applications, where fully resolved unsteady separated-flow calculations were too expensive for routine aeroelastic analysis.

The Leishman-Beddoes model may be implemented either as time-marching recurrence relations or in state-space form. The state-space formulation is especially useful for aeroelastic stability analysis because the aerodynamic states can be appended directly to the structural equations of motion.^[14]

The important point is that the Leishman-Beddoes model is not simply a corrected lift curve. It is a dynamic model with internal aerodynamic states. These states represent the lagged response of the attached flow, the delay in the onset of flow separation, the growth and convection of the dynamic-stall vortex, and the flow recovery after vortex shedding. In this sense, it extends the same general idea used in indicial-response modeling, i.e., the aerodynamic loads depend not only on the instantaneous angle of attack but also on the airfoil’s previous motion history.

In time-marching calculations, the model is often implemented using recurrence relations. At each time step, the effective angle of attack and the attached-flow response are updated, separation-point delays are computed, vortex-induced lift and pitching moment increments are added when appropriate, and the normal-force, chord-force, and pitching-moment coefficients are assembled from these components, which can be written as

(213)  

and

(214)  

where denotes the aerodynamic hereditary effects at the current time step and represents the resulting aerodynamic coefficients, such as normal force, chord force, and pitching moment. This recurrence form is convenient for time-domain simulations of rotor blades, wind-turbine blades, or airfoils undergoing large-amplitude motion.

The same model can also be cast in state-space form. In that case, the aerodynamic lag variables are collected into a state vector and advanced through first-order evolution equations, e.g.,

(215)  

with output equations of the form

(216)  

This form is especially useful in aeroelastic analysis because the aerodynamic states can be appended to the structural states, just as with the linear indicial-response states described earlier. The resulting coupled system can then be integrated forward in time, linearized locally for stability analysis, or used to predict the onset of stall flutter and limit-cycle oscillations.

The Leishman-Beddoes model reproduces the hysteresis in the lift and pitching moment that occurs during dynamic stall. Figure 26 compares measured data with model predictions and static airfoil data for two oscillatory pitching cases. The loops in the pitching-moment coefficient are particularly important for aeroelasticity because their orientation indicates whether the unsteady aerodynamic moment extracts energy from the motion or feeds energy into it. Portions of the loop may correspond to positive aerodynamic damping, while others may correspond to negative aerodynamic damping. If the net aerodynamic work over a cycle becomes destabilizing and exceeds the available structural damping, stall flutter or a limit-cycle oscillation may occur.

Therefore, the Leishman-Beddoes model bridges the linear unsteady aerodynamic models discussed earlier and the nonlinear separated-flow behavior associated with dynamic stall. Theodorsen’s theory, Wagner’s function, and related indicial models describe attached-flow unsteady aerodynamics, whereas the Leishman-Beddoes model adds semi-empirical state variables for flow separation, vortex shedding, and reattachment. This capability is why it has been widely used in rotorcraft, wind-turbine, and other aeroelastic calculations involving dynamic stall.

Stall Flutter

Stall flutter is a nonlinear aeroelastic phenomenon that occurs when an airfoil or structure experiences flow separation, such as dynamic stall, leading to oscillatory instability. This separation produces unsteady aerodynamic forces that can destabilize the structural response. Unlike classical flutter, which occurs in a linear aerodynamic regime before stall, stall flutter happens post-stall, when the flow has already separated from the surface. Because it depends on unsteady, nonlinear aerodynamics, stall flutter is more complex to analyze and predict.

The primary cause of stall flutter is flow separation, which occurs when the angle of attack exceeds the critical stall limit. Dynamic stall further intensifies this phenomenon, as lags in periodic flow separation and reattachment reduce or even reverse aerodynamic damping, as shown in Figure 27 for the pitching moment. If a structure has limited structural stiffness and/or damping, it can lead to stall flutter.

From a dynamical standpoint, stall flutter can be interpreted as a loss of net structural damping. Consider, for simplicity, a single generalized structural degree of freedom governed by

(217)  

where is the generalized mass, is the structural damping coefficient, is the structural stiffness, and represents the unsteady aerodynamic generalized force. Although is inherently nonlinear in stalled flow, it may be locally linearized about a mean post-stall condition or a small oscillatory motion as

(218)  

where and are effective aerodynamic damping and stiffness coefficients, respectively. The resulting equation of motion becomes

(219)  

In attached flow, is typically stabilizing, so that the effective damping remains positive. Near and beyond stall, however, phase lags associated with flow separation and reattachment can cause the aerodynamic contribution to become destabilizing, reducing or even reversing the net damping. Stall flutter corresponds to the condition

(220)  

at which point small perturbations grow rapidly until nonlinear aerodynamic effects limit the motion, establishing a finite-amplitude oscillation (a limit cycle).

Stall flutter can be an issue for wings and control surfaces, such as flaps and ailerons, at high angles of attack. Helicopter rotor blades have been known to encounter stall flutter in the retreating blade region, where local angles of attack can become high during forward flight. Stall flutter can also affect compressor and turbine blades in jet engines, where unsteady flow interactions are crucial in their performance. Civil engineering structures, such as bridges and tall buildings, may exhibit vortex-induced oscillations that resemble stall flutter.

Mathematically, stall flutter is challenging to model because of its nonlinearity. Traditional linear aeroelastic models, such as Theodorsen’s theory or the linear indicial response method, are insufficient because they assume attached flow. Instead, nonlinear aerodynamic models are required to capture the complex flow physics. Semi-empirical models, such as the Leishman-Beddoes dynamic-stall model, provide approximations of the unsteady aerodynamic forces encountered in stall flutter. Coupled fluid-structure interaction simulations are employed to investigate the impact of separated flow on structural vibrations. Mitigating stall flutter involves several approaches, including modifying the airfoil’s geometry to delay flow separation and reduce vortex shedding. Increasing structural stiffness and damping can also help dissipate energy and suppress the structural response.

Buffeting & Buzzing

While buffeting, buzzing, and stall flutter all involve unsteady aerodynamic forces acting on flexible structures, they differ fundamentally in their underlying mechanisms. Buffeting is a forced structural response to aerodynamic disturbances, whereas buzzing and stall flutter are self-excited aeroelastic instabilities. Buzzing is a self-excited, high-frequency oscillation of control surfaces, often because of aeroelastic feedback and insufficient damping.

Buffeting

Buffeting is a random, low-frequency, forced vibration caused by unsteady aerodynamic forces acting on a structure. It typically occurs when a lifting surface, such as an aircraft wing or tailplane, is subjected to turbulence or wake disturbances originating from another part of the airplane. One of the most common causes of buffeting is flow separation and vortex shedding behind a stalled airfoil. When a surface encounters separated flow, turbulent eddies form and impinge on downstream airframe structures, generating fluctuating forces, as shown in the schematic in Figure 28.

For example, an aircraft’s horizontal stabilizer may experience buffeting when positioned in the wake of the wing or fuselage, as shown in this video of MD-11 airliner stall tests. This behavior may occur at higher angles of attack, such as during takeoff or landing, especially when flaps are deployed. Similarly, shock waves interacting with boundary layers in transonic or supersonic flight can induce separated flow regions that produce turbulence and cause severe buffeting on downstream airframe structures.

Buffeting is a crucial consideration in airplane design because it can increase aerodynamic drag, lead to control difficulties, and cause structural fatigue. Engineers analyze buffeting using wind tunnel testing, CFD simulations, and flight testing, and develop mitigation strategies based on these results. Techniques to reduce buffeting include optimizing the shape and size of aerodynamic surfaces, adding vortex generators to reduce or control flow separation, and designing control systems to minimize buffeting response.

Buzzing

Buzzing is a high-frequency, self-excited aeroelastic behavior that occurs primarily in control surfaces such as ailerons, rudders, and elevators. It is caused by the interaction between the aerodynamic forces acting on a control surface and the structural dynamics of the hinge mechanism. Unlike buffeting, which is an externally forced vibration, buzzing is a self-excited oscillation resulting from unsteady aerodynamic loading on a lightly damped control surface. A common cause of buzzing is aeroelastic feedback from hinged surfaces, as shown in Figure 29.

If a control surface is not sufficiently stiff or damped, minor disturbances in the airflow can cause rapid surface oscillations that are then amplified by aerodynamic forces. In transonic or supersonic flight, buzzing can also result from shock-induced separation and reattachment over a control surface. This creates fluctuating pressure loads that drive oscillations in the strength and location of high-frequency shock waves. Engineers can mitigate buzzing by adding stiffness, mass-balancing weights, or dampers to suppress rapid, undesirable control-surface movements.

For further clarity, the table below compares buffeting, buzzing, dynamic stall, stall flutter, and classical flutter with respect to their excitation mechanisms, flow regimes, and aeroelastic coupling.

Rotating-System Aeroelasticity

Rotating systems add a distinct class of aeroelastic effects to aircraft design. Propellers, helicopter rotors, tiltrotors, ducted fans, and the many small rotors used on UAVs and eVTOL aircraft all feature rotating blades. These blades are connected to hubs, shafts, bearings, motors, gearboxes, pylons, wings, booms, or fuselages, so their aeroelastic behavior involves the complete rotor-support-airframe system. A rotating blade may be viewed, in a first approximation, as a slender beam subjected to radial tensile loading caused by its own rotation. Consider a small blade element of length at radius . If the blade has mass per unit length , then the radial force required to keep this element moving in a circular path is

(221)  

The blade section at radius must carry the radial force from all of the blade material outboard of that station. Therefore, the centrifugal force at radius is

(222)  

where is the blade radius and is the rotor angular speed. For a uniform blade, this expression becomes

(223)  

so the centrifugal force increases in proportion to .This force changes the blade bending dynamics. For a simple blade-bending model, the transverse deflection may be represented by a modal coordinate and a mode shape , i.e.,

(224)  

The corresponding modal equation may be written in the form

(225)  

where , , and are the modal mass, damping, and elastic bending stiffness, respectively. The term is the additional bending stiffness caused by the centrifugal force. A simple expression for this centrifugal-stiffness contribution is

(226)  

Therefore, the blade bending frequency is approximately

(227)  

or, in simplified form,

(228)  

where is the nonrotating bending frequency and is a nondimensional coefficient that depends on the blade mass distribution, stiffness distribution, and mode shape.

This equation is the key point for the fan plot, i.e., rotating-blade bending modes curve upward as rotor speed increases, while wing bending, wing torsion, and nacelle pitch/yaw modes appear as nearly horizontal support-structure frequencies. The periodic excitation of a rotor system occurs at integer multiples of the rotor speed. Therefore, the forcing frequencies are

(229)  

where . These are called , , , and higher harmonic excitations. The corresponding frequencies in cycles per second are

(230)  

A resonance condition is approached when a blade, hub, nacelle, pylon, wing, boom, or fuselage natural frequency is close to one of these harmonic excitation frequencies, i.e.,

(231)  

This condition is the mathematical basis of a rotor frequency diagram, commonly called a fan plot. As shown schematically in Figure 30, many rotating-blade bending frequencies increase with rotor speed because of centrifugal stiffening. Torsional modes and coupled rotating-system modes may follow more complicated trends because of centrifugal, gyroscopic, elastic, and aerodynamic coupling, so their frequency variation should be obtained from an appropriate rotating-system modal analysis. The wing bending, wing torsion, and nacelle pitch/yaw frequencies are shown as nearly horizontal support-structure frequencies. The diagonal straight lines are the rotor harmonics , , , and . The figure shows the main design concern, i.e., if a modal-frequency curve approaches or crosses one of the harmonic lines within the operating rotor-speed range, then resonant vibration may occur.

Fan plots are useful because they show whether important blade or support-structure modes are safely separated from the main rotor harmonics over the operating range. Frequency separation is only the first check. Damping, aerodynamic phase lag, mode coupling, control-system dynamics, and structural nonlinearities also affect the response.

Whirl flutter is one important rotating-system instability. It is especially important for wing-mounted propellers, turboprops, and tiltrotors, where the propeller or proprotor is mounted on a flexible nacelle, pylon, or wing. Tiltrotor aircraft such as the Bell-Boeing V-22 Osprey are particularly sensitive to this phenomenon because the large proprotors and tilting nacelles can couple with pylon motion, wing bending, and wing torsion. The instability occurs when the aerodynamic and gyroscopic moments from the propeller or rotor couple with the pitch and yaw motion of the supporting structure. If this coupling transfers energy into the nacelle, pylon, and wing faster than structural damping can dissipate it, the oscillation grows.

Consider a propeller mounted on a nacelle, pylon, or other support that can pitch and yaw. Let represent the pitch motion of the nacelle and represent its yaw motion. A simple linearized structural model may be written as

(232)  

where and are the pitch and yaw mass moments of inertia, and are damping coefficients, and are structural stiffnesses, and and are the aerodynamic and gyroscopic moments produced by the propeller. The gyroscopic coupling of a spinning propeller introduces cross-coupled pitch and yaw moments of the form

(233)  

where is the polar mass moment of inertia of the rotating propeller. This equation shows that a yaw rate can produce a pitch moment and a pitch rate can produce a yaw moment. The aerodynamic moments from the propeller also depend on nacelle attitude and angular rate, so they may be represented in linearized form as

(234)  

Therefore, the complete system contains structural stiffness, structural damping, gyroscopic coupling, and aerodynamic coupling. In matrix form, it may be written as

(235)  

where

(236)  

The onset of whirl flutter occurs when one of the system modes loses damping. In eigenvalue form, the solution may be written as , where . The motion is stable when all modes have . Flutter occurs when at least one mode reaches neutral stability, i.e., , and becomes unstable when . This form of the equations shows why propeller whirl flutter is a system-level aeroelastic instability. The instability depends on the propeller, engine or motor, nacelle, pylon, wing, stiffness, damping, mass distribution, gyroscopic coupling, and aerodynamic forces. The complete installation determines the stability boundary.

Recall from earlier in this chapter that the Lockheed Electra accident history is a classic example of this type of instability. The propeller, engine, nacelle, and wing formed a coupled aeroelastic system. Under unfavorable conditions, the propeller’s gyroscopic and aerodynamic effects, coupled with nacelle and wing flexibility, produced an unstable whirl motion. The important lesson is that rotating-system aeroelasticity must be assessed at the aircraft design level.

The same concepts are increasingly relevant to modern UAVs and eVTOL aircraft. Many of these aircraft use multiple propellers or rotors mounted on lightweight arms, booms, pylons, or wings. The rotor speeds may vary widely, especially when electric motors are used for thrust control. Therefore, the structure may pass through several possible excitation conditions during startup, shutdown, maneuvering, gust response, or transition flight. If rotor harmonics align with a blade, boom, wing, motor-mount, or fuselage mode, then vibration levels may become large enough to affect payload performance, flight-control sensors, autopilot behavior, or structural integrity.

For multirotor UAVs, an additional source of excitation can occur when the rotors operate at slightly different speeds. This effect is less relevant to conventional helicopter rotors, which usually have one main rotor speed, but it can be important for aircraft with several independently controlled electric rotors. If rotor has angular speed and rotor has angular speed , then their harmonic excitations may combine to produce beat frequencies. A general beat frequency between the th harmonic of one rotor and the th harmonic of another rotor is

(237)  

The simplest case is the difference between the two frequencies, i.e.,

(238)  

These beat frequencies usually appear as low-frequency excitation lines on a fan plot. They can be important because low-frequency beats may approach wing-bending, boom-bending, fuselage, or motor-mount modes even when the main rotor harmonics are well separated from those modes.

For this reason, rotating-system aeroelasticity requires attention to blade stiffness, rotor speed, hub design, motor and gearbox support stiffness, pylon or boom flexibility, damping, and control-system interactions. For multirotor UAVs, the assessment should also include rotor-to-rotor speed differences and the beat frequencies they may produce. Fan plots, ground vibration tests, rotor and structural dynamic models, and flight testing all contribute to the assessment. The objective is to avoid dangerous combinations of rotor speed, structural frequency, aerodynamic forcing, gyroscopic coupling, weak damping, and inter-rotor beat excitation.

Summary & Closure

Aeroelasticity concerns the coupling between aerodynamic forces, elastic deformation, and structural dynamics. These interactions can produce static instabilities, such as divergence and control reversal, as well as dynamic instabilities, such as classical flutter, whirl flutter, buzzing, and stall flutter. The common thread is that aerodynamic loading and structural motion cannot always be treated separately; under some conditions, they couple in ways that reduce stiffness or damping, or feed energy into the structure’s motion.

A central theme in aeroelastic analysis is the role of unsteady aerodynamics. Theodorsen’s and Sears’s functions describe harmonic unsteady responses in the frequency domain, while Wagner’s and Küssner’s functions describe corresponding indicial responses in the time domain. These classical incompressible theories clarify wake hereditary effects, phase lag, amplitude attenuation, and the distinction between circulatory and noncirculatory loading. However, they are baseline models; compressibility, separated flow, dynamic stall, and nonlinear structural effects may require semi-empirical models, state-space formulations, or higher-fidelity aerodynamic methods.

Modern aeroelastic analysis uses a hierarchy of models, from simple quasi-steady and two-degree-of-freedom section models to finite-element structural models coupled with panel, doublet-lattice, and Euler methods, or Reynolds-averaged Navier-Stokes CFD. Wind-tunnel testing, ground vibration testing, and flight flutter testing remain essential for validation. The practical goal is to ensure that aircraft, rotor systems, and other aerospace vehicles remain free from destructive aeroelastic behavior throughout their operating envelopes while meeting structural, aerodynamic, performance, and certification requirements.