Chapter 11

Aerodynamics
Antony Jameson
Stanford University, Stanford, USA
1 Focus and Historical Background 1
2 Mathematical Models of Fluid Flow 6
3 Potential Flow Methods 10
4 Shock-capturing Algorithms for the Euler and
NavierStokes Equations 25
5 Discretization Scheme for Flows in
Complex Multidimensional Domains 36
6 Time-stepping Schemes 42
7 Aerodynamic Shape Optimization 57
Acknowledgment 76
References 76
1 FOCUS AND HISTORICAL
BACKGROUND
1.1 Classical aerodynamics
This article surveys some of the principal developments of
computational aerodynamics, with a focus on aeronautical
applications. It is written with the perspective that com-
putational mathematics is a natural extension of classical
methods of applied mathematics, which has enabled the
treatment of more complex, in particular nonlinear, math-
ematical models, and also the calculation of solutions in
very complex geometric domains, not amenable to classical
techniques such as the separation of variables.
Encyclopedia of Computational Mechanics, Edited by Erwin
Stein, Ren e de Borst and Thomas J.R. Hughes. Volume 3: Com-
putational Fluid Dynamics. 2004 John Wiley & Sons, Ltd.
ISBN: 0-470-84699-2.
This is particularly true for aerodynamics. Efcient ight
can be achieved only by establishing highly coherent ows.
Consequently, there are many important applications where
it is not necessary to solve the full NavierStokes equations
in order to gain an insight into the nature of the ow, and
useful predictions can be made with simplied mathemati-
cal models. It was already recognized by Prandtl (1904),
and Schlichting and Gersten (1999), essentially contem-
poraneous with the rst successful ights of the Wright
brothers, that in ows at the large Reynolds numbers typi-
cal of powered ight, viscous effects are important chiey
in thin shear layers adjacent to the surface. While these
boundary layers play a critical role in determining whether
the ow will separate and how much circulation will be
generated around a lifting surface, the equations of inviscid
ow are a good approximation in the bulk of the ow eld
external to the boundary layer. In the absence of separation,
a rst estimate of the effect of the boundary layer is pro-
vided by regarding it as increasing the effective thickness
of the body. This procedure can be justied by asymptotic
analysis (Van Dyke, 1964; Ashley and Landahl, 1965).
The classical treatment of the external inviscid ow is
based on Kelvins theorem that in the absence of discontinu-
ities the circulation around a material loop remains constant.
Consequently, an initially irrotational ow remains irrota-
tional. This allows us to simplify the equations further by
representing the velocity as the gradient of a potential. If
the ow is also regarded as incompressible, the governing
equation reduces to Laplaces equation. These simplica-
tions provided the basis for the classical airfoil theory of
Glauert (1926) and Prandtls wing theory (Ashley and Lan-
dahl, 1965; Prandtl and Tietjens, 1934). Supersonic ow
over slender bodies at Mach numbers greater than two is
also well represented by the linearized equations. Tech-
niques for the solution of linearized ow were perfected in
2 Aerodynamics
the period 19351950, particularly by Hayes, who derived
the supersonic area rule (Hayes, 1947).
Classical aerodynamic theory provided engineers with a
good insight into the nature of the ow phenomena, and
a fairly good estimate of the force on simple congura-
tions such as an isolated wing, but could not predict the
details of the ow over the complex conguration of a
complete aircraft. Consequently, the primary tool for the
development of aerodynamic congurations was the wind
tunnel. Shapes were tested and modications selected in
the light of pressure and force measurements together with
ow visualization techniques. In much the same way that
Michelangelo, della Porta, and Fontana could design the
dome of St. Peters through a good physical understanding
of stress paths, so could experienced aerodynamicists arrive
at efcient shapes through testing guided by good physi-
cal insight. Notable examples of the power of this method
include the achievement of the Wright brothers in leaving
the ground (after rst building a wind tunnel), and more
recently Whitcombs discovery of the area rule for transonic
ow, followed by his development of aft-loaded super-
critical airfoils and winglets (Whitcomb, 1956; Whitcomb,
1974; Whitcomb, 1976). The process was expensive. More
than 20 000 hours of wind-tunnel testing were expended
in the development of some modern designs, such as the
Boeing 747.
1.2 The emergence of computational
aerodynamics and its application to
transonic ow
Prior to 1960, computational methods were hardly used in
aerodynamic analysis, although they were already widely
used for structural analysis. The NACA 6 series of air-
foils had been developed during the forties, using hand
computation to implement the Theodorsen method for con-
formal mapping (Theodorsen, 1931). The rst major suc-
cess in computational aerodynamics was the introduction
of boundary integral methods by Hess and Smith (1962)
to calculate potential ow over an arbitrary congura-
tion. Generally known in the aeronautical community as
panel methods, these continue to be used to the present
day to make initial predictions of low speed aerodynamic
characteristics of preliminary designs. It was the com-
pelling need, however, both to predict transonic ow and
to gain a better understanding of its properties and char-
acter that was a driving force for the development of
computational aerodynamics through the period 1970 to
1990.
In the case of military aircraft capable of supersonic
ight, the high drag associated with high g maneuvers
forces them to be performed in the transonic regime. In
the case of commercial aircraft, the importance of transonic
ow stems from the Breguet range equation. This provides
a good rst estimate of range as
R =
V
sfc
L
D
log
W
0
+W
f
W
0
(1)
Here V is the speed, L/D is the lift to drag ratio, sfc is
the specic fuel consumption of the engines, W
0
is the
landing weight, and W
f
is the weight of the fuel burnt.
The Breguet equation clearly exposes the multidisciplinary
nature of the design problem. A lightweight structure is
needed to minimize W
0
. The specic fuel consumption is
mainly the province of the engine manufacturers, and in
fact, the largest advances during the last 30 years have been
in engine efciency. The aerodynamic designer should try
to maximize VL/D. This means that the cruising speed
should be increased until the onset of drag rise due to the
formation of shock waves. Consequently, the best cruis-
ing speed is the transonic regime. The typical pattern
of transonic ow over a wing section is illustrated in
Figure 1.
Transonic ow had proved essentially intractable to ana-
lytic methods. Garabedian and Korn had demonstrated
the feasibility of designing airfoils for shock-free ow in
the transonic regime numerically by the method of com-
plex characteristics (Bauer, Garabedian and Korn, 1972).
Their method was formulated in the hodograph plane,
and it required great skill to obtain solutions correspond-
ing to physically realizable shapes. It was also known
from Morawetzs theorem (Morawetz, 1956) that shock-
free transonic solutions are isolated points.
A major breakthrough was accomplished by Murman and
Cole (1971) with their development of type-dependent dif-
ferencing in 1970. They obtained stable solutions by simply
switching from central differencing in the subsonic zone
to upwind differencing in the supersonic zone and using
Sonic line
Shock wave
Boundary layer
M < 1 M > 1
Figure 1. Transonic ow past an airfoil.
Aerodynamics 3
4
2
1 0
K
s
= 1.3
1
0
2
C
p
C
p
*
X
Figure 2. Scaled pressure coefcient on surface of a thin,
circular-arc airfoil in transonic ow, compared with experimental
data; solid line represents computational result.
a line-implicit relaxation scheme. Their discovery provided
major impetus for the further development of computational
uid dynamics (CFD) by demonstrating that solutions for
steady transonic ows could be computed economically.
Figure 2 taken from their landmark paper illustrates the
scaled pressure distribution on the surface of a symmetric
airfoil. Efforts were soon underway to extend their ideas to
more general transonic ows.
Numerical methods to solve transonic potential ow over
complex congurations were essentially perfected during
the period 1970 to 1982. The AIAA First Computational
Fluid Dynamics Conference, held in Palm Springs in July
1973, signied the emergence of CFD as an accepted tool
for airplane design, and seems to mark the rst use of the
name CFD. The rotated difference scheme for transonic
potential ow, rst introduced by the author at this confer-
ence, proved to be a very robust method, and it provided
the basis for the computer program o22, developed with
David Caughey during 1974 to 1975 to predict transonic
ow past swept wings. At the time we were using the
CDC 6600, which had been designed by Seymour Cray
and was the worlds fastest computer at its introduction,
but had only 131 000 words of memory. This forced the
calculation to be performed one plane at a time, with mul-
tiple transfers from the disk. Flo22 was immediately put
into use at McDonnell Douglas. A simplied in-core ver-
sion of o22 is still in use at Boeing Long Beach today.
Figure 3, shows the result of a recent calculation, using
o22, of transonic ow over the wing of a proposed air-
craft to y in the Martian atmosphere. The result was
obtained with 100 iterations on a 192 32 32 mesh in 7
seconds, using a typical modern workstation. When o22
was rst introduced at Long Beach, the calculations cost
$3000 each. Nevertheless, they found it worthwhile to use it
extensively for the aerodynamic design of the C17 military
cargo aircraft.
In order to treat complete congurations, it was neces-
sary to develop discretization formulas for arbitrary grids.
An approach that proved successful (Jameson and Caughey,
1977), is to derive the discretization formulas from the
Bateman variational principle that the integral of the pres-
sure over the domain,
I =
_
D
p d
is stationary (Jameson, 1978). The resulting scheme is
essentially a nite element scheme using trilinear isopara-
metric elements. It can be stabilized in the supersonic
zone by the introduction of articial viscosity to produce
an upwind bias. The hour-glass instability that results
from the use of one point integration scheme is suppressed
by the introduction of higher-order coupling terms based
on mixed derivatives. The ow solvers (o27-30) based
on this approach were subsequently incorporated in Boe-
ings A488 software, which was used in the aerodynamic
design of Boeing commercial aircraft throughout the eight-
ies (Rubbert, 1994).
In the same period, Perrier was focusing the research
efforts at Dassault on the development of nite element
methods using triangular and tetrahedral meshes, because
he believed that if CFD software was to be really use-
ful for aircraft design, it must be able to treat complete
congurations. Although nite element methods were more
computationally expensive, and mesh generation continued
to present difculties, nite element methods offered a route
toward the achievement of this goal. The Dassault/INRIA
group was ultimately successful, and they performed tran-
sonic potential ow calculations for complete aircraft such
as the Falcon 50 in the early eighties (Bristeau et al.,
1985).
1.3 The development of methods for the Euler
and NavierStokes equations
By the eighties, advances in computer hardware had made
it feasible to solve the full Euler equations using software
that could be cost-effective in industrial use. The idea of
directly discretizing the conservation laws to produce a
nite volume scheme had been introduced by MacCormack
and Paullay (1972). Most of the early ow solvers tended
to exhibit strong pre- or post-shock oscillations. Also, in a
workshop held in Stockholm in 1979, (Rizzi and Viviand,
1979) it was apparent that none of the existing schemes
converged to a steady state. These difculties were resolved
during the following decade.
4 Aerodynamics


S
y
m
b
o
l
S
o
u
r
c
e

F
L
O

2
2

+

L
/
N
M

+

S

A
l
p
h
a



6
.
7
0
0





C
D




.
0
3
1
9



C
M

0
.
0
1
2
2
5
C
o
m
p
a
r
i
s
o
n

o
f

c
h
o
r
d
w
i
s
e

p
r
e
s
s
u
r
e

d
i
s
t
r
i
b
u
t
i
o
n
s
b
a
s
e
l
i
n
e

M
A
R
S
0
0

f
l
y
i
n
g

w
i
n
g

c
o
n
f
i
g
u
r
a
t
i
o
n
M
a
c
h

=

0
.
6
5
0
,


C
L

=

0
.
6
1
5
S
o
l
u
t
i
o
n


1
u
p
p
e
r
-
s
u
r
f
a
c
e

i
s
o
b
a
r
s
(
c
o
n
t
o
u
r
s

a
t

0
.
0
5

C
p
)
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0

2
.
0

1
.
5

1
.
0

0
.
5
0
.
0
0
.
5
1
.
0
0
.
0
0
.
5
1
.
0
C
p

2
.
0

1
.
5

1
.
0

0
.
5
0
.
0
0
.
5
1
.
0
C
p

2
.
0

1
.
5

1
.
0

0
.
5
C
p
X

/

C
X

/

C
X

/

C
1
4
.
6
%

s
p
a
n
3
9
.
0
%

s
p
a
n
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0

2
.
0

1
.
5

1
.
0

0
.
5
0
.
0
0
.
5
1
.
0
C
p
X

/

C
X

/

C
X

/

C

9
2
.
7
%

s
p
a
n
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0

2
.
0

1
.
5

1
.
0

0
.
5
0
.
0
0
.
5
1
.
0
C
p

7
8
.
0
%

s
p
a
n
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0
2
4
.
4
%

s
p
a
n
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0

2
.
0

1
.
5

1
.
0

0
.
5
0
.
0
0
.
5
1
.
0
C
p
X

/

C
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0
0
.
0
%

s
p
a
n
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0

2
.
0

1
.
5

1
.
0

0
.
5
0
.
0
0
.
5
1
.
0
C
p
6
3
.
4

s
p
a
n
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0

2
.
0

1
.
5

1
.
0

0
.
5
0
.
0
0
.
5
1
.
0
C
p
X

/

C
5
3
.
7
%


s
p
a
n
F
i
g
u
r
e
3
.
P
r
e
s
s
u
r
e
d
i
s
t
r
i
b
u
t
i
o
n
o
v
e
r
t
h
e
w
i
n
g
o
f
a
M
a
r
s
L
a
n
d
e
r
u
s
i
n
g

o
2
2
(
s
u
p
p
l
i
e
d
b
y
J
o
h
n
V
a
s
s
b
e
r
g
)
.
Aerodynamics 5
The JamesonSchmidt Turkel (JST) scheme (Jameson,
Schmidt and Turkel, 1981a), which used RungeKutta time
stepping and a blend of second- and fourth-differences
(both to control oscillations and to provide background
dissipation), consistently demonstrated convergence to a
steady state, with the consequence that it has remained one
of the most widely used methods to the present day.
A fairly complete understanding of shock-capturing algo-
rithms was achieved, stemming from the ideas of Godunov,
Van Leer, Harten, and Roe. The issue of oscillation
control and positivity had already been addressed by
Godunov (1959) in his pioneering work in the 1950s
(translated into English in 1959). He had introduced the
concept of representing the ow as piecewise constant
in each computational cell, and solving a Riemann prob-
lem at each interface, thus obtaining a rst-order accu-
rate solution that avoids nonphysical features such as
expansion shocks. When this work was eventually rec-
ognized in the West, it became very inuential. It was
also widely recognized that numerical schemes might ben-
et from distinguishing the various wave speeds, and
this motivated the development of characteristics-based
schemes.
The earliest higher-order characteristics-based methods
used ux vector splitting (Steger and Warming, 1981),
but suffered from oscillations near discontinuities similar
to those of central-difference schemes in the absence of
numerical dissipation. The monotone upwind scheme for
conservation laws (MUSCL) of Van Leer (1974) extended
the monotonicity-preserving behavior of Godunovs scheme
to higher order through the use of limiters. The use of
limiters dates back to the ux-corrected transport (FCT)
scheme of Boris and Book (1973). A general framework for
oscillation control in the solution of nonlinear problems was
provided by Hartens concept of total variation diminishing
(TVD) schemes. It nally proved possible to give a rigorous
justication of the JST scheme (Jameson, 1995a; Jameson,
1995b).
Roes introduction of the concept of locally lineariz-
ing the equations through a mean value Jacobian (Roe,
1981) had a major impact. It provided valuable insight
into the nature of the wave motions and also enabled the
efcient implementation of Godunov-type schemes using
approximate Riemann solutions. Roes ux-difference split-
ting scheme has the additional benet that it yields a
single-point numerical shock structure for stationary nor-
mal shocks. Roes and other approximate Riemann solu-
tions, such as that due to Osher, have been incorporated
in a variety of schemes of Godunov type, including the
essentially nonoscillatory (ENO) schemes of Harten et al.
(1987).
Solution methods for the Reynolds-averaged Navier
Stokes (RANS) equations had been pioneered in the seven-
ties by MacCormack and others, but at that time they were
extremely expensive. By the nineties, computer technol-
ogy had progressed to the point where RANS simulations
could be performed with manageable costs, and they began
to be fairly widely used by the aircraft industry. The
need for robust and reliable methods to predict hypersonic
ows, which contain both very strong shock wave and near
vacuum regions, gave a further impetus to the development
of advanced shock-capturing algorithms for compressible
viscous ow.
1.4 Overview of the article
The development of software for aerodynamic simulation
can be broken down into ve main steps.
1. The choice of a mathematical model that represents
the physical phenomena that are important for the
application;
2. mathematical analysis of the model to ensure existence
and uniqueness of the solutions;
3. formulation of a stable and convergent discretization
scheme;
4. implementation in software;
5. validation.
Thorough validation is difcult and time consuming. It
should include verication procedures for program cor-
rectness and consistency checks. For example, does the
numerical solution of a symmetric prole at zero angle of
attack preserve the symmetry, with no lift? It should also
include mesh renement studies to verify convergence and,
ideally, comparisons with the results of other computer pro-
grams that purport to solve the same equations. Finally, if it
is sufciently well established that the software provides an
accurate solution of the chosen mathematical model, com-
parisons with experimental data should show whether the
model adequately represents the true physics or establish
its range of applicability.
This article is primarily focused on the third step, dis-
cretization. The complexity of predicting highly nonlinear
transonic and hypersonic ows has forced the emergence
of an entirely new class of numerical algorithms and a
supporting body of theory, which is reviewed in this arti-
cle. Section 2 presents a brief survey of the mathematical
models of uid ow that are relevant to different ight
regimes. Section 3 surveys potential ow methods, which
continue to be useful for preliminary design because of their
low computational costs and rapid turn around. Section 4
focuses on the formulation of shock-capturing methods
6 Aerodynamics
for the Euler and RANS equations. Section 5 discusses
alternative ways to discretize the equations in complex
geometric domains using either structured or unstructured
meshes. Section 6 discusses time-stepping schemes, includ-
ing convergence acceleration techniques for steady ows
and the formulation of accurate and efcient time-stepping
techniques for unsteady ows. The article concludes with
a discussion of methods to solve inverse and optimum
shape-design problems.
2 MATHEMATICAL MODELS OF FLUID
FLOW
The NavierStokes equations state the laws of conservation
of mass, momentum, and energy for the ow of a gas in
thermodynamic equilibrium. In the Cartesian tensor nota-
tion, let x
i
be the coordinates, p, , T , and E the pressure,
density, temperature, and total energy, and u
i
the velocity
components. Each conservation equation has the form
w
t
+
f
j
x
j
= 0 (2)
For the mass equation
w = , f
j
= u
j
(3)
For the i momentum equation
w
i
= u
i
, f
ij
= u
i
u
j
+p
ij

ij
(4)
where
ij
is the viscous stress tensor, which for a Newto-
nian uid is proportional to the rate of strain tensor and the
bulk dilatation. If and are the coefcients of viscosity
and bulk viscosity, then

ij
=
_
u
i
x
j
+
u
j
x
i
_
+
ij
_
u
k
x
k
_
(5)
Typically = 2/3. For the energy equation
w = E, f
j
= Hu
j

jk
u
k

T
x
j
(6)
where is the coefcient of heat conduction and H is the
total enthalpy,
H = E +
p

In the case of a perfect gas, the pressure is related to the

density and energy by the equation of state
p = ( 1)
_
E
1
2
q
2
_
(7)
where
q
2
= u
i
u
i
and is the ratio of specic heats. The coefcient of thermal
conductivity and the temperature satisfy the relations
k =
c
p

Pr
, T =
p
R
(8)
where c
p
is the specic heat at constant pressure, R is the
gas constant, and Pr is the Prandtl number. Also the speed
of sound c is given by the ratio
c
2
=
p

(9)
and a key dimensionless parameter governing the effects of
compressibility is the Mach number
M =
q
c
where q is the magnitude of the velocity.
If the ow is inviscid, the boundary condition that must
be satised at a solid wall is
u n = u
i
n
i
= 0 (10)
where n denotes the normal to the surface. Viscous ows
must satisfy the no-slip condition
u = 0 (11)
Viscous solutions also require a boundary condition for the
energy equation. The usual practice in pure aerodynamic
simulations is either to specify the isothermal condition
T = T
0
(12)
or to specify the adiabatic condition
T
n
= 0 (13)
corresponding to zero heat transfer. The calculation of heat
transfer requires an appropriate coupling to a model of the
structure.
Aerodynamics 7
For an external ow, the ow variables should approach
free-stream values
p = p

, =

, T = T

, u = u

for upstream at the inow boundary. If any entropy is

generated, the density for downstream at the outow bound-
ary cannot recover to

if the pressure recovers to p

.
In fact, if trailing vortices persist downstream, the pres-
sure does not recover to p

. In general, it is necessary
to examine the incoming and outgoing waves at the outer
boundaries of the ow domain. Boundary values should
then only be imposed for quantities transported by the
incoming waves. In a subsonic ow, there are four incoming
waves at the inow boundary, and one escaping acoustic
wave. Correspondingly, four quantities should be speci-
ed. At the outow boundary, there are four outgoing
waves, so one quantity should be specied. One way to
do this is to introduce Riemann invariants corresponding to
a one-dimensional ow normal to the boundary, as will be
discussed in Section 5.4. In a supersonic ow, all quantities
should be xed at the inow boundary, while they should
all be extrapolated at the outow boundary. The proper
specication of inow and outow boundary conditions is
particularly important in the calculation of internal ows.
Otherwise spurious wave reections may severely corrupt
the solution.
In smooth regions of the ow, the inviscid equations can
be written in quasilinear form as
w
t
+A
i
w
x
i
= 0 (14)
where A
i
are the Jacobians f
i
/w. By transforming to the
symmetrizing variables, which may be written in differen-
tial form as
d w =
_
dp
c
, du
1
, du
2
, du
3
, dp c
2
d
_
T
(15)
the Jacobians assume the symmetric form

A
i
=
_
_
_
_
_
_
u
i

i1
c
i2
c
i3
c 0

i1
c u
i
0 0 0

i2
c 0 u
i
0 0

i3
c 0 0 u
i
0
0 0 0 0 u
i
_

_
(16)
where
ij
are the Kronecker deltas. Equation (14) becomes
w
t
+

A
i
w
x
i
= 0 (17)
The Jacobians for the conservative variables may now be
expressed as
A
i
= T

A
i
T
1
(18)
where
T
1
=
w
w
=
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
(1)
q
2
2c
(1)
u
1
c
( 1)
u
2
c
(1)
u
3
c
1
c

u
1

0 0 0

u
2

0
1

0 0

u
3

0 0
1

0
(1)
q
2
2
c
2
(1)u
1
(1)u
2
(1)u
3
1
_

_
and
T =
w
w
=
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_

c
0 0 0
1
c
2
u
1
c
0 0
u
1
c
2
u
2
c
0 0
u
2
c
2
u
3
c
0 0
u
3
c
2
H
c
u
1
u
2
u
3

q
2
2c
2
_

_
(19)
The decomposition (17) clearly exposes the wave struc-
ture in solutions of the gas-dynamic equations. The wave
speeds appear as the eigenvalues of the linear combination

A = n
i

A
i
(20)
where n is a unit direction vector. They are
_
q
n
+c, q
n
c, q
n
, q
n
, q
n
_
T
(21)
where q
n
= q n. Corresponding to the fact that

A is
symmetric, one can nd a set of orthogonal eigenvectors,
which may be normalized to unit length. Then one can
express

A =

M

M
1
(22)
where is diagonal, with the eigenvalues as its elements.
The modal matrix

M containing the eigenvectors as its
8 Aerodynamics
columns is

M =
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
1

2

1

2
0 0 0
n
1

2
n
1

2
0 n
3
n
2
n
2

2
n
2

2
n
3
0 n
1
n
3

2
n
3

2
n
2
n
1
0
0 0 n
1
n
2
n
3
_

_
(23)
and

M
1
=

M
T
. The Jacobian matrix A = n
i
A
i
now takes
the form
A = MM
1
(24)
where
M = T

M, M
1
=

M
T
T
1
(25)
In the design of difference schemes, it proves useful to
introduce the absolute Jacobian matrix |A|, in which the
eigenvalues are replaced by their absolute values, as will
be discussed in Section 4.4.
Corresponding to the thermodynamic relation
dp

= dh T dS
where S is the entropy log(p/
1
), the last variable of d w
corresponds to p dS, since c
2
= (dp/d). It follows that in
the absence of shock waves S is constant along streamlines.
If the ow is isentropic, then (dp/d)
1
, and the rst
variable can be integrated to give 2c/( 1). Then we may
take the transformed variables as
w =
_
2c
1
, u
1
, u
2
, u
3
, S
_
T
(26)
In the case of a one-dimensional ow, the equations for the
Riemann invariants are recovered by adding and subtracting
the equations for 2c/( 1) and u
1
.
In order to calculate solutions for ows in complex geo-
metric domains, it is often useful to introduce body-tted
coordinates through global, or, as in the case of isopara-
metric elements, local transformations. With the body now
coinciding with a coordinate surface, it is much easier to
enforce the boundary conditions accurately. Suppose that
the mapping to computational coordinates (
1
,
2
,
3
) is
dened by the transformation matrices
K
ij
=
x
i

j
, K
1
ij
=

i
x
j
, J = det(K) (27)
The NavierStokes equations (26) become

t
(Jw) +

i
F
i
(w) = 0 (28)
Here the transformed uxes are
F
i
= S
ij
f
j
(29)
where
S = JK
1
(30)
The elements of S are the cofactors of K, and in a
nite volume discretization, they are just the face areas
of the computational cells projected in the x
1
, x
2
, and x
3
directions. Using the permutation tensor
ijk
we can express
the elements of S as
S
ij
=
1
2

jpq

irs
x
p

r
x
q

s
(31)
Then

i
S
ij
=
1
2

jpq

irs
_

2
x
p

i
x
q

s
+
x
p

2
x
q

i
_
= 0 (32)
Dening scaled contravariant velocity components as
U
i
= S
ij
u
j
(33)
the ux formulas may be expanded as
F
i
=
_

_
U
i
U
i
u
1
+S
i1
p
U
i
u
2
+S
i2
p
U
i
u
3
+S
i3
p
U
i
H
_

_
(34)
If we choose a coordinate system so that the boundary is
at
l
= 0, the wall boundary condition for inviscid ow is
now
U
l
= 0 (35)
An indication of the relative magnitude of the inertial
and viscous terms is given by the Reynolds number
Re =
UL

(36)
where U is a characteristic velocity and L a representa-
tive length. The viscosity of air is very small, and typical
Aerodynamics 9
Reynolds numbers for the ow past a component of an
aircraft such as a wing are of the order of 10
7
or more,
depending on the size and speed of the aircraft. In this sit-
uation, the viscous effects are essentially conned to thin
boundary layers covering the surface. Boundary layers may
nevertheless have a global impact on the ow by caus-
ing separation. Unfortunately, unless they are controlled
by active means such as suction through a porous sur-
face, boundary layers are unstable and generally become
turbulent.
Using dimensional analysis, Kolmogorovs theory of tur-
bulence (Kolmogorov, 1941) estimates the length scales of
the smallest persisting eddies to be of order (1/Re
3/4
) in
comparison with the macroscopic length scale of the ow.
Accordingly the computational requirements for the full
simulation of all scales of turbulence can be estimated as
growing proportionally to Re
9/4
, and are clearly beyond the
reach of current computers. Turbulent ows may be sim-
ulated by the RANS equations, where statistical averages
are taken of rapidly uctuating components. Denoting uc-
tuating parts by primes and averaging by an overbar, this
leads to the appearance of Reynolds stress terms of the
form u

i
u

j
, which cannot be determined from the mean
values of the velocity and density. Estimates of these addi-
tional terms must be provided by a turbulence model. The
simplest turbulence models augment the molecular viscos-
ity by an eddy viscosity that crudely represents the effects
of turbulent mixing, and is estimated with some character-
istic length scale such as the boundary layer thickness. A
rather more elaborate class of models introduces two addi-
tional equations for the turbulent kinetic energy and the rate
of dissipation. Existing turbulence models are adequate for
particular classes of ow for which empirical correlations
are available, but they are generally not capable of reli-
ably predicting more complex phenomena, such as shock
waveboundary layer interaction. The current status of tur-
bulence modeling is reviewed by Wilcox (1998), Haase
et al. (1997), Leschziner (2003), and Durbin in an article
in this Encyclopedia.
Outside the boundary layer, excellent predictions can be
made by treating the ow as inviscid. Setting
ij
= 0 and
eliminating heat conduction from equations (3, 4 and 6)
yields the Euler equations for inviscid ow. These are
a very useful model for predicting ows over aircraft.
According to Kelvins theorem, a smooth inviscid ow that
is initially irrotational remains irrotational. This allows one
to introduce a velocity potential such that u
i
= /x
i
.
The Euler equations for a steady ow now reduce to

x
i
_

x
i
_
= 0 (37)
In a steady inviscid ow, it follows from the energy
equation (6) and the continuity equation (3) that the total
enthalpy is constant
H =
c
2
1
+
1
2
u
i
u
i
= H

(38)
where the subscript is used to denote the value in the
far eld. According to Croccos theorem, vorticity in a
steady ow is associated with entropy production through
the relation
u +T S = H = 0
where u and are the velocity and vorticity vectors, T is Q1
the temperature, and S is the entropy. Thus, the introduction
of a velocity potential is consistent with the assumption of
isentropic ow.
Substituting the isentropic relationship p/

= constant,
and the formula for the speed of sound, equation (38) can
be solved for the density as

=
_
1 +
1
2
M
2

_
1
u
i
u
i
u
2

__
1/(1)
(39)
It can be seen from this equation that

u
i
=
u
i
c
2
(40)
and correspondingly in isentropic ow
p
u
i
=
dp
d

u
i
= u
i
(41)
Substituting (/x
j
) = (/u
i
)(u
i
/x
j
), the potential
ow equation (37) can be expanded in quasilinear form as
c
2

x
2
i
u
i
u
j

x
i
x
j
= 0 (42)
If the ow is locally aligned, say, with the x
1
axis, equa-
tion (42) reads as
(1 M
2
)

x
2
1
+

2

x
2
2
+

2

x
2
3
= 0 (43)
where M is the Mach number u
1
/c. The change from an
elliptic to a hyperbolic partial differential equation as the
ow becomes supersonic is evident.
The potential ow equation (42) also corresponds to the
Bateman variational principle that the integral over the
10 Aerodynamics
domain of the pressure
I =
_
D
p d (44)
is stationary. Here d denotes the volume element. Using
the relation (41), a variation p results in a variation
I =
_
D
p
u
i
u
i
d =
_
u
i

x
i
d
or, on integrating by parts with appropriate boundary con-
ditions
I =
_
D

x
i
_

x
i
_
d
Then I = 0 for an arbitrary variation if equation (37)
holds.
The equations of inviscid supersonic ow admit discon-
tinuous solutions, both shock waves and contact discon-
tinuities, which satisfy the Rankine Hugoniot jump con-
ditions (Liepmann and Roshko, 1957). Only compression
shock waves are admissible, corresponding to the pro-
duction of entropy. Expansion shock waves cannot occur
because they would correspond to a decrease in entropy.
Because shock waves generate entropy, they cannot
be exactly modeled by the potential ow equations. The
amount of entropy generated is proportional to (M 1)
3
where M is the Mach number upstream of the shock.
Accordingly, weak solutions admitting isentropic jumps
that conserve mass but not momentum are a good approxi-
mation to shock waves, as long as the shock waves are quite
weak (with a Mach number < 1.3 for the normal velocity
component upstream of the shock wave). Stronger shock
waves tend to separate the ow, with the result that the
inviscid approximation is no longer adequate. Thus this
model is well balanced, and it has proved extremely useful
for the prediction of the cruising performance of transport
aircraft. An estimate of the pressure drag arising from shock
waves is obtained because of the momentum decit through
an isentropic jump.
If one assumes small disturbances about a free stream in
the x
i
direction, and a Mach number close to unity, equa-
tion (43) can be reduced to the transonic small disturbance
equation in which M
2
is estimated as
M
2

_
1 ( +1)

x
1
_
This is the simplest nonlinear model of compressible ow.
The nal level of approximation is to linearize equa-
tion (43) by replacing M
2
by its free-stream value M
2

. In
the subsonic case, the resulting Prandtl Glauert equation
can be reduced to Laplaces equation by scaling the x
i
coor-
dinate by (1 M
2

)
1/2
. Irrotational incompressible ow
satises the Laplaces equation, as can be seen by setting
= constant, in equation (37). The relationships between
some of the widely used mathematical models is illustrated
in Figure 4. With limits on the available computing power,
and the cost of the calculations, one has to make a trade-off
between the complexity of the mathematical model and the
complexity of the geometric conguration to be treated.
The computational requirements for aerodynamic simu-
lation are a function of the number of operations required
per mesh point, the number of cycles or time steps needed
to reach a solution, and the number of mesh points needed
to resolve the important features of the ow. Algorithms
for the three-dimensional transonic potential ow equation
require about 500 oating-point operations per mesh point
per cycle. The number of operations required for an Euler
simulation is in the range of 1000 to 5000 per time step,
depending on the complexity of the algorithm. The number
of mesh intervals required to provide an accurate repre-
sentation of a two-dimensional inviscid transonic ow is
of the order of 160 wrapping around the prole, and 32
normal to the airfoil. Correspondingly, about 200 000 mesh
cells are sufcient to provide adequate resolution of three-
dimensional inviscid transonic ow past a swept wing, and
this number needs to be increased to provide a good simu-
lation of a more complex conguration such as a complete
aircraft. The requirements for viscous simulations by means
of turbulence models are much more severe. Good resolu-
tion of a turbulent boundary layer needs about 32 intervals
inside the boundary layer, with the result that a typical mesh
for a two-dimensional NavierStokes calculation contains
512 intervals wrapping around the prole, and 64 intervals
in the normal direction. A corresponding mesh for a swept
wing would have, say, 512 64 256 8 388 608 cells,
leading to a calculation at the outer limits of current com-
puting capabilities. The hierarchy of mathematical models
is illustrated in Figure 5, while Figure 6 gives an indication
of the boundaries of the complexity of problems which can
be treated with different levels of computing power. The
vertical axis indicates the geometric complexity, and the
horizontal axis the equation complexity.
3 POTENTIAL FLOW METHODS
3.1 Boundary integral methods
The rst major success in computational aerodynamics
was the development of boundary integral methods for the
solution of the subsonic linearized potential ow equation
(1 M
2

)
xx
+
yy
= 0 (45)
Aerodynamics 11
p
n

___
t
= 0
= 0

___
Viscosity = 0
Linearize

No
streamwise
viscous terms
Vorticity = 0 Density = const. Density = const.
Boundary
layer eqs.
Thinlayer
NS eqs.
NavierStokes
eqs.
Unsteady viscous
compressible flow
Euler eqs. Laplace eq.
Parabolized
NS eqs.
PrandtlGlauert
eq.
Potential eq.
Small
perturb.
Transonic small
perturb. eq.
Figure 4. Equations of uid dynamics for mathematical models of varying complexity. (Supplied by Luis Miranda, Lockheed
Corporation).
IV. RANS (1990s)
I
n
c
r
e
a
s
i
n
g

c
o
m
p
l
e
x
i
t
y
m
o
r
e

a
c
c
u
r
a
t
e

f
l
o
w

p
h
y
s
i
c
s
III. Euler (1980s)
D
e
c
r
e
a
s
i
n
g

c
o
m
p
u
t
a
t
i
o
n
a
l

c
o
s
t
s
II. Nonlinear potential (1970s)
I. Linear potential (1960s)
Inviscid, irrotational
linear
+ Viscous
+ Rotation
+ Nonlinear
Figure 5. Hierarchy of uid ow models.
This can be reduced to Laplaces equation by stretching the
x coordinate by the factor

(1 M
2

). Then, according
to potential theory, the general solution can be repre-
sented in terms of a distribution of sources or doublets, or
both sources and doublets, over the boundary surface. The
boundary condition is that the velocity component normal
to the surface is zero. Assuming, for example, a source dis-
tribution of strength (Q) at the point Q of a surface S,
this leads to the integral equation
2
p

__
S
(Q)n
p

_
1
r
_
= 0 (46)
where P is the point of evaluation, and r is the distance
from P to Q. A similar equation can be found for a doublet
distribution, and it usually pays to use a combination.
12 Aerodynamics
Linear
1 Mflops
CDC 6600
100 Gflops 100 Mflops
CRAY XMP
10 Mflops
CONVEX
Inviscid
Euler
Nonlinear
potential
flow
NavierStokes Reynolds
averaged
2-D airfoil
3-D wing
Aircraft
Figure 6. Complexity of the problems that can be treated with different classes of computer (1 op = 1 oating-point operation per sec-
ond; 1 Mop = 10
6
ops; 1 Gop = 10
9
ops). A color version of this image is available at http://www.mrw.interscience.wiley.com/ecm
Wing surface pressure distributions
Experiment
Theory
2.0
1.0
1.0
0
Sta 2
C
p
x/c
2.0
1.0
1.0
1.0
0
x/c
Sta 6
C
p
2.0
1.0
1.0
0
C
p
1.0
Sta 4
1.0
1.0
0
2.0
C
p
1.0
Sta 8
(c)
0.3
2.0 0.0 2.0
()
Lift variation with angle of attack
4.0 6.0
C
L
0.6
(b)
1.0
x/c
x/c
Surface panel
representation
S
t
a

2
S
t
a

4
S
t
a

6
S
t
a

8
(a)
Figure 7. Panel method applied to Boeing 747. (Supplied by Paul Rubbert, the Boeing Company.)
Aerodynamics 13
Equation (46) can be reduced to a set of algebraic equations
by dividing the surface into quadrilateral panels, assuming
a constant source strength on each panel, and satisfying
the condition of zero normal velocity at the center of each
panel. This leads to N equations for the source strengths on
N panels.
The rst such method was introduced by Hess and Smith
(1962). The method was extended to lifting ows, together
with the inclusion of doublet distributions, by Rubbert and
Saaris (1968). Subsequently higher-order panel methods (as
these methods are generally called in the aircraft industry)
have been introduced. A review has been given by Hunt
(1978). An example of a calculation by a panel method is
shown in Figure 7. The results are displayed in terms of
the pressure coefcient dened as
c
p
=
p p

1
2

q
2

Figure 8 illustrates the kind of geometric conguration that

can be treated by panel methods.
In comparison with eld methods, which solve for the
unknowns in the entire domain, panel methods have the
advantage that the dimensionality is reduced. Consider a
three-dimensional ow eld on an n n n grid. This
would be reduced to the solution of the source or doublet
strengths on N = O(n
2
) panels. Since, however, every
panel inuences every other panel, the resulting equations
have a dense matrix. The complexity of calculating the
N N inuence coefcients is then O(n
4
). Also, O(N
3
) =
O(n
6
) operations are required for an exact solution. If one
directly discretizes the equations for the three-dimensional
domain, the number of unknowns is n
3
, but the equations
are sparse and can be solved with O(n) iterations or even
Figure 8. Panel method applied to ow around Boeing 747 and
space shuttle. (Supplied by Allen Chen, the Boeing Company.)
with a number of iterations independent of n if a multigrid
method is used.
Although the eld methods appear to be potentially more
efcient, the boundary integral method has the advan-
tage that it is comparatively easy to divide a complex
surface into panels, whereas the problem of dividing a
three-dimensional domain into hexahedral or tetrahedral
cells remains a source of extreme difculty. Moreover
the operation count for the solution can be reduced by
iterative methods, while the complexity of calculating
the inuence coefcients can be reduced by agglomera-
tion (Vassberg, 1997). Panel methods thus continue to be
widely used both for the solution of ows at low Mach
numbers for which compressibility effects are unimportant,
and also to calculate supersonic ows at high Mach num-
bers, for which the linearized equation (45) is again a good
approximation.
3.2 Formulation of the numerical method for
transonic potential ow
The case of two-dimensional ow serves to illustrate the
formulation of a numerical method for solving the transonic
potential ow equation. With velocity components u, v and
coordinates x, y equation (37) takes the form

x
(u) +

y
(v) = 0 (47)
The desired solution should have the property that is
continuous, and the velocity components are piecewise
continuous, satisfying equation (47) at points where the
ow is smooth, together with the jump condition,
[v]
dy
dx
[v] = 0 (48)
across a shock wave, where [ ] denotes the jump, and
(dy/dx) is the slope of the discontinuity. That is to say,
should be a weak solution of the conservation law (47),
satisfying the condition,
__
(u
x
+v
y
) dx dy = 0 (49)
for any smooth test function , which vanishes in the far
eld.
The general method to be described stems from the
idea introduced by Murman and Cole (1971), and sub-
sequently improved by Murman (1974), of using type-
dependent differencing, with central-difference formulas in
the subsonic zone, where the governing equation is ellip-
tic, and upwind difference formulas in the supersonic zone,
14 Aerodynamics
where it is hyperbolic. The resulting directional bias in
the numerical scheme corresponds to the upwind region
of dependence of the ow in the supersonic zone. If we
consider the transonic ow past a prole with fore-and-
aft symmetry such as an ellipse, the desired solution of
the potential ow equation is not symmetric. Instead it
exhibits a smooth acceleration over the front half of the
prole, followed by a discontinuous compression through
a shock wave. However, the solution of the potential ow
equation (42) is invariant under a reversal of the velocity
vector, u
i
=
x
i
. Corresponding to the solution with a
compression shock, there is a reverse ow solution with an
expansion shock, as illustrated in Figure 9. In the absence
of a directional bias in the numerical scheme, the fore-
and-aft symmetry would be preserved in any solution that
could be obtained, resulting in the appearance of improper
discontinuities.
Since the quasilinear form does not distinguish between
conservation of mass and momentum, difference approxi-
mations to it will not necessarily yield solutions that satisfy
the jump condition unless shock waves are detected and
special difference formulas are used in their vicinity. If
we treat the conservation law (47), on the other hand, and
preserve the conservation form in the difference approxi-
mation, we can ensure that the proper jump condition is
satised. Similarly, we can obtain proper solutions of the
small-disturbance equation by treating it in the conservation
form.
The general method of constructing a difference approx-
imation to a conservation law of the form
f
x
+g
y
= 0
is to preserve the ux balance in each cell, as illustrated in
Figure 10. This leads to a scheme of the form
F
i+
1
2
,j
F
i
1
2
,j
x
+
G
i,j+
1
2
G
i,j
1
2
y
= 0 (50)
where F and G should converge to f and g in the limit as
the mesh width tends to zero. Suppose, for example, that
(a) (b) (c)
Figure 9. Alternative solutions for an ellipse. (a) Compression
shock, (b) expansion shock, (c) symmetric shock.
g
i, j +1/2
g
i, j 1/2
f
i +1/2, j
f
i 1/2, j
Figure 10. Flux balance of difference scheme in conservation
form. A color version of this image is available at http://www.
mrw.interscience.wiley.com/ecm
equation (50) represents the conservation law (47). Then
on multiplying by a test function
ij
and summing by
parts, there results an approximation to the integral (49).
Thus, the condition for a proper weak solution is satised.
Some latitude is allowed in the denitions of F and G,
since it is only necessary that F = f +O(x) and G =
g +O(x). In constructing a difference approximation,
we can therefore introduce an articial viscosity of the
form
P
x
+
Q
y
provided that P and Q are of order x. Then, the difference
scheme is an approximation to the modied conservation
law

x
(f +P) +

y
(g +Q) = 0
which reduces to the original conservation law in the limit
as the mesh width tends to zero.
This formulation provides a guideline for constructing
type-dependent difference schemes in conservation form.
The dominant term in the discretization error introduced
by the upwind differencing can be regarded as an articial
viscosity. We can, however, turn this idea around. Instead
of using a switch in the difference scheme to introduce
an articial viscosity, we can explicitly add an articial
viscosity, which produces an upwind bias in the difference
Aerodynamics 15
scheme at supersonic points. Suppose that we have a
central-difference approximation to the differential equation
in conservation form. Then the conservation form will
be preserved as long as the added viscosity is also in
conservation form. The effect of the viscosity is simply
to alter the conserved quantities by terms proportional to
the mesh width x, which vanish in the limit as the
mesh width approaches zero, with the result that the proper
jump conditions must be satised. By including a switching
function in the viscosity to make it vanish in the subsonic
zone, we can continue to obtain the sharp representation
of shock waves that results from switching the difference
scheme.
There remains the problem of nding a convergent iter-
ative scheme for solving the nonlinear difference equations
that result from the discretization. Suppose that in the
(n +1)st cycle the residual R
ij
at the point ix, jy is
evaluated by inserting the result
(n)
ij
of the nth cycle in
the difference approximation. Then, the correction C
ij
=

(n+1)
ij

(n)
ij
is to be calculated by solving an equation of
the form
NC +R = 0 (51)
where N is a discrete linear operator and is a scaling
function. In a relaxation method, N is restricted to a lower
triangular or block triangular form so that the elements
of C can be determined sequentially. In the analysis of
such a scheme, it is helpful to introduce a time-dependent
analogy. The residual R is an approximation to L, where
L is the operator appearing in the differential equation.
If we consider C as representing t
t
, where t is an
articial time coordinate, and Nt is an approximation
to a differential operator D, then equation (51) is an
approximation to
D
t
+L = 0 (52)
Thus, we should choose N so that this is a convergent
time-dependent process.
With this approach, the formulation of a relaxation
method for solving a transonic ow is reduced to three
main steps.
Construct a central-difference approximation to the
differential equation.
Add a numerical viscosity to produce the desired
directional bias in the hyperbolic region.
Add time-dependent terms to embed the steady state
equation in a convergent time-dependent process.
Methods constructed along these lines have proved
extremely reliable. Their main shortcoming is a rather slow
rate of convergence. In order to speed up the convergence,
we can extend the class of permissible operators N.
3.3 Solution of the transonic small-disturbance
equation
3.3.1 Murman difference scheme
The basic ideas can conveniently be illustrated by consider-
ing the solution of the transonic small-disturbance equation
(Ashley and Landahl, 1965)
_
1 M
2

( +1)M
2

x
_

xx
+
yy
= 0 (53)
The treatment of the small-disturbance equation is simpli-
ed by the fact that the characteristics are locally symmetric
about the x direction. Thus, the desired directional bias can
be introduced simply by switching to upwind differenc-
ing in the x direction at all supersonic points. To preserve
the conservation form, some care must be exercised in the
method of switching as illustrated in Figure 11. Let p
ij
be
A > 0:
Central
differencing
A < 0:
Upwind
differencing
Figure 11. MurmanCole difference scheme: A
xx
+
yy
= 0.
A color version of this image is available at http://www.mrw.
interscience.wiley.com/ecm
16 Aerodynamics
a central-difference approximation to the x derivatives at
the point ix, jy:
p
ij
= (1 M
2

i+1,j

ij
(
ij

i1,j
)
x
2
( +1)M
2

(
i+1,j

ij
)
2
(
ij

i1,j
)
2
2x
3
= A
ij

i+1,j
2
ij
+
i1,j
x
2
(54)
where
A
ij
= (1 M
2

) ( +1)M
2

i+1,j

i1,j
2x
(55)
Also, let q
ij
be a central-difference approximation to
yy
:
q
ij
=

i,j+1
2
ij
+
i,j1
y
2
(56)
Dene a switching function with the value unity at
supersonic points and zero at subsonic points:

ij
= 0 if A
ij
> 0;
ij
= 1 if A
ij
< 0 (57)
Then, the original scheme of Murman and Cole (1971) can
be written as
p
ij
+q
ij

ij
(p
ij
p
i1,j
) = 0 (58)
Let
P = x

x
_
(1 M
2

)
x

+1
2
M
2

x
2
_
= xA
xx
where A is the nonlinear coefcient dened by equa-
tion (55). Then, the added terms are an approximation to

P
x
= xA
xxx
(59)
This may be regarded as an articial viscosity of order
x, which is added at all points of the supersonic zone.
Since the coefcient A of
xxx
= u
xx
is positive in the
supersonic zone, it can be seen that the articial viscos-
ity includes a term similar to the viscous terms in the
NavierStokes equation.
Since is not constant, the articial viscosity is not
in conservation form, with the result that the difference
scheme does not satisfy the conditions stated in the previous
section for the discrete approximation to converge to a weak
solution satisfying the proper jump conditions. To correct
this, all that is required is to recast the articial viscosity in
a divergence form as (/x)(P). This leads to Murmans
fully conservative scheme (Murman, 1974)
p
ij
+q
ij

ij
p
ij
+
i1,j
p
i1,j
= 0 (60)
At points where the ow enters and leaves the supersonic
zone,
ij
and
i1,j
have different values, leading to
special parabolic and shock-point equations
q
ij
= 0
and
p
ij
+p
i1,j
+q
ij
= 0
With the introduction of these special operators, it can
be veried by directly summing the difference equations at
all points of the ow eld that the correct jump conditions
are satised across an oblique shock wave.
3.3.2 Solution of the difference equations by
relaxation
The nonlinear difference equations (5457, and 58 or 60)
may be solved by a generalization of the line relaxation
method for elliptic equations. At each point we calculate
the coefcient A
ij
and the residual R
ij
by substituting the
result
ij
of the previous cycle in the difference equations.
Then we set
(n+1)
ij
=
(n)
ij
+C
ij
, where the correction C
ij
is determined by solving the linear equations
C
i,j+1
2C
i,j
+C
i,j1
y
2
+(1
i,j
)A
i,j
(2/)C
i,j
+C
i1,j
x
2
+
i1,j
A
i1,j
C
i,j
2C
i1,j
+C
i2,j
x
2
+R
i,j
= 0
(61)
on each successive vertical line. In these equations, is
the overrelaxation factor for subsonic points, with a value
in the range 1 to 2. In a typical line relaxation scheme for
an elliptic equation, provisional values

ij
are determined
on the line x = ix by solving the difference equations
with the latest available values
(n+1)
i1,j
and
(n)
i+1,j
inserted
at points on the adjacent lines. Then, new values
(n+1)
i,j
are
determined by the formula

(n+1)
ij
=
(n)
ij
+(

ij

(n)
ij
)
By eliminating

ij
, we can write the difference equations in
terms of
(n+1)
ij
and
(n)
ij
. Then, it can be seen that
yy
would
Aerodynamics 17
be represented by (1/)
2
y

(n+1)
+[1 (1/)]
2
y

(n)
in
such a process, where
2
y
denotes the second central-
difference operator. The appropriate procedure for treating
the upwind difference formulas in the supersonic zone,
however, is to march in the ow direction, so that the values

(n+1)
ij
on each new column can be calculated from the val-
ues
(n+1)
i2,j
and
(n+1)
i1,j
already determined on the previous
columns. This implies that
yy
should be represented by

2
y

(n+1)
in the supersonic zone, leading to a discontinuity
at the sonic line. The correction formula (61) is derived by
modifying this process to remove this discontinuity. New
values
(n+1)
ij
are used instead of provisional values

ij
to
evaluate
yy
, at both supersonic and subsonic points. At
supersonic points,
xx
is also evaluated using new values.
At subsonic points,
xx
is evaluated from
(n+1)
i1,j
,
(n)
i+1,j
and
a linear combination of
(n+1)
ij
and
(n)
ij
. In the subsonic
zone, the scheme acts like a line relaxation scheme, with a
comparable rate of convergence. In the supersonic zone, it
is equivalent to a marching scheme, once the coefcients
A
ij
have been evaluated. Since the supersonic difference
scheme is implicit, no limit is imposed on the step length
x as A
ij
approaches zero near the sonic line.
3.3.3 Nonunique solutions of the difference
equations for one-dimensional ow
Some of the properties of the Murman difference formulas
are claried by considering a uniform ow in a parallel
channel. Then
yy
= 0, and with a suitable normalization
of the potential, the equation reduces to

x
_

2
x
2
_
= 0 (62)
with and
x
given at x = 0, and given at x = L.
The supersonic zone corresponds to
x
> 0. Since
2
x
is
constant,
x
simply reverses sign at a jump. Provided we
enforce the entropy condition that
x
decreases through
a jump, there is a unique solution with a single jump
whenever
x
(0) > 0 and (0) +L
x
(0) (L) (0)
L
x
(0).
Let u
i+1/2
= (
i+1

i
)/x and u
i
= (u
i+1/2
+
u
i1/2
)/2. Then, the fully conservative difference equations
can be written as
Elliptic:
u
2
i+
1
2
= u
2
i
1
2
when u
i
0 u
i1
0 (a)
Hyperbolic:
u
2
i
1
2
= u
2
i
3
2
when u
i
> 0 u
i1
> 0 (b)
Shock Point:
u
2
i+
1
2
= u
2
i
3
2
when u
i
0 u
i1
> 0 (c)
Parabolic:
0 = 0 when u
i
> 0 u
i1
< 0 (d)
These admit the correct solution, illustrated in Figure 12(a)
with a constant slope on the two sides of the shock.
The shock-point operator allows a single link with an
intermediate slope, corresponding to the shock lying in the
middle of a mesh cell.
The nonconservative difference scheme omits the shock-
point operator, with the result that it admits solutions of the
type illustrated in Figure 12(b), with the shock too far for-
ward and the downstream velocity too close to the sonic
speed (zero with the present normalization). The direct
switch in the difference scheme from (b) to (a) allows a
break in the slope as long as the downstream slope is nega-
tive. The magnitude of the downstream slope cannot exceed
the magnitude of the upstream slope, however, because
then u
i1
< 0, and accordingly the elliptic operator would
be used at the point (i 1)x. Thus, the nonconservative
(a)

Shock point
x = 0 x = L
(b)
x = 0 x = L
Figure 12. One-dimensional ow in a channel. value of
at node i. A color version of this image is available at
http://www.mrw.interscience.wiley.com/ecm
18 Aerodynamics
scheme enforces the weakened shock condition,

i

i2
>
i

i+2
> 0
which allows solutions ranging from the point at which the
downstream velocity is barely subsonic up to the point at
which the shock strength is correct. When the downstream
velocity is too close to sonic speed, there is an increase
in the mass ow. Thus, the nonconservative scheme may
introduce a source at the shock wave.
The fully conservative difference equations also admit,
however, various improper solutions. Figure 13(a) illus-
trates a sawtooth solution with u
2
constant everywhere
except in one cell ahead of a shock point. Figure 13(b)
illustrates another improper solution in which the shock is
too far forward. At the last interior point, there is then an
expansion shock that is admitted by the parabolic operator.
Since the difference equations have more than one root, we
must depend on the iterative scheme to nd the desired root.
The scheme should ideally be designed so that the correct
solution is stable under a small perturbation and improper
solutions are unstable. Using a scheme similar to equa-
tion (61), the instability of the sawtooth solution has been
conrmed in numerical experiments. The solutions with an

Shock point too far forward

Parabolic point
Shock point
(a)
x = 0 x = L
(b)
x = 0 x = L
Figure 13. One-dimensional ow in a channel (a) sawtooth
solution and (b) solution with downstream parabolic point. A color
version of this image is available at http://www.mrw.interscience.
wiley.com/ecm
expansion shock at the downstream boundary are stable,
on the other hand, if the compression shock is too far for-
ward by more than the width of a mesh cell. Thus there is
a continuous range of stable improper solutions, while the
correct solution is an isolated stable equilibrium point.
3.4 Solution of the exact potential ow equation
3.4.1 Difference schemes for the exact potential ow
equation in quasilinear form
It is less easy to construct difference approximations to
the potential ow equation with a correct directional bias,
because the upwind direction is not known in advance. Fol-
lowing Jameson (1974), the required rotation of the upwind
differencing at any particular point can be accomplished by
introducing an auxiliary Cartesian coordinate system that
is locally aligned with the ow at that point. If s and n
denote the local stream-wise and normal directions, then
the transonic potential ow equation becomes
(c
2
q
2
)
ss
+c
2

nn
= 0 (63)
Since u/q and v/q are the local direction cosines,
ss
and
nn
can be expressed in the original coordinate system
as

ss
=
1
q
2
_
u
2

xx
+2uv
xy
+v
2

yy
_
(64)
and

nn
=
1
q
2
_
v
2

xx
2uv
xy
+u
2

yy
_
(65)
Then, at subsonic points, central-difference formulas are
used for both
ss
and
nn
. At supersonic points, central-
difference formulas are used for
nn
, but upwind difference
formulas are used for the second derivatives contributing to

ss
, as illustrated in Figure 14.
At a supersonic point at which u > 0 and v > 0, for
example,
ss
is constructed from the formulas

xx
=

ij
2
i1,j
+
i2,j
x
2

xy
=

ij

i1,j

i,j1
+
i1,j1
xy

yy
=

ij
2
i,j1
+
i,j2
y
2
(66)
It can be seen that the rotated scheme reduces to a form
similar to the scheme of Murman and Cole for the small-
disturbance equation if either u = 0 or v = 0. The upwind
difference formulas can be regarded as approximations
Aerodynamics 19
Characteristic
q

nn

ss
Figure 14. Rotated difference scheme.
to
xx
x
xxx
,
xy
(x/2)
xxy
(y/2)
xyy
, and

yy
y
yyy
. Thus at supersonic points, the scheme intro-
duces an effective articial viscosity
_
1
c
2
q
2
_
_
x
_
u
2
u
xx
+uvv
xx
_
+y
_
uvu
yy
+v
2
v
yy
__
(67)
which is symmetric in x and y.
3.4.2 Difference schemes for the exact potential ow
equation in conservation form
In the construction of a discrete approximation to the
conservation form of the potential ow equation, it is
convenient to accomplish the switch to upwind differencing
by the explicit addition of an articial viscosity. Thus, we
solve an equation of the form
S
ij
+T
ij
= 0 (68)
where T
ij
is the articial viscosity, which is constructed
as an approximation to an expression in divergence form
P/x +Q/y, where P and Q are appropriate quanti-
ties with a magnitude proportional to the mesh width. The
central-difference approximation is constructed in the nat-
ural manner as
S
ij
=
(u)
i+
1
2
,j
(u)
i
1
2
,j
x
+
(v)
i,j+
1
2
(v)
i,j
1
2
y
(69)
Consider rst the case in which the ow in the supersonic
zone is aligned with the x coordinate, so that it is sufcient
to restrict the upwind differencing to the x derivatives. In
a smooth region of the ow, the rst term of S
ij
is an
approximation to

x
(u) =
_
1
u
2
c
2
_

xx

uv
c
2

xy
We wish to construct T
ij
so that
xx
is effectively
represented by an upwind difference formula when u > c.
Dene the switching function
= min
_
0,
_
1
u
2
c
2
__
(70)
Then set
T
ij
=
P
i+
1
2
,j
P
i
1
2
,j
x
(71)
where
P
i+
1
2
,j
=

ij
x
[
i+1,j
2
ij
+
i1,j
(
ij
2
i1,j
+
i2,j
)] (72)
The added terms are an approximation to P/x, where
P = [(1 )x
xx
+x
2

xxx
]
Thus, if = 0, the scheme is rst-order accurate; but if
= 1 x and is a constant, the scheme is second-
order accurate. Also, when = 0 the viscosity cancels the
term (1 u
2
/c
2
)
xx
and replaces it by its value at the
adjacent upwind point.
In this scheme, the switch to upwind differencing is
introduced smoothly because the coefcient 0 as u
c. If the rst term in S
ij
were simply replaced by the upwind
difference formula
(u)
i
1
2
,j
(u)
i
3
2
,j
x
the switch would be less smooth because there would
also be a sudden change in the representation of the term
(uv/c
2
)
xy
, which does not necessarily vanish when u =
c. A scheme of this type proved to be unstable in numerical
tests.
The treatment of ows that are not well aligned with the
coordinate system requires the use of a difference scheme in
which the upwind bias conforms to the local ow direction.
The desired bias can be obtained by modeling the added
terms T
ij
on the articial viscosity of the rotated difference
20 Aerodynamics
scheme for the quasilinear form described in the previous
section. Since equation (47) is equivalent to equation (63)
multiplied by /c
2
, P and Q should be chosen so that
P/x +Q/y contains terms similar to equation (67)
multiplied by /c
2
. The following scheme has proved
successful. Let be a switching function that vanishes in
the subsonic zone:
= max
_
0,
_
1
c
2
q
2
__
(73)
Then, P and Q are dened as approximations to

_
(1 )ux
x
+ux
2

xx
_
and
[(1 )vy
y
+vy
2

yy
]
where the parameter controls the accuracy in the same
way as in the simple scheme. If = 0, the scheme is rst-
order accurate, and at a supersonic point where u > 0 and
v > 0, P then approximates
x
_
1
c
2
q
2
_
u
x
= x

c
2
_
1
c
2
q
2
_
(u
2
u
x
+uvv
x
)
When this formula and the corresponding formula for
Q are inserted in P/x +Q/y, it can be veried
that the terms containing the highest derivatives of are
the same as those in equation (67) multiplied by /c
2
.
In the construction of P and Q, the derivatives of P
are represented by upwind difference formulas. Thus, the
formula for the viscosity nally becomes
T
ij
=
P
i+
1
2
,j
P
i
1
2
,j
x
+
Q
i,j+
1
2
Q
i,j
1
2
y
(74)
where if u
i+1/2,j
> 0, then
P
i+
1
2
,j
= u
i+
1
2
,j

ij
_

i+
1
2
,j

i
1
2
,j

i
1
2
,j

i
3
2
,j
__
and if u
i+1/2,j
< 0, then
P
i+
1
2
,j
= u
i+
1
2
,j

i+1,j
_

i+
1
2
,j

i+
3
2
,j

i+
3
2
,j

i+
5
2
,j
__
while Q
i,j+1/2
is dened by a similar formula.
3.4.3 Analysis of the relaxation method
Both the nonconservative rotated difference scheme and
the difference schemes in conservation form lead to dif-
ference equations that are not amenable to solution by
marching in the supersonic zone, and a rather careful anal-
ysis is needed to ensure the convergence of the iterative
scheme. For this purpose, it is convenient to introduce
the time-dependent analogy proposed in Section 3.2. Thus,
we regard the iterative scheme as an approximation to the
articial time-dependent equation (52). It was shown by
Garabedian (1956) that this method can be used to estimate
the optimum relaxation factor for an elliptic problem.
To illustrate the application of the method, consider
the standard difference scheme for Laplaces equation.
Typically, in a point overrelaxation scheme, a provisional
value

ij
is obtained by solving

(n+1)
i1,j
2

ij
+
(n)
i+1,j
x
2
+

(n+1)
i,j1
2

ij
+
(n)
i,j+1
y
2
= 0
Then the new value
(n+1)
ij
is determined by the formula

(n+1)
ij
=
(n)
ij
+
_

ij

(n)
ij
_
where is the overrelaxation factor. Eliminating

ij
, this is
equivalent to calculating the correction C
ij
=
(n+1)
ij

(n)
ij
by solving

1
(C
ij
C
i1,j
) +
2
(C
ij
C
i,j1
) +
3
C
i,j
= R
ij
(75)
where R
ij
is the residual, and

1
=
1
x
2

2
=
1
y
2

3
=
_
2

1
__
1
x
2
+
1
y
2
_
Equation (75) is an approximation to the wave equation

1
t x
xt
+
2
t y
yt
+
3
t
t
=
xx
+
yy
This is damped if
3
> 0, and to maximize the rate of
convergence, the relaxation factor should be chosen to
give an optimal amount of damping.
If we consider the potential ow equation (63) at a
subsonic point, these considerations suggest that the scheme
(75), where the residual R
ij
is evaluated from the difference
Aerodynamics 21
approximation described in Section 3.4.1, will converge if

1

c
2
u
2
x
2
,
2

c
2
v
2
y
2
,
3
> 0
Similarly, the scheme

1
(C
ij
C
i1,j
) +
2
(C
i,j+1
2C
ij
+C
i,j1
)
+
3
C
i,j
= R
ij
(76)
which requires the simultaneous solution of the corrections
on each vertical line, can be expected to converge if

1

c
2
u
2
x
2
,
2
=
c
2
v
2
y
2
,
3
> 0
At supersonic points, schemes similar to (75) or (76) are
not necessarily convergent (Jameson, 1974). If we intro-
duce a locally aligned Cartesian coordinate system and
divide through by c
2
, the general form of the equivalent
time-dependent equation is
_
M
2
1
_

ss

nn
+2
st
+2
nt
+
t
= 0 (77)
where M is the local Mach number, and s and n are the
stream-wise and normal directions. The coefcients , ,
and depend on the coefcients of the elements of C on
the left-hand side of (75) and (76). The substitution
T = t
s
M
2
1
+n
reduces this equation to the diagonal form
(M
2
1)
ss

nn

_

2
M
2
1

2
_

T T
+
T
= 0
Since the coefcients of
nn
and
ss
have opposite signs
when M > 1, T cannot be the time-like direction at a super-
sonic point. Instead, either s or n is time-like, depending
on the sign of the coefcient of
T T
. Since s is the time-
like direction of the steady state problem, it ought also to
be the time-like direction of the unsteady problem. Thus,
when M > 1, the relaxation scheme should be designed so
that and satisfy the compatibility condition
>
_
M
2
1 (78)
The characteristics of the unsteady equation (77) satisfy
(M
2
1)(t
2
+2nt ) 2st (s n)
2
= 0
Thus, the characteristic cone touches the s n plane. As
long as condition (78) holds with > 0 and > 0, it
slants upstream in the reverse time direction, as illustrated
in Figure 15. To ensure that the iterative scheme has the
proper region of dependence, the ow eld should be swept
in a direction such that the updated region always includes
the upwind line of tangency between the characteristic cone
and the s n plane.
A von Neumann analysis (Jameson, 1974) indicates that
the coefcient of
t
should be zero at supersonic points,
reecting the fact that t is not a time-like direction. The
mechanism of convergence in the supersonic zone can be
inferred from Figure 15. An equation of the form of (78)
with constant coefcients reaches a steady state because
with advancing time the cone of dependence ceases to
intersect the initial time plane. Instead, it intersects a surface
containing the Cauchy data of the steady state problem.
The rate of convergence is determined by the backward
inclination of the most retarded characteristic
t =
2s
M
2
1
, n =

s
and is maximized by using the smallest permissible coef-
cient for the term in
st
. In the subsonic zone, on the
other hand, the cone of dependence contains the t axis, and
it is important to introduce damping to remove the inuence
of the initial data.
sn plane
Initial guess
Cauchy
data
t
t
Initial guess
sn plane
(a)
(b)
Figure 15. Characteristic cone of equivalent time-dependent
equation. (a) Supersonic, (b) subsonic.
22 Aerodynamics
3.5 Treatment of complex geometric
congurations
An effective approach to the treatment of two-dimensional
ows over complex proles is to map the exterior domain
conformally onto the unit disk (Jameson, 1974). Equa-
tion (47) is then written in polar coordinates as

_
+

r
(r
r
) = 0 (79)
where the modulus h of the mapping function enters only
in the calculation of the density from the velocity
q =

h
(80)
The Kutta condition is enforced by adding circulation
such that = 0 at the trailing edge. This procedure is
very accurate. Figure 16 shows a numerical verication of
Morawetzs theorem that a shock-free transonic ow is an
isolated point, and that arbitrary small changes in boundary
conditions will lead to the appearance of shock waves
(Morawetz, 1956).These calculations were performed by
the authors program o6.
Applications to complex three-dimensional congura-
tions require a more exible method of discretization, such
as that provided by the nite element method. Jameson and
Caughey proposed a scheme using isoparametric bilinear
or trilinear elements (Jameson and Caughey, 1977; Jame-
son, 1978). The discrete equations can most conveniently
be derived from the Bateman variational principle. In the
scheme of Jameson and Caughey, I is approximated as
I =

p
k
V
k
where p
k
is the pressure at the center of the kth cell and
V
k
is its area (or volume), and the discrete equations are
obtained by setting the derivative of I with respect to the
nodal values of potential to zero. Articial viscosity is
added to give an upwind bias in the supersonic zone, and an
iterative scheme is derived by embedding the steady state
equation in an articial time-dependent equation. Several
widely used codes (o27, o28, o30) have been developed
using this scheme. Figure 17 shows a result for a swept
wing.
An alternative approach to the treatment of complex con-
gurations has been developed by Bristeau et al. (1980a);
Bristeau et al. (1980b). Their method uses a least squares
formulation of the problem, together with an iterative
scheme derived with the aid of optimal control theory. The
method could be used in conjunction with a subdivision
into either quadrilaterals or triangles, but in practice trian-
gulations have been used.
The simplest conceivable least squares formulation calls
for the minimization of the objective function
I =
_
S

2
dS
KORN AIRFOIL
MACH .752 ALPHA 0.000
CL .6367 CD .0001 CM .1482
1.60
1.20
.40
.40
1.20
.80
.00
.80
Cp
KORN AIRFOIL
MACH .745 ALPHA 0.000
CL .6196 CD .0003 CM .1452
.00
.40
1.60
1.20
.80
.80
.40
1.20
Cp
KORN AIRFOIL
MACH .750 ALPHA 0.000
CL .6259 CD .0000 CM .1458
.40
1.60
1.20
.00
.80
.40
1.20
.80
Cp
Figure 16. Sensitivity of a shock-free solution.
Aerodynamics 23
2.00
1.60
1.20
0.80
0.40
0.00
0.40
0.80
1.20
C
p
2.00
1.60
1.20
0.80
0.40
0.00
0.40
0.80
1.20
C
p
Onera wing M6 Upper surface pressure Lower surface pressure
Onera wing M6
Mach .840 Yaw 0.000 Alpha 3.060
Onera wing M6
Mach 0.840 Yaw 0.000 Alpha 3.060
Z 0.20 C
L
0.2733 C
D
0.0151
Onera wing M6
Mach 0.840 Yaw 0.000 Alpha 3.060
Z 0.65 C
L
0.2936 C
D
0.0006
. Experiment
. Theory
. Experiment
. Theory + +
Figure 17. Swept wing.
where is the residual of equation (47) and S is the domain
of the calculation. The resulting minimization problem
could be solved by a steepest descent method in which
the potential is repeatedly modied by a correction
proportional to (I/). Such a method would be very
slow. In fact, it simulates a time-dependent equation of
the form

t
= L

L
where L is the differential operator in equation (47), and L

is its adjoint. Much faster convergence can be obtained by

the introduction of a more sophisticated objective function
I =
_
S

2
dS
where the auxiliary function is calculated from

2
= ()
24 Aerodynamics
Let g be the value of (/n) specied on the boundary C
of the domain. Then, this equation can be replaced by the
corresponding variational form
_
S
v dS =
_
S
v dS
_
C
gv dS
which must be satised by for all differentiable test
functions v. This formulation, which is equivalent to the use
of an H
1
norm in Sobolev space, reduces the calculation
of the potential to the solution of an optimal control
problem, with as the control function and as the
state function. It leads to an iterative scheme that calls
for solutions of Poisson equations twice in each cycle.
A further improvement can be realized by the use of a
conjugate gradient method instead of a simple steepest
descent method.
The least squares method in its basic form allows expan-
sion shocks. In early formulations, these were eliminated
by penalty functions. Subsequently, it was found best to
use upwind biasing of the density. The method has been
extended at Avions Marcel Dassault to the treatment of
extremely complex three-dimensional congurations, using
a subdivision of the domain into tetrahedra (Bristeau et al.,
1985).