AE4ASM506 Notes 2012

AE4930 Aeroelasticity
60
50
40
30
20
10
0
0 50 100 150 200
S. J. Hulshoff
HSL 0.36, Aerodynamics Group
Faculty of Aerospace Engineering, TU Delft
S.J.Hulshoff@TUDelft.NL
Version 11.1
2
Contents
1 Introduction 9
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2 Course Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3 Additional Course Information . . . . . . . . . . . . . . . . . . . 13
1.3.1 Course format: . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.2 Prerequisite Concepts . . . . . . . . . . . . . . . . . . . . 14
1.3.3 Main References . . . . . . . . . . . . . . . . . . . . . . . 14
1.3.4 Instructor Contact . . . . . . . . . . . . . . . . . . . . . . 14
1.4 Characteristics of Aeroelastic Problems . . . . . . . . . . . . . . 15
2 Simplified Aeroelastic Problems 21

2.1 The Typical Section . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1.1 Simple Aerodynamic Models . . . . . . . . . . . . . . . . 24
2.1.2 Equations of motion for the typical section . . . . . . . . 25
2.1.3 Torsional Divergence . . . . . . . . . . . . . . . . . . . . . 25
2.1.4 Control Reversal . . . . . . . . . . . . . . . . . . . . . . . 27
2.1.5 Flutter of a Typical Section . . . . . . . . . . . . . . . . 28
2.1.6 Solution of the equations of motion . . . . . . . . . . . . . 30
2.1.7 Mode Shapes . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.1.8 Conditions which determine the types of motion . . . . . 34
2.1.9 Pines’ Conditions for flutter . . . . . . . . . . . . . . . . . 36
2.1.10 The effect of a low-frequency aerodynamic model . . . . . 37
2.2 The Semi-Rigid Wing . . . . . . . . . . . . . . . . . . . . . . . . 39
2.2.1 Equations of motion for the Semi-rigid wing . . . . . . . . 40
2.3 Propeller Whirl Flutter . . . . . . . . . . . . . . . . . . . . . . . 41
2.3.1 A Simplified model . . . . . . . . . . . . . . . . . . . . . . 42
2.3.2 Gyroscopic Coupling . . . . . . . . . . . . . . . . . . . . . 43
2.3.3 Aerodynamic Coupling . . . . . . . . . . . . . . . . . . . . 43
2.4 Practice Problems . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3 Unsteady Aerodynamics 47
3.1 Vorticity Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2 Compressibility effects . . . . . . . . . . . . . . . . . . . . . . . . 52
3.3 Viscous effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3
4 CONTENTS
3.4 Model Equations for Unsteady Flows . . . . . . . . . . . . . . . . 59

3.4.1 Navier-Stokes Equations . . . . . . . . . . . . . . . . . . . 59
3.4.2 Boundary-Layer Equations . . . . . . . . . . . . . . . . . 62
3.4.3 Euler Equations . . . . . . . . . . . . . . . . . . . . . . . 63
3.4.4 Full Potential Equations . . . . . . . . . . . . . . . . . . . 63
3.4.5 Small-Disturbance Potential Equation . . . . . . . . . . . 65
3.4.6 Linear Potential Equation . . . . . . . . . . . . . . . . . . 66
3.5 Analytic Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.5.1 Incompressible Flows . . . . . . . . . . . . . . . . . . . . . 67
3.5.2 Subsonic Compressible Flows . . . . . . . . . . . . . . . . 67
3.5.3 Transonic Flows . . . . . . . . . . . . . . . . . . . . . . . 67
3.5.4 Supersonic Flows . . . . . . . . . . . . . . . . . . . . . . . 67
4 Analytic Solutions for Incompressible Flows 71

4.1 Problem definition . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.2 Far-Field Perturbation . . . . . . . . . . . . . . . . . . . . . . . 73
4.3 Velocities on the z = 0 plane . . . . . . . . . . . . . . . . . . . . 73
4.4 Flow Tangency . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.5 Pressure Jump Across the Wake . . . . . . . . . . . . . . . . . . 75
4.6 Kutta Condition using Kelvin’s Theorem . . . . . . . . . . . . . 76
4.7 Harmonic Solutions . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.8 Harmonic Pitch and Plunge solutions . . . . . . . . . . . . . . . 79
4.9 Indicial Response (Wagner’s Function) . . . . . . . . . . . . . . . 82
4.10 Vertical Step Gust Response (Küssner’s Function) . . . . . . . . 83
5 Analytic Solutions for Subsonic Compressible Flows 87

5.1 Steady Flow: The Prandtl-Glauert Correction . . . . . . . . . . 87
5.2 Unsteady Harmonic Solutions for M < 1 . . . . . . . . . . . . . 88
5.3 Piston theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.4 Indicial Responses for M < 1 . . . . . . . . . . . . . . . . . . . . 92
6 Analytic Solutions for Supersonic Flows 95

6.1 Problem definition . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.1.1 General approach: . . . . . . . . . . . . . . . . . . . . . . 97
6.2 Time Integration Limits . . . . . . . . . . . . . . . . . . . . . . . 98
6.3 Spatial Integration Limits . . . . . . . . . . . . . . . . . . . . . . 99
6.4 Potential induced by the airfoil . . . . . . . . . . . . . . . . . . . 99
6.5 Normal velocity at the Airfoil . . . . . . . . . . . . . . . . . . . 100
6.6 Flow-Tangency Condition . . . . . . . . . . . . . . . . . . . . . . 101
6.7 Solutions for Harmonic Motion . . . . . . . . . . . . . . . . . . . 101
6.8 Finite wings, Multiple Surfaces . . . . . . . . . . . . . . . . . . . 102
CONTENTS 5
7 Arbitrary Unsteady Aerodynamic Responses 107

7.1 Indicial and Impulse Functions . . . . . . . . . . . . . . . . . . . 108
7.2 Fourier Transforms . . . . . . . . . . . . . . . . . . . . . . . . . 109
7.2.1 Computation of an arbitrary response . . . . . . . . . . . 109
7.2.2 Re-expression of the incompressible results . . . . . . . . 109
7.2.3 Conversion to an Indicial response . . . . . . . . . . . . 110
7.2.4 Conversion from an Indicial response . . . . . . . . . . . 111
7.3 State-Space Representations . . . . . . . . . . . . . . . . . . . . 111
8 Compact Models of Structural Systems 113

8.1 Generalised Coordinates . . . . . . . . . . . . . . . . . . . . . . 114
8.2 Lagrange’s Equations of Motion . . . . . . . . . . . . . . . . . . 115
8.3 Lagrange’s Equations for a Flexible Aircraft . . . . . . . . . . . 117
8.3.1 Options for efficient treatment . . . . . . . . . . . . . . . 118
8.4 Methods for Constructing Compact Models . . . . . . . . . . . . 119
8.5 Lumped Parameter Method . . . . . . . . . . . . . . . . . . . . . 120
8.5.1 Example Application . . . . . . . . . . . . . . . . . . . . . 120
8.6 Rayleigh-Ritz (Assumed Mode) Method . . . . . . . . . . . . . . 122
8.6.1 Example Application . . . . . . . . . . . . . . . . . . . . . 123
9 Flutter Calculation Methods 127

9.1 Non-dimensionalisation . . . . . . . . . . . . . . . . . . . . . . . 128
9.2 k method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
9.2.1 Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
9.3 p-k method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
9.3.1 Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
9.4 p method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
10 Dynamic Responses 135

10.1 Types of Dynamic Response problems . . . . . . . . . . . . . . . 136
10.2 Dynamic Response Equations . . . . . . . . . . . . . . . . . . . . 136
10.3 Solution in the frequency domain: Fourier Transform . . . . . . . 138
10.4 Solution in the time domain: Numerical Integration . . . . . . . 139
10.5 Stochastic Response Analysis - PSD Technique . . . . . . . . . . 141
11 Introduction to Computational Aeroelasticity 147

11.1 Treatment of CA . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
12 Computational structural dynamics 151

12.1 Rayleigh-Ritz finite-element formulation . . . . . . . . . . . . . . 152
12.2 Example: Formulation for an axially-loaded rod . . . . . . . . . 153
12.3 Solution of the equations of motion . . . . . . . . . . . . . . . . . 155
6 CONTENTS
12.4 Numerical Integration multi-DOF systems . . . . . . . . . . . . . 157

12.5 Newmark method . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
12.6 Efficient integration using modal superposition . . . . . . . . . . 158
13 Computational Methods for Unsteady Flows 161

13.1 Linearised potential equation . . . . . . . . . . . . . . . . . . . . 163
13.1.1 Kernel-function method . . . . . . . . . . . . . . . . . . . 163
13.1.2 Panel Methods . . . . . . . . . . . . . . . . . . . . . . . . 165
13.2 Small-disturbance potential equation . . . . . . . . . . . . . . . 167
13.3 Euler and Navier-Stokes equations . . . . . . . . . . . . . . . . . 168
13.4 Discretisation of the Euler and Navier Stokes equations . . . . . 170
13.4.1 Semi-discrete approach . . . . . . . . . . . . . . . . . . . . 170
13.4.2 Fully-discrete approach . . . . . . . . . . . . . . . . . . . 171
13.5 Dealing with moving meshes . . . . . . . . . . . . . . . . . . . . . 171
13.5.1 The Geometric Conservation Law (GCL) . . . . . . . . . 172
13.5.2 The Discrete Geometric Conservation Law (DGCL) . . . 173
13.6 Example: A semi-discrete finite-difference method . . . . . . . . 173
13.6.1 Enforcing the DGCL (optional section) . . . . . . . . . . 175
13.7 Example: A semi-discrete finite-volume method . . . . . . . . . 175
13.7.1 Enforcing the DGCL (optional section) . . . . . . . . . . 176
13.8 Example: A semi-discrete finite-element method
(optional section) . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
13.9 Example: A space-time finite-element method
(optional section) . . . . . . . . . . . . . . . . . . . . . . . . . . . 178
13.10Discretisation Error Analysis . . . . . . . . . . . . . . . . . . . . 180
13.10.1 Grid convergence studies . . . . . . . . . . . . . . . . . . . 181
13.11Phase and Amplitude Error Analysis . . . . . . . . . . . . . . . . 182
13.11.1 Adjoint-based error estimation . . . . . . . . . . . . . . . 185
13.12 Example Unsteady Flow Solutions . . . . . . . . . . . . . . . . . 186
13.13Practice Problems . . . . . . . . . . . . . . . . . . . . . . . . . . 190
14 Computing Fluid-Structure interactions 191

14.1 Generalised force models from CFD computations . . . . . . . . 191
14.2 Coupling of fully-discretised fluids and structures . . . . . . . . . 194
14.3 Partitioned Time Integration Techniques . . . . . . . . . . . . . 195
14.3.1 Volume-continuous methods . . . . . . . . . . . . . . . . . 196
14.3.2 Volume-discontinuous methods . . . . . . . . . . . . . . . 196
14.3.3 Analysis of errors due to partitioning . . . . . . . . . . . . 197
14.3.4 Improved accuracy by extrapolation . . . . . . . . . . . . 201
14.4 Monolithic Time Integration Techniques . . . . . . . . . . . . . . 201
14.5 Efficiency of Monolithic and Partitioned Techniques . . . . . . . 202
14.7 Handy Formulae . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
CONTENTS 7
15 Aeroelasticity in Aircraft Design 207

15.1 Requirements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207
15.2 Design Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
15.3 Wind Tunnel Testing Techniques . . . . . . . . . . . . . . . . . . 209
15.4 Ground Vibration Testing . . . . . . . . . . . . . . . . . . . . . . 210
15.5 Flight testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211
15.5.1 Excitation and Measurement . . . . . . . . . . . . . . . . 211
15.5.2 Ground data processing and damping prediction . . . . . 213
8 CONTENTS
Chapter 1
Introduction
1.1 Background
The field of aeroelasticity considers phenomena in which interactions occur be-
tween aerodynamic flows and elastic structures. In this definition, the word
“interaction” is significant. To be precise, we mean phenomena in which the
aerodynamic forces are a function of the structural displacement, and the struc-
tural displacement is a function of the aerodynamic forces. This is easily demon-
strated by an example, such as the glider shown below:
Figure 1.1: A flexible glider
If the glider’s wing is flexible, then we can expect it to deform. A change

in shape, however, implies a change in the aerodynamic forces acting on it. In
turn, these new forces will impose a new structural deflection. Therefore if we
wish to predict the glider’s response, we must treat interactions by considering
the complete fluid-structure system simultaneously.
The simultaneous analysis of fluids and structures can be complex. It is
worthwhile from the scientific point of view, as such combined systems are often
rich in unexpected and interesting behaviours. It can be essential from the
engineering point of view, as unexpected behaviours are potentially undesirable
ones.
Aeroelasticity is intimately connected with the field of aerospace engineering,
where it has a major impact on configuration design. Its importance has grown
9
10 CHAPTER 1. INTRODUCTION
Modern Aircraft
Maximum Efficiency/Performance
@
@
R
@
Slender Lightweight
@
@
R
@
Flexible !
Figure 1.2: Trends in high-performance aircraft
over the last few decades as the performance of aircraft has increased. High-
performance aircraft tend to be light and slender, which implies reduced induced
aerodynamic drag, but also increased flexibility. In addition, high operating
speeds mean that the changes in aerodynamic force for a given amplitude of
deflection are relatively large, leading to large-amplitude aeroelastic responses.
In the early days of aerospace engineering, the significance of aeroelastic

effects was not completely recognised. Although the Wright brothers did in-
vestigate losses in thrust due to the flexibility of their propellers, the design
of aircraft was mostly based on decoupled analysis until world war I. During
the war, higher aircraft operating speeds brought several unexpected types of
aeroelastic phenomena to light. These included those producing large static
deflections, and several forms of self-sustained oscillations whose rate of growth
could be dramatic. Often such phenomena lead to structural failure. This
meant that they were difficult to diagnose (due to a lack of returning aircraft)
and difficult to research (due to a lack of willing volunteers).
After world war I, the process of carrying out static and dynamic aeroelastic
analysis using simplified models was established. These provided insight into
why aeroelastic phenomena occurred, and allowed for the solution of problems
such as aileron flutter through the use of mass balancing. At that point struc-
tural and aerodynamic models were too crude to be qualitatively accurate, how-
ever, especially for the more complex dynamic responses. In the 1930’s it was
recognised that unsteady aerodynamic effects played a substantial role dynamic
aeroelastic response, which drove major developments in unsteady aerodynamic
theory. In the same decade the basic techniques of flight flutter testing were
established, although safety remained a substantial problem.
Before and during world war II aircraft entered new speed ranges, resulting
in increased aerodynamic modelling complexity due to the effects of compress-
ibility. They also employed more complex structural forms whose dynamics were
difficult to predict. It was during this period that wind-tunnel testing evolved
to the point where it became essential tool for aeroelastic analysis. With the
invention of the jet engine, aircraft began operating in the transonic aerody-
1.2. COURSE OUTLINE 11
namic regime. Here they experienced shock-induced aerodynamic phenomena

which could promote aeroelastic instability. The aerodynamics of the transonic
regime were found to be very difficult to predict accurately.
Since the 1950’s digital computers have been used to support aeroelastic
analysis. During the 1960s the first methods for predicting three-dimensional in-
compressible unsteady flows appeared, and during the 1970’s, discrete structural
prediction techniques advanced to the point where they could begin to predict
the dynamics of realistic aircraft configurations. In the 1980’s and 1990’s, meth-
ods for computing inviscid transonic compressible flows on general geometries
became widely available. Unfortunately, the accurate prediction of flows involv-
ing separation or shock boundary-layer interaction, which are of considerable
interest in aeroelasticity, remains computationally intractable.
Aeroelasticity is a complex research field which there is active participation
by both scientists and engineers. Their goal is not only to predict aeroelastic
phenomena, but also to understand aeroelastic mechanisms well enough to con-
trol them. Therefore, a combination of techniques is presently in use. Where
possible, simplified models are employed to develop maximum insight. Where
practical, numerical computations are used to verify assumptions and make de-
tailed changes in design. Wind-tunnel testing is employed to reduce the risk
associated with new design concepts, while flight testing is used to deliver the
final word on the characteristics of complete configurations.
Until recently, the role of aeroelasticity in aircraft design has been primarily
directed towards analysing and correcting the aeroelastic deficiencies of existing
design concepts. The push towards lighter construction techniques using active
structures, however, is changing this role. As aeroelasticity becomes incorpo-
rated into multi-disciplinary design optimisation (MDO), the design problem
is evolving from “how do we prevent this configuration from exhibiting insta-
bility” to “how can aerodynamics, structures and active controls be optimally
combined to produce a configuration with maximum performance”. Such chal-
lenges, along with those presented by the complex physical behaviours exhibited
in many conditions of interest, will ensure that aeroelasticity remains a vibrant
field for years to come.
1.2 Course Outline

This course will focus on the underlying analytical questions of aeroelasticity:
how does one predict and understand the behaviour of an aeroelastic system?
These two questions are not the same. One can perform a full aeroelastic simu-
lation using hundreds of thousands of elements to discretise the structure, and
tens of millions of elements for the fluid, and perhaps receive a reasonably ac-
curate prediction of a particular aeroelastic response. But which of the factors
taken into account are responsible for the observed response? More importantly,
how can the results be translated in to a design action which ensures the desired
behaviour over the complete range of operating conditions? In order develop
the insight necessary to answer such questions, it is useful to simultaneously
At the end of this course you should:
• Understand the physical processes which drive aeroelastic phenomena

• Be able to formulate and solve simple aeroelastic response and instability
problems
• Be able to identify strengths and weaknesses of different aerodynamic and
structural models for the analysis of a given aeroelastic condition
• Understand the basic design of computational aeroelastic solution tech-
niques
• Be familiar with the role of aeroelasticity in aircraft design
Figure 1.3: Course Goals
consider simplified (preferably analytic) models which condense complex physi-

cal interactions into conceptually manageable forms. This course with therefore
examine the problem of aeroelasticity using models with varying levels of com-
plexity. Familiarity with and comparison of these models will do much to achieve
the course goals, which are summarized in figure 1.3.
The path we will take through the material of the course is shown in fig-
ure 1.4 (chapter numbers are indicated in brackets). We will start by considering
simplified aeroelastic systems, specifically one or two degree-of-freedom struc-
tures with steady and quasi-steady aerodynamic models. These clearly do not
contain all relevant physicals effects, but are still general enough to exhibit
some classic aeroelastic phenomena. Furthermore, the basic dynamic analysis
techniques which will be presented are similar to those used for more complex
physical models.
We will then increase the fidelity of our modelling by considering unsteady
aerodynamics in detail. We will examine a range of theoretical aerodynamic
models, each of which is representative an unsteady physical phenomena which
contributes to aeroelastic response.
Next, we will investigate how general structures can be represented with
using systems with modest numbers of degrees of freedom. Then we will re-
consider dynamic response prediction in the context of these more general de-
scriptions.
Following that, we will briefly review a number of high-fidelity computa-
tional techniques for computing the unsteady behaviour of both fluids and
solids, including boundary-element, finite-element and finite-volume methods.
We will then consider problems which arise when performing fully-coupled fluid-
structure simulations. Finally, we will end with a summary of how the concepts
described in the course are employed in the process of aircraft design.
1.3. ADDITIONAL COURSE INFORMATION 13
Simplified aeroelastic problems (cp 1-2)

|
Unsteady aerodynamics (cp 3-7)
|
Compact models of structural systems (cp 8)
|
General dynamic response prediction (cp 9-10)
|
Computational aeroelasticity (cp 11-14)
|
Aeroelasticity in aircraft design (cp 15)
Figure 1.4: Course Content
1.3 Additional Course Information

1.3.1 Course format:
The course will consist of 14 lectures, an experimental demonstration, two
projects and the final exam. In the lectures, course notes for the relevant sec-
tions will be passed out, and the key associated concepts and problems will be
described and worked out. The experimental demonstration is short, and not
mandatory. However, it is a chance to see some classic aeroelastic phenomena
first hand.
Two projects will be given out. These primarily designed to reinforce the
concepts treated in the course. They are not mandatory, but if completed and
handed in they can count for up to 30% of your mark. You may also write a
100% final exam. Students who hand in projects but achieve a higher average
mark on the exam will receive the exam mark instead. The marking scheme is
summarised below:
• Project part 1: (15%)
Project part 2: (15%)
• Final examination:
(70%) - if project marks are highest
(100%) - if exam mark is highest
Students may work together on a project, but each student is expected to pro-
duce a unique report. (Near-) Identical reports will have their marks divided
by the number of copies. Note that projects handed in late will not be graded.
Finally, there is a blackboard site where you can find announcements, in-
formation about the projects, practice exams, and some additional links and
background information.
1.3.2 Prerequisite Concepts

Aeroelasticity is a truly multi-disciplinary topic, and thus relies on several pre-
requisite concepts. From aerodynamics, these include potential flow, the basics
of compressible flow, and familiarity with thin airfoil theory. From structural
mechanics, one should be comfortable with the concept of virtual work, La-
grange’s equations of motion, and simple beam theory. You should also know
how to analyse multiple degree-of-freedom vibrations, predict free and forced
responses, and be familiar with Fourier analysis and the convolution of step and
impulse responses.
1.3.3 Main References

The primary source of information for the material on the final exam are the
course notes. The references listed below are useful for providing background
information. Reference 1 highlights much of the material in the course, and
contains several modern applications. Reference 2 is less up to date, but com-
prehensive, and perhaps easier to read than reference 1. Reference 3 is the
easiest to read, and gives a good foundation for the course, but is less detailed
than the other two. Do not confuse Reference 2 with another Dover book by
Bisplinghoff and Ashley, titled “Principles of Aeroelasticity”.
• Introduction to Aircraft Aeroelasticity and Loads, J.R. Wright and
J.E. Cooper, Wiley (2007)
• A Modern Course in Aeroelasticity, 4th edition, E.H. Dowell (Edi-
tor), Kluwer Academic Publishers (2004)
• Aeroelasticity, C. S. Bisplinghoff, H. Ashley, and R. L. Halfman, Dover,
(1955)
• An introduction to the theory of Aeroelasticity, Y. C. Fung, Dover
(1955)
Many of the concepts presented in these notes appear in the previous course
notes of Professor Zwaan, last published in 1999. Much of the material is also
different, however. Please note that the exam and assignments are based on the
current course notes.
1.3.4 Instructor Contact

I will be happy to help you with any material from the course, but please make
an appointment before dropping by my office. You can contact me by telephone
at +31-15-278-1538, or by email at S.J.Hulshoff@TUDelft.NL . I can also answer
your questions directly by email.
1.4. CHARACTERISTICS OF AEROELASTIC PROBLEMS 15
1.4 Characteristics of Aeroelastic Problems

Before we continue with the rest of the material it is useful to first introduce
precise terminology to describe different classes of aeroelastic problems.
Basically we will consider four types of problems, in two groups. The first
group are called “Instability Boundary problems” (figure 1.5). These are con-
cerned with determining the boundary which separates a stable fluid-structure
interaction from an unstable one, and can be static or dynamic. Restricting
intrest to an instability boundary often allows simplifications in the analysis
procedure. Not only can the solution search space be limited, but it is also often
possible to use local linearizations of the governing equations. Static problems
are a subset of instability boundary problems for which it is possible to ignore
inertia effects. Dynamic problems are the more general type in which inertia
effects must be included.
The second group of problems are called “Response problems” (figure 1.6).
In this case one computes the response of the system to a given excitation.
These problems can also be either static or dynamic. The last type includes
every type of possible response, and is thus perfectly general. Therefore, a
method developed for computing dynamic response problems can also be used
to find instability boundaries. It just won’t be as efficient as a method which
exploits the fact that it is looking for a boundary.
Determination of Static and Dynamic Instability Boundaries
? ?
Static Dynamic
- negligible time gradients - oscillatory
- Known as “divergence” - Known as “flutter”
∆L ∆L
∆α ∆α
t t
Figure 1.5: Instability problem types
Another way of visualising aeroelastic problem types is the famous Collar’s

triangle, shown in figure 1.7. This shows the positioning of Aeroelastic problem
types in relation to the forces considered within the model. Beside the aeroe-
lastic problems, it is possible to locate several other classic aerospace problems
on the diagram.
Determination of Static and Dynamic Response
? ?
Static Dynamic
Determine the steady Determine the unsteady
operating condition response of an aircraft
of an aircraft taking to an applied excitation:
structural deformation
into account: - pilot input
- atmospheric turbulence
- Elevator angle to - landing loads
trim for a flexible fuse-
lage
δe
Figure 1.6: Response problem types

Forces:
SS A: Aerodynamic
A E: Elastic
CE I: Inertial
D DS
Problem Type:
F DR
F: Flutter
DR: Dynamic Response
D: Divergence
E I CE:
SS:
Control Effectiveness
Static Aircraft Stability
DS: Dynamic Aircraft Stability
MV MV: Mechanical Vibrations
Figure 1.7: Collar’s Triangle (1946)
Effects:
I A: Aerodynamic
E: Elastic
E I: Inertial
A C: Control
T: Thermodynamic
Figure 1.8: Friedmann’s Hexahedron (1999)
The date on figure 1.7 gives you some idea of the context. Recently it has be-
come fashionable to generalise this a bit, as shown in figure 1.8. This is because
high-speed designs may also experience thermal coupling, and there is consider-
able ambition to optimise the response of multi-physics systems using advanced
control concepts. These newer problems fall under the areas of Aeroservoelas-
ticity (Aeroelasticity + high-gain control systems), and Aerothermoelasticity
(Aeroelasticity + thermal effects).
A nice example of Aeroservoelastic engineering is the X-53, an aircraft which
makes use of active leading and trailing edge controls to achieve aeroelastically
deformed wing shapes with optimal performance (figure 1.9). This is a big
contrast to the old way of doing things, (i.e. make the wing stiff enough to
resist the moments incurred by control motion) and can lead to much lighter
structures and increased control effectiveness. The latter is demonstrated figure
Figure 1.9: X53: control surfaces (upper) and in flight (lower)
1.10, which shows that for the X-53, the largest roll control authority is obtained
using leading-edge outboard devices on a lighter, flexible wing.
Figure 1.10: X53 control effectiveness for leading and trailing edge (LE, TE)
inboard and outboard (I,O) devices
Chapter 2
Simplified Aeroelastic
Problems
We will begin our look at aeroelastic analysis techniques using structural models
which have only one or two degrees of freedom (DOF), and using only quasi-
steady aerodynamics. This allows all the basic concepts of aeroelastic analysis
to be introduced, however, and provides a first look at some of the physical phe-
nomena that play a role in realistic systems. The models presented here may be
too crude for use in detailed design, but they can still be applied in preliminary
design stages as their simplicity allows the influence of design variables to be
clearly identified.
We will start by describing simple structural and aerodynmaic models for
a typical airfoil section. Then approaches to deriving both static and dynamic
aeroelastic boundaries will given for the problems of divergence, control reversal
and flutter. Afterwards we will briefly introduce two other simplified systems
which can be used for preliminary analysis, namely the semi-rigid wing and a
model for propeller whirl flutter.
2.1 The Typical Section

For a first demonstration of static and dynamic aeroelastic phenomena it is
sufficient to consider a 2D airfoil section suspended using two springs as shown
in figure 2.1 (right). This is sometimes used to approximate the characterstics of
a complete wing by using its equivalent aerodynamic and structural properties
at the 75% span section. The structural and geometric parameters of the typical
section are given in figure 2.2. The center of rotation is referred to with EA,
short for elastic axis. This terminology comes from three-dimensional wings,
for which it may be thought of as the spanwise-running line which connects the
rotation centres of each section. Note that vertical displacement, h is defined to
be positive downwards. Also, in aerodynamics the chord c is used as a length
scale. In aeroelasticity, the half-chord b is normally used instead (be careful not
21
22 CHAPTER 2. SIMPLIFIED AEROELASTIC PROBLEMS
Typical
@@@@
Section
@@@@ Mean
@@@@
y
Position
z
U Linear Spring
Torsional Spring
Figure 2.1: The Typical Section
@@@@@@
@@@@@@
@@@@@@
Mean
h
@@@@@@ Position
EA
h
Kh
ab
Kθ
EA θ b b
c
Structural Parameters Geometric Parameters

EA: Elastic Axis
b: Half chord
h: Bending deflection (at EA)
c: Chord
θ: Torsional deflection (about EA)
a: Displacement of EA from half chord
Kh: Bending stiffness
S: Reference area
Kθ: Torsional stiffness
Figure 2.2: Structural and geometric parameters
to confuse b with wing span when working out problems).

We will use force and mass quantities such as L, MAC , H, m, Sθ , and Iθ ,
defined in figures 2.3 and 2.3. These are sometimes considerd to be local values
on a three-dimensional wing, measured per unit span. This requires that the
reference area, S, also be defined per unit span, i.e. S = 2b. We will occasionally
work with
p derived quantities such as the uncoupled natural bending frequency,
p
ωh = Kh /m and the uncoupled natural torsional frequency: ωθ = Kθ /Iθ .
Be careful not to confuse these with two the natural frequencies of the coupled
bending-torsion system. For problems involving control surfaces, we consider
an alternate degree of freedom, δ, the control deflection, as shown in figure 2.4.
2.1. THE TYPICAL SECTION 23
EA CG L
h
θ AC
α
U EA
MAC
xθ b
ec
b b
Mass Parameters:
Aerodynamic Parameters
CG: Center of Gravity
xθ: Displacement of CG from EA U: Freestream velocity
m: Mass AC: Aerodynamic center
Sθ: Static moment related to EA (mxθb) e: eccentricity factor
rθ b: Radius of gyration about EA α: Angle of attack
L: Lift
Iθ: Mass moment of inertia related to EA MAC: Moment about AC (independent of α)
(m rθ2 b2) Dynamic pressure ( ρU /2)
2
q:
Figure 2.3: Inertial and aerodynamic parameters
h
EA H
θ
xδ b
δ
Control Surface Parameters
δ: Control deflection
xδ b: Hinge axis position
H: Hinge moment
Sδ: Static moment related to hinge axis
I δ: Mass moment of inertia related
to hinge axis
Figure 2.4: Control surface parameters

2.1.1 Simple Aerodynamic Models

In the following sections we will only consider aerodynamic forces which are an
instantaneous function of the local angle of attack. Although the forces may
change in time due to changes in the angle of attack, they are not influenced by
the time history of the flow.
We will use α to denote the instantaneous angle of attack measured from the
zero lift line of the airfoil. Note that in general this not the same as θ, the pitch
deflection angle. Be aware that in much of the aeroelastic literature, α is also
used for the pitch deflection angle. This is due to the prevalence of linearised
perturbation problems considered by early aeroelasticians.
Steady Aerodynamic Model

For the steady aerodynamic model, we assume that the lift and moment for the
typical section in a freestream with density ρ and speed U is given by:
L(t) = qSCLα α(t) (2.1)

MEA (t) = 2bqS(CMAC + eCLα α(t)) (2.2)
α
EA
U
Figure 2.5: Steady angle of attack
Low-Frequency Aerodynamic Model

The low-frequency aerodynamic model is a little more advanced than the steady
aerodynamic model, as it includes the vertical translation speed in angle of
attack:
ḣ(t)
αLF (t) = α(t) + (Small-angle approximation) (2.3)
U
Thus:
!
ḣ(t)
L(t) = qSCLα α(t) + (2.4)
U
!
ḣ(t)
MEA (t) = qSecCLα α(t) + (2.5)
U
α
ULF EA
h
U
Figure 2.6: Quasi-steady Angle of attack
2.1.2 Equations of motion for the typical section

The equations of motion for the typical section can be derived using small-angle
approximations and Lagrange’s equations or more intuitively from a lumped-
mass point of view. They express equilibrium in the vertical direction, and
about the axis of rotation of the airfoil and control hinge:
X P X
Fz = 0, MElasticAxis = 0, MHinge = 0 (2.6)
In terms of the previously-derived parameters, these can be expressed:
mḧ + Sθ θ̈ + Sδ δ̈ + Kh h + L = 0 (2.7)
Sθ ḧ + Iθ θ̈ + (Iδ + xδ bSδ )δ̈ + Kθ θ − MEA = 0 (2.8)
Sδ ḧ + (Iδ + xδ bSδ )θ̈ + Iδ δ̈ + Kδ δ − H = 0 (2.9)
where θ and δ are small and the aerodynamic moment about the elastic axis is
defined by:
MEA = L(ec) + MAC (2.10)
Often simplified forms of the equations are used in order to clarify the roles
of the structural and aerodynamic properties for certain classes of problems.
For the static problem classes introduced in in chapter 1, for example, one may
ignore the terms with time derivatives.
2.1.3 Torsional Divergence

Divergence was the first type of aeroelastic instability that was recognised and
understood. If one considers the section shown below, it is obvious that a
perturbation in angle of attack will tend to increase θ and thus also α and L.
If the spring is strong enough to resist, there are no consequences. But as the
dynamic pressure q increases, the aerodynamic forces become larger while the
spring stiffness is unchanged. At a certain value of q, the aerodynamic moment
balances the spring moment. The lowest speed at which this happens for a given
density is the torsional divergence speed. Beyond this speed θ increases without
bound in response to a perturbation.
α θ AC
MAC
Kθ@@@@
@@@@
αo
@@@@
@@@@
Zero−Lift
@@@@
EA
Line
U ec
Figure 2.7: Typical section for the divergence problem
The divergence speed can be found by considering the equation of motion

for a single degree of freedom. We will ignore inertial effects for this example,
so it becomes a static instability boundary problem. Summing moments about
the EA (positive nose-up):
X
M = [CLα (αo + θ)e + CMAC ] cqS − Kθ θ (2.11)
P
For equilibrium, M = 0, or:
cqS [CLα αo e + CMAC ]

θ= C ecqS
(2.12)
Kθ 1 − LαKθ
examining the denominator, θ is seen to be unbounded when CLα ecqS = Kθ .

This gives a certain value for qdivergence . It can be shown that this is an upper
limit for q if one considers the stability of the system to small perturbations, for
which we require:
∂M
≤0 (2.13)
∂θ
Equation (2.13) simply says that if there is a positive perturbation in θ, then a
negative (nose-down) opposing moment must be produced. Applying it to the
moment equation leads immediately to:
ecqS
1− CLα ≥ 0 (2.14)
Kθ
from which can be concluded that the configuration is unconditionally stable
if e ≤ 0 (AC downstream of EA), but only conditionally stable if e > 0. The
former is the case for weather vanes, while the latter is often true for realistic
wings. We can also re-arrange (2.14) to obtain:
Kθ
q≤ (2.15)
CLα ecS
Which explicitly shows the configuration is limited to a maximum operating

dynamic pressure as a function of aerodynamic and structural parameters.
2.1.4 Control Reversal

The actuation of trailing-edge control devices normally changes the pressure
distribution so that the center of pressure shifts aft. Examining figure 2.8, it is
clear that a positive control deflection will therefore give a nose-down pitching
moment that tends to reduce the angle of attack of the configuration. A speed
will be reached, qreversal , where a positive control deflection decreases the angle
of attack sufficiently to have no change in force. At speeds higher than this
speed the control will act in the sense opposite to normal.
L
AC @@@@
Kθ@@@@
θ @@@@
@@@@
MAC
@@@@
U EA
Figure 2.8: Typical section for the control reversal problem
At equilibrium:
X
M = [(CLα θ + CLδ δ) e + CMacδ δ] cqS − Kθ θ = 0 (2.16)
solving for the ratio of structural twist to control deflection:

θ cqS CLδ e + CMacδ
= (2.17)
δ Kθ 1 − CLα ecqS
Kθ
Control effectiveness is indicated by lift/(rigid-wing lift) :

C C cqS
L (CLα θ + CLδ δ) 1 + LαCLM K acδ
θ
(R)
= = C
δ
ecqS
(2.18)
L CLδ δ L
1 − Kθ α
Noting CLα , CLδ , S, c, Kθ > 0, and CMacδ < 0, there will be a q > 0 at which
L
L(R)
will switch from positive to negative, denoted as qreversal :
CLδ Kθ
qreversal = − (2.19)
CLα CMacδ cS
Above this speed, the loss in lift due to (nose-down) structural deflection
will be greater than the lift gained by (positive) control deflection
Figure 2.9: B-47
One famous example of unexpected control reversal cropping up in a design

is that of the B-47 (figure 2.9). The torsion-strip-theory originally used did not
predict aileron reversal within the envelope. It turned out that bending effects
were significant as well, which is perhaps not surprising when you consider
that the difference in wing-tip deflections between the maximum positive and
negative loads is 11m for this aircraft. The early solution to this problem was to
use roll spoilers, which have a relatively low CMacδ . This is one of the reasons
that roll spoilers are used modern airliners, along with the fact that they are
also quite effective when placed ahead of a flap.
Note that flight beyond qreversal can be made by actuating the controls
in the opposite sense. To do so, however, the pilot would first need to spend a
significant amount of time with little or no control while flying near the relatively
high speed of qreversal .
2.1.5 Flutter of a Typical Section

Flutter is a word used to describe self-sustained oscillations arising from fluid-
structure interactions. It is normally implied that under certain conditions, the
oscillations can grow without bound. The flutter point is defined as the point
where a damped oscillation transitions to one growing in amplitude. A group of
flutter points therefore forms a boundary between a region where one can expect
to operate a vehicle “normally” from one where structural failure is likely. The
region within the boundary is referred to as a flutter envelope (figure 2.10).
Determination of the flutter envelope forms a major and mandatory part of
aircraft certification, as it must be shown to be 15% larger than the aircraft’s
flight envelope (10% if the aircraft has sustained damage).
The fact that the flutter point is a transition between damped and growing
oscillations implies that the response at the flutter point is of constant ampli-
tude. This idea is used repeatedly in analysis techniques, and often makes the
problem of determining a flutter point easier than computing the general re-
sponse of a configuration. In spite of this, determining the flutter boundary of
an aircraft is a major job. One reason for this is the complexity of the aero-
dynamic phenomena which tend to exist in flutter conditions. Not only must
unsteady aerodynamic effects be included, but the fact that the configuration
is far from its normal operating point can lead to the appearance of difficult-
to-deal-with features such as strong shocks and boundary-layer separation. In
Flutter Point
Altitude
Decaying Growing
Oscillation Oscillation
ary
nd
ou
rB
tte
Flu
Speed
Figure 2.10: A flutter boundary
addition, one must normally consider several types of flutter modes, for which
variations with fuel, payload and aircraft configuration must be established.
The combination of these factors translates into a large number of flight cases
which must be analysed.
To keep things simple, we start our analysis with the case of “Classical
Flutter” which refers to a class of two-degree of freedom oscillations which are
characteristic of 2D airfoil sections. Our approach will be to assume a solution of
an amplified harmonic form, substitute the assumed solution into the equations
of motion, and then solve for the unknown frequencies and amplification factors.
Our task will then be to determine under which conditions the amplification
factors transition from negative (damped case) to positive (growing case). The
analysis is precise (within the bounds of the assumed simple aerodynamic and
dynamic behaviours) but does not necessarily give a phenomological answer
to the question “What is flutter?”. For that, it is useful to appeal to energy
concepts, which at the very least can account for the existence of the phenomena.
When an aircraft travels in relative motion to air, it is in fact contact with
an energy reservoir. More succinctly, certain motions of the aircraft can extract
work from the airstream. If the speed (and structural integrity) of the aircraft
is maintained, such energy extractions will ultimately show up in the aircraft’s
fuel bill. So how is energy transferred to the structure? In our analysis we
will ignore the drag force, so that the transfer of energy must occur when the
airfoil’s vertical speed vector is aligned with the lift vector. When the lift and
vertical speed have the same sign, work is done on the structure (figure 2.11).
When the lift opposes the motion, energy is extracted (this is similar to what
occurs when pushing in and out of phase on a child’s swing).
If we consider an airfoil which can plunge but not pitch (h 6= 0, θ = 0), the
low-frequency aerodynamic model predicts that the motion is always damped
(check this yourself). The key to classical flutter, therefore, is that the timing
of the motion in θ and h must be such that over one cycle, the lift is on av-
Lift
(L) L(t)
Vertical
Speed (z) z(t)
−
+ +
+ t
− −
Figure 2.11: Timing of vertical motion and aerodynamic force
erage
R aligned with the airfoil motion, rather than opposing it (in other words
L(t)ż(t)dt > 0). Such a situation is shown in figure 2.11.
Classical flutter with a growing amplitude can therefore occur with the
proper timing of θ and h, which implies the proper timing of aerodynamic forces
with structural motion. However, non-classical flutter modes are also possible.
For example, single-degree-of-freedom flutter can occur when unsteady aerody-
namic forces have a phase different from that applied by the steady aerodynamic
model. A famous example is single-DOF control-surface flutter, which can occur
in transonic conditions. In this case, shock movement as a function of δ is lagged
such that the net aerodynamic damping is negative. Another possibility is the
pitching airfoil at very low speeds, for which wake effects can also introduce
negative damping.
Both of the single-degree-of-freedom cases above highlight the importance of
unsteady aerodynamic effects in determining the timing of the forces and overall
energy balance. Although in many cases the steady aerodynamic model we are
about to use may perform reasonably, it can be insufficient for the conditions of
most relevance for modern aircraft, the transonic regime. We will examine the
reasons behind the limits of the steady aerodynamic model in Chapters 3-6.
2.1.6 Solution of the equations of motion

First we will consider the various motions possible for a 2-DOF (pitch and
plunge) airfoil section with a steady aerodynamic model. To keep things sim-
ple, we will consider uncambered airfoils only (CMac = 0) so that the lift and
moment about EA is:
L(t) = qSCLα θ(t) (2.20)

MEA (t) = e2bqSCLα θ(t) (2.21)
The equations of motion can then be written:
mḧ + Sθ θ̈ + Kh h + qSCLα θ(t) = 0 (2.22)

L
ec
AC h
EA CG
MAC θ
ab
xθ b
b b
Figure 2.12: The typical section
Sθ ḧ + Iθ θ̈ + Kθ θ − 2ebqSCLα θ(t) = 0 (2.23)
or in matrix form (with {x} = [h, θ]t ):
[M ]{ẍ} + ([K] − q[A0 ]) {x} = {0} (2.24)
where

m Sθ Kh 0 0 −SCLα
[M ] = [K] = [A0 ] = (2.25)
Sθ Iθ 0 Kθ 0 2SebCLα
This form is typical for the larger systems we will consider later in the course.
To find solutions to the system, we begin by assuming that the 2 DOF response
ˆ pt so that:
of the section, {x} = [h, θ]T , is of the form {x} = {x}e
h = ĥept and θ = θ̂ept where p = (σ + iω)
This assumption therefore includes diverging, damped and neutral responses,

both periodic and aperiodic. The equations of motion are then:
(mp2 + Kh )ĥ + (Sθ p2 + qSCLα )θ̂ = 0 (2.26)

2 2
Sθ p ĥ + (Iθ p + Kθ − 2ebqSCLα )θ̂ = 0 (2.27)
or
[M ]p2 + ([K] − q[A0 ]) {x} = 0

(2.28)
where for the latter q, the dynamic pressure, has been separated out to indicate
its importance as a parameter. For non-trivial {x}, the determinant of the
qf1 qf2 qd
ω σ
1 4 1 3
q q
Figure 2.13: Flutter diagram for 2DOF section with Kθ = 0.05, xθ = 0.2
coefficient matrix must be zero. This requirement leads to the characteristic

equation, which can be written:
a4 p 4 + a2 p 2 + a0 = 0
a4 = mIθ − Sθ2
a2 = mKθ + Iθ Kh − (2meb + Sθ )qSCLα
a0 = Kh (Kθ − 2ebqSCLα )
This can be solved for its roots:

r
1
q
p1..4 = (σ + iω)1..4 = ± (−a2 ± a22 − 4a4 a0 ) (2.29)
2a4
so that the solution to the equations of motion can be written:
{xj } = {x̂j }e(σj +iωj )t ; j = 1..4 (2.30)
To get an idea of how the motion changes over the operating range of a given
configuration, we can plot the components of p as a function of the dynamic
pressure. Such diagrams are known as “flutter diagrams” and are frequently
used to express the overall characteristics of a system (they can also be plotted
with speed as the ordinate). Note that the imaginary component of p corre-
sponds to the frequency of oscillation, while the real component corresponds
to the growth rate. While plotting we reject p’s with negative frequencies as
unphysical solutions.
qf qd
ω σ
1 4 2 3
q q
Figure 2.14: Flutter diagram for 2DOF section with Kθ = 1.0, xθ = 0.2
qd
ω σ
1 3
q q
Figure 2.15: Flutter diagram for 2DOF section with Kθ = 1.0, xθ = −0.1
A number of distinct behaviours can be identified in figures 2.13 to 2.15. In figure

2.13, there are initially two constant-amplitude (neutral) modes with frequen-
cies that decrease with increasing q (motion type I). Then the two frequency
branches merge, and one decaying and one growing oscillatory mode appear
(motion type II). The latter is a flutter mode, so that the dynamic pressure
qf 1 thus corresponds to a flutter boundary. As q is further increased, the mo-
tion becomes neutral again ( qf 2 ). Increasing q further, the previous aperiodic
(ω = 0) modes, suddenly cease to be neutral (at qd ) , and one starts experience
positive growth rates while the other decays. The one with the positive growth
rate corresponds to the dynamic version of the divergence mode we studied
previously.
For figure 2.14, the stiffness of the configuration has been increased by a
factor of 20. Initially the mode frequencies are therefore much larger, but these
too decrease with increasing q. A flutter point is again encountered, but in this
case it merges into the divergence mode without an intermediate stable region
as q is increased.
For figure 2.15, the center of mass has been moved forward. The frequency
branches remain well separated, and flutter does not occur. The divergence
mode again appears at higher q.
The coalescing of frequency branches on flutter diagrams is often used as
an indication of potential flutter problems. This is not proof that flutter is
occurring, however, as σ must also be positive for one of the modes. Furthermore
flutter can occur when the branches are separated, as we will see later in this
chapter. For a full interpretation, both parts of the flutter diagram are required.
2.1.7 Mode Shapes

Associated with each value of p is a relative amplitude of pitch and plunge. It
can be found by substituting the value of p into either equation (2.26) or (2.27)
and solving for ĥθ̂ . This ratio is sometimes referred to as the mode “shape”,
for reasons that will become obvious when we consider continuous structural
models. In fact if p is considered to be an eigenvalue of the system, then
the mode shape is its corresponding eigenvector, where one component can
be assumed to be unity and the other obtained through the ratio derived above.
Note that the eigenvector, or mode shape, changes with dynamic pressure, as p
itself is a function of dynamic pressure.
2.1.8 Conditions which determine the types of motion

In this section we will examine in more detail how the parameters of the problem
effect the resulting motion types. We will assume that for a given configuration,
the parameters are fixed, except for the dynamic pressure q. By definition, we
would like to be sure that any undesirable motions occur at values of q which
are outside the configuration’s normal range of operation.
We will focus our analysis on equation (2.29), as the form of p, (real, imagi-
nary or complex), determines the type of motion present in the airfoil’s response.
Doing so we first note that a4 must always be positive, as rθ ≥ xθ . As a result,

we need only consider four sets of conditions.
Set 1: a2 2 − 4a4 a0 > 0, a0 > 0, a2 > 0 .

In this case:

1
q
p21..2 = −a2 ± (< a2 ) = −ω12 , −ω22
2
2a4
so that p1..4 = ±iω1 , ±iω2 . We reject negative frequencies as unphysical solu-
tions, and are left with two neutral harmonic motions.
Set 2: a2 2 − 4a4 a0 > 0, a0 > 0, a2 < 0 .

This is similar to set 1, but now the negative a2 coefficient produces two
positive p2 , so that the p are pure real: p1..4 = ±σ1 , ±σ2 . Here we have
two diverging and two converging aperiodic motions.
Set 3: a2 2 − 4a4 a0 > 0, a0 < 0, a2 ± .

This set produces both of the motions defined above.

1
q
2 2
p1..2 = −a2 ± (> a2 ) = r1 (−ve), r2 (+ve) (2.31)
2a4
so that p1..4 = ±σ1 , ±iω2 In this case, rejecting the negative frequency, we have
one converging and one diverging aperiodic motion, and
one neutral harmonic motion.
Set 4: a2 2 − 4a4 a0 < 0, a0 ±, a2 ± .

Finally, we have more general motions, defined by:
1
p21..2 = (−a2 ± iK) = −m + in (2.32)
2a4
This produces the following possible motions:


 σ1 + iω1 oscillatory diverging
−σ2 + iω2 oscillatory converging

p=

 σ3 − iω3 negative frequency - rejected
−σ4 − iω4 negative frequency - rejected

The above results are summarised in table 2.1.

a22 − 4a4 a0 >0 <0

a0 >0 <0 ·
a2 >0 <0 · ·
p2 −ω12 , −ω22 σ12 , σ22 σ2 , −ω 2 −m ± in
p ±iω1 , ±iω2 ±σ1 , ±σ2 ±σ, ±iω ±σ ± iω
Harmonic: Aperiodic: Aperiodic: Oscillatory:

2 pos. freq. 2 diverging 1 diverging 1 div. pos. freq.
Type of Motion 2 neg. freq. 2 converging 1 converging 1 conv. pos. freq.
Harmonic: 1 div. neg. freq.
1 pos. freq. 1 conv. neg. freq.
1 neg. freq.
Type of motion Neutral Divergence Diver./Neutral Flutter

Set 1 2 3 4
Table 2.1: Responses of a typical section with a steady aerodynamic model
2.1.9 Pines’ Conditions for flutter

Flutter can be defined as a self-sustained oscillation which, under certain con-
ditions can grow in time. The flutter point is defined as the point where the
self-sustained oscillation transitions from convergent to divergent motion. For
a given typical section combined with a steady aerodynamic model, the motion
will be a function of q (note how a2 and a0 are functions of q). For certain
groups of parameters, there will be a flutter point, qf , beyond which an oscil-
latory motion diverges. Given the local density (i.e. altitude and temperature)
this can be converted into a flutter speed, Vf .
From the previous section, we can see that the type of motion we are inter-
ested in occurs under condition set 4. It can therefore be useful to re-express
the condition a2 2 − 4a4 a0 ≤ 0 in terms of the section’s geometric and mass
properties, with the assumption that the only flutter solutions of interest are
for real and positive qf . We define:
qSbCLα
Q = (2.33)
Kθ
and require it to be real and positive. We can then re-express a2 2 − 4a4 a0 < 0
as:
C2 Q2 + C1 Q + C0 < 0 (2.34)
with:
C2 = (xθ + 2e)2 (2.35)

" 2 #
x2

wh
C1 = −2 xθ + 2e + (xθ − 2e + 4e 2θ ) (2.36)
wθ rθ
" 2 #2 2
x2

wh wh
C0 = 1− + 4 2θ (2.37)
wθ rθ wθ
We can find Q at the the flutter boundary using the quadratic formula:
p
−C1 ± C12 − 4C2 C0
Q1 , Q2 = (2.38)
2C2
Requiring Q to be real leads to:
C12 − 4C2 C0 > 0 (2.39)
which can be expressed:

" 2 #
wh xθ
xθ xθ + 2e − 2e 1 + 2e 2 >0 (2.40)
wθ rθ
In order to derive a second condition we note that:
Q1 + Q2 = −C1 /C2 (2.41)
We consider the case where both Q1 , Q2 are positive. Since C2 is necessarily

positive, C1 should be negative. This leads to :
2
x2

wh
xθ + 2e + xθ − 2e + 4e 2θ > 0 (2.42)
wθ rθ
Conditions (2.40) and (2.42) are know as Pine’s conditions. They can be used
to investigate the possibility of flutter for various combinations of geometric,
stiffness, and inertial parameters.
2.1.10 The effect of a low-frequency aerodynamic model

Changing the aerodynamic model can have a profound effect on the observed
flutter behaviour. This can be demonstrated just by considering the low-
frequency aerodynamic model, which introduces an additional damping term,
[A1 ], into the equations:
q
[M ]{ẍ} − [A1 ]{ẋ} + ([K] − q[A0 ]) {x} = {0} (2.43)
U
where

−SCLα 0
[A1 ] = (2.44)
SecCLα 0
To solve the system we once again assume a solution of the form:
{x} = {x̂}ept (2.45)
Which leads to the characteristic equation:
a4 p 4 + a3 p 3 + a2 p 2 + a1 p + a0 = 0 (2.46)
(2.47)
σ
ω qf
q q
Figure 2.16: Flutter diagrams with the steady and low-frequency aerodynamic
models (∗) - steady aerodynamic model (x,+) - low-frequency aerodynamic
model
where
a4 = mIθ − Sθ2 (2.48)

q
a3 = SCLα (2ebSθ + Iθ ) (2.49)
U
a2 = mKθ + Iθ Kh − (2meb + Sθ )qsCLα (2.50)
q
a1 = SCLα Kθ (2.51)
U
a0 = Kh (Kθ − 2ebqSCLα ) (2.52)
As before, the roots of the characteristic equation define the motion types.
Unfortunately, it is now a full-forth order equation so it can not be so easily
manipulated. We therefore resort to its numerical solution. Some example
results are shown compared to their steady aerodynamic model counterparts in
figure 2.16.
The low-frequency aerodynamic model results show the same initial behaviour
of the mode frequency branches, but these become substantially different as q is
increased. Note that there is no longer an intersection of the frequency branches.
Examining the σ plot, we observe that at small q the motion is damped, as we
might have anticipated from simple reasoning, but then for one oscillatory mode
σ grows rapidly, so that the critical dynamic pressure is more than 30% less than
that computed with the steady aerodynamic model.
2.2. THE SEMI-RIGID WING 39
2.2 The Semi-Rigid Wing
Figure 2.17: The Semi-Rigid Wing
Sweep is often used in high-speed aircraft, not only to raise the critical Mach
number of a wing, but also to favourably affect its aeroelastic properties. A very
simple model which can be used to illustrate this is the semi-rigid wing. This
model combines bending and torsion with wing sweep as a parameter.
7777
7777 Kγ
7777
7777 U
Λ: Sweep angle
7777
7777 θ: Twist
Kθ7777
7777 γ: Bending angle
7777
7777 Kθ : Torsional stiffness
7777 ycp
7777
7777 Λ
Kγ : Bending stiffness
7777
7777
ecp Iθ : Torsional inertia
7777
7777
Iγ : Bending inertia
7777 Iθγ : Product of inertia
7777 y S: Total wing area
γ θ b: Semi-span
x b ȳ: Coordinate along
elastic axis
x̄: Coordinate normal
to elastic axis
ȳcp , ēcp : Center of pressure
location
Figure 2.18: Semi-rigid wing parameters

Figure 2.19: X29 forward-swept aircraft
2.2.1 Equations of motion for the Semi-rigid wing

In this case we assume that the system is in equilibrium about both axes of
rotation, so that:
X P
Mx̄ = 0, Mȳ = 0 (2.53)
Iθ θ̈ − Iθγ γ̈ + Kθ θ − L ēcp = 0 (2.54)

Iγ γ̈ − Iθγ θ̈ + Kγ γ − L ȳcp = 0 (2.55)
The lift is usually assumed to be a function of the time history of the local
effective angle of attack:
L = L(α(t)) (2.56)
here α = αo + αd , where αo is the initial angle of attack, and αd is the angle of

attack due to deflection:
αd = θcos(Λ) − γsin(Λ) (2.57)
Following a procedure similar to that used for the typical section (see the prob-
lems at the end of this section), the divergence speed for the semi-rigid wing
can be derived. Doing so reveals that forward-swept wings (e.g. figure 2.19)
have lower divergence speeds. If aft sweep is used, there will be a point where
the effect of bending cancels the effect of torsion and no divergence is obtained.
This is referred to as the aeroisoclinic condition.
2.3. PROPELLER WHIRL FLUTTER 41
2.3 Propeller Whirl Flutter

Propeller whirl flutter arises with the correct sychronisation of aerodynamic
forces with gyroscopic-coupled structural modes. It was first analysed as early
as the late 1930’s , but was not brought to the forefront until 1960, after two
crashes of the Lockheed Electra 188 turboprop, a modified version of which is
shown in figure 2.20. Both crashes occurred during passenger-carrying service.
Figure 2.20: Lockheed Electra turboprop
For these cases, it is believed that the whirl flutter mode also induced flutter
in the main wing, leading to the wing’s separation from the aircraft. Photos
of a model of the aircraft before and after wind-tunnel testing are shown in
figure 2.21.
Figure 2.21: L188 model in the NASA Langley Transonic Dynamics Tunnel
We will examine the basics dynamics of whirl flutter by considering a simple

two degree of freedom system, with appropriate simplifications. Although the
following discussion is focused on the propeller whirl case, similar analyses could
be carried out for other interesting configurations, including wind turbines and
helicopter rotors.
2.3.1 A Simplified model

We consider a propellor rotating with a shaft at a very high angular velocity Ω
(figure 2.22). The shaft is stiff, but is mounted such that it can deflect with the
small angles ψ and θ. These motions are resisted by the structural stiffnesses
Kψ and Kθ . An additional aerodynamic force FA operates at the propeller
end of the shaft. The propeller is assumed to have 3 or more blades, so it can
reasonably be approximated as a disk. Because Ω is large, we assume that the
Κψ
Κθ
θ y
x ψ Ω
FA
b
x’
Figure 2.22: 2-DOF model of whirl flutter
angular momentum in rotor-fixed coordinates about the axis x′ , HΩ , is much

larger than that of the other two directions, i.e.:
HΩ >> Hψ̇ , Hθ̇ (2.58)
as a consequence:
dHΩ
=0 (2.59)
dt
Conservation of angular momentum requires:
dH~
~ = ~
X
+ ~ω × H M (2.60)
dt
where the vectors are measured with respect to a frame moving with ω ~ =
(0, θ̇, ψ̇). For ~
ω × (HΩ , 0, 0) >> ω
~ × (0, Hθ̇ , Hψ̇ ) (2.60) becomes:
Iy θ̈ + Ix Ωψ̇ = −Kθ θ − FAz b + MAy (2.61)

Iz ψ̈ − Ix Ωθ̇ = −Kψ ψ + FAy b + MAz (2.62)
where F~A and M~A are the aerodynamic forces and moments, and b is the shaft
length. These equations are coupled both gyroscopically and aerodynamically.
2.3. PROPELLER WHIRL FLUTTER 43
2.3.2 Gyroscopic Coupling

If the rotor operates in a vacuum (FA = MA = 0), the gyroscopic moment terms
Ix Ωψ̇ and −Ix Ωθ̇ will act to produce two coupled motions:
1. a forward whirl mode (frequency higher than uncoupled frequencies)
2. a backward whirl mode (frequency lower than uncoupled frequencies)
z z
y y
Ω Ω
x x
Figure 2.23: Forward (left) and backward (right) whirl modes
Both motions will be circular if Iy = Iz
2.3.3 Aerodynamic Coupling

Aerodynamic forcing can act to produce unstable (flutter) modes. How this
occurs can be visualised using figure 2.24. Consider first the angle of attack
θ Ω
x U
U θ
r −Ωr
x’
Figure 2.24: Aerodynamic Coupling
of a section on the up-going blade. If θ > 0 is positive, the angle of attack

will be increased. Conversely, the section angle of attack on down-going blade
decreased for θ > 0. This results in MAz < 0 when ψ = 0 and θ = θmax . The
result is that the aerodynamic coupling tends to drive the backward whirl mode,
and oppose the forward whirl mode.
2.4 Practice Problems
p2.4.1 Derive and expression for the qdivergence of the semi-rigid wing (hint:
first use the equations of motion and angle of attack equation to derive an
explicit expression for αd ). Plot the variation of the divergence speed with wing
sweep. Derive an expression for the aeroisoclinic condition.
2.4. PRACTICE PROBLEMS 45
p2.4.2 An experimental setup for measuring control-surface flutter is shown

in left part of the figure below. The airfoil section can translate in the vertical
direction, but cannot rotate. Motion in the vertical direction is resisted by a
spring, Kh . The control surface is free to rotate about its hinge (Kδ = 0). The
mass of the complete configuration is m, while the static moment and moment
of inertia of the control surface at the hinge axis are Sδ and Iδ , respectively.
555555555
555555555
555555555
555555555
K
@@@ h
@@@
@@@
@@@
555555555
555555555
@@@
555555555
555555555
@@@
555555555
555555555
L @@@
@@@ h
@@@
@@@
@@@
H
@@@
@@@
@@@
@@@
@@@
@@@ δ
!!!!
!!!!
@@@
555555555
555555555
@@@
555555555
555555555 @@@@
!!!!
@@@
555555555
555555555
@@@
555555555
555555555
@@@
@@@
@@@@
!!!!
@@@@
@@@
@@@@
a) Derive the characteristic equation for this configuration. Assume the lift
is given by L = [CLδ · δ]qS and the aerodynamic hinge moment by
H = [CHδ · δ]cqS (i.e. a steady aerodynamic model).
b) What types of motion are possible if the experimenter forgets to hook

up the spring (Kh = 0)?
c) Assume an aerodynamic balance (see figure above, right) has been added
to the control surface, such that the net aerodynamic hinge moment is
zero (CHδ = 0). Draw flutter diagrams for the configuration for both
positive and negative values of Sδ (assume Kh 6= 0 and Sδ2 < mIδ ).
p2.4.3 For the typical section with two degrees of freedom (θ, h), under what
conditions is flutter possible if the AC lies ahead of the EA and the EA lies
ahead of the CG?
p2.4.4
z
Kθ
x
θ ψ
y
Kψ
The motion of a ducted fan mounted on a wind-tunnel support stand is

described by:
Iy θ̈ + Ix Ωψ̇ + Kθ θ − a b qθ = 0
Iz ψ̈ − Ix Ωθ̇ + Kψ ψ − a b qψ = 0
where a is the lift-curve slope of the engine nacelle, q is the dynamic pressure,
and Ω is the engine rotation speed.
a) Derive the characteristic equation for this system.
b) Assume Iz = Iy = 1, Ix = √12 , Kθ = Kψ = 1, ab = 2 and Ω = 2. Would

you permit this setup to be run at the test condition of q = 1?
Chapter 3
Unsteady Aerodynamics
In the previous chapter, it was demonstrated how energy can be added to the
structural system if we have the correct phase between structural motion and
aerodynamic forces. But what if the aerodynamic forces are themselves not in
phase with the instantaneous position of the aerodynamic surfaces? This is pre-
cisely what occurs when the effects of unsteady aerodynamics are significant. In
fact, not only may the relative phases be large, but the magnitudes of unsteady
forces can also differ considerably from their quasi-steady counterparts. As a
result, neglecting unsteady aerodynamic effects can result in substantial errors
in predictions of aeroelastic behaviour.
There are several physical processes which contribute to unsteady aerody-
namic effects. For airfoils, important examples include the influence of vortical
convection, compressible effects such as time delays due to the finite speed of
sound and shock movement, and viscous effects including flow separation and
boundary-layer dynamics. We will begin this chapter with some qualitative
descriptions of these phenomena. Then we will examine a hierarchy of model
equations which can be used to predict all or some of their behaviour. We will
finish by listing some of the analytic solutions which have been derived for these
equations to date.
3.1 Vorticity Dynamics

If viscosity is neglected, one can show that the model equations for fluid flow
express (among other things) the convection of vorticity along streamlines. The
dynamics of such convecting vorticity can profoundly effect the forces gener-
ated on nearby solid bodies. To see why this is so, we will consider a thought
experiment based on the theory of incompressible potential flow for an airfoil.
You should be familiar with the idea of the lift force on an airfoil being
related with its value of “bound circulation”, often associated with a point
vortex at the airfoil 1/4 point. In steady aerodynamic theory, the value of
the bound circulation is determined by the Kutta condition, a condition which
47
48 CHAPTER 3. UNSTEADY AERODYNAMICS
Figure 3.1: A flow with compressible effects (e.g. shocks), viscous effects
(e.g. shock-boundary layer interaction, separation), and vorticity dynamics (e.g.
wake and turbulent regions).
prevents a singularity from appearing at an airfoil’s sharp trailing edge. It is

useful to consider the Kutta condition in more detail. To start with, note that
there is nothing in the potential flow equations themselves which dictate that
no singularities can be present in the flow (in fact singularities are often used
as building blocks to construct potential flow solutions). Thus the zero bound
circulation solution shown in figure 3.2 is a perfectly valid one.
Figure 3.2: Zero-circulation potential flow solution
Now consider the impulsive start of an airfoil from rest. Initially, the flow
pattern will indeed resemble figure 3.2. Zooming into the trailing edge, however,
reveals some basic problems for a real viscous flow (figure 3.3). At the trailing
edge, the sharp geometry will lead to velocities which tend towards infinity,
implying very low pressures. Meanwhile, at the stagnation point on the rear
3.1. VORTICITY DYNAMICS 49
+
+
−
−
− −
Figure 3.3: Pressure field of the zero-circulation potential flow solution
upper surface, the pressure reaches its highest possible value. Therefore the flow
which tracks around the trailing edge towards the stagnation point experiences a
severe adverse (positive) pressure gradient. From boundary layer theory, adverse
pressure gradients in wall-bounded flows tend to induce separation. This can
be seen by considering the boundary-layer streamwise momentum equation (
see equation (3.7) later in this chapter) for a constant viscosity evaluated at the
airfoil surface, where the velocities are zero due to the no-slip condition:
∂ 2 ui ∂pe
µ = (3.1)
∂x21 ∂xi
Where x1 is the wall-normal direction. Thus in a positive pressure gradient,
the second derivative of the profile in the normal direction also tends to be
positive at the wall. This is characteristic of profiles near separation, as shown
in figure 3.4.
Thus there is separation from the sharp trailing edge, leading to the forma-
tion of a vortex, known as the “starting vortex” above the airfoil (figure 3.5).
This is in turn carried downstream with the flow (figure 3.6). Kelvin’s theorem
states that the total vorticity within a closed contour which convects with the
flow remains constant. If we draw this contour around the airfoil alone we re-
cover the circulation of the bound vortex. If we draw it large enough to enclose
both the bound and starting vortices, we see that the total vorticity must be
zero (its value before the impulsive start) and that the starting vortex is equal
and opposite in strength to the bound vortex.
While the starting vortex is in the immediate vicinity of the airfoil, its in-
duced velocity field tends to lower the local angle of attack. The result is a lift
value which is lower that of the final steady-state value (in an incompressible
inviscid flow, the initial value of airfoil lift is roughly half that of its final value).
As the starting vortex convects downstream, its influence wanes so that the final
steady lift value is approached asymptotically.
x1
u( x1 )
∂ 2 ui
Figure 3.4: A profile with ∂x21
> 0 near the wall
Figure 3.5: Formation of the starting vortex
Any change in circulation of an airfoil is associated with the shedding of

vorticity with equal magnitude and opposite sign. If we consider an oscillating
profile, it is clear that it must be continuously shedding vorticity increments
(figure 3.7). The sign of these increments is determined by the position within
the oscillation cycle. The resulting oscillatory wake can induce a significant
phase lag and magnitude change of the aerodynamic forces relative to those
computed with a purely steady theory.
The wake vorticity also induces a velocity field which influences the shape
of the wake itself (remember vorticity is carried with the local velocity). This
leads to complex wake shapes as time evolves, a factor which can considerably
complicate numerical methods which rely on wake tracking.
3.1. VORTICITY DYNAMICS 51
Γ
−Γ
Airfoil Starting vortex
Figure 3.6: Circulation within an initial contour
Figure 3.7: Wake vorticity due to oscillating airfoil

3.2 Compressibility effects

The inviscid compressible equations can also be used to show that acoustic
waves propagate with the finite speeds u + c and u − c, where u is the local
flow velocity and c is the local value for the speed of sound. Steady supersonic
airfoil theory demonstrates that for Mach numbers greater than one, this results
in finite zones of influence and dependence in space. For unsteady flows, there
are finite zones of of influence and dependence in space-time, even in subsonic
conditions.
Finite acoustic wave speeds imply that the state of a given point is influenced
by its neighbours at different times. Consider the airfoil shown in figure 3.8. If
(u+c)local (u−c)local
Figure 3.8: Acoustic wave propagation in a subsonic flow
the airfoil is moved suddenly, the three points shown in the figure will immedi-
ately begin to adjust their local states to the new boundary conditions. Their
transients won’t be complete, however, until the information reaching them from
all other points on the airfoil remains fixed. The speed at which the information
travels depends on the wave direction (±c) and on the local velocity produced
by the flow u. If u is of a comparable size to c, then the point in the middle
will hear from the downstream point at a much later time, resulting in a pro-
longed transient, and a complex local behaviour in time. Note that information
travelling from the airfoil’s wake can experience significant delays, resulting in
significant response time lags.
As the Mach number is increased into the transonic region, shocks will begin
to appear. A typical transonic pressure flow topology, and its associated pressure
distribution is shown in figure 9. It is clear from the pressure distribution that
the position of the shock has a strong influence on the the forces and moments,
and thus the aeroelastic response. At very high excitation frequencies it is
possible to show that the amplitude of shock movement is small, albeit with
significant phase lag. At very low excitation frequencies one can expect a quasi-
static transonic response. In between this range, where the majority of cases
lie, shock dynamics are hard to predict. This is because they are determined
by the time history of the acoustic wave pattern, which is even more complex
when the flow contains embedded regions of supersonic flow. Figure 3.10 shows
the time it takes for a signal travelling forward from the trailing edge of an
airfoil to reach a given streamwise position. In the transonic case (right), the
delay is significantly increased by the presence of the shock, as the information
3.2. COMPRESSIBILITY EFFECTS 53
−Cp
x/c
Shock Wave
Sonic Line
M >1 Vortical and entropic wake
Figure 3.9: A typical transonic flow topology (below) and its associated pressure
distribution (above)
must first travel around the supersonic region before reaching the points ahead
of the shock. The time required for this is a strong function of the shock
geometry, which is a function of the airfoil geometry (e.g. compare supercritical
profiles to conventional ones). If the shock is moving, the transient becomes even
more complex. The difficulty of predicting responses in the transonic regime is
further highlighted by figures 11 and 12. which compare linearised theory to
experiment for an oscillating flap as a function of the flow Mach number. The
′
mean measured pressure coefficient, Cp , along with the in-phase Cp and out-of-
′′
phase Cp measured and computed unsteady pressure perturbations are shown.
Clearly the linearized theory does fairly well until transonic effects begin to
appear M = 0.85, after which it becomes increasingly inaccurate.
Figure 3.10: Time delays due to the presence of an airfoil (From Tijdeman, NLR
TR 77090U)
3.3 Viscous effects

When analysing aerodynamic flows it is normally possible to ignore viscous
effects except near the solid boundaries, where the complex dynamics of either
the laminar or turbulent boundary layers ultimately determine the transient
behaviour of the viscous drag. For thin wings, however, viscous drag may be only
a minor contributor to the aeroelastic response. Still the effects of neglecting
viscosity can be large, since they can also significantly affect the unsteady lift,
as discussed in the examples below.
Returning first to the discussion of the Kutta condition in section 3.1, we
note that the dynamics of the trailing separation can influenced by the frequency
of oscillation. At very high frequencies, it is possible to have a reduction in the
effect of the Kutta condition via the appearance of a counter vortex near the
trailing edge (figure 3.13). In this case the change in forces can be substantially
reduced. The exact effect observed, however, will be a function of the trailing
edge geometry and sharpness.
Assuming the Kutta condition is satisfied, the unsteady lift can still affected
by the transient response of the attached boundary layer over the wing. In
3.3. VISCOUS EFFECTS 55
Figure 3.11: Pressure distributions for a NACA 64006 with an oscillating flap.
Steady (upper) and real and imaginary components (lower) comparison of exper-
iment (symbols) with linear theory (lines). (From Tijdeman, NLR TR 77090U)
steady flows, the upper surface boundary layer tends to thicken more rapidly
than lower surface boundary layer as the angle of attack is increased. This
results in a reduction in the effective camber of the section. (figure 3.14). This
effect is further enhanced if a control surface is deflected. In unsteady flows, the
decambering effect becomes a function of the boundary-layer dynamics, which
is itself different from the main flow dynamics due to the relatively low-inertia
near-wall flow.
If separation occurs, the forces and moments are much more directly af-
fected. Conventional separation patterns include the leading and trailing edge
separations shown in the left part of figure 3.15, and the transonic ones shown
at right. The latter occur due to the strong adverse gradients encountered when
passing through the shock wave, and can manifest themselves as either bubbles
or complete separated regions. What makes it interesting is that the shock po-
sition is also a function of the separation pattern. These types of problems pose
difficult challenges to even the most advanced computational techniques.
Under certain circumstances, it is possible to have a self-sustained oscillation
which is driven by cycles of separation and re-attachment. A prominent example
Figure 3.12: Pressure distributions for a NACA 64006 with an oscillating flap.
Steady (upper) and real and imaginary components (lower) comparison of exper-
iment (symbols) with linear theory (lines). (From Tijdeman, NLR TR 77090U)
of this is dynamic stall, where leading edge separation, vortex convection and re-
attachment can provide significant pitching moment excitation. If the structure
is flexible, its motion can in turn drive the dynamic stall cycle. The basic
pattern is illustrated in figures 16 and 17. (Note that the letter indications do
not correspond between the two figures).
Another type of self-sustained oscillation, known as aerodynamic resonance,

can occur in transonic conditions. The basic cycle is shown in figure 18. As the
shock moves aft, it strengthens, leading to boundary-layer separation. The re-
sulting loss in circulation encourages the shock to move forward again, while the
pressure field encourages aft movement of the lower shock. The cycle then alter-
nates on both surfaces. What is remarkable about this self-sustained oscillation
is that it can even happen when the structure is fixed.
3.3. VISCOUS EFFECTS 57
Counter vortex
@@@@@
@@@@@
@@@@@
@@@@@@@
@@@@@
@@@@@@@
@@@@@
@@@@@@@
@@@@@@@
Trailing edge vortex
Figure 3.13: Reduction in the effectiveness of the Kutta condition
δ∗U
δ∗L
Figure 3.14: Viscous de-cambering effect
Shock
Separation bubble
@@@@@
Leading edge separation
@@@@@
Separated wake
@@@@@@@@@@@@@@@@@@@
@@@@@@@@@@@@@@@@@@@
@@@@@@@@@@@@@@@@@@@
Trailing edge separation
@@@@@@@@@@@@@@@@@@@
@@@@@@@@@@@@@@@@@@@
Figure 3.15: Conventional and transonic separation patterns
Figure 3.16: Dynamic Stall (1) From McCorskey, Journal of Fluid Engineering
March 1977, 8-39
Figure 3.17: Dynamic Stall (2)

3.4. MODEL EQUATIONS FOR UNSTEADY FLOWS 59
Figure 3.18: Aerodynamic resonance
3.4 Model Equations for Unsteady Flows

Here we will review the various sets of model equations which can be used
to predict unsteady aerodynamic flows. Starting with the Navier-Stokes equa-
tions, models with varying levels of approximation will be described, and their
advantages and disadvantages discussed. Figure 19 illustrates the hierarchical
relationship between the model equations presented here.
3.4.1 Navier-Stokes Equations

The Navier-Stokes equations may be used to model all effects in homogeneous
flows which may be described as a continuum. They express the conservation
of mass, momentum, and energy for a viscous compressible fluid, and may be
written in integral form for a moving control volume (figure 20) as:
d
Z Z
W dV + (F~ i − F~ v − ẋW ) · ~n dS = 0 (3.2)
dt
V S
ρ~u 0
     
ρ
= =
W =  ρ~u  , F~ i =  ρ~u~u + p I  , F~ v =  τ 
=
ρE ρ~uH τ ·~u − ~q
Here W is known as the state vector, and contains the unknown conservative
variables density (ρ), momentum (ρ~u), and total specific energy (ρE), where ~u
is the local fluid velocity and E = ǫ + 21 |~u|2 , the sum of internal and kinetic
energies. The total enthalpy H, is given by H = E + p/ρ, where p is the local
value of pressure.
The inviscid flux vector, F~ i , in combination with the term ẋW , expresses
the transport of the conserved variables across the control surface. It also in-
cludes the effect of the pressure acting on the control surface on the balance of
momentum within the control volume. The viscous flux vector, F~ v , contains
the effects of viscous stresses. These are described by the viscous stress tensor,
Navier−Stokes Equations
Boundary−Layer Equations Euler Equations
Full Potential Equations
Small−Disturbance Potential
Equation
Linearized Potential Equation
Figure 3.19: Models for unsteady aerodynamic flows
=
τ , which has components given by:

∂ui ∂uj ∂uk
τij = µ + + λδij (3.3)
∂xj ∂xi ∂xk
where µ and λ are the viscosity coefficients. Normally λ is taken to be 2µ/3. The
viscosity coefficients are functions of temperature, and for most circumstances
their variation is determined experimentally.
To complete the system, a thermodynamic relation is required which relates
the pressure to the remaining flow variables. For many aerodynamic flows, an
ideal gas model may be used. In some specialized cases (e.g. re-entry vehicles)
this assumption is inappropriate, and more elaborate thermodynamic relations
must be employed.
Although a few specialized analytic solutions to the Navier-Stokes equations
exist, numerical techniques must usually be used to solve flow problems of en-
gineering interest. The computational cost of numerical simulations, however,
is strongly affected by the enormous range of length scales present at the high
Reynolds numbers typical of aerodynamic flows. This arises from the need to
represent the effects of turbulence, for which relevant length scales can be several
n
V
V : Control volume
S: Surface of control volume
x
n: Local surface normal vector
x: Local surface velocity vector
dS
Figure 3.20: Nomenclature for a moving control volume
orders of magnitude smaller than the global scales of interest for an aeroelastic
problem (airfoil chord, wing span etc.) Various approaches for dealing with this
difficulty are summarized below:
DNS
DNS stands for direct numerical simulation, and refers broadly to computations
which attempt to resolve all the turbulent lengths scales of a flow. These com-
putations must be performed in three spatial dimensions in order to properly
capture the dynamics of turbulence. Although feasible for low-Reynolds-number
flows, the expense of DNS increases rapidly with Reynolds number. For the
foreseeable future, DNS techniques will be impractical for direct application to
aeroelastic problems. They will continue to play an important indirect role,
however, in the the formulation of models of turbulent behaviour, which are
required to exploit the techniques described in the next two sections.
LES
Large Eddy Simulation (LES) techniques make use of a sub-grid scale model
(SGS) to eliminate the smallest turbulence length scales from consideration,
while fully simulating the dynamics of the larger structures. This technique
works well for free turbulent flows, where the dynamics of the smallest scales can
be sufficiently well represented using simplified models. For wall-bounded flows,
however, it can be difficult to model smaller scales accurately, due to the highly
anisotropic nature of turbulent boundary layers. Although the development of
accurate near-wall SGS models is an area of active research, presently DNS-
like resolution is required near the wall in order to obtain accurate results.
Consequently, LES computations remain too expensive to be routinely used for
aeroelastic computations.
F/RANS
The classical approach to dealing with turbulence is based on the concept of
averaging, as introduced by Reynolds in 1895. In summary, one may decompose
the quantities of interest in the flow, for example the pressure, p(x, t), into mean
and fluctuating components:
′
p(x, t) = p̄(x, t) + p (x, t) (3.4)
where the mean value is given by averaging over a time scale, T :
1 t+T
Z
p̄(x, t) = p(x, t)dt (3.5)
T t
where T must be chosen to be large compared to the time scale of turbulent
fluctuations. Substitution of the mean and fluctuating components of each
variable into the Navier-Stokes equations and averaging over time leads to the
Reynolds-averaged Navier-Stokes equations (RANS) for incompressible flows, or
the Favre-averaged Navier-Stokes equations for compressible flows. These equa-
tions contain additional terms arising from the correlation of the fluctuating
solution components. The correlations represent new unknowns, the behavior
of which must be described by an additional turbulence model.
There are numerous approaches to turbulence modelling, each of which has
different ranges of applicability. Although reasonable results can be obtained for
certain classes of flows, the development of turbulence models remains an active
research area, as obtaining accurate results for arbitrary flow conditions can be
difficult. When applied to unsteady flows, an additional problem arises in the
choice of the averaging period, T . In practice, it may be difficult to separate
the time scales of slowest turbulent fluctuations from those of the aeroelastic-
induced flow unsteadiness. In such cases, the use of the time-averaging approach
becomes questionable, and other modelling procedures must be applied.
In spite of the difficulties mentioned above, the application of the Reynolds
or Favre-averaged Navier-Stokes equations to aeroelastic problems can produce
useful results. However, for realistic three-dimentional aircraft geometries sim-
ulation costs remain relatively high, limiting their repeated use for aeroelastic
design procedures.
3.4.2 Boundary-Layer Equations

The boundary-layer equations are an approximation of the Navier-Stokes equa-
tions based on the observation that in high Reynolds-number flows, viscous
effects are usually confined to thin regions adjacent to solid boundaries. When
the radius of curvature of such regions is large compared to their thickness, the
local streamwise flow gradients can be considered to be much smaller than the
gradients normal to the surface. This results in a considerable simplification of
the momentum equation written for the direction normal to the surface:
∂p
=0 (3.6)
∂x1
while the remaining momentum equations collapse to:

∂ρui ∂ ∂pe ∂ ∂ui
+ (ρui uj ) = − + µ (3.7)
∂t ∂xj ∂xi ∂x1 ∂x1
where x1 is the direction normal to the surface, and pe is the pressure of the
flow at the edge of the boundary layer, which is that of the external inviscid
flow. The resulting system, which includes the continuity and energy equations,
is solved only within the thin viscous layer, with boundary conditions given by
the no-slip condition at the wall, and the external flow velocities and pressures
on the boundary layer edge. In many cases, the solution to the boundary layer
equations may be found using marching techniques which exploit the parabolic
character of the equations. This makes such techniques very efficient.
3.4.3 Euler Equations

The Euler equations can be obtained from the Navier-Stokes equations by di-
rectly eliminating the viscous terms (i.e. F~ v = 0 in (3.2)). From the point
of view of aeroelasticity, this limits their application to problems where vis-
cous effects do not significantly affect the forces and moments transferred to
the aircraft structure. In practice, this requires attached flow, and shock Mach
numbers below about 1.4 (in order to avoid significant shock/boundary-layer
interactions).
Alternatively, the influence of viscous effects can be accounted for by solving
the boundary-layer equations near solid walls based on an Euler solution for
the external flow, and in turn correcting the Euler solution using the computed
displacement thickness of the boundary layer. This approach is known as the
viscous-inviscid interaction method, and can be used to treat attached flows, or
those with small regions of separation.
The Euler equations are powerful in that they fully represent all of the in-
viscid effects associated with unsteady flow, including the convection of entropy
and vorticity, and the propagation of information via acoustic waves. In the
context of aeroelastic applications, the fact that vorticity dynamics are auto-
matically included in the solution is a great advantage when considering complex
wakes, or multiple-wake interactions. The solution of the unsteady Euler equa-
tions in three dimensions is still relatively expensive, but is certainly feasible for
exploratory computations.
3.4.4 Full Potential Equations

The full potential equations can be derived from the Euler equations by assum-
ing the flow to be irrotational and isentropic. To do so, we must first convert
the integral form (3.2) with F~ v = 0 to a differential form. Leaving out the
motion of the control volume for now, we can use Gauss’s formula to convert
the surface integral term to a volume integral:
Z Z
F~ i · ~n dS = ∇ · F~ dV (3.8)
S V
As the size of the control volume is arbitrary, the complete volume integral
(including the first term of (3.2)) can only be true if it is true for all points
within the volume, allowing us to write the differential form of the equations as:
∂W
+ ∇ · F~ i = 0 . (3.9)
∂t
Consider for a moment the continuity and momentum equations:
∂ρ
+ ∇ · (ρ~u) = 0 (3.10)
∂t
∂ρ~u
+ ∇ · (ρ~u ⊗ ~u) + ∇p = 0 (3.11)
∂t
By subtracting the continuity equation from the momentum equation, and using
the identity:
1
~u · ∇~u = ∇|~u|2 + (ω × ~u) (3.12)
2
where ω is the local vorticity vector, ∇ × ~u, we can also also express the mo-
mentum equation as:
∂~u 1 1
+ ∇|~u|2 + (ω × ~u) + ∇p = 0 (3.13)
∂t 2 ρ
Recalling the second law of thermodynamics, which can be expressed in terms
of the enthalpy, h = ǫ + p/ρ, as:
1
dh = T ds + dp (3.14)
ρ
we then can express the momentum in terms of the total enthalpy, H, as:
∂~u
+ ∇H = −(ω × ~u) + T ∇s = 0 (3.15)
∂t
where T is the local temperature, and s is the local entropy. Equation (3.15) is
known as Crocco’s relation and is still exact within the framework of the Euler
equations. Now we assume irrotational, isentropic flow, allowing us to replace
the representation of the velocity field with the gradient of a scalar potential
function, ~u = ∇Φ. Equation (3.15) then becomes:
∂∇Φ
+ ∇H = 0 (3.16)
∂t
or
∂Φ
+ H = const ≡ H∞ (3.17)
∂t
We can use this equation to construct an explicit relationship between density
and the potential function by using the isentropic formula for a perfect gas:
1 1
∂Φ |∇Φ|2
γ−1
ρ h γ−1 1
= = H∞ − − (3.18)
ρ∞ h∞ h∞ ∂t 2
√
Noting that h = Cp T and a = γRT , this can be written as:
1
γ−1
ρ γ−1 2 ∂Φ 2
= 1+ U ∞ − 2 − |∇Φ| (3.19)
ρ∞ 2a2∞ ∂t
Finally, we can also relate the potential function to the density via the continuity
equation:
∂ρ
+ ∇ · (ρ∇Φ) = 0 (3.20)
∂t
Equations (3.19) and (3.20) form a coupled non-linear system for only two un-
knowns, even for three-dimensional problems. These are the full potential equa-
tions, and represent a considerable simplification of the original Euler system.
As a result, numerical schemes for the full potential equations are relatively
inexpensive.
The assumptions used in deriving the full potential equations have two con-
sequences. Firstly, the assumption of isentropic flow limits consideration to
weak shocks, i. e. those with Mach numbers less than 1.3. Secondly, assuming
the flow to be irrotational requires that vortical wakes be explicitly modelled,
as opposed to being “captured” as part of the solution. This can be difficult to
do in three dimensions, especially for situations which involve strong wake-wake
interactions.
3.4.5 Small-Disturbance Potential Equation

To derive the small-disturbance potential equations, we consider the perturba-
tions to the full potential equations for the case of a uniform freestream aligned
with the x-axis:
Φ = Φ∞ + φ = U∞ x + φ (3.21)
During the derivation, however, we note that that small terms which accompany
(1 − M 2 ) must be retained, as (1 − M 2 ) is itself small in the transonic regime.
This leads to the following non-linear equation for the disturbance potential:
2
∂2φ ∂2φ ∂2φ ∂ 2φ

2 2 ∂φ ∂ φ 1
1 − M∞ − (γ + 1)M∞ b + 2 + 2 = 2 2U∞ + 2 (3.22)
∂x ∂x2 ∂y ∂z a∞ ∂x∂t ∂t
Equation (3.22) is still capable of representing shocks, although the assump-

tion of small disturbances further limits its accuracy relative to the basic isen-
tropic, irrotational flow restrictions. A significant advantage, however, is that
the boundary conditions can now be linearised and applied on a reference plane
instead of the true wing and wake surfaces. This eliminates the problem of
dealing with moving grids for wings which change their shape.
In order to complete the description of the flow, we need to derive a equation
relating the local value of the pressure to the small-disturbance potential. This
can be done by considering the momentum equation written using (3.21), and
eliminating higher-order terms:

∂φ ∂φ
p − p∞ = −ρ∞ + U∞ (3.23)
∂t ∂x
Note that the equation for pressure is decoupled from (3.22), resulting in further
reductions in computational work. The pressure may be computed as a post-
processing step.
Equation (3.22), also known as the transonic small-disturbance (TSD) equa-
tion, has been widely used in aerodynamic analysis for transonic aeroelasticity,
due to its relatively low computational cost.
3.4.6 Linear Potential Equation

If we assume (1 − M 2 ) is large, we are justified in a full linearisation of the
potential equations for small disturbances. Neglecting all higher-order terms
leads us to:
2
1 ∂2φ ∂2φ

2 ∂ φ
∇2 φ − M∞ − + 2U ∞ =0 (3.24)
∂x2 a2∞ ∂t2 ∂x∂t
Two limits of this equation are often employed:
Steady flow
Removing the time derivatives, we obtain:
2 ∂2φ ∂2φ ∂2φ

(1 − M∞ ) + 2 + 2 =0 (3.25)
∂x2 ∂y ∂z
which is the classical Prandtl-Glauert equation for subsonic and supersonic

steady flows.
Incompressible flow
Taking the limit of a∞ → ∞ gives us Laplace’s equation:
∇2 φ = 0 (3.26)
In this case, the time derivatives are lost from the equation for the perturbation
potential. The problem remains unsteady, however, as the boundary conditions
on the airfoil surface will be changing in time. Equation (3.26) is linear, which
allows us to use superposition techniques in its solution. If the changing wake
shape is not specified, however, determining its position may still require a
non-linear iteration process.
3.5. ANALYTIC SOLUTIONS 67
3.5 Analytic Solutions

Below is a short list of analytic solutions for unsteady flows which are of rele-
vance to aeroelasticity.
3.5.1 Incompressible Flows

• Birnbaum (1924), Theodorsen (1934); linearised incompressible potential
flow for oscillating airfoils
• Wagner (1925); linearised incompressible potential flow for a step change

in angle of attack
3.5.2 Subsonic Compressible Flows

• Possio (1938); Series solution of Integral equation for acceleration potential
for a flat plate.
• Lomax (1953); Piston theory
• Balakrishnan (1998); Lin and Iliff (2000); Closed-form solutions of approx-

imations to the Possio equation
3.5.3 Transonic Flows

• Stahara and Spreiter (1973), Isogai (1974), Dowel (1975); Linearisation
of potential equation for M ≈ 1. Good agreement for thin airfoils in
shock-free flows.
• Tijdeman (1977) Local analysis at shock for frozen-Mach flows. Can use
to investigate effects of shock movement
• Nixon (1979) Linearisation of small-disturbance equation using indicial

function approach.
3.5.4 Supersonic Flows

• von Borbely (1942); Garrick and Rubinow (1946) Flat plates in Supersonic
flow - integral equation solution.
In the chapters following this one we will examine selected analytic solutions
for the simplest of the model equations. These have been chosen in order to
consider the most important of the unsteady effects separately. Each solution
can be used to provide reasonable predictions under a restricted set of conditions,
but their chief advantage is that they help to develop insight into the relevant
physical phenomena. The discussion of prediction methods for combined or
more complex phenomena is deferred to later in the notes, where we will examine
numerical solution methods for unsteady flows.
The analytic solutions we will consider are representative of some of the

major results from unsteady aerodynamic theory over the past 80 years. In
most cases, the full details of the solutions are too lengthy to be presented in
their entirety, however, so we will omit some of the complex details where they
do not interfere with understanding. For the exam, you will only be responsible
for what is presented in class:
• Basic assumptions
• Governing equations
• Geometric simplifications
• Boundary conditions
• Outline of the solution procedure
• Interpretation of results
In addition, you should be able to solve simple problems employing the

information above (see the practice problems in the chapaters following this
one).
p3.6.5 Compare the advantages and disadvantages of applying the Euler equa-
tions versus the linear potential equation to the computation of whirl flutter in
the propeller-wing-nacelle configuration shown below. Assume the configuration
is operating a with a freestream Mach number of 0.6.
M = 0.6
p3.6.6 Consider a 2D thin airfoil moving in a uniform flow in the x direction,

U∞ , with a known perturbation potential, φ:
φ
z
Using the momentum equation:
∂~u ∇p
+ (~u · ∇)~u + =0
∂t ρ
derive the linearised expression for pressure:

∂φ ∂φ
p = p ∞ − ρ∞ + U∞
∂t ∂x
Start by considering the x and z momentum equations for ~u = (U∞ , 0)+(u′ , w′ ),

ρ = ρ∞ +ρ′ , p = p∞ +p′ and eliminating the higher-order terms (Note that there
are no spatial gradients in the mean-flow quantities). Then write the remaining
equations as a vector equation in terms of ∇φ and ∇p′ and integrate it along a
streamline.
Chapter 4
Analytic Solutions for

Incompressible Flows
In this chapter we will examine potential-flow solutions for unsteady incom-

pressible flows around thin airfoils. We will begin by considering a closed-form
solution procedure for the problem of a harmonically oscillating airfoil, then
review some solutions for step changes in angle of attack and step gusts.
The first solutions for the forces on a moving airfoil were produced in the the
1920’s by aerodynamicists such as W. Birnbaum, H. Wagner, and H. Glauert.
Their analyses established a clear strategy for solving unsteady aerodynamic
problems, although their solution methods sometimes relied on series approx-
imation techniques which did not always converge. In 1934, T. Theodorsen
was able to express the solution to the harmonically oscillating airfoil problem
explicitly in terms of Hankel functions. He also developed some of the first
flutter analysis techniques incorporating the effects of unsteady aerodynamics.
His solution to the oscillating airfoil problem has been widely used in the field
of aeroelasticity.
Figure 4.1: An unsteady airfoil and its wake
In the following sections we will consider solutions based on a vortex-sheet

model for an moving airfoil and its wake. Note that due to the influence of
convecting shed vorticity, we cannot omit the wake in an unsteady analysis.
This contrasts with the 2D steady airfoil case, where the wake geometry is
irrelevant so long as the Kutta condition is satisfied.
71
72CHAPTER 4. ANALYTIC SOLUTIONS FOR INCOMPRESSIBLE FLOWS
Zt(x) Zc(x) Za(x,t)

Z(x,t)
t
= + +
x
Figure 4.2: Decomposition of an unsteady airfoil problem
4.1 Problem definition

We will consider the linearized equations for the potential of an unsteady flow,
which can be expressed as (see chapter 3):
∇2 φ = 0 (4.1)

∂φ ∂φ
p − p∞ = −ρ∞ + U∞ (4.2)
∂t ∂x
Where φ is the perturbation potential and the ∞ denotes freestream conditions.
To be consistent with the linearization of the potential, we will consider only
small angle excursions, and a linearized airfoil geometry.
If we linearize the boundary conditions on the airfoil, the solution of an un-
steady thin-airfoil problem can be broken up into its thickness (zt ), camber(zc ),
and flat-plate (za ) components, as shown in figure 4.2. As thickness and camber
are not usually considered to be functions of time, their contributions can be
obtained from standard steady thin-airfoil theory. Therefore, for the remainder
of this analysis, we will focus on the solution for a moving flat plate. In doing
so we will solve (4.1) and (4.2) subject to the following conditions:
1. Zero perturbation velocity in the far-field
2. Flow tangency on the airfoil surface
3. No pressure jump across the wake
4. Kutta condition using Kelvin’s Theorem
The first of these conditions is simple enough, we require that the perturbation
vanish when considering points far from the airfoil. The second is the standard
no-flow-through-the-airfoil condition familiar from steady aerodynamic theory.
The third is a new condition which arises from having to consider the wake in
the problem, which of course cannot maintain a difference in pressure across
it. The final condition represents an interpretation of the Kutta condition for
unsteady flow, that the vorticity shed at the trailing edge must be equal and
opposite in strength to the change in bound vorticity of the airfoil (see the
previous chapter for a discussion).
In order to solve the problem we will represent the airfoil and its wake by
a sheet of vortices of unknown strength, γ(x, t), lying on the z = 0 plane, as
4.2. FAR-FIELD PERTURBATION 73
x
b
−b
γa(x,t) γw(x,t)
Figure 4.3: Replacement of airfoil and wake by a vortex sheet
shown in figure 4.3. Each of these vorticies is a elementary solution of equation

(4.1). This results in the following expression for the perturbation potential
induced by the complete vortex sheet.
Z ∞
γ(ξ) z
φ(x, z) = − tan−1 dξ (4.3)
−b 2π x−ξ
The use of a vortex sheet exploits the linear nature of the governing equations
to simplify the problem. As each of the infinitesimal point vortices in the vortex
sheet satisfies the governing equations, then the entire vortex sheet satisfies the
governing equations, even for arbitrary local strengths. To find the solution,
it thus remains only to find the distribution of vortex strengths which satisfies
the boundary conditions 1-4 above. We will derive mathematical expressions of
these conditions in the following sections.
4.2 Far-Field Perturbation

To examine if our vortex sheet representation does indeed have a zero pertur-
bation far from the airfoil, we start by deriving expressions for the perturbation
velocities anywhere in the domain by differentiating the potential function for
the vortex sheet:
Z ∞
γ(ξ) z
u(x, z) = φx (x, z) = dξ (4.4)
−b 2π (x − ξ)2 + z 2
∞
γ(ξ) x−ξ
Z
w(x, z) = φz (x, z) = − dξ (4.5)
−b 2π (x − ξ)2 + z 2
Both of these components tend to zero in the far-field, as one would expect
considering the behaviour of a single potential-flow vortex.
4.3 Velocities on the z = 0 plane

Before proceeding to the flow tangency and zero-wake pressure jump conditions,
it is worthwhile first deriving the expressions for the local velocity at the surface
Figure 4.4: Integration in the region near x
of the vortex sheet. Taking the limit of the expressions for velocity (4.4) and
(4.5) as z → 0, we have:
Z ∞
γ(ξ) z
u(x, 0) = lim dξ (4.6)
z→0 −b 2π (x − ξ)2 + z 2
and
∞
1 γ(ξ)
Z
w(x, 0) = − dξ (4.7)
2π −b (x − ξ)
Note that the integral for u(x, 0) tends to vanish except where x ≈ ξ.
In order to find an expression for u(x, 0), we restrict the range of integration
to small region near x of width 2ǫ (figure 4.4). Introducing the intermediate
variable ξ ′ = ξ − x:
Z +ǫ
γ(ξ ′ )

z
u(x, 0) = Limz→0 ′ )2 + z 2
dξ ′ (4.8)
−ǫ 2π (ξ
A Taylor series approximation for γ(ξ ′ ) is:
γ(x − ǫ → x + ǫ) = γ(x) + O(ǫ) (4.9)
Integrating and eliminating higher-order terms:

γ(x) + O(ǫ) −1 ǫ

−1 −ǫ
u(x, 0) = Limǫ→0 Limz→0 tan − tan (4.10)
2π z z
We assume ǫ is small, but is still large compared with z → 0. This leads to the
following expression for the velocity above (+) and below (−) the vortex sheet:
γ(x)
u(x, ±0) = ± (4.11)
2
Alternatively we can express the jump across the vortex sheet directly as:
∆u(x, 0) = γ(x) (4.12)
4.4. FLOW TANGENCY 75
4.4 Flow Tangency

Now that we have derived expressions for the velocity at z = 0, we can employ
them directly in the application of the flow tangency condition. This requires
~ · ~n = Vn = 0, where
that at the surface of the airfoil V
Vn = Vnf reestream + Vnmotion + Vnγ (4.13)
Here Vnf reestream is the sheet-normal component of the freestream vector due
to the instantaneous pitch angle of the airfoil, Vnmotion is the normal velocity
induced by the vertical motion of the airfoil, and Vnγ is the local value of induced
normal velocity due to the vortex sheet. Recalling that we are using small-angle
approximations, we may replace Vnf reestream and Vnmotion with space and time
derivatives of the local z coordinate of the sheet, as shown in figure 4.5. If
n dz(x,t)
dt
n
U
|Vn| − dz(x,t)
motion
|Vn| − U dz(x,t) dt
freestream
dx
Figure 4.5: Geometric linearizations for the flow-tangency condition
we also note that for small angles w ≈ Vn , the total vertical velocity due to
freestream and motion, wF M , may be expressed:
dza (x, t) dza (x, t)
wF M (x, t) ≈ −U∞ − (4.14)
dx dt
This allows us to express the flow tangency condition as:
wF M (x, t) + wγ (x, t) = 0 (4.15)
or, by substituting our previous expressions, as:
Z b Z ∞
dz(x, t) dz(x, t) 1 γa (ξ, t) 1 γw (ξ, t)
−U∞ − − dξ − dξ = 0 (4.16)
dx dt 2π −b x − ξ 2π b x−ξ
where za (x, t) is the specified motion of the airfoil. Note that we have separated
the integration of the velocity induced by the vortex sheet into components
associated with the airfoil and wake.
4.5 Pressure Jump Across the Wake

In general, the wake is a free surface which will move to ensure that a pressure
jump does not occur across it. In our model, we will assume that the excursions
of the wake in the vertical direction are so small that they can be neglected while
enforcing the zero-pressure-jump condition. This assumption is quite reasonable
when considering the relatively high forward velocities characteristic of most
aircraft wings.
We can construct an expression for the pressure jump across the wake using
(4.2):

∂∆φ ∂∆φ
pu − pl = −ρ∞ + U∞ (4.17)
∂t ∂x
(4.18)
where ∆φ = φu − φl . This equation can be re-expressed in terms of γ if we note
that:
∂∆φ(x)
= ∆u(x) = γ(x) (4.19)
∂x
Also, for x < −b, ∆φ = 0, so that for x > −b:
Z x Z x
∂∆φ(ξ)
∆φ(x) = dξ = γ(ξ)dξ (4.20)
−b ∂x −b
This leads to:

Z x
∂
pu − pl = −ρ∞ U∞ γ + γ(ξ)dξ (4.21)
∂t −b
The wake pressure condition then requires pu − pl = 0 for x = b → ∞
4.6 Kutta Condition using Kelvin’s Theorem

For our analysis, we assume that the Kutta condition is satisfied at each instant
in time. As we saw in chapter 3, this implies that changes in the airfoil’s
condition will result in the shedding of vorticity from the trailing edge. In order
to comply with Kelvin’s theorem, we require that the shed vorticity is equal
and opposite to the change in bound circulation over the airfoil. To express this
mathematically, we start with an equation for the airfoil’s bound circulation:
Z b
Γ(t) = γa (ξ, t)dξ (4.22)
−b
Consider an increment in vorticity shed from the trailing edge, as shown in

figure 4.6. In terms of Γ(t), Kelvin’s theorem can be now expressed as:
dΓ(t)
γw (b, t)dx = − dt (4.23)
dt
If we consider infinitesimal time and space intervals, we can get the final ex-
pression for the magnitude of the vorticity shed instantaneously at the trailing
edge:
1 dΓ(t)
γw (b, t) = − (4.24)
U∞ dt
4.7. HARMONIC SOLUTIONS 77
γw(b,t)
CCCCCC
CCCCCC
CCCCCC
CCCCCC
dx = U dt
Figure 4.6:
4.7 Harmonic Solutions

We have already seen that the far-field condition is automatically satisfied
though our choice of potential-flow singularities used to represent the airfoil
and wake. For the case of harmonic motion, however, we can also show that if
we assume that wake vorticity is convected with the freestream velocity, U∞ ,
then the wake pressure condition is also automatically satisfied.
Let’s begin by assuming a harmonic forms for the airfoil excitation and
resulting solution:
wF M (x, t) = ŵF M (x)eiωt ; γa (x, t) = γ̂a (x)eiωt ; Γ(t) = Γ̂eiωt (4.25)
We also assume a harmonic response in the wake, so that:
γw (x, t) = γ̂w (x)eiωt (4.26)
allowing Kelvin’s theorem to be expressed as:
iω
γ̂w (b) = − Γ̂ (4.27)
U∞
Assuming downstream convection of vorticity with U∞ , the complex coefficient

of the wake circulation distribution can then be expressed:
iω iω
γ̂w (x) = − Γ̂e U∞ (b−x) (4.28)
U∞
Now we substitute these expressions into the wake pressure condition pu −pl = 0.
The second term in (4.21) can be expressed as:
x Z x
∂ dΓ ∂
Z
γ(ξ)dξ = + γw (ξ)dξ (4.29)
∂t −b dt ∂t b
Z x
= iω Γ̂ + iω γ̂w (ξ)dξ eiωt (4.30)
b
G(k)
0.5 1.0
0.0 F(k)
k= 1.0 k=0
0.6
0.05
−0.2 0.4 0.1
0.2
Figure 4.7: Components of Theodorsen’s function C(k) = F (k) + iG(k)
Using the expression for γ̂w (ξ) given by Kelvin’s theorem:

x x
iω Γ̂
Z Z
iω
γ̂w (ξ)dξ = − e U∞ (b−ξ) dξ (4.31)
b U∞ b
Integration and substitution into the expression for the pressure change across
the wake shows that the wake pressure condition is automatically satisfied.
Now we are left with just the expression for flow tangency (4.16), augmented
by an explicit expression for the wake vorticity distribution (4.28) for harmonic
flows. (4.16) evaluated for harmonic solutions can be combined with (4.28) to
obtain the following expression for flow tangency on a harmonically oscillating
airfoil:
b−ξ
1 b
γ̂a (ξ) iω eiω( U∞ )
∞
Z Z
− dξ + Γ̂ dξ (4.32)
2π −b x−ξ 2πU∞ b x−ξ
∂ ẑa (x)
= iω ẑa (x) + U∞
∂x
where za (x, t) = ẑa (x)eiωt . This integral equation has been solved analytically
for γ̂a (x) Theodorsen (1934) and Schwarz (1940), with ẑa (x) defined for a har-
monically plunging and pitching flat plate. Solutions of the integral equation
are usually expressed in terms of Theodorsen’s function, C(k):
(2)
H1 (k)
C(k) = F (k) + iG(k) = (2) (2)
(4.33)
H1 (k) + iH0 (k)
(2)
where Hn are nth-order Hankel functions of the second kind. Theodorsen’s
function is also known as the circulation function. The parameter k ≡ (ωb)/U∞ ,
is referred to as the reduced frequency. It expresses the ratio of periods for
oscillation and (wake) convection. The variation of the components of C(k) as
a function of k is shown diagrammatically in figure 4.7.
4.8. HARMONIC PITCH AND PLUNGE SOLUTIONS 79
−b b
x
ab
θ h
Figure 4.8: Flat plate undergoing combined pitch and plunge
4.8 Harmonic Pitch and Plunge solutions

Although their detailed derivation is beyond the scope of this course, it is useful
to examine solutions to (32) for some basic cases. We first consider the solution
for a flat plate undergoing combined harmonic pitching and plunging with a
reduced frequency k, as shown in figure 4.8. For a center of pitch positioned at
ab, the solution for the airfoil’s lift and nose-up moment about the rotation for
combined harmonic pitch can be written as:
L = L̂eiωt = πρ∞ b2 [ḧ + U∞ θ̇ − baθ̈] (4.34)

1
+ 2πρ∞ U∞ bC(k)[ḣ + U∞ θ + b( − a)θ̇]
2
and :
1 1
M = M̂ eiωt = πρ∞ b2 [baḧ − U∞ b( − a)θ̇ − b2 ( + a2 )θ̈] (4.35)
2 8
1 1
+ 2πρ∞ U∞ b2 (a + )C(k)[ḣ + U∞ θ + b( − a)θ̇]
2 2
The first group of terms in both the lift and moment equations is known as the
non-circulatory component of the response, while the second group of terms is
known as the circulatory component, as it arises from the satisfaction of the
Kutta condition. The non-circulant components are sometimes referred to as
apparent mass terms, as they are proportional to accelerations felt by a fluid
cylinder with radius b. Consider also the low k limit of the lift and moment
expressions. At k = 0, all the rates dissappear, C(k) = 1 (see figure 4.7), and
the standard flat-plate solution is recovered (note b = c/2). We will often use
an alternative notation for the lift and moment expressions, in terms of their
unsteady (complex) coefficients:
" #
2 ĥ
L̂ = πρ∞ U∞ b Lh (k) + Lθ (k)θ̂ (4.36)
b
" #
2 2 ĥ
M̂ = πρ∞ U∞ b Mh (k) + Mθ (k)θ̂ (4.37)
b
3 2
Lh’ Mh’
Lh’’ Mh’’
2
1.5
1
0
1
-1
-2 0.5
-3
0
-4
0 0.5 1 1.5 2 0 0.5 1 1.5 2
k k
Figure 4.9: Unsteady lift and moment coefficients for a flat plate airfoil under-
going pure harmonic plunge
4 1.5
Lθ’ Mθ’
Lθ’’ Mθ’’
1
3
0.5
2 0
1 -0.5
-1
0
-1.5
-1 -2
0 0.5 1 1.5 2 0 0.5 1 1.5 2
k k
Figure 4.10: Unsteady lift and moment coefficients for a flat plate airfoil under-
going pure harmonic pitch
with M̂ measured nose-up positive relative to x = ab, and
Lh = −k 2 + C(k) [2ik] (4.38)

Lθ = ak 2 + ik + C(k) [2 + ik(1 − 2a)] (4.39)
Mh = −ak 2 + C(k) [ik(1 + 2a)] (4.40)

1 1 1
Mθ = (a − )ik + (a2 + )k 2 + C(k) (2a + 1) + ik( − 2a2 ) (4.41)
2 8 2
The variation of the total lift and moment coefficients with k is shown in figures
4.9 and 4.10. In figure 4.10, the lift curve slope is seen to drop from its standard
value of 2 fairly quickly with increasing k.
It is also possible to find expressions for the chordwise perturbations to
the steady pressure distributions of the airfoil. For harmonic flows, these are
4.8. HARMONIC PITCH AND PLUNGE SOLUTIONS 81
6 10
∆Cp’ ∆Cp’’
4 8
h/b h/b
2 k=0.2 6
0 1 x
0 b 4 k=1.0
−2 0.6
0.6 2 0.2
−4 1.0 x
−1 0 1 b
Figure 4.11: Pressure perturbations for a flat plate airfoil undergoing pure har-
monic plunge
k=1.0
8 8
∆Cp’ ∆Cp’’
6 6
θ θ 0.6
4 4
2 0.6 0.2 2 0.2

k=1.0
x 0 x
0 b
0 0 1b
−2 −2
Figure 4.12: Pressure perturbations for a flat plate airfoil undergoing pure har-
monic pitch
expressed as:
∆Cp (x, t) = ∆Cp′ (x) + i∆Cp′′ (x) eiωt

(4.42)
The real and imaginary perturbations to the change in pressure across the plate,
∆Cp′ and ∆Cp′′ are shown plotted as a function of chordwise position in figures
4.11 and 4.12. For the plunge case, real component ∆Cp′ corresponds to
the perturbation obtained when the plate is at its most downward position,
while ∆Cp′′ corresponds to the mid-position while moving downwards. Each
distribution shows the characteristic leading-edge singularity of flat plate steady
solutions. Examination of the pure pitch perturbations reveals that they only
resemble the steady flat plate pressure distribution at low k.
t* < 0 t* > 0
1.0 φ(t*)
0.8
0.6
0.4
0.2
t*
0 4 8 12 16 20
Figure 4.13: Wagner’s indicial lift function
4.9 Indicial Response (Wagner’s Function)

Although harmonic responses are clearly useful for analysing flutter instability
boundaries, indicial (step) responses are often more convenient when consider-
ing arbitrary excitations. The solution for the indicial plunge response of a flat
plate airfoil was first determined by Wagner in 1925. It can also be derived di-
rectly from the previously-presented harmonic plunge results using the following
excitation:
w 34 c = 0 t<0 (4.43)
−U∞ αo t>0 (4.44)
The harmonic solution’s non-circulatory component will provide a singular pulse

at t=0, but otherwise it will have no contribution. The remainder of the response
can be obtained by using a Fourier transform of the circulatory terms, resulting
in:
2
L = 2πρ∞ U∞ bαo φ(t∗ ) (4.45)
where
∞
1 C(k) ikt∗ U∞ t
Z
∗
φ(t ) = e dk ; t∗ = (4.46)
2π −∞ ik b
In incompressible flow, the lift perturbation due to an indicial plunge acts at the
quarter chord, so there is no perturbation to the moment about c/4. Wagner’s
indicial lift function can be approximated using:
∗ ∗
φ(t∗ ) ≈ 1 − 0.165e−0.0455t − 0.355e−0.3t (4.47)
4.10. VERTICAL STEP GUST RESPONSE (KÜSSNER’S FUNCTION) 83
4.10 Vertical Step Gust Response (Küssner’s Func-

tion)
In an analogous manner, the response to a step gust is also very useful. Note
that unlike the previous indicial plunge case, in this case the local angle of attack
is not constant along the airfoil chord until the gust has past. This excitation
can be expressed as (see figure 4.14).
wG = 0 x > U∞ t − b (4.48)
wo x < U∞ t − b (4.49)
The response to this excitation can be expressed in terms of Küssner’s Function,
wG
Figure 4.14: Step gust excitation
ψ(t∗ ):
L = 2πρ∞ U∞ bwo ψ(t∗ ) (4.50)
As for the indicial plunge case, ψ(t∗ ) can be derived via a Fourier transform of
the solution to a sinusoidal gust problem. It has the form shown in figure 4.15.
Once again in incompressible flow the lift response acts at the quarter chord, so
there is no perturbation to the moment about c/4. Küssner’s step gust response
function can be approximated by:
∗ ∗
ψ(t∗ ) ≈ 1 − 0.5e−0.13t − 0.5e−t (4.51)
Note that the non-dimensional time, t∗ , can also be interpreted as the number
of half-chords travelled after the leading edge encounters the gust.
1.0 φ(t*)
0.8
0.6 ψ(t*)
0.4
0.2
t*
0 4 8 12 16 20
Figure 4.15: Küssner’s step gust response function
p4.11.7 Sketch the pressure distributions for a harmonically plunging airfoil

at 8 points in the cycle.
p4.11.8 Plot the quasi-steady aerodynamic model response in figure 4.9.
p4.11.9 Show that the wake pressure condition is automatically satisfied in

the harmonic motion case using the definitions from sections 4.5 to 4.7.
p4.11.10 Consider a thin airfoil in an incompressible flow, as shown below:

γa(x,t) γw(x,t)
x
−b b
x=b
e
γ (b,t)
CCCCCC
w
CCCCCC
CCCCCC
CCCCCC P
dx = U dt
If the airfoil deforms harmonically, so that the time-variant vorticity distribution

on the airfoil surface is given by:
γa (x, t) = (x − b)2 (1 + 2i)eiωt ; −b ≤ x ≤ b
What is the vorticity as as a function of time at the point P located a distance

e behind the trailing edge [ i.e. γw (b + e, t) ]? Use assumptions consistent with
the theory presented in class.
Chapter 5

Subsonic Compressible
Flows
5.1 Steady Flow: The Prandtl-Glauert Correc-

tion
The simplest equation describing compressible flows is the steady linear poten-
tial equation (see chapter 3), written below for two dimensions:
2 ∂2φ ∂2φ
(1 − M∞ ) + 2 =0 (5.1)
∂x2 ∂y
Comparing (5.1) to the linear potential equation for incompressible flow reveals
that solutions to the former may be obtained from those of the latter if a simple
coordinate transformation is applied. Specifically, The solution of (5.1) can be
related to that of Laplace’s equation for φ̄(ξ,
pη) = βφ(x, y) in the transformed
coordinates ξ = x , η = βy; where β = 1 − M∞ 2 . Using this concept, it
is possible to obtain forces and moments for a steady compressible flow by

applying corrections to the forces and moments for the corresponding steady
incompressible flow. For thin airfoils, the linearised, steady pressure coefficient
can be shown to be (see an introductory aerodynamics text, such as Anderson):
2u 2 ∂ φ̄ Cpo
Cp = − =− = (5.2)
U∞ U∞ β ∂x β
where Cpo is the pressure coefficient from the solution to Laplace’s equation
(i.e. the incompressible problem). The behaviour of the corrected lift curve
slope over the complete Mach number range is illustrated in figure 5.1. The
correction was first derived by Prandtl in in the early 1920’s and more formally
87
88CHAPTER 5. ANALYTIC SOLUTIONS FOR SUBSONIC COMPRESSIBLE FLOWS
2π
Clα 1−Μ2
4
2π Μ2−1
0 1.0 2.0 3.0

M
Figure 5.1: Prandtl-Glauert correction to the lift-curve slope of a flat-plate

airfoil
by Glauert in 1928. It is important to note that since the linear potential

equation for compressible flow (5.1) is obtained by ignoring time derivatives,
the Prandtl-Glauert correction applies to steady flows only. Therefore, it can
only be reasonably applied to aeroelastic problems if the frequencies of structural
motion are small enough to allow the use of a steady aerodynamic model.
5.2 Unsteady Harmonic Solutions for M < 1

In 1938, Possio rewrote the unsteady linearized potential equation for compress-
ible flows (equation 24 in chapter 3) in terms of the acceleration potential, and
formulated the problem for the unsteady thin airfoil. This approach is more
general than that described in chapter 4 as it includes the effects of time delays
in addition to those of vorticity convection. Unfortunately, it has not been pos-
sible to derive closed form solutions to Possio’s equation, so traditionally it has
been solved using series-approximation techniques or other numerical methods.
Recently, however, Balakrishnan (1998) Lin and Iliff (2000) have derived closed-
form solutions to approximations to Poisso’s equation. The Lin and Iliff (2000)
approximation is given below (moment positive clockwise):
k4

Lh = −k 2 + C(k) [2ik] + M 2 logM − 2ik 3 C(k) − 2k 2 C(k)2 (5.3)
2
Lθ = ak 2 + ik + C(k) [2 + ik(1 − 2a)]

ik 3 k4

2 2 1 3
+M logM − − a + C(k) −2k − − 2a ik
2 2 2
+M 2 logM C(k)2 2ik − (1 − 2a)k 2

(5.4)
5.3. PISTON THEORY 89
Mh = −ak 2 + C(k) [ik(1 + 2a)]

k4

2 2 2 1 3
+M logM −(1 + 2a)k C(k) − + 2a ik C(k) + a (5.5)
2 2

1 1 1
Mθ = (a − )ik + (a2 + )k 2 + C(k) (2a + 1) + ik( − 2a2 )
2 8 2
3 2 4

aik a k 1
+M 2 logM − − + C(k) (1 + 2a)ik − ( − 2a2 )k 2
2 2 2

1
−M 2 logM C(k) + 2a k 2 − 2a2 ik 3 (5.6)
2
Some results this approximation at a Mach number or 0.7 are shown compared
to incompressible results in figures 5.2 and 5.3. The compressible effects can be
seen to increase as a function of k. You can further investigate the influence of
Mach number and the centre of rotation on your own by visiting:
https://aerodynamics.lr.tudelft.nl/cgi-bin/lpfp
5.3 Piston theory

If there are time delays, then information transfer is small for very fast motions.
This concept allows a particularly simple aerodynamic theory to be derived.
Consider the finite wing shown in figure 14.6. When the wing is impulsively
plunged, then the interior of the wing initially acts as if the span and chord were
infinite, since waves coming from the edges require time to propagate inwards.
As a result, a simple 1D compression wave is generated on the lower surface of
the wing, and a simple expansion wave is generated on the upper surface. From
linearised acoustic theory, the initial pressure perturbation on the surface of a
piston which is impulsively accelerated is:
p′ = ρ∞ a∞ wo (5.7)
where wo is the speed to which the piston is accelerated to, and a∞ is the
speed of sound in the fluid. This perturbation appears anti-symmetrically on
the upper and lower surfaces so that ∆p = 2ρ∞ a∞ wo . The initial lift coefficient
of the wing is then:
(p′lower − p′upper )S 4wo

CL (0) = = (5.8)
1/2ρ∞ U 2 S M∞ U
These expressions provide a useful model for oscillations with very high frequen-
cies at both subsonic and supersonic Mach numbers. Formally, the theory is
accurate at large values of the product of Mach number and reduced frequency,
M k.
0.5 3
0
2.5
-0.5
2
-1
Lh’’
Lh’
-1.5 1.5
-2
1
-2.5
0.5
-3
-3.5 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
k k
2 0
1.8
1.6 -0.05
1.4
1.2 -0.1
Mh’’
Mh’
1
0.8 -0.15
0.6
0.4 -0.2
0.2
0 -0.25
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
k k
Figure 5.2: Plunging case - comparison of M=0.7 (solid) and M=0 (dashed)
responses (a=-1/2).
5.3. PISTON THEORY 91
3 5
4.5
2.5 4
3.5
2 3
2.5
Lθ’’
Lθ’
1.5
2
1 1.5
1
0.5 0.5
0
0 -0.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
k k
1.2 0
1 -0.5
0.8 -1
Mθ’’
Mθ’
0.6 -1.5
0.4 -2
0.2 -2.5
0 -3
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
k k
Figure 5.3: Pitching case - comparison of M=0.7 (solid) and M=0 (dashed)
responses (a=-1/2).
@@@@@
@@@@@
@@@@@
@@@@@
Expansion Wave
Compression Wave
wo
Figure 5.4: Piston theory
5.4 Indicial Responses for M < 1

It is also possible to compute indicial responses using unsteady linearized po-
tential theory. These solutions are quite different from their incompressible
counterparts, in that:
• The non-circulatory components are active after t=0, due to the finite
times required for wave propagation.
• The solution is no longer singular at t=0; its initial peak can be predicted
using “piston theory”.
• The circulatory response is no longer determined by w 43 c
• The change in moment about the quarter chord is no longer zero.
Example responses for an indicial plunge at different Mach numbers are shown
in figures 5.5 and 5.6, corresponding following definitions of lift and moment:
2 ḣo
L = πρU∞ Φc (t∗ ) (5.9)
U∞
2 ḣo
M = πρU∞ 2b ΦcM (t∗ ) (5.10)
U∞
After the initial piston theory response, there is an increasingly slow climb
to the steady state value as the Mach number is increased. This is due to
the increasing upstream propagation times of information describing the wake’s
dynamics. For cases where linearised theory is adequate, the final large-time
asymptotic value can be estimated using the Prandtl-Glauert factor √ 1 2 .
1−M∞
5.4. INDICIAL RESPONSES FOR M < 1 93
M = 0.7
1.2
φc(t*) M = 0.5
1.0
M=0
0.8
0.6
0.4 t*
0 4 8 12 16 20 24 28
Figure 5.5: Indicial lift responses in compressible flow
(φcM)c/4(t*)
−0.3
M = 0.5
−0.2
−0.1
M = 0.7
M=0
0
t*
0 4 8 12
Figure 5.6: Indicial moment responses in compressible flow

p5.5.11 Which of the real or imaginary components of Lh , Lt , Mh or Mt has

the largest relative deviation from the incompressible unsteady flat-plate model
corrected by the Prandtl-Glauert factor?
p5.5.12 Using piston theory and the Prandtl-Glauert correction, sketch the
M=0.5 indicial plunge function for a finite wing with an incompressible lift-curve
slope of 4. Give values for t = 0 and t = ∞.
Chapter 6

Supersonic Flows
In this chapter we will consider solutions to the linear potential equation for
supersonic flows. We will start by considering an isolated airfoil (figure 6.1).
This case is quite different from the solutions we have considered up to now, as
convecting vorticity in the wake has no influence on the airfoil’s forces due to
the former’s restricted domain of influence. This allows us to concentrate on the
analysis of time-delay effects. Note that in more general supersonic problems,
such as those with highly swept or multiple lifting surfaces, convecting vorticity
must again be taken into account. This will be discussed in more detail at the
end of the chapter.
Although the unsteady supersonic airfoil problem was first solved by von
Borbely in 1942, we will consider the later solution by Garrick and Rubinow
(1946), because of its more physical interpretation. As for the incompressible
case, we will exploit the linearity of the governing equations and solve the prob-
lem using a distribution of potential flow singularities. Once again, since each
of the singularities is a solution to the governing equations, to solve the prob-
lem we only need to find the combination of their strengths which satisfies the
boundary conditions.
The boundary conditions we need to consider are simpler than those of the
incompressible case. Firstly, as the wake has no influence on the isolated airfoil,
we may ignore Kelvin’s theorem and the zero-pressure-change-across-the-wake
condition. Secondly, although in the far-field the disturbance is finite and a
function of the solution, we need not take it explicitly into account as it will be
automatically provided once the singularity strengths are determined (this was
also true for the incompressible case). This leaves us with only the flow-tangency
condition to enforce.
Note that although the perturbation in the far-field is non-zero, for an earth-
fixed observer, it is only non-zero for a finite time. This can be appreciated
by imagining the airfoil in figure 6.1 passing over you. Initially you will feel
95
96 CHAPTER 6. ANALYTIC SOLUTIONS FOR SUPERSONIC FLOWS
µ
x
−b b
Figure 6.1: Thin airfoil in supersonic flow
the compression after you are passed by the wave emanating from the leading
edge, and then expansion back to ambient pressure after you are passed by the
wave emanating from the trailing edge. If you are sufficiently far away, the
perturbation you feel due to the wake will be negligible.
6.1 Problem definition

As given in chapter 3, the linearised potential equations for unsteady compress-
ible flows are:
∂2φ ∂2φ 1 ∂2φ ∂2φ

2
(1 − M ) + = + 2U (6.1)
∂x2 ∂z 2 a2 ∂t2 ∂t∂x

∂φ ∂φ
p − p∞ = −ρ∞ + U∞ (6.2)
∂t ∂x
As for the incompressible case, the airfoil will be located in the z = 0 plane, so
that the linearized flow tangency condition can be expressed:

∂φ
= wF M (x, t) (6.3)
∂z z=0
Our intent is to represent the airfoil using a distribution of solutions to (6.1).

A particularly useful one for this purpose is the source-pulse solution:
A(ξ, ζ, T )
φs (x, z, t) = p (6.4)
a (t − T ) − [(x − ξ) − U (t − T )]2 − (z − ζ)2
2 2
φs (x, z, t) is the potential induced by a pulse of strength A at (ξ, ζ) and time

T , measured at a point located at (x, z) at time t. Its effect can be visualised
by considering the expanding sphere of the pulse disturbance front in the airfoil
frame of reference, as shown in figure 6.2.
At time T the pulse is made, but due to the finite time of information prop-
agation, it takes time for the surrounding regions to become aware of it. We use
6.1. PROBLEM DEFINITION 97
Disturbed
Region
(x,z)
a(t−T)
(ξ,ζ)
U(t−T)
(x−ξ)
Figure 6.2: Expansion of a travelling pulse
the term “pulse disturbance front” to refer to the suface which separates points
which are aware of the pulse from those who are not. The pulse disturbance
front is a sphere (a cylinder in 2D problems) which expands with the speed of
sound a. Since the airfoil is moving, however, the sphere is also convected down-
stream with a speed U . At a certain time, the disturbance front will expand
and convect enough to reach (x, z). After this time (x, z) is influenced by the
pulse. Since U > a, however, after a certain time the sphere will be swept past
(x, z). After this time the effect of the pulse at (x, z) disappears. Note that
points above the dashed line in figure 6.2 never feel the effect of the pulse. The
dashed line therefore defines the limit of the pulse’s zone of influence. The cone
obtained by revolving the dashed line around the x-axis is called the Mach cone
(this becomes a wedge in 2D problems).
6.1.1 General approach:

We will represent the airfoil using a superposition of source pulses located on
z = 0 between x = −b and x = +b (see figure 6.3). Then we will:
1. Derive a general expression for the potential at (x, z) due to source-pulses

distributed along the airfoil
∂φ
2. Derive ∂z , and take the limit as z → 0±
3. Apply the flow-tangency condition to determine the source-pulse strength

distribution
∂φ ∂φ
4. Evaluate ∂t and ∂x from the determined source-pulse potentials
5. Compute the pressure distribution, forces, and moments

As explained above, the the effects of our component solutions are only felt over
limited extents of time and space. Thus, in order to complete step 1, we must
first establish formal limits of integration for the potential at (x, z).
φ (x,0,t)
x
−b b
Figure 6.3: Source-pulse representation of a thin airfoil
6.2 Time Integration Limits

Consider the pulse disturbance front a two different times, τ1 and τ2 , as shown
in figure 6.4. The point (x, z) first feels the potential of the pulse at τ1 , when the
expanding cylinder of the disturbance first makes contact. (x, z) is then kept
within the cylinder as it expands and convects until τ2 , when the expanding
disturbance front is finally swept past (x, z).
Cylinder position Cylinder position when

when first felt at (x,z) swept downstream past (x,z)
(x,z)
aτ2
a τ1
(ξ,ζ)
Uτ1
U τ2
Figure 6.4: Times of influence for a travelling pulse
Mathematically, the interval within which the influence of the pulse is felt
at (x, z) corresponds to when the expression for φs (x, z, t) is real. Defining
τ = t − T , we can determine the times τ1 , τ2 between which the pulse is felt at
6.3. SPATIAL INTEGRATION LIMITS 99
(x, z) by setting the denominator of the source function to zero:

p
M (x − ξ) ∓ (x − ξ)2 − (M 2 − 1)(z − ζ)2
τ1 , τ2 = (6.5)
a(M 2 − 1)
Thus when deriving an expression for the potential at (x, z), we need only to
integrate between times τ1 and τ2 .
6.3 Spatial Integration Limits

In the same way that we defined a rearward-looking Mach cone as containing
the region of influence of a pulse, we can also construct a forward-looking Mach
cone which defines the region of influence for a point. Such a forward-looking
Mach cone is shown in figure 6.5. Note that the half-angle of the cone, µ (the
Mach angle), is the same as that of the rearward-looking Mach cone.
2
µ= tan−1 [1/(M −1) ]
1/2
Forward−looking
Mach cone
(x,z)
−b ξ1
Figure 6.5: Spatial limits of integration for point (x, z)
With a little bit of geometry, we can establish the spatial limits of integration
for sources on the airfoil which can influence the point (x, z). If (x, z) is above
the airfoil, then the limits are:
p
ξ = −b to ξ1 = x − z M 2 − 1 (6.6)
if (x, z) is below the airfoil, then the limits are:
p
ξ = −b to ξ1 = x + z M 2 − 1 (6.7)
6.4 Potential induced by the airfoil

Assuming we have a continuous distribution of sources on the airfoil running
from −b to b, we can then express the potential induced at (x, z) as:
Z ξ1 Z τ2
A(ξ, 0, T ) dτ dξ
φ(x, z, t) = p (6.8)
−b τ1 a (t − T ) − [(x − ξ) − U (t − T )]2 − (z − ζ)2
2 2
which, noting that τ1 , τ2 are the roots of the denominator of the source function,
can also be expressed:
Z ξ1 Z τ2
1 Â(ξ, 0, t − τ )
φ(x, z, t) = √ p dτ dξ (6.9)
2 2
U − a −b τ1 (τ − τ1 )(τ2 − τ )
√
where we have chosen to re-express the integral in terms of Â = U 2 − a2 A. If
we use an angle transformation defined in terms of an angle γ as :
(τ2 − τ1 ) (τ2 + τ1 )
τ= cosγ + (6.10)
2 2
then after changing variables from τ to γ and substituting for dτ the denomi-
nator cancels perfectly, allowing us to simplify the expression to:
φ(x, z, t)
ξ1 π
1 (τ2 − τ1 ) (τ2 + τ1 )
Z Z
= √ Â ξ, 0, t − cos γ − dγdξ(6.11)
U 2 − a2 −b 0 2 2
ξ1 Z π
1
Z
= √ Â(ξ, 0, T ) dγdξ (6.12)
U − a2
2
−b 0
where
p
M (x − ξ) (x − ξ)2 − (M 2 − 1)z 2
T =t− − cosγ (6.13)
a(M 2 − 1) a(M 2 − 1)
6.5 Normal velocity at the Airfoil

To evaluate w = ∂φ
∂z , we must consider the z-dependence of the limit ξ1 and the
time-delay variable, both of which are ultimately arguments to Â:
∂ξ1 π

∂φ 1
Z
= √ Â ( ξ1 (z), T (t, ξ1 (z), z, γ) ) dγ
∂z U 2 − a2 ∂z 0
Z ξ1 Z π #
∂ Â ( ξ, T (t, ξ, z, γ) )
+ dγdξ (6.14)
−b 0 ∂z
which leads to:

∂φ 1 ∂ξ1 Mz
=√ π Â ξ1 , t − √
∂z U 2 − a2 ∂z a M2 − 1
Z ξ1 Z π #
∂ Â ((x − ξ)2 − (M 2 − 1)z 2 )−1/2
+ (zcosγ)dγdξ (6.15)
−b 0 ∂T a(M 2 − 1)
Taking the limit as you approach the airfoil from the upper side:

∂φ π
= − Â(x, 0, t) (6.16)
∂z + az=0
while the negative of this expression is true when approaching from the lower
side (see the definition of ξ1 ).
6.6. FLOW-TANGENCY CONDITION 101
xo
Figure 6.6: A pitching and plunging flat plate
6.6 Flow-Tangency Condition

From (6.16) we see that flow tangency can be expressed:
Â(x, 0, t) = −wF M (x, t)a/π (upper side) (6.17)

= +wF M (x, t)a/π (lower side) (6.18)
The source-pulse distribution is therefore directly related to the local vertical

velocity! The resulting potential function at the airfoil is:
Z x Z τ2
1 wF M (ξ, t − τ )
φ(x, t) = ∓ √ p dτ dξ (6.19)
2
π M − 1 −b τ1 (τ − τ1 )(τ2 − τ )
This expression can be used to find the force response of the airfoil to any
given wF M , by integration and substitution into the linearised pressure relation.
In general it is integrated numerically, but it is also possible to integrate it
analytically for the case of harmonic motion, as discussed in the next section.
6.7 Solutions for Harmonic Motion

To construct a harmonic solution we assume the following forms for the input
and response:
θ = θ̂eiωt ; h = ĥeiωt ; wF M (ξ, t − τ ) = ŵF M eiω(t−τ ) (6.20)
After some manipulation the potential at the airfoil can then be expressed as:
Z xh
2b i
φ(x, t) = ± √ U θ + ḣ + 2b(ξ − xo )θ̇ I(ξ, x, M, ω̄)dξ (6.21)
M2 − 1 0
h ω̄ i ωb 2M 2
where: I(ξ, x, M, ω̄) = e−iω̄(x−ξ) Jo (x − ξ) and ω̄ = (6.22)
M U M2 − 1
The details of the analytic procedure can be found in NACA report R846 by
Garrick and Rubinow, or in the text by Bisplinghoff. The harmonic case can
also be derived very efficiently using a Laplace transform, as shown in the text
by Fung.
The explicit expression for the potential can be differentiated with respect to
time and space, and the results substituted into the linearised pressure relation
to determine the pressure perturbation. This can be integrated along the chord
to arrive at expressions for the airfoil’s lift and moment, which are usually
written compactly as:
" #
2 2 iωt ĥ
L = 4ρbU k e (L1 + iL2 ) + θ̂(L3 + iL4 ) (6.23)
b
" #
2 2 2 iωt ĥ
M = −4ρb U k e (M1 + iM2 ) + θ̂(M3 + iM4 ) (6.24)
b
(6.25)
where L1..4 = L1..4 (M, ω̄) and M1..4 = M1..4 (M, ω̄) are coefficients dependent
on ω̄ and M . Alternatively, the lift and moment can be expressed in terms
of the Lh , Mh , Lθ , and Mθ coefficients introduced previously, which are then
functions of k and M (note that ω̄ is essentially a modified reduced-frequency
parameter). Figures 6.7 and 6.8 show the variation of these coefficients with
Mach number and reduced frequency. Both the real and imaginary components
of the response can be dramatically affected by these parameters. You are free
to experiment for yourself using the program at:
https://aerodynamics.lr.tudelft.nl/cgi-bin/lpfp
6.8 Finite wings, Multiple Surfaces

The 2D method described in the previous sections can readily be extended to 3D
wings, provided that the leading edges of the 3D wing are “supersonic”. A wing
with supersonic leading edges is shown in the left part of figure 6.9. A point
P lying on the wing surface is only aware of pulses from the upper surface,
as defined by its forward-looking Mach cone. There is no way a disturbance
from the wake or the lower surface can reach it. The situation is different for a
highly-swept wing operating at a low Mach number, as shown on the right. In
this case, a disturbance on the lower surface near the forward edge of the wing
will disturb the inflow region to point P (draw a rearward-looking Mach cone to
convince yourself of this). This makes the integration procedure more complex,
so it is normally necessary to resort to numerical methods. In the case shown,
P is also influenced by a portion of the wake. Therefore, it is also necessary to
impose the unsteady Kutta condition and the zero-pressure-change-across-the-
wake condition in order to obtain the correct solution.
When multiple lifting surfaces are present, it can easily be the case that the
wake of the upstream surface influences the forces and moments developed on
downstream surfaces (see, for example, problem 6.4). Where this is the case, all
6.8. FINITE WINGS, MULTIPLE SURFACES 103
0.5 3
0 2.5
-0.5 2
Lh’’
Lh’
-1 1.5
-1.5 1
-2 0.5
-2.5 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
k k
1.2 1.4
1 1.2
0.8 1
0.6
0.8
Mh’’
Mh’
0.4
0.6
0.2
0 0.4
-0.2 0.2
-0.4 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
k k
Figure 6.7: Plunging case - comparison of M=2.0 (solid), M=4.0 (dotted) and
M=0 (dashed) responses (a=-1/2). (Supersonic values multiplied by π)
2.5 3
2.5
2
2
1.5 1.5
Lt’’
Lt’
1 1
0.5
0.5
0
0 -0.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
k k
1.2 2
1 1.5
1
0.8
0.5
Mt’’
Mt’
0.6
0
0.4
-0.5
0.2 -1
0 -1.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
k k
Figure 6.8: Pitching case - comparison of M=2.0 (solid), M=4.0 (dotted) and
M=0 (dashed) responses (a=-1/2). (Supersonic values multiplied by π)
P P
Supersonic Leading Subsonic Leading

and Trailing Edges and Trailing Edges
Figure 6.9: Zones of influence for finite wings

appropriate boundary conditions must be applied to the wake of the upstream

surface.
p6.9.13 How long does a point (x, z) feel a pulse at (ξ, ζ, T ) if (ξ, ζ) is on the
forward-looking Mach cone of (x, z) ?
p6.9.14 Briefly outline how you could use linearised theory to determine the
lift response of an airfoil to a sharp-edged gust in supersonic flow. Use equations
from the notes.
p6.9.15 After how many half-chords travelled, s, (s = U ∆t/b) should an

indicially-plunged supersonic airfoil reach its steady-state lift?
√
p6.9.16 A wing-tail combination is travelling at a Mach number of M = 5,
when it is subjected to a step plunge with amplitude ḣo at t = 0.
e
d
c
P
U f
CG
c/ 4
x
a) Estimate the aerodynamic moment, P, about the CG at the initiation of

the plunge.
b) Assume the step plunge response is to be computed using a linear po-

tential method. Calculate the length of the wing wake which must be
included in the computation. What condition should be applied on the
wake?
Chapter 7
Arbitrary Unsteady
Aerodynamic Responses
Most of the analytic solutions we have considered up to now exploit mathe-

matical simplifications inherent with the use of harmonic functions. Unsteady
experimental data is also most often available in harmonic form, since it avoids
the application of sudden or irregular forces, and allows the certainty of the
data to be increased by averaging over several cycles.
If one is confronted with a problem where a non-harmonic response is re-
quired, one could in principle just use more complex versions of the same ana-
lytic or experimental procedures discussed up to now. This can be inconvenient,
however, if many different types of responses must be considered, or if the exact
form of the motion is itself determined by the solution of an aeroelastic interac-
tion. It is often far smarter to directly use the harmonic motion results in a way
which can provide predictions for aerodynamic responses to arbitrary motions.
A different reasoning might apply if one is using expensive numerical compu-
tations, or expensive flight tests to determine unsteady aerodynamic responses.
In such cases one may wish to limit measurement times, and therefore limit
the excitations considered to a few which contain a range of frequencies. One
would still like to use the measured results to predict arbitrary responses, how-
ever, including pure harmonic ones if necessary for flutter computations or for
comparison with harmonic experiments.
This section will describe how arbitrary responses can be computed from
indicial, impulse, harmonic, and state-space responses. For the moment we will
assume that the responses are linear, allowing us to use simple integral relations.
It also is possible to treat fully non-linear responses using similar techniques,
as we will see in a later chapter. In practice, however, the linear assumption
is often used even when non-linear conditions are considered, as linearisations
about non-linear base flows are often reasonably accurate for practical ranges
of perturbations.
107
108CHAPTER 7. ARBITRARY UNSTEADY AERODYNAMIC RESPONSES
∆τ
α(t)
α(t) ∆τ ::::
::::::::::::: ::::
::::::::::::: ::::
::::::::::::: ::::
::::::::::::: ::::
:::: t
t ::::
::::
τ ::::
∆α τ δα :::: α(t)
∆α = dα ∆τ ::::
Step
:::::::::::::
dτ Impulse
::::
Input
:::::::::::::
Input
::::
::::::::::::: t :::: t
:::::::::::::
dα A(t−τ) ∆τ
∆L dτ δL
Impulse α(t) h(t−τ) ∆τ
Step
Response Response
t t
Figure 7.1: Convolution of step (left) and impulse (right)responses
7.1 Indicial and Impulse Functions

The response to an arbitrary input can be obtained by summing indicial (step)
or impulse (spike) responses of a system. If the system has multiple inputs, one
indicial or impulse function per input degree of freedom is required. Indicial and
impulse functions are convenient for time-domain simulations where the final
form of the response is unknown. When identifying responses using numerical
methods (CFD), indicial inputs tend to be preferred over impulse responses due
to their less singular nature. In some cases smoothed indicial functions are also
used in order to limit the excitation of high-frequency spurious numerical modes
within the solution.
Given an indicial response, one can obtain an arbitrary response by convolu-
tion using Duhamel’s integral. Referring to figure 13.13(left), for a step change
in angle of attack, ∆α, we may write:
L(t) = A(t)∆α (7.1)
where A(t) is the indicial function. For a continuously changing α(t):
Z t
dα
L(t) = A(t)α(0) + (τ ) A(t − τ ) dτ (7.2)
0 dt
A similar approach can be followed to compute arbitrary responses using impulse
responses (figure 13.13(right)):
Z t
L(t) = α(t)h(t − τ )dτ (7.3)
0
where h(t) is the aerodynamic impulse response.

7.2. FOURIER TRANSFORMS 109
7.2 Fourier Transforms

If one has a description of the aerodynamic forces in the frequency domain,
Fourier transforms can be used construct responses to arbitrary inputs. For
reference, the definition of the Fourier transform pair is:
Z ∞
F (ω) = f (t)e−iωt dt (Direct Fourier Transform) (7.4)
−∞
Z ∞
1
f (t) = F (ω)eiωt dω (Inverse Fourier Transform) (7.5)
2π −∞
7.2.1 Computation of an arbitrary response

Assume we have an arbitrary input signal, for example an airfoil’s instantaneous
angle of attack, α(t), and unsteady results for the lift coefficient for a range of
harmonic excitations, CLα (ω). To compute the arbitrary response, we first
transform α(t) to the frequency domain using the direct transform to obtain
α(ω). Then we use the inverse transform to compute the time-domain lift
response as:
ρU 2 S ∞
Z
L(t) = α(ω)CLα (ω)eiωt dω (7.6)
4π −∞
We can think of taking the direct transform of α(t) as determining the harmonic
components of α, and the inverse transform as a summation of each component
multiplied by the CLα at the corresponding frequency. Although we will consider
analytic expressions here, in practice the procedure is often performed using fast
Fourier transform (FFT) techniques.
7.2.2 Re-expression of the incompressible results

Using the above ideas we can re-express our results for the 2D flat plate in
incompressible harmonic motion in order to compute arbitrary responses. Recall
that the component of the lift response dependent on the flow time history can
be expressed as a function of the motion of the three-quarter chord point:
L(t) = πρb2 [ḧ(t) + U∞ θ̇(t) − baθ̈(t)] + 2πρU∞ bC(k)[w 34 c (t)] (7.7)
where k = ωb/U∞, and w 43 c (t) can be computed for any combination of pitch,
θ, and plunge, h . If we focus on the second term, the frequency-domain repre-
sentation of its input is:
Z ∞
∗
W (k) = w 34 c (t)e−ikt dt∗ (7.8)
−∞
where t∗ = tUb∞ . W (k) can be seen as the magnitude of a group of sinusoids

with frequencies ranging from −∞ to ∞ which add up to w 43 c (t).
Since our aerodynamic theory is linear, we can replace the second term in
the lift expression with the sum of the responses to each one of the sinusoids:
Z ∞
∗
L(t) = πρb2 [ḧ(t) + U∞ θ̇(t) − baθ̈(t)] + ρU∞ b C(k) W (k) eikt dk (7.9)
−∞
This allows us to use the analytic expressions for C(k) in the computation of
arbitrary lift responses.
7.2.3 Conversion to an Indicial response

If we want to convert our frequency-domain representations to an indicial re-
sponse, we need to compute the Fourier transform of a step input. This can be
done by considering the limit as a → 0 for a one-sided decaying exponential:
g(t) = 0 t<0 (7.10)

−at
e t≥0 (7.11)
which for a > 0 has the Fourier transform:

1 a iω
G(w) = = 2 − 2 (7.12)
a + iω a + ω2 a + ω2
Since for a → 0
∞
a
Z
dω = π (7.13)
−∞ a2 + ω2
the Fourier transform for a unit step is finally:
1
πδ(ω) + (7.14)
iω
To speed manipulation, however, we will consider the step function to be defined
by (7.11) where a is very small but finite, so that its approximate transform is:
1
(7.15)
iω
This approximation allows us to express the indicial plunge lift response of a
flat plate in incompressible flow as:
∞
1 C(k) ikt∗
Z
φ(t∗ ) = e dk (7.16)
2π −∞ ik
where we have changed to the non-dimensional variables k and t∗ , and have also
ignored the initial infinite response peak arising from the apparent mass terms.
7.3. STATE-SPACE REPRESENTATIONS 111
7.2.4 Conversion from an Indicial response

For the analysis of aeroelastic instability boundaries, it can be convenient to
have the aerodynamic forces for harmonic motion. We can obtain these from
an indicial response using the direct Fourier transform.
For example, we can derive an approximate relation for Theodorsen’s func-
tion in the frequency domain using the approximation to the indicial function
in the time domain using:
Z ∞
C(k) ∗
= φ(t∗ )e−ikt dt∗ (7.17)
ik −∞
0.165k 0.335k
C(k) ≈ 1 − − (7.18)
k − 0.0455i k − 0.3i
Note that this approximation contains spurious poles at k = 0.0455i and k =
0.3i.
7.3 State-Space Representations

Control system designs often based on a “state space” approach, where the
system is represented in terms of first-order differential equations in time:
d
{x} = [A] {x} (7.19)
dt
For steady aerodynamic models, one can easily cast the dynamic equations in
this form by introducing additional variables for the derivatives of the structural
DOFs (degrees of freedom). In the case of unsteady aerodynamics, however, a
theoretically infinite number of additional variables would be required in order
to express the applied forces in terms of multi-variate Taylor expansions.
Using sophisticated identification techniques, it is possible to produce rea-
sonable approximations for the unsteady aerodynamic operators while keeping
the number of additional variables relatively low, typically on the order of 10
times the number of structural DOFs. (see Silva and Raveh, AIAA Paper 2001-
1213 for an example).
p7.4.17 The Sears function expresses the response of an airfoil to a sinu-

soidal gust. Derive an approximation for the Sear’s function, S(k), using the
approximation to Küssner’s function given in chapter 4.
Chapter 8
Compact Models of
Structural Systems
With modern computing facilities it is possible to resolve most of the detailed

dynamics of structural systems. Yet compact (low-DOF) structural models
are often used in practical aeroelastic computations, particularly during the
preliminary design stage. To understand why this is so, let’s consider the flutter
analysis technique we discussed in chapter 2. A complex representation of a
structure, such as that obtained by finite element analysis, will have many
degrees of freedom, or a large vector of unknowns x. In principal such a system
can also be treated with the standard flutter analysis technique, although it
may be necessary to use specialised procedures for finding the p roots of large
systems.
The problem arises when one tries to interpret the flutter diagrams. A
very large number of modes will be generated, and correspondingly a large
number of interactions will occur. A detailed analysis will be necessary to
determine which modes are physically relevant. This will involve, among other
procedures, rejecting modes of high spatial wave number due to their large
associated numerical errors. Depending on how complex the retained modes
are, however, it can still be a challenge deciding what is the most appropriate
action to change the behaviour of the system. In contrast, if one can choose
a reasonably-accurate compact representation of the system, the chances of
developing design insight are far higher.
This chapter will address the question of how to construct compact represen-
tations of complex structural systems, such as those present in aircraft. We will
start by discussing the concept of generalised coordinates, leading to a deriva-
tion of Lagrange’s equations. These provide a straightforward way to determine
the equations of motion in terms of compact sets of variables that are only
indirectly related to local structural displacements. We will then specialise La-
grange’s equations to the case of a free-flying vehicle, and discuss opportunities
for their simplification.
113
114 CHAPTER 8. COMPACT MODELS OF STRUCTURAL SYSTEMS
Finally, two approaches to deriving compact structural models using La-

grange’s equations will be demonstrated. The first of these, known as the
lumped-parameter approach, makes intuitive use of local mass-spring systems.
The second, known as the Rayleigh-Ritz approach, instead makes use of global
modes which are known to be representative of the structure’s motion. These
methods are discussed in an introductory manner. For further information the
reader is referred to one of many textbooks on the suject of structural dynamics.
8.1 Generalised Coordinates

The use of generalised rather than particle coordinates for representing the
motion of a body can lead to substantial savings in the numbers of degrees of
freedom required. For example, consider two representations of a rigid triangle
moving in the x-y plane, as shown in figure 8.1.
θ
3
L2
L3
2
(Xo,Yo)
L1
1
Figure 8.1: Particle (left) and generalised (right) coordinates
In figure 8.1 (left) the triangle is visualised as nodes connected by rigid bars.
In order to describe its motion, we use the particle coordinates of the nodes as
unknowns, and Newton’s second law to describe their evolution. This leads to
six equations of motion, as there are six particle coordinates:
rp = (xp , yp ) p = 1..3 (8.1)
Since the bars are rigid, we also need to apply three constraints between the
nodes of the form:
(x2 − x1 )2 + (y2 − y1 )2 = L21 , ... (8.2)
However, had we made a choice of non-particulate, generalised coordinates,

such as the displacement of the centroid and the rotation angle of the triangle,
q1 = Xo , q2 = Yo , q3 = θ we would immediately have arrived at only three
equations of motion, without the need for constraints.
In the particle coordinate case, the additional constraints are holonomic, in
that they can be written as explicit equations which express the behaviour of
the particle coordinates in time. In theory these can also be used to reduce
8.2. LAGRANGE’S EQUATIONS OF MOTION 115
the equations of motion to three. In many cases, however, constraints are non-
holonomic and can be difficult to apply (e.g. inequalities forcing a fuel particle
to remain in the fuel tank). By choosing our coordinates carefully, we can
apply the constraints implicitly, thereby reducing the degrees of freedom while
avoiding the need to deal with difficult constraints directly.
8.2 Lagrange’s Equations of Motion

Once appropriate generalised coordinates have been chosen, how does one de-
termine their associated equations of motion? This can be done by starting by
defining a transformation between generalised and particle coordinates:
rp = rp (q1 , q2 , q3 ... qN ) (8.3)
Using the chain rule, one may express a small change in a particle coordinate,
δrp , as:
N
X ∂rp
δrp = δqi (8.4)
i=1
∂qi
Similarly:
N
X ∂rp ∂ ṙp ∂rp
ṙp = q̇i from which: = (8.5)
i=1
∂qi ∂ q̇i ∂qi
Consider Newton’s second law for a system of M particles:
mp r̈p − Fp = 0 p = 1...M (8.6)
We may write the virtual work expression for the system of particles as:
M
X
(mp r̈p − Fp ) · δrp = 0 (8.7)
p=1
The first term of this equation may be expressed:

M M N
!
X X X ∂rp
(mp r̈p ) · δrp = mp r̈p · δqi (8.8)
p=1 p=1 i=1
∂qi
N X
M
X d ∂rp d ∂rp
= (mp ṙp · ) − mp ṙp · ( ) δqi (8.9)
i=1 p=1
dt ∂qi dt ∂qi
N X M
X d ∂ ṙp ∂ ṙp
= (mp ṙp · ) − mp ṙp · ( ) δqi (8.10)
i=1 p=1
dt ∂ q̇i ∂qi
N X M
X d ∂ ∂ mp
= − ṙp · ṙp δqi (8.11)
i=1 p=1
dt ∂ q̇i ∂qi 2
N
X d ∂T ∂T
= − δqi (8.12)
i=1
dt ∂ q̇i ∂qi
where
M
X mp
T = ṙp · ṙp (8.13)
p=1
2
is the total kinetic energy of the system. So the first term of the virtual work
expression for the system can now be expressed entirely in terms of the gener-
alised coordinates as the kinetic energy of the system is a scalar that does not
depend on what coordinates it is expressed in. The second term of the virtual
work expression for the system can be written:
M M N N
X X X ∂rp X
Fp · δrp = Fp · δqi = Q̃δqi ≡ δW (8.14)
p=1 p=1 i=1
∂qi i=1
where Q̃ is the net generalised force and δW is the virtual work.

It is ultimately required to express the behaviour of generalised forces di-
rectly in terms of generalised coordinates. In some cases this can be done
directly, such as when they explicit change in potential (strain) energy: U =
U (q1 , q2 , ..qN ) or dissipation D = D(q̇1 , q̇2 , ...q̇N ). For these cases the virtual
work can be expressed:
N
X
δW = −δU − δD + Qδqi (8.15)
i=1
N
X ∂U ∂D
= Q− − δqi (8.16)
i=1
∂qi ∂ q̇i
Substituting into the virtual work expression, and noting that by definition δqi
are arbitrary and independent, we can write:

d ∂T ∂T ∂U ∂D
− + + = Qi , i = 1, N (8.17)
dt ∂ q̇i ∂qi ∂qi ∂ q̇i
where:
T : Kinetic energy
U: Potential (strain) energy
D: Rayleigh’s dissipation function
qi : Generalised coordinate
Q: Remaining generalised forces
Equations (8.17) are known as Lagrange’s equations of motion.
8.3. LAGRANGE’S EQUATIONS FOR A FLEXIBLE AIRCRAFT 117
8.3 Lagrange’s Equations for a Flexible Aircraft

If one considers a flexible aircraft with deflections measured in an aircraft-fixed
frame which is rotating relative to the inertial frame of reference, Lagrange’s
equations can take on a relatively complex form. Consider the following defini-
tions, as shown in figure 8.2:
r=s+d (8.18)
r′ = r′o + r (8.19)
(8.20)
from which
dr
ṙ′ = ṙ′o + ω × r + (8.21)
dt
x, y, z : coordinates in local frame
x′ , y ′ , z ′ : coordinates in space-fixed frame
r′o : vector giving position of the local coordinate origin
r′ : vector giving position of an element in space-fixed frame
r: vector giving position of an element in local frame
s: undeflected position vector of an element in local frame
d: deformation vector of an element in local frame
ω: angular velocity of local frame
r
x
y
z
r’ r’
o
x’
y’
z’
Figure 8.2: Local and space-fixed frames

The kinetic energy terms of Lagrange’s equations then expand to:

1
Z
T = ṙ′ · ṙ′ ρdV (8.22)
2 V
1
Z
= [ṙ′ · ṙ′ + 2ṙ′o · (ω × r) + (ω × r) · (ω × r) (8.23)
2 V o o

dr dr dr dr
+2ṙ′o + 2(ω × r) + · ρdV (8.24)
dt dt dt dt
d
where dt denotes the time derivative in the local frame system and the dot
denotes the time derivative in the space-fixed frame. (recall ṙ = dr
dt + (ω × r)).
The virtual work terms are:
Z Z Z t Z
δW = F · δr′o dS + s×F·δ ωdt dS + F · δd dS (8.25)
S S 0 S
Where S is the surface of the aircraft and F is the force acting on the surface.
8.3.1 Options for efficient treatment

To keep things as simple as possible, a number of strategies are often applied:
1. Simplify the equations using orthogonal deformation modes as generalised

coordinates
2. Linearise the equations for small ω
3. Assume that rigid-body dynamics are decoupled from aeroelastic effects
(frequency separation) then:
• Solve the rigid-body dynamics problem neglecting flexibility
• Compute the resulting aeroelastic deformations
As well as reducing the number of terms in the equations, using orthogonal
modes such as the free vibration modes can be useful from the conceptual and
validation points of view. Linearizing the equations for small ω, on the other
hand can lead to a loss of accuracy, so it is only appropriate for small excursions.
Decoupling the rigid-body dynamics from aeroelastic motions is often done, but
is inappropriate for many large, flexible aircraft. For example, the B747 has
a short period mode ≈ 0.1 Hz while its fundamental symmetric wing bending
mode is ≈ 1 Hz. This is relatively little separation for analysis purposes.
8.4. METHODS FOR CONSTRUCTING COMPACT MODELS 119
8.4 Methods for Constructing Compact Models

There are systematic methods for choosing generalised coordinates and deriving
their equations of motion. We will consider two methods here, the lumped-
parameter method and the Rayleigh-Ritz method.
Lagrange’s Equations
Lumped-parameter Rayleigh-Ritz
method method
? ?
System of N exact System of N approximate
equations for the equations which converge to
lumped parameters an exact structural represen-
which approximate the tation as N → ∞
system
We will start by consider their application to beam-like wings. Before we do

so, it is worthwhile to review the following formulae for the kinetic and strain
energies of a beam of length l, mass distribution m(y), moment of inertia I(y),
and bending and torsional stiffness of EI(y) and GJ(y):
Kinetic and strain energies in bending
l 2

m dw
Z
T = dy (8.26)
0 2 dt
l 2
EI d2 w
Z
U = dy w = w(y, t) = bending deflection (8.27)
0 2 dy 2
Kinetic and strain energies in torsion
l 2
I dθ
Z
T = dy (8.28)
0 2 dt
l 2
GJ dθ
Z
U = dy θ = θ(y, t) = torsional deflection (8.29)
0 2 dy
8.5 Lumped Parameter Method

This straight-forward approach replaces the continuous system by a set of elas-
tically coupled rigid elements: The degrees of freedom associated with each
h
α
Figure 8.3: The lumped parameter method
element become the vector of generalised coordinates (Dimension N). To derive

the equations of motion:
1. Express the total kinetic, potential, and dissipative energies of the system
in terms of the generalised coordinates
2. Derive the generalised forces Qi by considering the work done on the

system by a virtual displacement of δqi
3. Express Lagrange’s equation for each generalised coordinate
(Note that for certain choices of generalised coordinates, the equations of

motion can also be derived in a straight-forward way by applying force and
moment equilibrium to each element.)
8.5.1 Example Application

Consider the two-degree-of-freedom system shown below, with masses m1 and
m2 , springs with constants k1 and k2 , dampers with constants c1 and c2 , and
applied force F2 (t) (figure 1). The first step is to define the generalised coor-
dinates. A good choice is q1 = x1 , and q2 = x2 . Next, we have to define each
component of the total energy:
1
T = [m1 ẋ21 + m2 ẋ22 ] (Kinetic) (8.30)
2
1
U = [k1 x21 + k2 (x1 − x2 )2 ] (Potential) (8.31)
2
1
D = [c1 ẋ21 + c2 (ẋ1 − ẋ2 )2 ] (Dissipative) (8.32)
2
8.5. LUMPED PARAMETER METHOD 121
9999 k1 k2
9999
9999
9999
9999
9999
F2(t)
m2
9999
m1
9999
9999
9999
9999 c1 c2
x1 x2
Figure 8.4: A 2 DOF Forced System
Then we need to consider the virtual work done by the force Fp . From the
expression for δW on slide 8-7, we can express the generalised force as:
M
X ∂rp
Qi = Fp (8.33)
p=1
∂qi
From our definition of the generalised coordinates, we find that Q1 = 0, and

Q2 = F2 (t). The final step is to substitute our expressions into Lagrange’s
equation:

d ∂T ∂T ∂U ∂D
− + + = Qi , i = 1, N (8.34)
dt ∂ q̇i ∂qi ∂qi ∂ q̇i
For i = 1, this produces:
m1 ẍ1 + k1 x1 + k2 (x1 − x2 ) + c1 ẋ1 + c2 (ẋ1 − x˙2 ) = 0 (8.35)
While for i = 2:
m2 ẍ2 + k2 (x2 − x1 ) + c2 (ẋ2 − x˙1 ) = F2 (t) (8.36)
You can verify these equations using a direct application of Newton’s second
law.
8.6 Rayleigh-Ritz (Assumed Mode) Method
γ1 γ2
1
0
0
1
0
1
0
1
0
1
0
1
0
1 w(y,t)
0
1
0
1 y y
0
1 y
0
1
0
1
0
1
Figure 8.5: The Rayleigh-Ritz method
In this case, the deflections are expressed in terms of a finite set of assumed
modes, the amplitudes of which are the generalised coordinates:
N
X
w(y, t) = γi (y)qi (t) (8.37)
i=1
where w(y, t) is the local deflection and γi (y) is the assumed mode. To derive
the equations of motion:
1. Assume N modes of unknown amplitude which satisfy the boundary con-

ditions and are linearly independent
2. Express the total kinetic, potential, and dissipative energies of the system
in terms of the assumed modes and generalised coordinates
3. Derive the generalised forces Qi by considering the work done on the

system by a virtual displacement of δqi
4. Express Lagrange’s equation for each generalised coordinate
In general, when applied to a system with zero damping, the Rayleigh-Ritz

method will produce a system of the form:
N
X N
X
Mij q̈j + Kij qj = Qi (t), i = 1, N (8.38)
j=1 j=1
This system can be considerably simplified if the chosen modes have the property
of orthogonality:
Z l
γi (y) γj (y) dy 6= 0 for i = j (8.39)
0
= 0 for i 6= j (8.40)
Orthogonal modes can be obtained by choosing γi (y) to be the undamped free-

vibration modes of the structure
8.6. RAYLEIGH-RITZ (ASSUMED MODE) METHOD 123
8.6.1 Example Application

Consider a bending wing of length l which has its root fixed to a wind-tunnel
wall (figure 2). The wing has a mass distribution m(y), and stiffness distribution
EI(y). The lift acting on the wing is L(y, t), and the resulting deflection is
w(y, t). We wish to derive the equations of motion of the wing in terms of N
assumed modes γi , where i = 1..N . We will assume the local deflection can
L(y,t)
9999
9999
9999
9999
9999
9999
w(y,t)
9999 y
9999
9999
9999
9999
l
9999
Figure 8.6: A Bending Wing
be written as a linear combination of assumed mode shapes, the amplitudes of

which vary in time:
N
X
w(y, t) = γi (y)qi (t) (8.41)
i=1
The kinetic and strain energy of the beam are given by:
1 l
Z
T = m(y)ẇ2 (y, t)dy (Kinetic) (8.42)
2 0
1 l
Z
U= EI(y)(w′′ )2 (y, t)dy (Strain) (8.43)
2 0
d
Here ˙ denotes the time derivative, and ′ denotes dy . The above equations can
be re-expressed in terms of the assumed modes as:
N N
1 l
Z X X
T = m(y)( γi (y)q̇i )( γi (y)q̇i )dy (8.44)
2 0 i=1 i=1
N N Z l
1 XX
= q̇i q̇j m(y)γi (y)γj (y)dy (8.45)
2 i=1 j=1 0
N N l
1 XX
Z
= Mij q̇i q̇j ; where Mij = m(y)γi (y)γj (y)dy (8.46)
2 i=1 j=1 0
Similarly for the strain energy:

N N l
1 XX
Z
U= Kij qi qj ; where Kij = EI(y)γi′′ (y)γj′′ (y)dy (8.47)
2 i=1 j=1 0
Finally, we define the ith generalised force by considering the work done by
L(y, t) for a generalised displacement δqi :
Z l
Qi (t)δqi (t) = (L(y, t)δw(t))dy (8.48)
0
Z l
= L(y, t)γi (y)δqi (t)dy (8.49)
0
or:
Z l
Qi (t) = L(y, t)γi (y)dy (8.50)
0
Substituting these expressions into Lagrange’s equation yields the following sys-
tem of N equations:
N
X N
X
Mij q̈j + Kij qj = Qi (t) (8.51)
j=1 j=1
As N → ∞, this system will approximate the dynamics of the continuous beam.

Often, however, only a few modes will suffice, especially if one is primarily
interested in the lower-frequency structural modes. For the latter case, a good
choice for the γi are the beam’s natural modes. As an alternative, Bisplinghoff
gives the following polynomial formula for a bending beam:
y y
γi (y) = 1 − (i + 3) − (1 − )i+3 (8.52)
l l
p8.7.18 Derive the equations of motion for the lumped-parameter model of a

cantilever beam subject to torsion as shown below. Find the natural modes of
the system and verify that they are orthogonal.
!!!!!!!!!!!
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
k1=k2
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
i1 = i2
::::
!!!!!!!!!!! k 1
::::
!!!!!!!!!!!
:::: !!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!!
:::: !!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
θ1
::::
!!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
::::
!!!!!!!!!!!
:::: !!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!!
:::: !!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!! i1
!!!!!!!!!!
!!!!!!!!!! k2
!!!!!!!!!!
!!!!!!!!!! !!!!!!!!!!θ2
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!! !!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!! !!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
i2
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
Figure 8.7: Lumped-parameter model
p8.7.19 Consider the continuous beam subjected to torsion shown below.

Apply the Rayleigh-Ritz method with two assumed modes to determine the
natural modes of the beam. How do these compare with those obtained from
the lumped-parameter model approach in question 1?
!!!!!!!!!!!
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!! Ιθ(y) = Constant
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
GJ(y) = Constant
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
y=0
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
θ(0) = 0
!!!!!!!!!!! θ(y)
::::
!!!!!!!!!!!
!!!
!!!
y = !!!
1
θ’(1) = 0
Figure 8.8: Continuous beam

Chapter 9
Flutter Calculation
Methods
In the last few chapters, we have considered more advanced representations

for unsteady aerodynamic forces and structural systems. Now we return to the
problem of flutter prediction. Like the simple two-degree-of-freedom examples of
chapter 2, our more advanced representations still lead to a system of second-
order differential equations, although now expressed in terms of generalised
coordinates:
[M ]{q̈} + [D]{q̇} + [K]{q} = {QA } (9.1)
where:
[M ] : Mass matrix
[D] : Structural damping matrix
[K] : Stiffness matrix
{q} : Vector of Generalised Coordinates
{QA } : Generalised aerodynamic forces
In chapter 2 we produced flutter diagrams by assuming a solution of the

form {q} = {q̂}ept , and solving the resulting characteristic equation for p. That
can still be done here, provided that the aerodynamic forces QA , can be ex-
plicitly expressed as algebraic functions of {q}. This is the case for the simple
aerodynamic models considered in chapter 2, as well as the state-space approx-
imations described in chapter 7. It is not the case when only pure harmonic
results are available, or when the analytic solution is expressed in terms of tran-
scendental functions, or when the data is available only in tabular form from
either experiments or CFD computations. In these circumstances, more general
flutter-calculation methods are required.
A brute-force approach is to add a harmonic forcing term F (t) to the equa-
tions and evaluate the system response as a function of operating conditions
127
128 CHAPTER 9. FLUTTER CALCULATION METHODS
U
Figure 9.1: Brute-force amplitude analysis
(e.g. the flight speed U ) and the frequency, ω, of F (t). Typically this results
in system response amplitude plots such as that shown in figure 9.1. Although
this procedure is analogous to those used in flight flutter testing, it is quite in-
efficient from the numerical point of view as it requires harmonic time-domain
simulations for each point in the U − ω plane.
Fortunately there are a number of more efficient methods, the most common
of which will be be described in this chapter. Each of these avoids the costs of
the brute-force approach by directly generating an approximate flutter diagram.
The first two methods considered, the k method and p-k method, assume that
only harmonic aerodynamic data are available. This is sufficient for the predic-
tion of flutter boundaries, as we know that at each flutter point, the response
will be purely harmonic. The last method we consider, the p method, is more
general in that it can also predict the response away from the flutter boundary,
but requires enough aerodynamic data to construct an approximate model for
QA as a function of p.
9.1 Non-dimensionalisation
For the remainder of this chapter, we will employ a non-dimensional expression
for p:
b
p= (σ + iω) = δ + ik (9.2)
U
9.2. K METHOD 129
Ut
and assume solutions of the form {q} = {q̂}ep b . The dynamic equations of
motion can then be written as:
2
U 2 U 1
[M ]p + [D]p + [K] {q̂} = ρU 2 [A(p)]{q̂} (9.3)
b2 b 2
Where A(p) is also a function of the Mach number M∞ and the Reynolds number
Re. In general, A(p) can be non-linear, and need not be expressed in terms of
explicit functions.
9.2 k method
The k method arose during the 1940’s as a technique to efficiently make use of
pure harmonic aerodynamic data during flutter analysis. The idea is simply to
add a ficticious structural damping term to the equations, and then to determine
its required value in order to keep the response purely harmonic. The flutter
boundary is then determined by finding the point at which the required ficticious
damping becomes zero. Note that since the response is always kept purely
harmonic, the use of harmonic aerodynamic data is justified throughout the
analysis.
9.2.1 Procedure
If we assume a ficticious structural damping, g, and a purely harmonic response
p = ik, the equations of motion become:
1
−ω 2 [M ] + (1 + ig)[K] {q̂} = ρU 2 [A(k)]{q̂}

(9.4)
2
We can then set up the following eigenvalue problem:

1 1 + ig
[K]−1 [M ] + ρb2 [A(k)/k 2 ] {q̂} = {q̂} (9.5)
2 ω2
or:
1 + ig
[B(k)]{q̂} = λ{q̂}, λ= (9.6)
ω2
where[B(k)] and λ = λ′ + iλ′′ are, in general, complex. Now to generate a point
on the k method’s “flutter diagram”, we follow the following procedure:
1. Choose an altitude, which results in a value for ρ
2. Choose k and calculate [B(k)], which will include the aerodynamic data
at frequency k
3. Compute the eigenvalues for each mode λ1 ...λN
p ′
4. Compute the frequency of each mode ωi = 1/ λi
5. Compute the flight speed Ui = ωi b/k, and fictitious damping

′′
gi = ωi2 λi corresponding to each mode
This procedure is repeated at the same altitude for different k’s , which generates
diagrams similar to that shown in figure 9.2 (note that a constant k is a straight
line on the U − ω plane). An estimate for the flutter speed, UF , can be obtained
by interpolation using two values near gi = 0.
k2
ω kF
k1
k2
k1
U
+
g
U
UF
Figure 9.2: k-method diagrams for a two-mode problem
The k method is efficient when the aerodynamic forces do not vary signifi-
cantly with M∞ or Re, as it then allows the flutter point to be found after a
small number of order N eigenvalue problems, where N is the number of gen-
eralised coordinates. It is worth noting that at the flutter point, no additional
sources of error are present, so the flutter solution is exact.
When the aerodynamic forces do vary significantly with M∞ or Re, addi-
tional iterations are required. This is because the flight velocity is determined
from the computation after A(k) has been specified. Additionally, when flight
control systems are involved for which the frequency parameter is ω instead of
k, an additional interpolation step has to be made to match both parameters.
This can make the application of the k method unwieldy in general design cases.
Outside of these efficiency issues, the k method suffers from practical problems
in implementation and interpretation. These include difficulties associated with
connecting the solution pairs in the flutter diagram to form sensible frequency
and damping curves, especially when the curves merge. In some cases it is
9.3. P-K METHOD 131
necessary to calculate the eigenvectors and to inspect their development along

the frequency and damping curves. Another problem is the appearance of more
frequency and damping solutions than the number of degrees of freedom. An
example of such anomalous results is shown in figure 9.3. This phenomenon
typically occurs when the aerodynamic forces are large with respect to inertia
and stiffness forces, for example in cases with free rotation of control surfaces.
+
g
Figure 9.3: Non-physical behaviour of k-method diagrams
9.3 p-k method

The p-k method emerged in the mid 1960’s as an alternative to the k method
with better properties for aerospace problems. It is based on the assumption
that we may approximate the aerodynamics of sinusoidal motions with slowly
increasing or decreasing amplitudes using purely harmonic aerodynamic results.
i.e. that:
[A(p)] = [A(δ + ik, M∞ , Re, ...)] ≈ [A(k, M∞ , Re, ...)] (9.7)
The full equations of motion are then used to determine the correct value of p
for each mode at a given flight condition.
9.3.1 Procedure
With our assumed form of the aerodynamic forces, the equations of motion can
be written:
2
U U 1
2
2
[M ]p + [D]p + [K] {q̂} = ρU 2 [A(k, M∞ , ...)]{q̂} (9.8)
b b 2
or
[F (p, k)] {q̂} = 0 (9.9)
where p = δ + ik. This equation can be used to determine values of p. Normally

an iterative solution method is used, with guesses for p determined from the
previously analysed flight condition. The most common technique is known as
“determinant iteration” (see Hassig, H, “An Approximate true damping solution
of the flutter equation by determinant iteration”, J. Aircraft, Vol. 8, No. 11,
Nov. 1971).
The basic procedure is as follows. First the mode which will be tracked is
chosen, along with a starting value for U near U = 0. The Mach and Reynolds
number are then computed for that U . A first guess, p1 , can be made using the
natural frequency of the mode in the absence of aerodynamic forces. Another
initial guess, p2 , can be made by adding a small value of damping to p1 . Then
the following sequence of steps is performed:
1. Compute or interpolate for [A(k1 , M∞ , ...)] and [A(k2 , M∞ , ...)].
2. Compute the determinants F1 =| [F (p1 , k1 )] | and F2 =| [F (p2 , k2 )] |
3. Update p using: p2 F1 − p1 F2 ;
p3 =
F1 − F2
4. Setp1 = p2 and p2 = p3 , then repeat from step 1 until converged
5. Choose the next flight condition (ρ, U , M∞ , Re, ...)
6. Choose two new values of p based on extrapolation from the previous flight
conditions
7. Return to step 1 and repeat to find the new p at the next flight condition
Once the mode is tracked through its desired range, a different mode is selected
and the entire procedure is repeated. The formula for p3 in step 3 can be
derived using the diagram shown in figure 9.4. Note that only one determinant
evaluation is required per iteration as after the procedure is started, as F1 can
be taken from the previous iteration.
Although the process of finding a values of p can be costly for complex F ,
this must be compared with the high cost of finding eigenvalues for large systems
as required by the k method. Additionally, no further iteration is required by
the p-k method to incorporate Mach or Reynolds number effects. An additional
9.4. P METHOD 133
F1
p3 p
0 2 p
p
1
F
2
Figure 9.4: The determinant iteration procedure
control transfer function [H(p)]{q̂}, can also be easily incorporated. In general,

the flutter diagrams from the p-k method are much closer to those of the actual
system than those obtained via the k method. They are still only exact, however,
at the flutter point.
9.4 p method
If we want flutter diagrams which are accurate over the complete range of condi-
tions, we will be forced to construct an approximate expression for [A(p)] which
is accurate for decaying or growing motions. One of the major complicating fac-
tors in doing so is representing the time lags inherent in unsteady aerodynamic
responses.
Once [A(p)] is specified, the resulting system can be solved using the same
determinant iteration procedure used in the p-k method. Alternatively, A[(p)]
can be expressed in the time domain using an inverse Laplace transform, and the
solution carried out as an eigenvalue problem in state space. The dimension of
state-space problem is usually higher, as new state variables have to be introduce
to properly represent the aerodynamic lag modes.
Approximate models for [A(p)] are typically of the form:
L
" #
2
X A(l+2)ij (M∞ )p
A0ij (M∞ ) + A1ij (M∞ )p + A2ij p + (9.10)
βlij + p
l=1
• The first three terms are instantaneous in nature and are often referred to
as aerodynamic stiffness, damping and mass terms
• Terms l = 1 to L approximate aerodynamic lag effects
• Coefficients Aij , βlij cab be determined by fitting to computed or experi-

mental results
• A2ij is usually zero, except for incompressible flow
It is instructive to compare the lag effect terms above to the approximation for
Theodorson’s function given in chapter 7 (Note that the latter approximation
accounts for the circulatory part of the response only).
The main advantage of the p method is that it can accurately represent non-
harmonic aerodynamic responses, including converging, diverging and aperiodic
motions. In its state-space form, it is also effective for control system design. Its
disadvantage is that it can require a complex procedure for accurately approx-
imating the aerodynamic forces, and its state-space approximations can lead
to very large eigenvalue problems. Also, the p method will find a slightly dif-
ferent value of flutter speed than determined using k and p-k methods, due to
additional interpolation errors in the aerodynamic forces.
p9.5.20 Consider the following single-DOF system:

1 2
[M ]q̈ + [K]q = ρU [A]q
2
with ρ = 1, b = 1, M = 2, and K = 1.
Given the following two experimental results:
k=1 A(k) = 4 + 2i
k=2 A(k) = 2 − 4i
compute the flutter speed of the system using the k method.
p9.5.21 For the single DOF system described above, use the p-k method to
estimate p = δ + ik at U=0.4. Use an assumed value of δ1 = −0.1 with the
experimental A(k) for k = 1, and δ2 = 0.1 with the experimental A(k) for k = 2.
Perform a single step of determinant iteration.
Chapter 10
Dynamic Responses
The need to ensure aeroelastic instabilities occur outside of an aircraft’s flight

envelop is self-evident, and is a natural concern for aeroelasticians. There are
also many important aeroelastic design problems, however, which occur well
away from instability boundaries. The most common of these are related to
achieving desired dynamic response characteristics. Such problems are becoming
more and more prevalent as we move towards lighter, more flexible structures,
particularly when active control systems are used.
The sources of excitation for dynamic responses are varied. They include at-
mospheric turbulence, manual or automated control deflections, water or ground
forces transmitted through the landing gear, or unsteady aerodynamic interfer-
ence occurring between aircraft components. In some cases failure modes asso-
ciated with dynamic responses can be critical. This can be due to the possibility
of structural failure, as is the case for tail surfaces subjected to strong atmo-
spheric gusts (figure 10.1). Alternatively, an undesired elastic response could
lead to loss of control.
Elastic Tail
Tail−
Root
Bending
Moment
Rigid T a i l
Half Chords Travelled
Figure 10.1: Response of a tail to sharp-edge gust
135
136 CHAPTER 10. DYNAMIC RESPONSES
Analysis of cases with less severe consequences can also be useful. Passenger
comfort, for example, is affected by the details of response to atmospheric turbu-
lence. On the other hand, non-critical dynamic responses can also be measured
and used by aeroelasticians in order to validate the numerical models they have
developed.
From the analysis point of view, dynamic response problems are more gen-
eral than dynamic instability problems, as we can no longer exploit the sim-
plifications inherent in assuming undamped harmonic solutions. In fact, the
complexity of modern aircraft structural and control systems, coupled with the
difficulties of accurately predicting unsteady aerodynamic flows, means that
dynamic response problems can pose a significant challenge.
10.1 Types of Dynamic Response problems

We can broadly divide dynamic response problems into two classes, determinis-
tic and stochastic. In the first class the excitation is know precisely. Examples
include specified control deflections, idealised test gusts or measured forces from
initial touchdown. Such problems may be treated directly by time-domain sim-
ulations. If the system is linear, we can also treat the problem using frequency-
domain techniques, providing another viewpoint which often aids interpretation.
The second class of problems are those with excitations which are only quan-
tified statistically. Examples include general atmospheric turbulence or landing
loads imparted by water or uneven terrain. In recent years a number of ad-
vanced techniques for analysing such stochastic systems have been developed.
We will limit our discussions to simple stochastic problems, however, which can
be treated using frequency-domain techniques. An overview of problems and
their analysis appears in figure 10.2.
10.2 Dynamic Response Equations

Most aeroelastic systems can be expressed as a system of equations of the form:
1 2
[M ]{q̈} + [D]{q̇} + [K]{q} = ρU {A(t, q, v, M∞ , Re)} + {E(t)} (10.1)
2
[M ] : Mass matrix
[D] : Structural damping matrix
[K] : Stiffness matrix
{A(t..)} : Aerodynamic forces
{q} : Vector of generalised deformation coordinates
{v} : Vector of external aerodynamic excitations
{E(t)} : External mechanical forces (e.g landing gear)
Here {A(t)} defines the combined unsteady aerodynamic forces arising from
both excitation (gusts) and deformation. For a truly general analysis, this must
10.2. DYNAMIC RESPONSE EQUATIONS 137
Dynamic Response Analysis
? ?
Deterministic
- Control deflections Stochastic
- Prescribed gusts - Atmospheric turbulence
- Initial landing impact - Ground roll, water waves
wg wg
? ?
Time Frequency ?
Domain Domain Frequency
- Numerical - Fourier Domain
integration Transform - PSD Technique
Figure 10.2: Dynamic Response Analysis
be determined together with the response. The most practical way to do so is

using time integration methods.
If the aerodynamic interaction between the excitation and the responding
configuration is low, we can split {A(t)} into excitation and response compo-
nents. This allows the use of simplified analysis techniques. The assumption
can be expressed as the the independence of the aerodynamic forces due to
excitation (v) from those arising due to deformation (q). This results in the
following equations of motion:
1 2
[M ]{q̈} + [D]{q̇} + [K]{q} = ρU ([A(t)]{q} + [Ae (t)]{v}) + {E(t)} (10.2)
2
[A(t)] : Aerodynamic forces due to {q}

[Ae (t)] : Aerodynamic forces due to {v}
A common example is the frozen gust assumption, where the gust is represented
by a vertical velocity distribution fixed in space, assumed to remain constant
during the passage of the wing (figure 10.3). For certification purposes, stan-
dardised gusts, such as the sharp-edged and (1-cosine) gusts may be used to
establish minimum response criteria (figure 10.4).
z
U x
y
U
wg αg
x
Figure 10.3: Frozen gust assumption
wg wg
sharp-edge gust 1-cosine gust
Figure 10.4: Standardised gusts
10.3 Solution in the frequency domain: Fourier

Transform
If we can separate the aerodynamic forces as described in the previous section,
then equation (10.2) can be expressed for a single harmonic frequency as:

1 2 1
2
−ω [M ] + iω[D] + [K] − ρU [A(k)] {q̂} = ρU 2 {Ae (k)}{v̂} (10.3)
2 2
or equivalently:
[Dq (ω)]{q̂} = [Dv (ω)]{v̂} (10.4)
The ratio of output to input is then given by:
[H(ω)] = [Dq ]−1 [Dv ] (10.5)
[H(ω)] is known as the mechanical admittance. It can be pre-computed before
considering a specific gust problem using harmonic experimental, analytic or
numeric data. Assuming [H(ω)] has been defined, we can compute the response
to an arbitrary gust using the following procedure:
Step 1: Transform {v} from time to frequency domain:

Z ∞
{v(ω)} = {v(t)}e−iωt dt (10.6)
−∞
10.4. SOLUTION IN THE TIME DOMAIN: NUMERICAL INTEGRATION139
Step 2: Compute response in frequency domain:
{q(ω)} = [H(ω)]{v(ω)} (10.7)
Step 3: Return response to time domain

Z ∞
1
{q(t)} = {q(ω)}eiωt dω (10.8)
2π −∞
In practice fast Fourier transforms are often used to perform these steps.
10.4 Solution in the time domain: Numerical

Integration
There is a wide range of numerical methods available for the integration of
systems of ODE’s. These are typically based on either finite-difference or vari-
ational techniques. Some of these methods require expressing the system of N
second-order equations as a system of 2N first-order ODE’s. This can be done
by introducing the state vector r = {q1 , q̇1 , ....., B1 , B2 ...}t , and the correspond-
ing excitation vector v = {v1 , v1′ , .....}t , where the extra state variables B1 , B2
are usually required in order to properly simulate aerodynamic lag modes. The
final system of first-order equations can be expressed:
ṙ = [A]r + [B]v (10.9)
A typical a numerical integration technique is the finite difference “θ-method”:
rn+1 = rn + ∆t (1 − θ)([A]r + [B]v)n + θ([A]r + [B]v)n+1

(10.10)
The above formula allows the solution at t = tn + 1 to be computed from that at

t = tn . This formula is applied recursively to advance the solution in time. For
θ = 0, the method reduces to Euler explicit integration, while for θ = 1, Euler
implicit integration is obtained. Both of these are O(∆t) accurate in time. The
trapezoidal method, which is O(∆t2 ) accurate in time, is obtained when θ = 12 .
Example: Numerical integration of a 2 DOF System:
The equations of motion for our standard 2 DOF (h,θ) representative section
with a low-frequency aerodynamic model (see chapter 2) can be written:
[M ]ṙ + [C]r = v (10.11)
where
r = {h, ḣ, θ, θ̇}t ,

v = {0, −CLα qS w(t)/U, 0, CLα qSec w(t)/U }t ,
and:
 
1 0 0 0
 0 m 0 Sθ 
[M ] = 
 0 0
 (10.12)
1 0 
0 Sθ 0 Iθ
 
0 −1 0 0
 Kh CLα qS/U CLα qS 0 
[C] = 
 0
 (10.13)
0 0 −1 
0 −CLα qSec/U Kθ − CLα qSec 0
We will consider some dynamic responses for a section with the following prop-
erties:
ρ = 0.53kg/m3
e = 0.2, c = 6 m
m = 400 kg/m
Sθ = 180 kg m/m
Iθ = 200 kg m2 /m
Kh = 1 × 105 N/m /m,
Kθ = 3 × 105 N m/m
For reference, the associated flutter diagram is shown in figure 10.5.
60
50
40
30
20
10
0
0 50 100 150 200
Figure 10.5: Flutter Diagram (low-frequency aerodynamic model)
Figures 10.6 and 10.7 show the response of the system to a sharp-edged gust with
a magnitude of 10m/s encountered at t = 0. Initially there is increased damping
10.5. STOCHASTIC RESPONSE ANALYSIS - PSD TECHNIQUE 141
0.08 0.2
0.06 theta theta
h 0.15 h
0.04
0.1
0.02
0 0.05
-0.02 0
-0.04 -0.05
-0.06 -0.1
-0.08
-0.15
-0.1
-0.12 -0.2
-0.14 -0.25
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Figure 10.6: Responses at U=60 m/s and U=100 m/s
as the velocity is increased, as can be anticipated from the flutter diagram. As
0.3 80
theta theta
h 60 h
0.2 40
0.1 20
0
0 -20
-40
-0.1 -60
-0.2 -80
-100
-0.3 -120
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Figure 10.7: Responses at U=110 m/s and U=120 m/s
the speed is increased further however, the damping decreases until the system
becomes unstable. Depending on the maximum allowed strain of the structure,
the 100 and 110 m/s cases may already be critical, as these produce relatively
large excursions.
10.5 Stochastic Response Analysis - PSD Tech-

nique
Many types of excitations, including atmospheric turbulence, can be appropri-
ately modelled as random processes. In such cases it is best to define their
properties in terms of average quantities such as:
the mean value of a variable:

T
1
Z
x̄ = lim x(t)dt (10.14)
T →∞ 2T −T
and the mean square value (or power) of a variable:

Z T
¯ 1
2
x = lim x2 (t)dt (10.15)
T →∞ 2T −T
The mean square value can also be expressed in terms of the Power Spectral
Density (PSD) function:
Z ∞
x¯2 = Φx (ω)dω (10.16)
0
where:
2πX(iω)X(−iω)
Φx (ω) = (10.17)
T
and
T
1
Z
X(iω) = x(t)e−iωt dt (10.18)
2π −T
The function Φx (ω) may be thought of the distribution of power across the
frequency domain. Expression (10.16) can be derived using the autocorrelation
function for x assuming it to be defined by a processes which is statistically
stationary (its statistical properties are independent of time) homogeneous (its
statistical properties are independent of space) and isotropic (its statistical prop-
erties are independent of direction) [see “Aeroelasticity” by Bisplinghoff et al,
appendix C ].
Power Spectral Density Technique
If we assume that the response to the excitation is also a random function, then
it can be shown that for a linear system:
Φr (ω) =| Hrw (ω) |2 Φx (ω) (10.19)
where:
Φr (ω) : response power spectral density

Φx (ω) : input (e.g gust velocity) power spectral density
Hrx (ω) mechanical admittance
This means that given the PSD of an input signal, we can quickly get an idea
of the power frequency distribution of the output if we have the system’s me-
chanical admittance. The complete analysis can be summarised as:
10.5. STOCHASTIC RESPONSE ANALYSIS - PSD TECHNIQUE 143
1. Obtain the input PSD , Φx (ω)

2. Determine Hrx (ω) for the system
3. Multiply | Hrx (ω) |2 by Φx (ω) to obtain response PSD, Φr (ω)
An example of two commonly-used input PSD functions for atmospheric turbu-
lence are shown in figure 10.8. Both of these are biased towards low frequencies,
Dryden Spectrum
Von Karman Spectrum
Gust PSD function
Frequency
Figure 10.8: Standard PSD functions for atmospheric turbulence
indicating atmospheric turbulence tends to excite low-frequency modes most

effectively. This is shown in figure 10.9, where the admittance function of the
2DOF representative section is plotted with an input PSD for atmospheric tur-
bulence and the corresponding output PSD. The input PSD is seen to act as a
filter which tends to enhance the relative power of the low-frequency mode.
1.2
|H|^2
Phi g
1 Phi r
0.8
PSD function
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Frequency
Figure 10.9: PSD response to atmospheric turbulence

p10.6.22 An wing with mass m enters a stochastic gust field with a power
spectral density, Φg (ω), approximated by:
A
Φg (ω) =
B + ω2
The wing moves with a single degree of freedom, h, and is supported by a spring
Kh . Its lift may be approximated by:
!
ḣ wg
L(t) = CLα qS +
U U
L(t)
wg
h K
h
111111111111
000000000000
000000000000
111111111111
Compute the mechanical admittance, H(ω), and the PSD of the response,
Φr (ω). If A and B are positive, will increasing Kh increase or decrease the
maximum response amplitude? Use sketches of H(ω) and Φg (ω) to illustrate
your answer.
p10.6.23 Consider the single DOF plunging airfoil of the previous problem.
Using the frozen gust assumption and results from unsteady incompressible flat-
plate theory, express the response PSD as a function of reduced frequency for
an arbitrary gust PSD function, Φ(k). Assume that the lift due to a harmonic
gust with amplitude ŵ can be expressed:
Lg = 2πρU b [S(k)] ŵ
Chapter 11
Introduction to
Computational
Aeroelasticity
In this short chapter we review some advantages and disadvantages of perform-

ing high-fidelity numerical simulations of the flow and structure in order to
estimate the aeroelastic performance of a configuration. Methods for doing so
will be considered in the following chapters.
With the widespread availability of computing facilities, one might ask with
justification, why would we not always approach aeroelastic problems using
highly-detailed numerical simulations? After all, these do not employ geometric
simplifications, and can automatically account for non-linear effects in both
the flow and the structure. Indeed these are powerful advantages which have
motivated the increasing use of numerical simulations in industry.
However it is worthwhile to acknowledge some of the limitations of the cur-
rent state-of-the art. A prime motivation for the use of high-fidelity numerical
simulations is the fact that many aeroelastic phenomena of interest occur in off-
design conditions (e.g. flutter), which requires the use of sophisticated modelling
techniques. This has two consequences:
High computational cost

The high cost of aeroelastic simulations is primarily associated with the need
to resolve unsteady flow phenomena. For models applicable to off-design con-
ditions a full volume mesh is normally required, which must in turn be suffi-
ciently resolved to track features which are moving relative to the mesh. Also
the simulation must be time accurate, and performed for the relatively long
time intervals associated with structural vibrations. This means that individual
aeroelastic simulations tend to be far more expensive than their steady design-
point counterparts.
Not only is the cost of individual simulations high, but many of them must be
147
148CHAPTER 11. INTRODUCTION TO COMPUTATIONAL AEROELASTICITY
performed. To establish a flutter envelope, for example, a wide range of airspeed-

altitude combinations must be considered, and at each of these, several inertial
loading scenarios (fuel, payload) must be simulated.
Difficulties with Validation

Each aeroelastic simulation is subject to:
• Flow / turbulence modelling errors
• Nonlinear structural modelling errors
• Grid and time step dependent errors
• Fluid-structure coupling errors
Of these, the flow and structural modelling errors are normally the most signif-
icant. In particular, off-design flow problems including separation and shock-
boundary layer interaction are particularly difficult to compute, even using the
most advanced analysis methods. This means that experimental validation is
often required.
Ground vibration tests can provide a reasonable means to validate the fi-
delity of a structural model from the aeroelastic viewpoint. Unsteady-flow
experiments, however, are complex and are carried out much less frequently.
Fully-coupled aeroelastic experiments are very difficult to perform reliably and
are relatively rare.
Given these difficulties, the objective of multidisciplinary design optimisa-
tion (MDO) including off-cruise conditions is indeed an ambitious one. As
design loops require a significantly higher number of computations than stan-
dard aeroelastic analyses, simplified models are still widely used. Adjoint-based
optimisation methods have potential for lowering this cost, so these will likely
see increasing application in the future.
11.1 Treatment of CA
Our look at computational aeroelasticity over the next three chapters (CA) is
summarised in the following diagram:
11.1. TREATMENT OF CA 149
Computational Methods for Structural Dynamics

?
Computational Methods for Unsteady Flows
?
? ?
CA methods using Fully-coupled CA solution
CFD models methods
150CHAPTER 11. INTRODUCTION TO COMPUTATIONAL AEROELASTICITY
Chapter 12
Computational structural
dynamics
The methods we considered in chapter 8 for deriving representations of continu-

ous structures can be awkward to apply to complex configurations which feature
multiple structural and inertial components (figure 12.1). Many of the concepts
we have used, however, can still be applied on a more local level. In this chapter
we will the main concepts of the finite element method when used for structural
dynamic problems. In order to build on the concepts we have already developed,
we will limit the discussion to to semi-discrete finite-element methods obtained
by application of the Rayleigh-Ritz method and Lagrange’s equation. Time dis-
cretisation will be treated with standard time-marching methods. We note that
what we consider forms a mere subset of available finite-element techniques.
Figure 12.1:
151
152 CHAPTER 12. COMPUTATIONAL STRUCTURAL DYNAMICS
12.1 Rayleigh-Ritz finite-element formulation

The defining characteristic of finite-element methods is their locality. In other
words, the solution within each element is defined purely in terms of unknowns
associated with the element or its immediate neighbours. This feature simplifies
integrations since only one element contribution need be considered at a time.
It also results in large but sparse matrices, for which a number of powerful
algebraic solution techniques have been developed.
An example finite-element method can be obtained by considering the Rayleigh-
Ritz method of chapter 8, but with “assumed modes” (referred to from now on
as shape functions) defined only in the immediate vicinity of an element 12.2.
This allows large local deviations in the solution to be represented accurately
without affecting the quality of the solution interpolation in the remaining re-
gions.
x
φ
φi
y
node i
Figure 12.2:
We can summarise the discretisation process using the following steps:

1. Discretize the domain into a number of elements
2. Choose the shape functions which will represent the local solution (struc-
tural deflection)
3. Find the contributions of an element to Lagrange’s equation expressed in
terms of the nodal unknowns
4. Transform the local contributions of each element to the global coordinate
system, and assemble the global mass and stiffness matrix
5. Apply the boundary conditions and generalized forces
This procedure will result a system of equations for the unknown nodal dis-
placements, u, of the form:
[M ]{ü} + [K]{u} = {f (t)} (12.1)
An illustrative example is given in the next section.

12.2. EXAMPLE: FORMULATION FOR AN AXIALLY-LOADED ROD153
12.2 Example: Formulation for an axially-loaded

rod
Steps 1 and 2: Introduce elements and shape functions

The bar is first divided into elements e, e + 1 etc. We will assume the solution
will be interpolated within element e by two “tent-like” linear shape functions,
one centered at node i on the left side of the element and one centered at node
i + i on the right.
φi
e−2 e−1 e+1
element e
ui ui+1 Nodal Displacements
φi = 1 − x/L φi+1 = x/L
x Shape Functions
0 L
Βi+1 = 1/L
x Derivatives of
0 Shape Functions
Βi = − 1/L
Figure 12.3:
Step 3: Compute the Element contribution to Lagrange’s equations

Assume that ρ (density) A (cross-sectional area) E (modulus of elasticity) con-
stant. The kinetic energy is then defined as:
L
1 ρA
Z Z
2
T = ρu̇ dV = (u̇i φi + u̇i+1 φi+1 )2 dx (12.2)
2 v 2 0
Recalling the first term from Lagrange’s equation, it is convenient to define an

element mass matrix as:
Z L
M eij = ρA φi φj dx ; i = i, i + 1; j = i, i + 1 (12.3)
0
which for this example is:

ρAL 2 1
Me = (symmetric) (12.4)
6 1 2
We can then define the strain energy as:

L
1 EA
Z Z
2
U= Eu′ dV = (ui Bi + ui+1 Bi+1 )2 dx (12.5)
2 v 2 0
and define an element stiffness matrix to be:

Z L
Keij = EA Bi Bj dx ; i = i, i + 1; j = i, i + 1 (12.6)
0
which for this example is:

EA 1 −1
Ke = (symmetric) (12.7)
L −1 1
Step 4: Assemble the global mass and stiffness matrix

Each equation in the global system represents a Lagrange equation for the gen-
eralised displacement ui φi . Since the φi are chosen with local support, the ith
equation can be obtained by simply adding (Assembling) the contributions of
the elements to the left and right of node i. This can be done by ”overlapping”
the element mass and stiffness matrices as shown in figure 12.4.
M = ρAL
6
Mei
4 1
1 4 1
1 4
Mei+1
Figure 12.4: Assembly of the global mass matrix
The final result is a produces a system of the form:
[M ]{ü} + [K]{u} = {F } (12.8)

12.3. SOLUTION OF THE EQUATIONS OF MOTION 155
where the ith row of the system expresses Lagrange’s equation for the generalised
coordinate ui :

d ∂T ∂T ∂U
− + = Qi (12.9)
dt ∂ u̇i ∂ui ∂ui
Step 5: Apply the B.C’s and generalised forces

Constraints such as u(0) = 0 can be introduced to eliminate variables from the
system. Due to the form chosen for φi , the generalised forces, Qi , are simply
the the external forces applied at each node.
12.3 Solution of the equations of motion

Once the basic semi-discrete system has been defined, we must define a time
integration method to advance the system in time. There are a large number
of available methods, each of which have different accuracy and stability prop-
erties. When performing computations with possible physical instabilities, it is
important to understand how introducing a discrete time march modifies the
solution in time. In a flutter simulation, for example, it can be difficult to tell
if the instability arises due to numerical or physical mode amplification.
It is therefore useful to gain insight into how such methods behave, by con-
sidering a simple single degree of freedom model systems. Consider the following
equation of motion:
mü + ku = 0 (12.10)
This can also be written as a system of two first-order ODE’s:

d u 0 1 u
= (12.11)
dt u̇ −k/m 0 u̇
The
p solution of this system is a constant-amplitude oscillation with frequency
k/m.
Now assume we will advance the system in time using the Explicit Euler
integration technique, which is defined by a simple first-order Taylor series ex-
pansion in time. In the following, we will assume that the solution is advanced
with a constant time-step ∆t. We denote the current time level by n so that
tn = n∆t. Considering a step from n to n + 1 with the Explicit Euler technique,
the discrete version of our system can be written:
un+1 = un + ∆tu̇n (12.12)

u̇n+1 k n
= u̇n − ∆t m u (12.13)
or:
n+1 n
u u 1 ∆t
= [C] where [C] = (12.14)
u̇ u̇ −∆t(k/m) 1
Using [X], the matrix of right eigenvectors of [C], we can also express this system
in “modal coordinates” w1 , w2 as:
n+1 n+1 n
w1 u u
= [X]−1 = [X]−1 [C][X][X]−1 (12.15)
w2 u̇ u̇
n
σ1 0 w1
= (12.16)
0 σ2 w2
where σ1 , σ2 are the eigenvalues of [C]. The transformation to diagonal form

using the matrix of eigenvalues is known as the similarity transform. For the
Euler explicit case:
p
σ1,2 = 1 ± ∆t −k/m (12.17)
n+1 n
Since w1,2 = (σ1,2 )w1,2 , and | σ1,2 |> 1 we can conclude that the numerically
integrated response is growing in time. This is in contrast to the original semi-
discrete one DOF system, which has constant amplitude and frequency. It can
also be seen that The error is first-order in ∆t, consistent with the order of the
truncation error in the Taylor series.
A similar analysis of the implicit Euler method produces:
p
1 ± ∆t −k/m
σ1,2 = (12.18)
1 + ∆t2 k/m
This method is also first-order accurate, but in this case, | σ1,2 |< 1 so that
the amplitude is decreased each time step. One has to be aware of this type of
error during a computation, as this damping introduced by the time march can
stabilise an otherwise unstable mode.
In general, all numerical integration methods produce either amplitude er-
rors (signal growth or dissipation) or phase errors (time leads or lags) or both,
depending on their design. We may express the amplitude error per time step
(local amplitude error) as:
eamp = 1− | σ | (12.19)
while the local phase error is:
epse = ω∆t − arg(σ) (12.20)

p
where for the system considered, ω = k/m. If the method is consistent, eamp
and epse will reduce with (∆t)m , where m is the order of the method.
Although the phase and amplitude errors of the Euler explicit and Euler
implicit methods are similar, in the case of the explicit method, the amplitudes
are growing with each time step. Repeated application of the explicit method
by taking multiple time steps will thus result in an exponential growth in the
amplitude. This phenomena is known as (linear) numerical instability. For the
Euler implicit method, the solution will remain bounded after multiple time
12.4. NUMERICAL INTEGRATION MULTI-DOF SYSTEMS 157
steps (with an amplitude tending to zero). It is said to be (linearly) stable for

this problem.
In general, numerical time-integration techniques applied to either analytic
or discretised systems with produce solutions which are either linearly stable,
unstable, or conditionally stable. In the latter case, there is normally a maxi-
mum time step with which the system can be advanced and remain stable. For
structural problems, an implicit method is usually required to ensure uncondi-
tional linear stability.
12.4 Numerical Integration multi-DOF systems

Now consider a system such as what would be obtained by application of the
finite-element method described in the previous section:
[M ]{ü} + [K]{u} = 0 (12.21)
for such a case it is clear that the eigenvalues which characterise the be-
haviour of the numerically time-integrated solution will be functions of [M ] and
[K]. Therefore, the type of spatial discretisation used also effects the result-
ing amplitude errors (stability) and phase errors of the solution. Multiple-DOF
systems produce multiple modes, each with an associated oscillation period.
If the oscillation periods are close to each other, then conditionally-stable
explicit methods may be preferred, since they require far less work per time
step then implicit methods. If the oscillation periods are far apart in magnitude
(i.e. the matrix resulting from the discretisation is stiff), then unconditionally
stable implicit methods are generally preferred. This is because typically only
the lower frequency modes are of interest, and implicit methods permit the use
of large time steps. In contrast, a conditionally-stable explicit method would be
forced to march at the time-step required to prevent instability of the highest-
frequency mode.
12.5 Newmark method

The Newmark time integration method commonly used in structural dynamic
applications. For the system:
[M ]{ü} + [C]{u̇} + [K]{u} = {f (t)} (12.22)
An implicit method can be defined based on the satisfaction of equilibrium at

tn+1 :
[M ]{ün+1 } + [C]{u̇n+1 } + [K]{un+1 } = {f (tn+1 )} (12.23)

and two Taylor series expansions in terms of the anticipated value of acceleration
at the next time step:
∆t2
un+1 = un + ∆tu̇n + [(1 − 2β)ün + 2β ün+1 ] (12.24)
2
u̇n+1 = u̇n + ∆t[(1 − γ)ün + γ ün+1 ] (12.25)
The Newmark method is unconditionally stable (no amplitude growth) when

2β ≥ γ ≥ 21 . It is similar to trapezoidal integration when β = 41 , and γ = 12 .
12.6 Efficient integration using modal superpo-

sition
If you will repeatedly be advancing a structural discretisation in time, it is
often worthwhile to first solve for the natural modes of the structure, and then
construct a decoupled dynamic system based on the natural modes.
To solve for the natural modes, we assume a harmonic solution: {u} =
{û}eiωt , Substituting into the equations of motion leads to:
([K] − ω 2 [M ]){û} = 0 (12.26)
This generalised eigenvalue problem can be solved for the eigenvalues, ω 2 , and
the eigenvectors (or natural modes) Φi . The latter are normalised so that
[Φ]t [M ][Φ] = I, where [Φ] is the matrix of eigenvectors (orthonormal modes).
We can then use orthogonal properties of the modes to decouple the dynamic
system, by noting the matrix of orthonormal modes, [Φ] has the following prop-
erty:
 2 
ω1 . . .
 . ω22 . . 
[Φ]t [K][Φ] = [Λ] =  (12.27)
 
. . 
 . . . . 
2
. . . ωN
(12.28)
while by definition:
[Φ]t [M ][Φ] = [I] (12.29)
so that if we define generalized coordinates {q} by:
{u} = [Φ]{q} (12.30)
the equations of motion can be expressed:
[I]{q̈} + [Λ]{q} = [Φ]t {f (t)} (12.31)

12.6. EFFICIENT INTEGRATION USING MODAL SUPERPOSITION 159
which is a decoupled system in modal coordinates.

A damping term [C]{u̇} may also be included in the system provided that it is
of the form:
[C] = α[M ] + β[K] (12.32)
Overall, using modal superposition can substantially reduce the cost of nu-
merical integration. Furthermore, we can now choose to retain only the low-
frequency modes. This approach to reduced-order modelling procedure allows us
to perform compact computations which retain the most important behaviours
present in a complex geometry.
p12.7.24 Derive the element mass and stiffness matrices for a bar element
loaded in torsion using linear shape functions.
p12.7.25 Compute the torsional normal modes of the cantilever bar shown in
the figure below. Use two finite elements with linear shape functions.
!!!!!!!!!!!
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!! Ιθ(y) = Constant
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
GJ(y) = Constant
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
y=0
!!!!!!!!!!!
::::
!!!!!!!!!!!
::::
θ(0) = 0
!!!!!!!!!!! θ(y)
::::
!!!!!!!!!!!
!!!
!!!
y=1!!!
θ’(1) = 0
p12.7.26 Compute the forced response θ(y, t), of the bar in question 2 to a
torque T = Asin(ωt) applied at y = l. Use the results of question 2 and the
modal superposition approach.
p12.7.27 Consider a single DOF system given by:
mü + ku = 0
and a trapezoidal direct integration method to advance this system in time

based on the equations:
un+1 = un + ∆t
2 (u̇
n+1
+ u̇n )
u̇n+1 = u̇n − ∆tk
2m (u
n+1
+ un )
Compute the phase and amplitude errors produced by this method per time
step when k/m = 1 and ∆t = 0.1.
Chapter 13
Computational Methods for

Unsteady Flows
This chapter will provide an overview of computational methods currently used

to predict unsteady aerodynamic forces for aeroelastic analyses. We will begin
by considering methods for the linearised potential equation. These are still
in wide use for preliminary design, as they can give a first estimation of the
critical points for a configuration, which can later be examined in detail using
more sophisticated methods. The typical level of fidelity achieved for classical
flutter problems is illustrated in figure 13.1, which shows a non-dimensional
flutter speed versus Mach number. Comparison of linear and non-linear results
Figure 13.1: Flutter speeds predicted with different aerodynamic models (from
Schuster et al, J. of Aircraft (40) 5 , 2003)
also gives an indication of the sensitivity of the design to non-linear effects.

In the transonic region, non-linear prediction methods are typically required
in order to predict the flow with acceptable accuracy. In addition to providing
161
162CHAPTER 13. COMPUTATIONAL METHODS FOR UNSTEADY FLOWS
better predictions of classical flutter modes, non-linear methods can also be used
to consider responses which are purely dependent on phenomena which cannot
be computed with linear methods. Examples include aileron buzz and stall
flutter, although the truly accurate computation of these phenomena remains a
considerable challenge.
In our description of non-linear methods, we will first consider the small-
disturbance potential equation, and then move on to the Euler and Navier-
Stokes equations, for which we will also discuss the analysis of numerical errors.
We will not discuss methods for the full-potential equations, as these are similar
in form to methods for the Euler equations. Figure 13.2 presents an overview
of the models, concepts, and methods considered in this chapter.
Aerodynamic Important Example

Model Concepts Methods
Superposition Kernel Function

Linearised
Potential Wake / time delay Panel
treatment
Small-Distubance Linearised
Finite Difference
Potential moving-wall B.C.s
Moving meshes Finite Difference

and the DGCL
Euler / Finite Volume
Navier-Stokes Methods for error
estimation Finite Element
Figure 13.2: Models, concepts and example methods considered in chapter 13

13.1. LINEARISED POTENTIAL EQUATION 163
13.1 Linearised potential equation

Methods for the linearised potential equation are not only attractive for their
computational speed, but also for the ease with which they can be validated
using exact solutions. In industry, these methods remain popular due to the
wealth of experience built up with their use. Linearised methods are also good
candidates for inverse design procedures, although the sensitive nature of the
flutter problem still makes the development of robust inverse algorithms a chal-
lenge.
Many solution methods for linear equations make use of superposition. For
the linearised potential equation, surface-singularity approaches are particularly
popular. These can be based on either the velocity- or acceleration-potential
form of the equations. The velocity-potential form of the equations was intro-
duced in chapter 3, and can be written:
2
1 ∂2φ ∂2φ

2 2 ∂ φ
∇ φ− M∞ − 2 + 2U∞ =0 (13.1)
∂x2 a∞ ∂t2 ∂x∂t
(13.2)

∂φ ∂φ
p − p∞ = −ρ + U∞ (13.3)
∂t ∂x
The acceleration-potential form can be obtained by defining ψ, the acceleration

potential, as:
∂φ ∂φ
ψ= + U∞ (13.4)
∂t ∂x
For a harmonic solution ψ = ψ̂(x, y, z)eiωt , the linearised potential equations

can be expressed in terms of acceleration potential as:
∂ 2 ψ̂ M∞ ∂ ψ̂ ω2
2
∇2 ψ̂ − M∞ 2
− 2iω + 2 ψ̂ = 0 (13.5)
∂x a∞ ∂x a∞
p − p∞ = −ρψ̂ (13.6)
The first equation is similar to the velocity-potential form for harmonic motions.
The second equation is now much more convenient for the application of wake
boundary conditions. In this case the ∆p = 0 wake condition can be simply
expressed as ψ = 0 on the wake surface.
13.1.1 Kernel-function method

Consider the isolated wing geometry of figure 13.3. If we assume the solution is
harmonic, it is possible to express the vertical velocity at all points on the wing
as:
w(x, y) = ŵ(x, y)eiωt (13.7)

z y,η
za
x,ξ
S
Figure 13.3: Basic wing geometry
where:
Z
ŵ(x, y) = ∆ψ̂(ξ, η)K(x − ξ, y − η, k, M∞ )dξdη (13.8)
S
and K, known as the kernel function, is a (rather complicated) analytic expres-

sion which can be derived for thin wings with planar wakes. In order to solve
the problem we may choose:
n X
X m
∆ψ̂(x, y) = âij fi (x)gj (y) (13.9)
i=1 j=1
where fi (x) and gj (y) are assumed continuous functions which satisfy the pres-
f(x)
g(y)
Figure 13.4: Assumed forms for the kernel function method
sure boundary conditions on the edges of the wing (singular at leading edge,
13.1. LINEARISED POTENTIAL EQUATION 165
zero at the tips and trailing edges). Example functions are illustrated in figure
13.4. The âij coefficients can be found by enforcing flow tangency condition at
n × m control points using the integral equation for w and the assumed form
for ∆ψ̂(x, y). This leads to a system of (n × m) equations for the unknown
coefficients. Note that the wake need not be included in the integration since
∆ψ̂wake = 0. This is because all of the assumed harmonic wake-lag effects are
already encapsulated in the definition of K.
13.1.2 Panel Methods

The kernel-function approach is extremely efficient, but it is restricted to sim-
ple, thin wings. Alternatively, panel methods can be used if more complex
geometries are to be considered. Panel methods can be formulated in the time
domain, using the velocity-potential form of the equations. Vorticity convection
in the wake is then treated explicitly. The general steps in the procedure are
Figure 13.5: Layout of a panel method
(see figure 13.5):

1. The surface of the wing or body and its associated wake is discretised
into a series of panels with surface singularity distributions of unknown
strength.
2. The potential induced at a point due a given panel is then expressed in
terms of the panel’s p unknown singularity-strength coefficients.
3. If there are N panels, a system of N ∗ p equations can be formed by
enforcing flow tangency at p control points on each panel. Since the PDE
is linear, the local normal velocity is simply the sum of the contributions
from all the panels.
4. The position of the wake is updated so that ∆pwake = 0
For incompressible flows, or harmonic compressible flows, the potential induced

by a panel can be expressed simply. For arbitrary compressible flows, time de-
lays must be taken into account, which complicates implementation. An equally
important problem is the determination of the wake position, as it can experi-
ence significant distortions in regions of high vorticity, such as near wing tips
(figure 13.6). It is possible to find the wake position using non-linear iteration
procedures, but these are generally unreliable if the wing is significantly loaded.
A robust alternative is to use a fixed reference wake, and solve for small normal
displacements from the wake as a function of time. This can be quite accurate
at the low reduced frequencies which are often crucial in aircraft configuration
analysis.
Figure 13.6: Concentration of wake vorticity
The singular behaviour of the potential solution in the immediate vicinity of

the panels means that wake-body interactions, such as produced by propellers
mounted ahead of wings, are difficult to treat accurately. Robust analysis of
such problems requires the solution of the Euler or Navier-Stokes equations,
which can capture arbitrary vortical phenomena within their solutions.
13.2. SMALL-DISTURBANCE POTENTIAL EQUATION 167
13.2 Small-disturbance potential equation

The earliest methods which could predict transonic flutter where based on the
small-disturbance potential equation, the simplest model which can include the
effects of shock motion. Since many companies have extensive experience with
its use, the small-disturbance potential equation (also known as transonic small
disturbance (TSD) equation), is still often applied for preliminary transonic
analysis. The small-disturbance potential equations can be written:
2
∂2φ ∂2φ

2 2 ∂φ ∂ φ
1 − M∞ − (γ + 1)M∞ b 2
+ 2 + 2 (13.10)
∂x ∂x ∂y ∂z
∂2φ ∂2φ

1
− 2 2U∞ + 2 =0
a∞ ∂x∂t ∂t

∂φ ∂φ
p − p∞ = −ρ + U∞ (13.11)
∂t ∂x
Since the main equation is non-linear, superposition cannot be used. This means
that the entire flow domain must be discretised (figure 13.7). As we are con-
sidering small disturbances, however, it is customary to linearise the boundary
conditions. The latter is an extremely extremely convenient simplification, as
then the application of flow tangency on the airfoil surface does not require dis-
tortion of the grid in order to account for structural deformations. A popular
Far−Field Boundary
j,k+1
j,k j+1,k
j−1,k
j,k−1
Airfoil Wake
Figure 13.7: Discretisation for the small-disturbance potential equation
discretisation option is to apply finite differences simultaneously in space and

time, where derivatives are replaced by approximate expressions derived from
Taylor series expansions such as:
φj−1,k − 2φj,k + φj+1,k
φxxij = + O(∆x2 ) (13.12)
∆x2
Cartesian meshes can be used for simplicity, although transformed meshes are
sometimes used to place the far-field boundary at infinity. Once discretised, the
following boundary conditions are applied:
• Linearised flow tangency on the airfoil surface

(φ±
z is prescribed by the airfoil motion)
• Linearised change in pressure across the wake surface is zero
• The perturbation is zero at the far field
An implicit discretisation is then obtained in terms of the unknowns at the next

time level. This is typically solved using a sub-iteration techniques to remove
linearisation errors in time.
13.3 Euler and Navier-Stokes equations

The 1980’s saw the introduction of the first reliable numerical methods for
the compressible Euler and Navier-Stokes equations. These have only recently
come into common use for aeroelastic analysis, however, due to their relatively
high cost when used for low-frequency unsteady simulations. The Euler equa-
tions represent a considerable step up from potential-based methods (including
full potential methods), as they admit vorticity directly into their solutions,
and can therefore capture complex wake interactions occurring between aerody-
namic components. In addition they automatically include the effects of time
delays, and include such features such as non-isentropic shock waves, non-linear
expansion waves, and contact surfaces. The Euler equations are still limited, of
Figure 13.8: An Euler mesh for unsteady aerodynamic computations
course, as they do not include the effects of viscosity. In theory this means that
13.3. EULER AND NAVIER-STOKES EQUATIONS 169
they also do not enforce the Kutta condition, but in practice this is enforced
automatically due to the necessary presence of finite amounts of numerical dis-
sipation.
When other viscous effects must be considered, the Reynolds-averaged Navier-
Stokes equations are normally used. This raises an important issue concerning
the fidelity of turbulence models used to predict unsteady flows. A lack of un-
steady data makes such models difficult to calibrate, and often models calibrated
using steady data are used. Another problem occurs if the turbulent fluctuations
of the largest eddies approaches the time scales of the structural vibration. In
this case the assumptions underlying the Reynolds-averaged approach are com-
promised, and one formally should switch to large-eddy simulation. Due to
resolution requirements, however accurate LES is currently not possible at the
Reynolds numbers of interest to aircraft applications.
In contrast to the small-disturbance equations, in most cases it is inap-
propriate to linearise the boundary conditions for the Euler and Navier-Stokes
equations. This adds considerable complication to the solution procedure, as the
discretisation must be general enough to treat deforming domains (figure 13.9).
Aside from the complication of implementation, it has been demonstrated that
Figure 13.9: A mesh with a deforming boundary
algorithms which discretely enforce conservation in the unsteady sense usually

deliver superior accuracy, making this an additional criteria in method design.
As far as wall boundary conditions are concerned, for the Euler equations
it is sufficient to apply the flow-tangency condition at the wall, while for the
Navier-Stokes equations no-slip conditions are required. Typically this is done
using methods similar to those used for steady flows, but accounting for wall
movement. At far field boundaries, however, the enforcement of boundary condi-
tions in unsteady flows can be far more complicated than for steady problems.
Integral or Partial Differential

Flow Equations
Mesh Generation
Fully−Discrete Approach
Semi−Discrete Approach Spatial Discretization
Semi−Discrete System Space-Time

(Ordinary Differential Discretization
Equations)
Time Integration
Fully−Discrete System
(Algebraic Equations)
Figure 13.10: Options for discretisation
Due to computational cost the domain of the problem must remain limited,
which means that important influences of physical phenomena (such as con-
vecting wake vorticity) can be lost as they leave the domain. This contrasts
with steady problems, which typically have uniform far-field conditions. The
latter can be be represented by relatively simple procedures, provided solution-
error transients can be properly expelled.
13.4 Discretisation of the Euler and Navier Stokes

equations
In general, either semi- or fully-discrete approaches can be applied using either
the Finite-difference, Finite-volume, or Finite-element method. An overview is
given in figure 13.10.
13.4.1 Semi-discrete approach

In the semi-discrete approach, a finite-difference, volume or element method is
first applied to the spatial part of the partial differential equations (PDE). This
results in a system of ordinary differential equations (ODE) for the unknown
13.5. DEALING WITH MOVING MESHES 171
discrete solution values as a function of time. These are integrated using either
explicit or implicit time-integration methods.
For unsteady problems the solution accuracy is a function of both the spatial
discretisation and time integration method. This is different from when steady-
flow problems are considered, for which the accuracy of the time integration
method normally used to solve the non-linear problem is unimportant.
In general, explicit time-integration methods (e.g. MacCormack, Runge-
Kutta) are restricted for stability reasons to small time steps, on the order of
the time necessary for an acoustic wave to traverse the smallest cell. Implicit
time-integration methods (e.g. Trapezoidal) are stable for much larger time
steps but require more work per time step than explicit methods. Implicit time-
integration methods are often preferred, as the time scales of fluid-structure
interactions are usually much greater than the acoustic scale of the smallest
mesh cell, particularly for viscous calculations.
13.4.2 Fully-discrete approach

In the fully-discrete approach a finite-difference, volume or element method is
applied simultaneously to the spatial and temporal parts of the PDE. This
leads directly to a (non-linear) algebraic system for the unknown nodal values.
Fully-discrete discretisations (also known as space-time discretisations) can be
advantageous from an implementation point of view, since they naturally in-
clude mesh motion. This is illustrated in figure 13.11, which shows space-time
elements for a one-dimensional domain moving to the right as time advances.
A moving 1D mesh in space−time
Figure 13.11: Fully-discrete approach
13.5 Dealing with moving meshes

Typically the equations of fluid mechanics are in Eulerian form, that is they are
derived for a control volume which is fixed in space. This is inconvenient if the
boundary of the domain is changing due to, for example, the movement of an
n
V
V : Control volume
S: Surface of control volume
x
n: Local surface normal vector
x: Local surface velocity vector
dS
Figure 13.12: A moving control volume
elastic wing. An alternative is to use a Lagrangian formulation, where the mesh

follows the fluid particles (and thus the fluid boundary). This is impractical
for most flows due to the large distortions associated with tracking fast-moving
particles. In fact, we need to be able to chose the mesh geometry arbitrarily,
in order to control mesh quality as the boundary deforms. This leads to the
concept of the Arbitrary Lagrangian-Eulerian (ALE) formulation. This is most
easily expressed by directly applying the Reynolds Transport theorem for a
moving control volume to the Euler or Navier-Stokes equations. Referring to
figure 13.13, this results in:
d
Z Z
W dV + (F~ i − F~ v − ~ẋW ) · ~n dS = 0 (13.13)
dt
V S
Here W represents the flow state vector:

 
ρ
W =  ρ~u  (13.14)
ρE
which contains components of density, specific momentum, and specific total

energy. The flux vector is defined by F~ i − F~ v , where F~ i and F~ v are the inviscid
and viscous contributions. Finally, ~ẋ denotes the local surface velocity of the
control volume. Note that (13.13) reverts to the Lagrangian description for
~ẋ = ~u and the Eulerian for ~ẋ = 0. In practice we will instead associate ~ẋ with
the movement of the mesh, as it changes position to represent aeroelastically-
deformed geometries.
13.5.1 The Geometric Conservation Law (GCL)

Consider the ALE formulation (13.13) for a constant W = W ∗ :
d
Z Z
∗
W dV + (F~ i (W ∗ ) − ~ẋW ∗ ) · ~n dS = 0 (13.15)
dt
V S
13.6. EXAMPLE: A SEMI-DISCRETE FINITE-DIFFERENCE METHOD173
In this case, the viscousRfluxes F~ v (W ∗ ) = 0, and the inviscid fluxes are constant.
Since for a closed cell, ~n dS = 0:
S
d
Z Z
dV = ẋ · ~n dS GCL (Integral form) (13.16)
dt
V S
This relation is called the Geometric Conservation Law (GCL), as it can be

interpreted physically as “The change in V must equal the volume swept out
by its surface”.
13.5.2 The Discrete Geometric Conservation Law (DGCL)

In the design of steady solution methods, it has long been known that methods
which are discretely conservative (those which conserve the total mass, momen-
tum, and energy exactly), tend to produce more accurate results. This becomes
crucial when discontinuous solutions are present, as their description is tied to
the integral rather than the differential form of the equations. In recent years,
it has also become clear that discrete enforcement of the geometric conserva-
tion law is desirable for methods which use moving meshes. Although the idea
of enforcing the DGCL dates back to Thomas and Lombard (AIAA J. 1979
18:1030-1037), comprehensive analyses of its advantages were first made in a se-
ries of results published by Farhat et al between 2000 and 2001. First, Guillard
and Farhat (Comput. Methods Appl. Mech Engrg. 2000 190:1496-1482) showed
that for a p-order accurate scheme, enforcing the p-discrete DGCL will ensure at
least 1st-order accuracy on a moving mesh. Then, Farhat and Geuzaine (AIAA
paper 2001-2607) showed that the second-order three-point backward-difference
time-integration scheme maintains its accuracy on a moving mesh if the DGCL
is enforced. Finally, in Farhat, Geuzaine and Grandmont (J. Comp. Physics
2001 174:669-694) it is proposed that the DGCL is a necessary and sufficient
condition for preserving the nonlinear stability of fixed-mesh methods when
generalised to moving-mesh problems.
In the following sections, some example discretisations for moving meshes
will be presented. The details of enforcement of the DGCL will also be discussed.
Note that the details of the enforcement of the DGCL are considered additional
material for this course, and are not required knowledge for the exam.
13.6 Example: A semi-discrete finite-difference

method
As the application of multidimensional finite-difference relations is inconvenient,
normally finite-difference techniques are implemented using a curvilinear trans-
form, which maps the mesh onto a Cartesian reference mesh. When moving
meshes are used, this transformation must be time-varying. The finite-difference
ξ=const
y η
η=const
x ξ
Figure 13.13: Curvilinear coordinate transform
discretisation is then applied to the transformed equations on the fixed reference

mesh. Consider the divergence form of the Euler equations:
dW ~ · F~ = 0
+∇ or in 2D: Wt + fx + gy = 0 (13.17)
dt
where fx are the flux components gy in two dimensions. The transformed equa-
tions can be obtained using the chain rule, and are thus expressed in terms of
the metrics of the transformation (ξt , ξx , ξy , ηt etc.):
Wt + Wξ ξt + Wη ηt + fξ ξx + fη ηx + gξ ξy + gη ηy = 0 (13.18)
(Note that dissipative terms which vanish in the limit of zero mesh size are often
added to (13.18) to improve the stability of the numerical solution procedure).
To evaluate the metrics of the transformation, note:
    
dt 1 ξt ηt dτ
 dx  =  0 ξx ηx   dξ  (13.19)
dy 0 ξy ηy dη
and:
    
dτ 1 xτ yτ dt
 dξ  =  0 xξ yξ   dx  (13.20)
dη 0 xη yη dy
By comparison, the second matrix of metrics is equivalent to the inverse of the

first. This leads to:
   
1 ξt ηt J (yτ xη − xτ yη ) (xτ yξ − yτ xξ )
 0 ξx ηx  = J −1  0 yη −yξ  (13.21)
0 ξy ηy 0 −xη xξ
where J = (xξ yη − yξ xη ) is the Jacobian of the transformation. The expressions

for xξ , yτ , yη etc. can easily be evaluated on the fixed (Cartesian) reference
mesh.
13.7. EXAMPLE: A SEMI-DISCRETE FINITE-VOLUME METHOD 175
F
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
S
CCCCCCCCCCCCCCC
Wi
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
Figure 13.14: A single cell of a finite-volume discretisation
13.6.1 Enforcing the DGCL (optional section)

In order to satisfy the DGCL, the Jacobian is updated separately by integrating
a transformed version of the geometric conservation law:
Jt + (Jξt )ξ + (Jηt )η = 0 (13.22)
Ensuring the time and space discretisation for this relation is identical to that
of the flow equations enforces the DGCL (see Thomas and Lombard AIAA J.
1979 18:1030-1037 for more details).
13.7 Example: A semi-discrete finite-volume

method
A moving-mesh finite-volume method can be derived directly from (13.13), that
is the conservation laws written for an arbitrarily-moving control volume:
d
Z Z
W dV + (F~ − ~ẋW ) · ~n dS = 0 (13.23)
dt
V S
By defining cell-averaged solution values:

1
Z
Wi = W dV (13.24)
Vi
Vi
The conservation laws can applied directly to each cell (figure 13.14), resulting
in a system of ODEs for the time evolution of the average values:
M
d
(F~m − Wm~ẋm ) · S
~m
X
(Vi Wi ) = − (13.25)
dt m=1
where S~m = ~nm · | Sm | These can be integrated using standard time-integration

techniques, although some care must be taken to satisfy the DGCL.
n+1
n+1 Xb
Xa
t
Xa(t) Xb(t)
x n
n Xb
Xa
Figure 13.15: Movement of the side of a control volume in space-time
13.7.1 Enforcing the DGCL (optional section)

The finite-volume version of the DGCL can be expressed:
Z tn+1 X
M
(Vin+1 n
− Vi ) = ~m = 0
x~˙m · S (13.26)
tn m=1
Normally it is chosen to compute Vin+1 and Vin exactly from the given mesh
deformation data, and then formulate the right-hand side in a way which sat-
isfies the DGCL. How this is done depends on both the spatial and temporal
discretisation. Consider a single (flat) side of a 2D cell as shown in figure 13.15,
combined with a first-order accurate time-integration scheme, and define:
Z tn+1 Z b
Iab = ẋ · ~n ds dt (13.27)
tn a
A linear parameterisation for the mesh movement can be made as:

x(t) = αxa (t) + (1 − α)xb (t) α ∈ [0, 1] (13.28)
ẋ(t) = αẋa (t) + (1 − α)ẋb (t) (13.29)
xa (t) = δ(t)xn+1
a + (1 − δ(t))xna δ ∈ [0, 1] (13.30)
xb (t) = δ(t)xn+1
b + (1 − δ(t))xnb (13.31)
Then, noting ~n is a function of time but not α, and ds = ldα:
Z tn+1 Z 1
Iab = (αẋa + (1 − α)ẋb ) ~nl dα dt (13.32)
tn 0
Z tn+1
1
= (ẋa + ẋb )H(xa − xb ) dt (13.33)
tn 2

0 −1
where H = converts the vector between the two endpoints to the
1 0
normal area, ~nl. Substituting our parameterisations, and noting dt = ∆tdδ,
1 1
Z
(xn+1 − xna ) + (xn+1 − xnb ) H

Iab = a b
2 0
δ(xn+1 − xn+1 ) + (1 − δ)(xna − xnb ) dδ

a b (13.34)
13.8. EXAMPLE: A SEMI-DISCRETE FINITE-ELEMENT METHOD(OPTIONAL SECTION)177
n+1 n+1
Xa Xb
Xa(t) Xb(t)
n
n Xb
Xa
Figure 13.16: Higher-order variation in time
n+1
Xb
n+1
y
n+1 Xc
x Xa
t Xa(t) Xc(t)
n n
Xa Xc
x
Figure 13.17: Three-dimensional DGCL
which has an integrand which is linear in δ. Therefore, Iab can be evaluated

exactly using the mid-point rule. The DGCL is thus satisfied in this case if the
1
configuration used for evaluating the flux integral is at xn+ 2 . A higher-order
time-integration scheme requires a higher-order parameterisation of the edge
displacement (figure 13.15). For quadratic variation, at least two integration
points are required in order to exactly evaluate Iab . As a consequence, the flux
integral must be computed on at least two mesh configurations.
If we instead consider the single side of a 3D tetrahedron in the first-
order case (figure 13.17), the integral Iabc contains a surface-normal expression
quadratic in δ, arising from the cross product of xab and xac . Since computing
Iabc exactly then requires two integration points, the discrete expression for the
flux integral must again be evaluated with two mesh geometries.
13.8 Example: A semi-discrete finite-element method

(optional section)
Semi-discrete finite-element methods are normally based on a modified ALE
form of the conservation laws, derived from considering a fixed reference volume
Vo , which is related to the moving volume V via a coordinate transform (xi , t) ↔
V x
Vo
Figure 13.18: Coordinate transform for semi-discrete FEM methods
(ξi , t) for which ∂V = J∂Vo (figure 13.18). The original ALE form then becomes:
d
Z Z
W J dVo + (F~ − ~ẋW ) · ~n dS = 0 (13.35)
dt
Vo S
d
Since Vo is fixed, we can bring dt inside the integral. By employing the diver-
gence theorem, the following weak form can be derived:

d(W J)
Z
w ~ ~
+ J∇ · (F − ẋW ) dVo = 0 (13.36)
dt
Vo
where w is the finite-element test function. Normally the second integral is

evaluated on V , leading to:
d(W J)
Z Z
w dVo + w∇ · (F~ − ~ẋW ) dV = 0 (13.37)
dt
Vo V
A semi-discrete ALE finite-element method can then be defined by considering

d(W J)
dt an unknown which is related to the flow solution via the time-marching
algorithm. Satisfaction of the DGCL using such methods requires an analysis
similar to that shown for finite-volume schemes in the previous section.
13.9 Example: A space-time finite-element method

(optional section)
Space-time (or fully-discrete) discretisations are advantageous as they provide
a natural setting for satisfying the DGCL on moving meshes. In fact finite-
element space-time discretisations automatically satisfy the DGCL for any order
of accuracy in space in time. In order to formulate a space time method, we
consider a fixed control volume in space-time, as shown in figure 13.19. If we
define a generalised divergence operator and flux vector:

˜ = ∂, ∂ , ∂ ;
∇ F̃ = (W, f, g, h) (13.38)
∂t ∂x ∂y
13.9. EXAMPLE: A SPACE-TIME FINITE-ELEMENT METHOD(OPTIONAL SECTION)179
t n
@@@@@@@@@
@@@@@@@@@
@@@@@@@@@
@@@@@@@@@
V
@@@@@@@@@
@@@@@@@@@
@@@@@@@@@
S
x
Figure 13.19: A fixed control volume in space-time coordinates
Q n+2
t
Q n+1
Qn
Q n−1
x
Ωn
Figure 13.20: The discrete space-time domain
where f, g, and h are the x, y and z components of the space flux vector, F~ , the
conservation laws can be written:
Z
˜ · F̃ dṼ = 0;
∇ (13.39)
Ṽ
or
Z
F̃ · ñ dS̃ = 0 where ñ = (nt , nx , ny ) (13.40)
S̃
To advance the solution, we note that due to causality, only one part of the
space-time domain needs to be considered at a time. The space-time domain
is therefore split into slabs Qn , n = 1..N with constant-time boundaries Ωn at
t = tn , as shown in figure 13.20. The solution can then be advanced in time
by solving the slabs sequentially. If the mesh topology does not change in time,
the number of unknowns in each slab will be comparable to those in standard
time-marching techniques.
If a doubling of unknowns can be afforded, one can also choose to make the
boundary between slabs discontinuous in time, as shown between the first and
Q
t y
t +∆t
t
R
Ωo
x
Figure 13.21: A two-dimensional space time mesh
second slabs of figure 13.20. The information from the previous slab is then
imposed as a weak boundary condition on Ωo , the boundary between Qn+1 and
Qn :
Z Z
˜ · F̃ dQ +
wh ∇ wh (uh+ − uh− ) = 0 (13.41)
Q Ωo
This gives enormous flexibility for the generation of moving meshes, as the mesh
in one time step can be completely unrelated to the one in the next time step
(the total number of elements, for example, can be arbitrarily changed from one
time step to another). In order to avoid evaluating the gradient of the flux,
(13.41) is usually expanded to:
Z Z
˜
−F̃ · ∇w dQ + wh F̃ · ñ dΓ
h
(13.42)
Q Γ
Z
+ wh (uh+ − uh− )dΩ = 0 (13.43)
Ωo
where Γ is the complete surface of Q. This expression maintains a uniform flow,

and thus satisfies the DGCL, provided:
Z Z
∇w dQ̃ − wh ñ dΓ = 0
h
(13.44)
Q̃ Γ
which in the absence of quadrature error, is evaluated exactly by the finite-

element procedure. Note that space-time discretisations are necessarily implicit
in time, although explicit sub-iteration techniques can be still be used as low-
memory solution options.
13.10 Discretisation Error Analysis

The large computational requirements of unsteady simulations means that they
are typically performed at the limit of computer resources, and often with
13.10. DISCRETISATION ERROR ANALYSIS 181
meshes which are coarser that what would be desirable. Consequently, the esti-
mation of discretisation errors due to finite space and time spacings is crucial,
especially when considering phenomena as sensitive as flutter.
13.10.1 Grid convergence studies

Space and time grid convergence studies are essential for verifying the quality of
unsteady computational solutions. These can also be used to provide direct error
estimates. In general, for a method which is nth-order accurate, the discrete
solution ud can be expressed in terms of the exact solution, ue as:
ud = ue + ghn + higher-order terms (13.45)
where g is a function of the gradients of the continuous solution, and h is the

step size (∆x or ∆t). This results in behaviour of the error ǫ =| ue − ud | similar
to that shown in figure 13.22. When h is small, the higher-order terms become
log(ε)
Asymptotic
Region
Non−Asymptotic
Region
n
1
log(∆x)
Figure 13.22: Standard behaviour of discretisation error
negligible and the slope on a log-log plot becomes n. The error is then said to
be in the asymptotic region.
Results from standard CFD calculations must be in the asymptotic region
to be reliable. This can depend on the definition of ǫ, as results for drag may
require a finer discretisation than those for lift. In order to have confidence in
your results, variations with both h = ∆x and h = ∆t must be checked.
In practice, one does not have ue . Therefore plots similar to figure 13.22 are
produced by comparing coarse-mesh solutions to those of the finest mesh in space
and time. A straight-forward way to estimate error can then be constructed by
considering solutions computed on a sequence of grids with different h.
Consider uh computed with element size h, and u2h computed with element
size 2h, for a method which is second-order accurate (n = 2). g can be eliminated
from the previous equation to obtain:

uh − u2h
ue − uh ≈ (13.46)
3
If this error is less than the required tolerance one can assume that the finest
numerical solution is adequate. The error estimate can also be used to correct
uh . This is known as Richardson extrapolation.
13.11 Phase and Amplitude Error Analysis

Finite-difference, finite-volume, and finite-element methods, and their associ-
ated time-integration schemes, all introduce errors in the representation of flow
features as they are convected through the domain. A truncation error analysis
of a method will yield its order of accuracy in space and time. It is also useful
to examine phase and amplitude errors, however, since methods of the same
order can differ significantly from each other in the representation of convective
phenomena. The Euler equations can be re-written in a form in which they
Exact
solution
eamp
t
Numerical
x1 x2 solution epse
Vortex travelling relative to mesh
Figure 13.23: Phase and amplitude errors
express:
• The convection of entropy with speed u
• The convection of vorticity with speed u
• The convection of acoustic perturbations with speeds u + a and u − a
where u is the local flow speed, and a is the local speed of sound. By analogy,
solutions for the linear convection equation:
∂u ∂u
+c =0 (13.47)
∂t ∂x
are often used to analyse and compare the performance of different numerical
methods for the Euler and Navier-Stokes equations.
13.11. PHASE AND AMPLITUDE ERROR ANALYSIS 183
m=6
m=5
∆x
m=4
m=1 m=2 m=M m=3
Normal 1D Domain m=2

m=M
m=1
Periodic 1D Domain
Figure 13.24: Normal and periodic 1D domains
If we consider a 1D periodic domain (figure 13.24), application of a discreti-

sation method to the linear convection equation will lead to a system of the
form:
~un+1 = [C]~un (13.48)
where ~un is the vector containing the unknown values of u at each node or
volume. This system can be diagonalised using the similarity transform:
~ n+1
w ≡ [X]−1 ~un+1 (13.49)
−1 −1 n
= [X] [C][X][X] ~u (13.50)
~n
= [σ]w (13.51)
where [X] is the matrix of eigenvectors of [C], and [σ] is the diagonal matrix of
eigenvalues of of [C]. For a periodic domain, [C] is circulant, and will always
have the same eigenvectors, independent of the discretisation method used:
Xjm = ei(2π(j−1)(m−1)/M) ; j = 1..M ; m = 1..M ; (13.52)
where j is the number of the eigenvector, m is the node index and M is the
total number of nodes.
The eigenvectors of the matrix [C] are a collection of sine waves of different
spatial frequencies, with the lowest being a constant, and the highest being a
signal which changes sign from node to node (figure 13.25) . In the exact solution
to the linear convection equation, each of these waves should be convected with
constant amplitude with a speed c. This leads to the following definitions of
phase and amplitude errors incurred per time step for each wave:
eamp = 1− | σj | (13.53)
c∆t
epse = βj − arg(σj ) (13.54)
∆x
1
0.8 j=2
j=4
0.6
0.4
Amplitude
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
1 2 3 4 5 6 7 8 9 10
Node index (M=10)
Figure 13.25: Eigenvectors for a 10 node mesh
Where in this case, βj = 2π(j − 1)/M is the frequency of the eigenvector (wave)
associated with the eigenvalue σj . For the linear convection equation, the phase
and amplitude errors produced by a numerical method are a function of both
the wave frequency and the Courant number = c∆t/∆x, where ∆t is the time
step and ∆x is the spacing between nodes. In figures 13.26 to 13.28, phase
and amplitude errors for three different numerical methods are given, along
with numerical solutions for the convection of a block wave. In the first
case, the combination of a second-order accurate finite-difference method with
first-order accurate implicit time integration produces a method which is high
in phase error at high frequencies (beta) but also has very little damping at
high frequencies. The result is a sawtooth effect due to the incorrect speeds of
the high-frequency portions of the solution. In the second case, a second-order
accurate explicit time integration method is used, with coefficients designed
to produced better high-frequency damping, resulting in an improved solution.
Finally, results from a third-order accurate space-time method are shown, to
which high-frequency damping has been added. The remaining oscillations near
the discontinuous part of the solution can be further controlled using operators
tuned for shock capturing.
The above observations will also apply to more general flow features con-
vected by similar schemes for the Euler and Navier-Stokes equations. If you
think of a flow feature (such as a vortex) as being made up from waves of dif-
ferent frequencies, you can anticipate that its high-frequency components will
normally exhibit large phase errors. For a well-designed scheme, these compo-
nents will also be strongly damped. A side effect will be the smoothing out of
the flow feature as it is convected. Marching with higher Courant numbers will
tend to produce greater errors. It can also be seen that higher-order schemes
can be really worthwhile, provided that their cost is not excessive.
13.11. PHASE AND AMPLITUDE ERROR ANALYSIS 185
1
Courant=1.5
Courant=1.0
Courant=0.5
0.8
0.6
e_amp
0.4 1.5
Initial
Exact
0.2 Numerical
1
0
0 0.5 1 1.5 2 2.5 3
0.5
u
beta
3 Courant=1.5
Courant=1.0 0
Courant=0.5
2.5
2 -0.5
0 0.2 0.4 0.6 0.8 1
e_pse
1.5 x
0.5
0
0 0.5 1 1.5 2 2.5 3
beta
Figure 13.26: Second-order accurate centred finite difference method with first-
order accurate implicit time integration: 100 nodes, 80 time steps.
13.11.1 Adjoint-based error estimation
Over the last few years methods for formally estimating the magnitude of local
error contributions to quantities of interest, such as lift and moment, have de-
veloped rapidly. These involve the solution of dual or adjoint problems which
express the sensitivity of the error in the desired quantity to local values of the
residual. There is great potential for such methods in the field of aeroelasticity,
but their exists two basic complicating issues. The first of these is that the dual
problem evolves in the reverse time direction, making its simultaneous compu-
tation with the main (primal) problem a bit awkward. The second more serious
problem arises from the moving interface between the flow and the structure,
which poses serious difficulties for the implementation of the dual problem. The
reader is referred to the recent TU Delft thesis of Kris van der Zee for more
information and successful strategies.
1
Courant=1.5
Courant=1.0
Courant=0.5
0.8
0.6
e_amp
0.4 1.5
Initial
Exact
Numerical
0.2
1
0
0 0.5 1 1.5 2 2.5 3
0.5
u
beta
3 Courant=1.5
Courant=1.0
Courant=0.5 0
2.5
2
-0.5
e_pse
0 0.2 0.4 0.6 0.8 1

1.5
x
1
0.5
0
0 0.5 1 1.5 2 2.5 3
beta
Figure 13.27: Second-order accurate centred finite difference method with ex-
plicit second-order accurate Runge-Kutta time integration, 100 nodes, 80 time
steps (additional dissipation added for high frequencies).
13.12 Example Unsteady Flow Solutions

In this section we present some example numerical solutions for the transonic
flow around the NLR F-5 wing shown in figure 13.29. In the experiment, the
wing was subjected to sinusoidal pitching about an axis at 50% of the root
chord in a transonic flow. Both experimental and computed results can be
found in RTO-TR-26 (October 2000). On initial inspection the TSD (small-
disturbance potential) results seem to be superior, but this cannot be true from
the theoretical point of view. A possible source of the discrepancy is the relative
coarseness of the Euler meshes used. In order to fully compute the effects of
convecting vorticity in the unsteady case, a fine mesh in the wake areas is
required. Also, since none of the inviscid methods take shock boundary-layer
interactions into account, it could also be that the errors present in the TSD
model fortuitously cancel errors incurred by neglecting viscosity. This seems to
be borne out by the steady results for Euler with boundary layer correction. In
summary, one must keep aware of all types of error when interpreting results,
13.12. EXAMPLE UNSTEADY FLOW SOLUTIONS 187
1
Courant=1.0
0.8
0.6
e_amp
0.4 1.5
Initial
Exact
Numerical
0.2
1
0
0 0.5 1 1.5 2 2.5 3
0.5
u
beta
3 Courant=1.0
0
2.5
2
-0.5
e_pse
0 0.2 0.4 0.6 0.8 1

1.5
x
1
0.5
0
0 0.5 1 1.5 2 2.5 3
beta
Figure 13.28: Third-order accurate (implicit) space-time discretisation, 50 ele-

ments, 2 nodes per element, 80 time steps (additional dissipation added for high
frequencies).
including both discretisation and modelling errors.

Figure 13.29: NLR F-5 wing
-0.6 -0.6
-0.4 -0.4
-0.2 -0.2
0 0
Cp
Cp
0.2 0.2
0.4 0.4
Exper_Upper Exper_Upper
0.6 Exper_Lower 0.6 Exper_Lower
TSD Euler
Full_Potential Euler + BL
0.8 0.8
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x/c x/c
Figure 13.30: Steady Flow M = 0.896, α = 0.5o , Re = 5.79 × 106

TSD: UTSPV21 (BAe) 96 000 nodes
Full Potential: TCITRON (Dassault) 85 470 nodes
Euler: EUGINE (Dassault) 294 851 cells
13.12. EXAMPLE UNSTEADY FLOW SOLUTIONS 189
-25 -25
-20 -20
-15 -15
-10 -10
Re(dCP/Theta)
Im(dCP/Theta)
-5 -5
0 0
5 5
10 Exper_Upper 10 Exper_Upper
15 Exper_Lower 15 Exper_Lower
TSD TSD
20 Full_Potential 20 Full_Potential
Euler Euler
25 25
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x/c x/c
Figure 13.31: Unsteady Flow M = 0.896, αo = 0.001o, θ = 0.111o, k = 0.275

Real (left) and imaginary (right) perturbation components
p13.13.28 Compare the application of a kernel function method with a panel

method for the computation of responses of the wing-tail configuration shown
below. What are their advantages and disadvantages?
p13.13.29 Given two solutions uh and u2h , estimate the error of uh if both
solutions were computed using a first-order accurate method.
p13.13.30 Describe the treatment of moving meshes in the example finite-

difference and finite-volume techniques presented in the notes.
p13.13.31 Is amplification error more important that phase error in the

prediction of flutter? Discuss in the context of a coarse-mesh Euler discretisation
to investigate tail flutter in the the wing-tail configuration above.
Chapter 14
Computing Fluid-Structure
interactions
As discussed in chapter 11, the repeated application of complete complex aeroe-

lastic simulations is limited in practice due to computational cost. Therefore,
most preliminary fully-coupled studies make use of reduced-order order models
for both the aerodynamic flow and structural response.
As we have already discussed how low-order approximations to structural
behaviour can be accomplished, we will start this chapter by describing how
compatible low-order aerodynamic approximations can be generated from com-
plex CFD simulations, and used with their structural counterparts to perform
aeroelastic simulations.
We will then consider optional methods for coupling fully discretised fluid
and structural systems. This is an interesting topic, as standard partitioned
approaches to coupling can lead to considerable errors in quantities of interest
for aeroelasticity. We will examine a technique for estimating these errors, and
then compare the efficiency of partitioned integration methods to that of their
strongly-coupled monolithic counterparts.
14.1 Generalised force models from CFD com-

putations
The equations of motion of a discrete fluid-structure system can typically be
expressed:
[M ]{q̈} + [D]{q̇} + [K]{q} = {Q(t)} (14.1)
If we choose the generalised coordinates, q, to be the amplitudes of the m most

dominant (natural) structural modes, we can construct an efficient reduced-
order representation of the fluid-structure system. In this case, the displacement
191
192 CHAPTER 14. COMPUTING FLUID-STRUCTURE INTERACTIONS
of the structure, u, can be expressed in terms of the dominant structural mode

shapes φ, as:
m
X
u= φi qi (14.2)
i=1
Examples of the dominant natural mode shapes for a finite wing are shown in
figure 14.1.
Figure 14.1: Examples of the first four mode shapes (φi , i = 1..4) for a wing.
(Taken from Raveh, Levy, and Karpel, AIAA paper 00-1325)
Simplified CFD methods are sometimes designed to produce solutions to

harmonic excitations. For such methods, generalised forces Q(t) may be con-
structed using inverse Fourier transforms, as described in previous chapters.
More sophisticated CFD methods usually produce time-domain results. For
these methods, the computation of harmonic responses is often too time con-
suming. Several periods of excitation are required to establish a harmonic result,
and several computations are required to obtain results over the frequency range
required for an inverse transform. Impulse and step responses, in contrast, pro-
vide information for a wide range of frequencies within one calculation. Impulse
responses require only short time intervals to be computed, but are difficult
to realise accurately using a fully discrete representation of the flow. Step re-
sponses, on the other hand, are only slightly longer to compute and are far easier
to realise computationally. If only a limited frequency content is required, it is
also possible to use a smoothed step input.
One of the motivations for performing a complex CFD simulation is to rep-
resent non-linear phenomena in the fluid. Depending on the unsteady behaviour
of these phenomena, it may also be necessary to represent the unsteady force
14.1. GENERALISED FORCE MODELS FROM CFD COMPUTATIONS193
response non-linearly. A general approach is to then use a sequence of impulse-

response problems to construct a Volterra series:
Z ∞
F (t) = h0 + h1 (t − τ )u(τ )dτ
0
Z ∞Z ∞
+ h2 (t − τ1 , t − τ2 )u(τ1 )u(τ2 )dτ1 dτ2
Z0 ∞ Z0 ∞ Z ∞
+ h3 (t − τ1 , t − τ2 , t − τ3 ). . . .
0 0 0
+ ... (14.3)
If sufficient terms are retained, this expression is capable of representing highly

non-linear phenomena. Experience has shown, however, that the third and
higher-order terms (h3 , h4 ,..) tend to be negligible for most aerodynamic flows.
In fact, it is often sufficient to retain only h0 and h1 , even when considering
responses of non-linear phenomena such as shock waves, provided the displace-
ments from the steady non-linear condition are reasonably small.
If is justified to retain only h0 and h1 , generalised forces can be quickly
constructed using a limited number of impulse or step response computations.
For example, the generalised force response to a step excitation with shape ǫφi
may be expressed as an integral over the surface of the geometry S as:
1
Z
Qis (t) = ǫφi · G(t)dS (14.4)
ǫ S
where G(t) is the computed response of the surface fluid forces due to a step
change in the structure geometry of shape φi and magnitude ǫ. An arbitrary
response can then be computed using discrete convolutions for the generalised
forces:
N =t/∆t
X
Qi (t) = Qi (0) + Qis (t) q˙i (t) ∆t i = 1, 2, ... (14.5)
n=1
If it is not sufficient to retain only h0 and h1 , then more general expressions for
Qi (t) must be constructed which include terms with multiple mode displace-
ments.
An alternative to the convolution approach is make a reduced- order model
of the aerodynamic operator. This can be done by computing the eigensystem
of the discrete CFD matrix which expresses the change in the fluid state over
a time step. The eigenvectors with largest magnitude are then used to make a
simplified model of the fluid response (see, for example Dowell, AIAA Journal
Vol 34, No 6, 1996). Although powerful in concept, this approach must be used
with care, as due to the multi-scale nature of fluid phenomena, mode truncations
can significantly reduce the range of behaviours which can be predicted.
14.2 Coupling of fully-discretised fluids and struc-

tures
We now consider the coupled integration of a fully discretised fluid and struc-
tural system. For non-linear fluid models, a mesh which fills the volume of
the fluid domain is normally required. If the boundary conditions are have not
been linearised, then the mesh must also move in response to the deflection of
structural components which border the fluid domain. As the movement of the
mesh must be determined during the solution procedure, the total problem is
often said to be a three-field problem (figure 14.2). The most common approach
F
M
S
Figure 14.2: Coupled integration as a three-field problem
to treating mesh movement is to use a structural analogy, that is to compute

the deflection of the fluid mesh nodes as if they were following the displacement
of a ficticious solid which fills the deforming fluid domain. The final system of
equations for the three-field problem can then be written:
1. Structural field: ∂2q ∂q

M +D + Kq = F (14.6)
∂t2 ∂t
d
Z Z
2. Fluid field: W dV + (F~ − ~ẋW ) · ~n dS = 0(14.7)
dt
V S
3. Fluid Mesh field: ∂2x ∂x

M̃ + D̃ + K̃x = Kc q (14.8)
∂t2 ∂t
Where q are the unknowns associated with structural deflections, W are the un-
knowns associated with the fluid state, and x are the unknowns associated with
the fluid mesh position. A visualisation of how the fields are coupled appears
14.3. PARTITIONED TIME INTEGRATION TECHNIQUES 195
in figure 14.3. Note that the choice of M̃ , K̃, and D̃, are in principal arbitrary,
Interface moving
displacement domain
Structure Fluid Fluid

Mesh
Interface forces
Figure 14.3: Coupling in the three-field representation
as they only affect the dynamics of the fluid-structure system indirectly by de-
termining the quality of the mesh on which the fluid solution is computed. One
can choose M̃ and D̃ to be zero, for example, to avoid introducing complex un-
steady motion of the mesh. K̃ is normally constructed so that fluid elements in
critical regions are stiffer, which prevents them from becoming overly distorted
in response to boundary deformations.
Assuming that the mesh responds instantaneously to the structural defor-
mation, we can make the following classification of time integration techniques
for the remaining fluid-structure system:
Partitioned
A partitioned method advances the solution in time using one (fully-converged)
fluid solution and one (fully-converged) structure solution per time step. Such
methods are asynchronous in nature (first fluid, then structure or first structure,
then fluid). Extrapolation of previous results is sometimes used to reduced the
errors incurred by asynchronicity.
Monolithic
A monolithic method either fully or partially converges the complete fluid-
structure interaction system at the end of each time step by using multiple
fluid and structure solutions. Depending on the level of convergence, errors in-
curred by asynchronicity are eliminated.
Mixed
Mixed methods use monolithic techniques on low-order representations of the
system (e.g. coarse meshes) and partitioned techniques on higher-order repre-
sentations (e.g fine meshes). These have low but significant errors incurred by
asynchronicity.
14.3 Partitioned Time Integration Techniques

In a partitioned time integration technique, the structure and fluid are each
advanced separately in a fixed sequence within the time step. The advantages
of this commonly-used approach are that it is straight-forward to implement,

and can be done with completely separate fluid and structural solvers. Its
disadvantage is that its inherent asynchronicity introduces new sources of error
which can significantly limit the maximum allowable time step for the fluid-
structure system. Partitioned solution techniques may be divided into volume-
continuous and volume-discontinuous methods, as described below.
14.3.1 Volume-continuous methods

Volume continuous methods are those which maintain a synchronous represen-
tation of the position of the structural boundary in the fluid and structure, but
as a consequence have an asynchronous representation the forces acting at the
boundary. Such a method is visualised in figure 14.4, which shows the sequence
of information transfer in a single discrete time step from time tn to time tn+1
. Denoting P n as the fluid pressure and viscous forces at time level n, and X n
as the structural interface displacement at time level n, the sequence is:
Fluid,
Fluid Mesh 4
Pn
Xn+1 3
1
Structure 2
tn tn+1
Figure 14.4: Volume-continuous method
1. Transfer P n to the structure

2. Advance the structure from tn to tn+1
3. Transfer the new structure location, X n+1 , to the fluid mesh
4. Redistribute the fluid mesh and advance the fluid from tn to tn+1
As a consequence, the fluid mesh interface matches structural mesh interface at
each time tn , but the fluid interface forces do not match those of the structural
interface forces at tn .
14.3.2 Volume-discontinuous methods

In a volume-discontinuous method, the synchronicity of the forces that is main-
tained in favour of that of the displacements. The sequence, as shown in figure
14.5, is then:
Fluid,
Fluid Mesh 2
Xn
1 3
Pn+1
Structure 4
tn tn+1
Figure 14.5: Volume-discontinuous method
1. Transfer X n to the fluid mesh
2. Redistribute the fluid mesh and advance the fluid from tn to tn+1
3. Transfer P n+1 to the from the fluid to the structure
4. Advance the structure from tn to tn+1
In this case, the forces in the fluid match those applied to the structure at each
time tn , while the displacements of boundary do not match at tn . In practice the
latter mismatch is not not necessarily as disturbing as it sounds, as the fluid and
structure interface meshes often differ in interface refinement, and thus already
do not match in the steady case.
14.3.3 Analysis of errors due to partitioning

As described above, partitioned schemes cannot simultaneously match both the
displacements and forces across the fluid-structure interface. This leads to errors
in the transfer of energy between fluid and structure. Assuming the partitioned
scheme is consistent, these errors can be reduced by decreasing the time step.
For many problems, however, including aeroelastic instability problems, the
energy balance is crucial. For such cases, certain partitioned schemes can require
very small time steps in order to achieve an acceptable level of accuracy.
In this section we describe a method to estimate the errors introduced parti-
tioning, using the simple one-dimensional piston problem shown in figure 14.6.
This can be used to design partitioned methods with relatively low errors in en-
ergy conservation. The procedure is described in detail in Piperno and Farhat,
Comput. Meths. Appl. Mech. Engrg., vol. 190, No. 24, pp. 3147-3170, (2001).
It can be summarised as:
1. Find the change in energy of the fluid, ∆EF , and the change in energy of
the structure ∆ES , during one time step
F S
CCCCCCCCCCCCCCCCCCCCCCC
P
X
Figure 14.6: A one-dimensional piston problem
2. Estimate the the changes in energy over one period of the oscillation,
(∆EF )T and (∆ES )T , using a continuous approximation for the n = T /∆t
discrete contributions
3. Determine how accurately the energy balance (∆EF )T + (∆ES )T = 0 is
satisfied by the partitioned scheme under consideration
For simple time-marching techniques, we may express the transfer of energy to
the fluid during a single time step ∆t as:
Z tn+1
∆EF = − ẊF APF dt = −(XFn+1 − XFn )AP̄F (14.9)
tn
where :
A = Surface area of the piston (14.10)

P̄F = P n (Euler explicit) (14.11)
P̄F = P n+1 (Euler implicit) (14.12)
(P n+1 + P n )
P̄F = (trapezoidal) (14.13)
2
Similarly, for the structure:
∆ES = +(XSn+1 − XSn )AP̄S (14.14)
If we consider a harmonic solution such as:
XS (t) = Xcos(ωt) (14.15)

PF (t) = P cos(ωt + φ) (14.16)
the energy change over a period can be expressed as the piecewise integration
of ∆E ∆E
∆t (figure 14.7). As ∆t represents an average value, This can be related to
the integration of a continuous function via the midpoint rule (figure 14.8). We
can therefore use the approximate identity:
2π
∆E
Z ω
(∆E)T ≈ dt (14.17)
0 ∆t
55555
55555
55555
∆E 55555
55555
55555
∆t 55555
55555
55555
55555
5555555555 5555555555 55555
5555555555 5555555555 55555
5555555555 5555555555 55555
55555
5555555555
55555 55555
5555555555
55555 55555
55555
5555555555 5555555555 55555
5555555555 5555555555 55555
5555555555 5555555555 55555
Figure 14.7: Energy change over a period of oscillation
∆E
∆t
Figure 14.8: Approximation of the discrete integral using a continuous function
in step 2 of the procedure.
Example Problem
Compute the energy-conservation accuracy of a coupled numerical solution pro-
cedure for the 1D model piston problem shown in figure 14.6. Consider the
case of a volume-discontinuous partitioned solution technique, combined with
trapezoidal time integration of both the fluid and the structure.
Solution
We will assume that the piston is oscillating harmonically, so that the discrete
solution can be described by:
XS (tn ) = Xcos(ωtn ) (14.18)

PF (tn ) = P cos(ωtn + φ) (14.19)
Step 1: For trapezoidal time integration, the change in the energy of the fluid
during one time step may be expressed:
(PFn+1 + PFn )
∆EF = −(XFn+1 − XFn )A (14.20)
2
Similarly, for trapezoidal integration of the structure:
(PSn+1 + PSn )
∆ES = (XSn+1 − XSn )A (14.21)
2
where X is the structural displacement, P is the fluid pressure, and A is the

surface area of the piston.
In order to express the changes energy in terms of the assumed solution’s
parameters, we need to consider the differences in the representation of X and
P for the fluid and the structure when a partitioned solution technique is used.
For a volume-discontinuous technique, the pressure is synchronous, so that:
PSn = PFn = P cos(ωtn + φ) (14.22)

PSn+1 = PFn+1 n
= P cos(ωt + h + φ) (14.23)
where h = ω∆t. The structural displacement is asynchronous, however, so that:
XFn+1 = XSn = Xcos(ωtn ) (14.24)

XFn = XSn−1 = Xcos(ωtn − h) (14.25)
Thus, the change in energy of the fluid and structure during a single time step
is:
P AX
∆EF = − [cos(ωtn ) − cos(ωtn − h)] [cos(ωtn + h + φ) + cos(ωtn + φ)]
(14.26)
2
P AX
∆ES = [cos(ωtn + h) − cos(ωtn )] [cos(ωtn + h + φ) + cos(ωtn + φ)]
(14.27)
2
Step 2: To compute the change in energy over one period of oscillation, we use
the approximate identity:
2π
∆E
Z ω
(∆E)T ≈ dt (14.28)
0 ∆t
For convenience, we define k = P AXcos(φ), d = P AXsin(φ). Then using

the trigonometric formulae from the appendix, the expression for the change in
energy in the fluid can simplified and integrated to yield:
kπ dπ
(∆EF )T ≈ − 2sin2 (h) − [2cos(h)sin(h)] (14.29)
2h 2h
Using the trigonometric series expansions, this is approximately:
h3 2h2

(∆EF )T ≈ kπ −h + + dπ −1 + (14.30)
3 3
Employing a similar procedure, we can also express the approximate change in

energy over a period for the structure as:
h2

dπ
(∆ES )T ≈ [2sin(h)] ≈ dπ 1 − (14.31)
2h 6
14.4. MONOLITHIC TIME INTEGRATION TECHNIQUES 201
Step 3: The energy balance over one period of oscillation is given by:
h3
2
h
(∆EF )T + (∆ES )T = kπ −h + + dπ (14.32)
3 2
The energy conservation accuracy is thus O(h) = O(∆t).
Note that although we have used (rather expensive) O(∆t2 ) accurate trape-
zoidal integration for both the fluid and the structure, the use of a partitioned
solution technique has reduced the time accuracy to O(∆t). For lower frequen-
cies, the displacement is out of phase with the change in pressure (pressure is
low for large X). Thus φ ≈ π and the change in energy of the complete system is
positive (note h > h2 > h3 ). As a result, this partitioned technique is formally
unstable for this undamped problem, and as a consequence the amplitude of the
oscillation will grow slowly in time. In the context of aeroelastic instabilities,
this means that the flutter speed of a more complex configuration might be
underpredicted by such a method.
14.3.4 Improved accuracy by extrapolation

If we are using a volume-discontinuous partitioned method, we can reduce the
discrepancy in the interface position as perceived by the fluid and structure
using an extrapolation procedure for X n+1 , such as:
X n+1 = X n + α1 Ẋ n + α2 Ẋ n−1 (14.33)
This value is then transferred to the fluid in the first part of the partitioned
sequence. Using the analysis presented in the previous section, α1 and α2 can
be chosen in a way which improves the energy-conservation accuracy to O(∆t)3 .
In practice, this can increase the allowable time step by a factor of 10-100.
A similar idea can also be applied to volume-continuous partitioned methods,
by transferring an extrapolated estimate for P n+1 to the structure in the first
step of the sequence. When using extrapolation for the forces, however, some
care must be taken near moving discontinuities in the flow.
14.4 Monolithic Time Integration Techniques

In monolithic time integration, the fluid and structure are advanced simultane-
ously within the time step. In theory, this would lead to a final system of the
form:
    
[F F ] [F M ] [0]  ∆W   
 [0] [M M ] [M S]  ∆x = RHS (14.34)
[SF ] [0] [SS] ∆q
   
Where F F , M S, etc denote sub-matrices coupling the fluid, fluid mesh, and
structural fields. The complete system is not normally solved simultaneously,
however, due to the requirement of having a single program describe both the
fluid and the structure, and due to the global stiffness introduced by the very
different natures of each sub-matrix. The development of optimal decomposition
and preconditioning procedures for the complete matrix remains an interesting
research topic, however.
Thus to avoid the issues mentioned above, monolithic methods are usually
implemented as an iteration loop involving sequential fluid and structure solu-
tion steps performed until the system reaches a specified level of convergence.
Since the fluid problem is normally solved by a subiteration procedure within
each time step, the communication with the structure can be included in each
fluid subiteration with little extra cost.
The prime advantage of monolithic methods is that energy conservation
errors at the fluid structure interface can be kept arbitrarily small, or eliminated.
Like partitioned methods, it is possible to use extrapolation at the interface
to speed the rate of the converge, but this does not influence the final energy
conservation property. In principle, energy conservation errors can be controlled
even if the accuracy of the solution representation is very different on the fluid
and structure sides of the interface (see for example, Van Brummelen, E.H. et
al.: Energy Conservation Under Incompatibility for Fluid Structure Interaction
Problems, Computer Methods in Applied Mechanics and Engineering, Volume
192, Number 25, 20 June 2003, pp. 2727-2748(22)).
14.5 Efficiency of Monolithic and Partitioned Tech-

niques
To put the application of partitioned and monolithic methods in context, in this
section we compare the performance of a volume-discontinuous partitioned with
structural extrapolation to that of an energy conservative monolithic method.
We consider the case of an oscillating 1 DOF piston connected to a cylinder
filled with a gas simulated using the Euler equations. It should be noted that
the partitioned method requires only one fluid-structure solution per time step,
while the (iterative) monolithic method uses several. Both methods use second-
order accurate discretisations for both the fluid and the structure.
The error in the piston’s displacement is shown as a function of time step in
figure 14.9. An additional non-conservative monolithic method is also shown for
comparison. As extrapolation is used in the partitioned case, each of the meth-
ods achieves their full second-order accuracy of convergence. The conservative
monolithic method is several orders of magnitude more accurate on the same
mesh. It is clear that it is the energy-conservation property that is contributing
to this result.
An alternate comparison is shown in figure 14.10, which plots the error in
the piston displacement versus the total number of fluid and structure solutions
14.5. EFFICIENCY OF MONOLITHIC AND PARTITIONED TECHNIQUES203
performed to compute one period of oscillation. In this case we (conserva-

tively) assume that the monolithic method’s fluid solutions are fully converged
before each fluid-structure communication. As a single period of oscillation is
considered, the time step of the methods is changing with the number of fluid-
structure solutions. It is clear that if higher levels of accuracy are desired, then
the monolithic method is considerably more efficient.
-2 partitioned with prediction

monolithic
-3
-4
log10(error)
-5
-6
-7
-8
-9
1 2
log10(fluid-structure solutions )
Figure 14.9: Error produced by different coupling methods as a function of time

step (plot from C Michler, T.U Delft - LR)
Although illustrative, this problem only considers a single mode of oscilla-

tion. If a partitioned method is used in a problem with multiple modal frequen-
cies, it may also be necessary to march at a time step small enough to keep
their highest-frequency mode from growing excessively. This can be a severe
restriction if only low-frequency accuracy is desired.
Finally, the errors committed due to fluid-structure coupling must be seen in
the context of errors introduced by other parts of the discretisation. A summary
of the various sources of error in a typical fluid-structure interaction computation
is shown in figure 14.11.
conservative monolithic
non-conservative monolithic
partitioned withprediction
-2
-3
-4
log10(error)
-5
-6
-7
-8 2
1
-9
-2 -1
log10(deltaT)
Figure 14.10: Error produced during a period of oscillation vs number of fluid

and structure solutions used (plot from C Michler, T.U Delft - LR)
14.5. EFFICIENCY OF MONOLITHIC AND PARTITIONED TECHNIQUES205
Interface energy conservation errors

(when partitioning is used)
influenced by:
− fluid and structure integration methods
− communication timing
− extrapolation procedures
Dissipation of fluid Numerical

transient energy: damping of structure
influenced by: influenced by:

−time march −time march
−time step −time step
−spatial discretisation −spatial discretisation
−mesh size and quality −mesh size and quality
(can be avoided with some

discretisations)
Figure 14.11: Sources of error in the energy balance

p14.6.32 Compute the order of accuracy of energy conservation for a volume-

discontinuous partitioned integration scheme with a first-order implicit time
march for the fluid, and a trapezoidal implicit integration of the structure.
14.7 Handy Formulae
Angle formulae
cos(α ± β) = cos(α)cos(β) ∓ sin(α)sin(β)

sin(α ± β) = sin(α)cos(β) ± cos(α)sin(β)
Integrals of sine and cosine products

2π
π
Z ω
cos2 (ωt)dt =
0 ω
2π
π
Z ω
sin2 (ωt)dt =
0 ω
Z 2π
ω
cos(ωt)sin(ωt)dt = 0
0
h4
Series expansions cos2 (h) ≈ 1 − h2 +
3
h2 h4 h4
cos(h) ≈ 1− + sin2 (h) ≈ h2 −
2 24 3
h3 h5 h3
sin(h) ≈ h− + cos(h)sin(h) ≈ h−2
6 120 3
(During an exam, all of these formulae will be made available when necessary
for completing a question)
Chapter 15
Aeroelasticity in Aircraft
Design
15.1 Requirements
There are a large number of design requirements associated with aeroelastic
phenomena. The basic bounds for these requirements are set by certification
standards. For large aircraft, these can be found in section 25 of the air regu-
lation documents of the relevant country or region (e.g. EASA CS-25, FAR-25,
CAR-525 etc). Aeroelasticity is mentioned most explicitly in section 629. Some
of the CS-25 aeroelasticity-related standards include:
• The aeroplane must be designed to be free from flutter and divergence

for all combinations of altitude and speed encompassed by the VD/MD
versus altitude envelope enlarged at all points by an increase of 15% in
equivalent airspeed at both constant Mach number and constant altitude,
except that the envelope may be limited to a maximum Mach number of
1.0.
• When MD < 1.0, A proper margin of damping exists at all speeds up

to MD, and there is no large and rapid reduction in damping as MD is
approached.
• If concentrated balance weights are used on control surfaces, their effec-

tiveness and strength, including supporting structure, must be substanti-
ated.
• The aeroplane must be designed to be free from control reversal and from
undue loss of longitudinal, lateral, and directional stability and control
within the previously described VD/MD envelope.
• The aeroplane must be free of flutter or divergence that would preclude

safe flight, at any speed up to VD after:
207
208 CHAPTER 15. AEROELASTICITY IN AIRCRAFT DESIGN
- Failure of any single element of the structure supporting any engine,

independently mounted propeller shaft, large auxiliary power unit, or large
externally mounted aerodynamic body.
- Absence of propeller aerodynamic forces resulting from the feathering of
any single propeller, and, for aeroplanes with four or more engines, the
feathering of the critical combination of two propellers.
- Any single propeller rotating at the highest likely overspeed.
- Failure of critical components identified in 25.571 (e.g. wing, empennage,
control surfaces fuselage, engine mounting landing gear).
- Any single failure in any flutter damping system.
where MD is the dive Mach number and VD is the dive speed. Normally,
compliance with the standards must be demonstrated by analysis, wind tunnel
tests, ground vibration tests, flight tests or any other means found necessary
by the certification authority. Aside from the certification standards, however,
an aircraft producer may have more stringent aeroelastic design requirements,
such as those relating to ride quality or control effectiveness.
15.2 Design Process

As aeroelastic design requirements are intimately connected to the operation
envelope of the aircraft, it is normally necessary to consider aeroelastic effects
at all stages in the design process. This helps to prevent later performance
reductions or costly re-designs. Consequently, a sequence of models is normally
employed. These become increasingly sophisticated as the details of the config-
uration are specified. A typically sequence of aeroelastic analysis procedures is
given below:
Pre-design stage
• Analytic modelling
• Wind tunnel testing to support modelling of new aeroelastic features
Design iteration stage

• Repeated aeroelastic calculations to verify constraints on or aeroelastically
optimize the evolving configuration
Frozen design stage

• Tests of structurally scaled models and comparison with analytic struc-
tural models
• Wind tunnel tests with aeroelastically scaled models to verify analytic
aerodynamic models and coupling
15.3. WIND TUNNEL TESTING TECHNIQUES 209
Prototype stage
• Ground vibration and flexibility tests to provide final validation of struc-
tural models
• Analytic calculations to prepare for flight testing
• Flight-test program
We have already discussed various levels of analytic modelling within the

course. The remaining sections in this chapter will therefore concentrate on
methods for determining aeroelastic characteristics experimentally.
15.3 Wind Tunnel Testing Techniques

Aeroelastic wind tunnel testing is particularly complex, due to the need to match
both flow and structural dynamic characteristics at a reduced scale.
In principal an accurate representation of the flow requires the simultane-
ous matching of both Mach and Reynolds numbers. This is not feasible at
the high-Mach and Reynolds numbers combinations at which many aeroelastic
instabilities occur. If the phenomena is dominated by either viscous or com-
pressible effects, it can be sufficient to approximate either the flight Reynolds
number or Mach number. However for fully-coupled problems, such those ac-
companied by strong shock boundary-layer interactions, combinations of tests
and modelling may be necessary.
The construction of structural models which approximate the dynamic char-
acteristics of the full-size configuration is also not a trivial matter. Scaled skin
thicknesses are usually far too small to allow use of materials and attachment
methods (e.g. riveting) similar to those on the full-size aircraft. Often plastics
are substituted due to their formability, but these typically result in exces-
sive internal damping. Normally the solution is not to attempt to construct a
scaled model of the structure, but instead to construct a simplified model which
closely reproduces the dynamic characteristics of the real structure. Figure 15.1
shows a full model of a tilt rotor in NASA’s Transonic Dynamics Tunnel, and
an dynamically-matched underlying structure being prepared for wind-tunnel
use. For a simple wing geometry, a dynamically-matched structure might sim-
ply match the correct mass, GJ and EI distributions as a function of spanwise
position. Normally finite-element methods are used to design and check the
structure, while final verification and fine-tuning might make use of vibration
test equipment.
There are also problems designing supports when there is significant cou-
pling between aircraft dynamic modes aeroelastic modes. One solution is to
support the model using a system of cables designed to allow the excitation of
the dynamic mode of interest, as shown in figure 15.2.
Once the configuration of the model has been specified, a careful testing
procedure must still be designed to accurately measure aeroelastic phenomena
Figure 15.1: Full aeroelastic model in tunnel (left) and a typical underlying
structural model (right)
near their points of instability. Normally this is accomplished by using harmonic

excitation devices to provide a forced response. These include rotating masses
attached to the structure, or oscillating aerodynamic vanes upstream of the
model. Frequency sweeps are carried out at multiple fixed operating conditions
near the point of instability. Values of damping are then estimated by examining
the response characteristics, and a more precise operating condition is chosen.
Special care must be taken, however, when approaching instability points
using such procedures. Aeroelastic damping is notoriously hard to extrapo-
late, as its behaviour is can be strongly non-linear. Often expert knowledge is
required, or use must be made of advanced computational models which incor-
porate information obtained in previous tests. As a last resort, a net is normally
placed downstream to protect the tunnel from large pieces of debris if the model
accidently crosses an instability boundary.
15.4 Ground Vibration Testing

In the final stages of the design process, full size structures become available,
permitting the accurate determination of dynamic characteristics. Normally
this is done by ground vibration testing. A typical setup is shown in figure 15.3.
The aircraft is placed on supports, and then shakers are applied at strategic
points determined by existing finite-element models or previous tests (i.e. away
from nodal lines). The supports can be soft if coupling with rigid-body modes
is required. The design of the shaker depends on the frequencies of interest and
amplitudes required for measurement. Normally mechanical shakers are used
for low frequencies( 1-100 Hz), electromagnetic for medium frequencies (10-1000
Hz) and acoustic for high frequencies (100-10000 Hz). Measurement is typically
carried out using accelerometers, which may be mechanical or electromagnetic.
Using advanced identification techniques, structural non-linearities (e.g control
surface play) can also be measured. For the ground vibration testing of large
15.5. FLIGHT TESTING 211
Figure 15.2: Model mounted on cables in NASA’s Transonic Dynamics Tunnel
aircraft, as shown in figure 15.4, hundreds of accelerometers may be used. When

testing a high-speed aircraft, it may also be necessary to place the structure in
a special enclosure to allow the heating to the conditions encountered in flight.
15.5 Flight testing

The first formal flutter tests were carried out in Germany in the mid 1930’s by
Von Schlippe. His strategy, which is still in use today, was to excite aeroelastic
modes of the aircraft in flight and measure the amplitude of the response while
gradually increasing the airspeed. An example result is shown in figure 15.5. It
can be seen that the flutter condition was quite closely approached in this case,
although fortunately the aircraft and crew survived. In general Von Schlippe’s
technique requires reliable excitation and measurement techniques, and suffers
from the complexity inherent in the extrapolation of damping characteristics.
In several cases, test aircraft have been lost due to deficiencies one or all of these
areas. At the end of the 1958 flight flutter test symposium it was stated that:
“It is hoped that ... tests will be considerably less hazardous than they are
today” Considerable improvements in measurement and prediction have been
made since then, but flight flutter testing remains one of the most hazardous
forms of experimental work encountered in aerospace engineering.
15.5.1 Excitation and Measurement

Due to the wide rage of flight and loading conditions, normally aircraft require
hundreds of test points in order to comply with certification standards. As
these often involve diving to achieve limit airspeeds, the tests must be of limited
Soft Supports
Shaker
Figure 15.3: Supports and shakers for a ground vibration test
duration to prevent significant changes in test altitude. Therefore the excitation

and measurement process must be as quick as possible.
The main requirement for excitation is that it has a sufficient amplitude
over the frequency range of interest. Insufficient excitation can lead to both
scattered results and inaccurate damping estimations. It is also important that
the excitation device alter the configuration as little as possible. Consequently,
several excitations techniques have developed, tailored to the phenomena under
consideration:
• Pilot control pulses
- Can provide excitation up to 10Hz
- Simple, but inconsistent and amplitude often too low
• Oscillating control surfaces
- Up to 50 Hz
- Can be semi-random
- No changes to structure or weight distribution
• Inertial shakers
- force from rotating-mass systems quadratic in rotation speed: limited at
low frequencies, excessive at high fequencies
- heavy, and often too large to be contained within original surface
(affects structural and aerodynamic representation)
- Not used much since the 1950’s, with some exceptions
• Thrusters (Single-use solid-propellant rockets)
- Run time 15-30ms, thrust up to 4000lbs
- lightweight
- limited frequency range
- difficult to time multiple burns for in or out of phase excitation
Figure 15.4: Large-scale ground vibration tests
• Aerodynamic vanes
- Good for low frequencies
- force varies with airspeed
- can distort aerodynamics
• Atmospheric turbulence
- No change to aircraft
- Excites both symmetric and asymmetric modes simultaneously
- Usually amplitude is too small to provide accurate damping estimate
- Flight time lost looking for it
A comparison of the measured response levels achieved via a rotating vane
and atmospheric turbulence is shown in figure 15.6. Strain gauges have been
traditionally used for the measurement of in flight responses, but piezoelectric
accelerometers now much more common as they are compact and light, and have
good frequency response. The current trend is towards self-contained “stick-on”
devices which can transmit data wirelessly to recorders on the aircraft or directly
to the ground station.
15.5.2 Ground data processing and damping prediction

Significant advances in air-to-ground telemetry over the last decades have rev-
olutionized flight flutter testing, resulting in a vastly improved capacity for the
prediction of instability points.
In the 1970’s the first high-speed fast Fourier transforms become possible,
allowing rapid estimation of damping ratios at the last measurement point. Data
Figure 15.5: Cessna AT8 and flight test result
from range of speeds was then used to estimate the zero-damping point before
the next condition was considered. This approach considerably improved safety,
although testing in areas with non-linear damping trends remained risky.
Modern on-line flutter prediction methods make use of robust analysis tech-
niques, where a model with uncertainties is constructed using complex numerical
aerodynamic and structural representations. As the flight progresses, the flight
data is used to conservatively identify the uncertainties in the model, thereby
improving the accuracy of the flutter prediction (see figure 15.7). The result is
then transmitted back to the aircraft, where it can be shown, for example, as a
continuously-adjusted instrument showing the distance from the flutter point.
The latter display device is referred to as a “flutterometer” (figure 15.8).
Rotating Vane Exciter Sweep Atmospheric Turbulence
Measured F16XL wing accelerometer responses (M=0.9, h=30 000 ft)

(taken from NASA TM −4720)
Figure 15.6: Responses to aerodynamic and atmospheric excitation
On−line damping estimate
FEM/CFD p−k analysis
Robust Analysis
Method
Flutter Prediction
Figure 15.7: Robust analysis for improved on-line flutter prediction (From Lind,
R.C, NASA TM-97-206220)
Figure 15.8: Flutterometer (From Lind, R.C, NASA TM-97-206220)

AE4ASM506 Notes 2012

Uploaded by

Copyright:

Available Formats

AE4ASM506 Notes 2012

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

AE4ASM506 Notes 2012

Uploaded by

Copyright:

Available Formats

AE4930 Aeroelasticity

2 Simplified Aeroelastic Problems 21

3.4 Model Equations for Unsteady Flows . . . . . . . . . . . . . . . . 59

4 Analytic Solutions for Incompressible Flows 71

5 Analytic Solutions for Subsonic Compressible Flows 87

6 Analytic Solutions for Supersonic Flows 95

7 Arbitrary Unsteady Aerodynamic Responses 107

8 Compact Models of Structural Systems 113

9 Flutter Calculation Methods 127

10 Dynamic Responses 135

11 Introduction to Computational Aeroelasticity 147

12 Computational structural dynamics 151

12.4 Numerical Integration multi-DOF systems . . . . . . . . . . . . . 157

13 Computational Methods for Unsteady Flows 161

14 Computing Fluid-Structure interactions 191

15 Aeroelasticity in Aircraft Design 207

Figure 1.1: A flexible glider

If the glider’s wing is flexible, then we can expect it to deform. A change

Figure 1.2: Trends in high-performance aircraft

In the early days of aerospace engineering, the significance of aeroelastic

namic regime. Here they experienced shock-induced aerodynamic phenomena

1.2 Course Outline

At the end of this course you should:

• Understand the physical processes which drive aeroelastic phenomena

Figure 1.3: Course Goals

consider simplified (preferably analytic) models which condense complex physi-

Simplified aeroelastic problems (cp 1-2)

Figure 1.4: Course Content

1.3 Additional Course Information

1.3.2 Prerequisite Concepts

1.3.3 Main References

1.3.4 Instructor Contact

1.4 Characteristics of Aeroelastic Problems

Determination of Static and Dynamic Instability Boundaries

Figure 1.5: Instability problem types

Another way of visualising aeroelastic problem types is the famous Collar’s

Determination of Static and Dynamic Response

Figure 1.6: Response problem types

Figure 1.7: Collar’s Triangle (1946)

Figure 1.8: Friedmann’s Hexahedron (1999)

Figure 1.9: X53: control surfaces (upper) and in flight (lower)

2.1 The Typical Section

Figure 2.1: The Typical Section

Structural Parameters Geometric Parameters

Figure 2.2: Structural and geometric parameters

to confuse b with wing span when working out problems).

Figure 2.3: Inertial and aerodynamic parameters

Control Surface Parameters

Figure 2.4: Control surface parameters

2.1.1 Simple Aerodynamic Models

Steady Aerodynamic Model

L(t) = qSCLα α(t) (2.1)

Figure 2.5: Steady angle of attack

Low-Frequency Aerodynamic Model

Figure 2.6: Quasi-steady Angle of attack

2.1.2 Equations of motion for the typical section

In terms of the previously-derived parameters, these can be expressed:

MEA = L(ec) + MAC (2.10)

2.1.3 Torsional Divergence

Figure 2.7: Typical section for the divergence problem