Lecture Notes D-42

Numerical Methods in Aircraft

Aerodynamics
Part I: Panel Methods TECHNISCHE UNIVERSITEIT
Laboratorium voor
Scheepshydromechan
Archief
Meketweg 2,2628 CD Dell!
ToL 015.786ß73 - Fac 015-781838
March 1993 Prof.ir. J.W. Slooff

- TUDeift Faculty of Aerospace Engineering

Deift University of Technology

Numerical Methods in Aircraft
Aerodynamics
Part I: Panel Methods

Prof.ir. J.W. Slooff

 NUMERICAL METHODS IN AIRCRAFT AERODYNAMICS

CONTENTh Page

O. Preface O-1

1. INTRODUCTION AND OVERVIEW 1-1

1.1 Numerical Aerodynamics; What, Why and How? l-1

1.2 More on 'How' 1-3

1.2.1 Sketch of the Physical Problem l-3
1.2.2 Mathematical-Physical Models i-k
1.2.3 Numerical Models 1-18
1.2.4 Relative Computational Efforts l-21
1.3 Aspects of Application 1-21
1.4 State-of-the-Art 1-23

2. PANEL METHODS 2-1

2.1 Mathematical Foundations 2-1

2.1.1 Basic Equations 2-1

2.1.2 Boundary Conditions and Additional Conditions; 2-5
Uniqueness of' the Boundary Value Problem
2.1.3 Modeling of Circulation/Lift 2-8
2.1.4 Summary of the Boundary Value Problem for Wing-Like bodies 2-20
2.1.5 Models for Other Types of Bodies 2-21
2.1.6 Elementary Solutions of Laplace's Equation 2-24
2.1.7 Green's identity 2-35
2.1.8 Boundary value problems and equivalent integral equations 2-37
2.1.9 General formulation of' boundary value problems in terms of 2-39
source and doublet distributions using Green's identity
2.1.10 Uniqueness of the integral representation derived from 2-42
the general formulation
2.1.11 Restrictions with respect to the choice of source and 2-44

doublet distributions
2.1.12 Summary of well-posed and some ill-posed boundary value 2-45
problems and integral representations
2.1.13 The application of the integral representation to the 2-47

boundary value problem for the flow around an aircraft

configuration

2.2 Numerical Aspects 2-55

2.2.1 General considerations on the discretisation of the 2-55

integral equations
2.2.2 Surface-grid generation ('paneling') 2-66

2.2.3 Local representation of the geometry of a panel 2-68

2.2.4 Representation of source and doublet distributions 2-72

over a panel
Aerodynamic influence coefficients 2-7'-t
2.2.5
2.2.6 Truncation errors; consistent approximations 2-86

2.2.7 Convergence and stability 2-101

2.2.8 Factors affecting the magnitude of the truncation error 2-107

2.2.9 Solution methods for the resulting system of linear 2-108

algebraic equations; computational effort

2.2.10 Post-processing; lift and drag 2-116

2.3 The Modeling of Compressibility (Subsonic) 2-119

2.3.1 The classical Göthert rule 2-120

2.3.2 Application of the Göthert rule in panel methods 2-122

Alternative treatments, mass flux boundary condition 2-124

2.3.3
2.3.4 Summary of linearized models for subsonic 2-125

compressible flow
2.4 Examples of Application (3D) of 'Ordinary Panel Methods 2-129

2.4.1 Status arid role of panel methods in the process of 2-129

aircraft aerodynamic design
2.4.2 Low speed flow 2-129
2.4.3 High speed, subsonic flow 2-131

2.5 Simulation of viscous effects 2-132

2.6 Special Methods 2-136

2.6.1 Methods based on thin-wing theory 2-136

2.6.2 Panel methods for supersonic flow 2-140

2.6.3 Methods for modeling vortex sheet roll-up 2-142

2.6.4 Other methods 2-146

2.7 Inverse Methods 2-146

 O-1

O. PREFACE

The purpose of these lecture notes is to make students in aeronautical

engineering of graduate or final-stage undergraduate level, specialising in
aerodynamics, acquainted with the (relatively young) discipline of NUMERICAL
AERODYNAMICS.
More specifically it is intended to provide the student with knowledge that will
enable, or help him/her to appreciate the possibilities and limitations of the
numerical methods that are currently utilised in aerodynamic design and research
environments. For this purpose the course attempts to provide a balanced
treatment of physical, mathematical, numerical and application aspects.
It is assumed that the student possesses undergraduate-level knowledge of
- physics in general and aerodynamics in particular
- methods of mathematical analysis and calculus, in particular linear algebra,
vector and tensor analysis, ordinary and partial differential equations
- methods of numerical (mathematical) analysis
- informatics.

The contents of the course is of course limited in width and depth, in the sense
that the emphasis is on methods for stationnary subsonic and transonic flows.
The student is encouraged to consult the following literature for further study.

General References

0.1 Moran, J. 'An introduction to theoretical and computational

aerodynamics'
John Wiley & Sons, 1984
ISBN O-471-87491-4

0.2 Katz, J. 'Low Speed Aerodynamics, from Wing Theory to

Plotkin, Panel Methods'
Mc Graw-Hill, Inc., 1991
ISBN O-O7-O5OL44614
 0-2

0.3 Hirsch, C. 'Numerical Computation of Internal and External

Flows'
Vol. 1, Fundamentals of Numerical Discretization
John Wiley & Sons, 1988
ISBN 0-471-91762-1

0.4 Hirsch, C. 'Numerical Computation of Internal and External

Flows'
Vol. 2, Computational Methods for Inviscid and
Viscous Flows
John Wiley & Sons, 1990
ISBN 0-471-92351-6

0.5 Hoffmann, K. 'Computational Fluid Dynamics for Engineers'

Engineering Education System, Austin, Texas, 1989
ISBN 09623731_i4l

0.6 Nixon, D. (editor) 'Transonic Aerodynamics'

Progress in Astronautics and Aeronautics, Vol. 81
AIAA, Wash, 1982
ISBN 0-915928-65-5

0.7 Habashi, W.G. (editor) 'Advances in Computational Transonics'

Recent Advances in Numerical Methods in Fluids,
Voi. 4; Pineridge Press, Swansea, UK; 1985
ISBN 0-906674-28-X

0.8 Henne, P.A. (editor) 'Applied Computational Aerodynamics'

Progress in Astronautics and Aeronautics,
Vol. 125
AIAA, Wash, 1990
ISBN 0-930403-69-X
 l-1

1. INTRODUCTION AND OVERVIEW

1.1. Numerical Aerodynamics; what, why and how?

Every student that is confronted with a new subject will probably immediately
feel the following questions rising within himself:
What is this about?
Why is it there?
How is it done?
To answer these questions completely probably takes a life-time but at least a
full course of lectures.

About 'WHAT?'
Numerical (or computational) Aerodynamics is a sub-discipline of Computational
Fluid Dynamics ('CFD' for short).
CFD is concerned with research, development and application of methods for
constructing numerical approximations to solutions ('numerical solutions') of
the partial differential equations that describe the motion of fluids and gasses
within or around bodies.
Numerical (aircraft) aerodynamics is concerned with the application of CFD in
aircraft aerodynamics.
Numerical Aerodynamics is also
- sitting at a graphics work-station to
generate and inspect geometries
generate and inspect meshes and grids (which form the basis of numerical
discretization) see figs. 1.1 - 1.5.
'post-process' (visualize) the results
- 'digging' in listings of codes to find the 'bug'
- sitting at your desk to
think about
how to attack a problem (make a plan!)
the error you made
call the guy that gave you the wrong data
call the guy that developed the code because it would not run (properly)
write your report
 l-2

Numerical Aerodynamics, as a separate sub-discipline, is generally recognized to

have started around 1970. It owes its existence in the first place to the
spectacular developments in computer technology. Processing speed and memory
increase by a factor 10 in roughly every seven years (fig. 1.9, ref. 1.1).

About 'WHY?'
Numerical Aerodynamics, as compared to windtunnel testing offers complementary
possibilities that. may improve the efficiency of the aerodynamic design process.
In particular it can lead to
- reduction of time and/or cost of design and development
- product/quality improvement
- improved accuracy of performance estimates in the early stages of the design
process (aircraft manufacturers must give performance guarantees to airlines
that 'buy from the drawing board'; hence reduced risk).

COMPETITIVE AIRCRAFT DESIGN IS NO LONGER POSSIBLE WITHOUT NUMERICAL AERO!

(everybody does it!).

About 'HOW?'
In the process of development and application of numerical methods for solving
physical problems we may distinguish the following steps
definition and description of the physical problem
formulation of the corresponding well-posed analytical-mathematical
problem (mathematical-physical model)
formulation of a numerical model that discretises (approximates) the
analytical-mathematical model and its solution which sufficient accuracy
and efficiency
formulation of algorithm(s) (set of calculation rules) for the efficient
execution of the solution process for the numerical model
y) coding of the algorithm(s)
('writing' the computer program (or 'code'))
verification (checking that the code is mathematically and numerically
sound and self-consistent)
validation (checking and establishing the boundaries of the area of
applicability of the code)
 l-3

calibration (tuning of free parameters in the code to improve the

efficiency and/or applicability)
application of the code

The mathematical-physical model plus the numerical model and algorithm(s)

constitute what is generally called a 'flow solver'.

In case (candidate) codes are already available steps iii) to vi) or vii) can be
replaced by 'selection of code'.

In the following paragraphs we will look at some of these steps in a little more
detail.

1.2 More on 'HOW?'

1.2.1 Sketch of Physical Problem

In numerical (aircraft) aerodynamics with emphasis on stationnary subsonic and

transonic flows we focus in particular on the following physical problem (fig.
1.10):

The (air)flow about an airplane configuration (or part of a configuration)

exhibiting, in general
- lift and drag

- friction effects that are limited to a relative thin layer (boundary layer)
adjacent to the surface of the configuration (high Reynolds numbers)
- 'smooth' flow 'separation' at the (sharp) trailing edge of wing-like
components and a thin 'trailing vortex wake'
- compressibility effects
- disturbances that, at least in subsonic flow (Mach number <1,) vanish at a
large distance from the configuration (except, possibly in the vicinity of the
trailing vortex wake).

Some of phenomena sketched above (thin boundary layers, separation at trailing

edge) are typical for attached flows. However, since the performance envelope of
an aircraft is usually determined by flow separation upstream of the trailing
edge we are also interested in the (much more complicated) problem of separated
flows. A particular kind of (controlled) separated flow characterised by the
presence of strong leading-edge vortices is typical for certain types of
military (fighter) aircraft.

1.2.2 Mathematical-Physical Models

Formulation of the mathematical(-physical) model implies

- choice of a system of conservation laws of physics, or Partial Differential
Equations (PDE's) derived there from, such that the physical phenomena that
are considered to be essential are modelled.
- formulation of boundary conditions and additional conditions;
such that a well-posed mathematical problem is obtained with a unique and
physically relevant solution.

The conservation laws or PDE's model the conservation of the flow quantities
mass
momentum (3 components!)
energy
supplemented with the
thermal and caloric equations of state (usually for a perfect gas)
and, for viscous, heat conducting fluids, expressions for
the stress tensor
viscosity coefficient
heat conduction coefficient

Conservation laws for flow quantities U within a volume with bounding surface

S can in the absence of internal and boundary surface sources, be written in the
integral form

f UdQ + d = O
Q S

or in the differential form

(1.1.2)
+ O
 l-5

In eqs. (1.1.1), (1.1.2) represents the amount of quantity crossing the

boundary surface. That is called the flux. In case U is a scalar quantity, such
as mass or energy, the flux f' is a vector. If U is a vector quantity, such as
momentum, the flux is a tensor.
Distinction can be made between fluxes due to transport of fluid (convective
fluxes) and flux due to molecular motion and thermal agitation (diffusive
fluxes). The diffusive flux can also be considered as the remaining flux when
the fluid is at rest; it appears only in the energy equation. Another type of
distinction is between 'inviscid' and 'viscous' fluxes. The viscous fluxes are
due, directly, to viscosity and heat conduction. The inviscid fluxes are those
present in the absence of viscosity and heat conduction.
Convective fluxes contain both inviscid and viscous parts, diffusive fluxes a
'viscous' (heat conduction) part only. The inviscid parts of convective fluxes
are sometimes referred to as advective fluxes.

In fluid mechanics systems of conservation laws or equivalent PDE's are usually

named after the form of the momentum equations (equations of motion) that is
utilized (e.q. Euler eqs., Navier-Stokes eqs.).

A hierarchy of mathematical-physical models cari be distinguished (fig. 1.11) in

which at lower 'sub-ordinate' levels the equations are simplified successively.
The process of simplification implies a loss of physically significant
information. This has two important consequences
Certain physical phenomena are no longer modelled
In a mathematical sense the problem may no longer have a unique solution.
If the problem does not have a unique solution it is generally not suitable for
numerical treatment (the algorithm may not be able to choose between multiple
solutions). Additional information will then have to be supplied in order to be
able to obtain a solution at all (preferably the physically relevant solution in
which we are interested).

Overview of' mathematical-physical models

The most complete description of the flow of a fluid continuum is given by the
Time-dependent Navier-Stokes equations. The table below lists the phenomena
contained by this flow model.
 1-6

Model eq(s) Phenomena that are modelled

Time-dependent pressure, inertia and friction forces

Navier-Stokes convection (advection), diffusion,

(N-S) dissipation, heat conduction

rotation (vorticity)
(Direct simulation) separation (vortex formation)
compressibility (shock waves)
transition
turbulence

The time-dependent Navier-Stokes equations can be written in the form (see eq.

(1.1.2))

au -
= o (1.1.3)
+ v.(
at mv + visc )

or, expressed in Cartesian coordinates, as

au a
p (ri. + H = O (1.1»)
+ - (F.
ax mv + F
visc. )
+ ay
- mv +
(G. G
visc
. ) + az iriv visc J
.

at

Here U represents the (5x1) column vector of conserved quantities (mass,

momentum (3x), energy) and a generalized (5x3) flux vector. F, G and H
'

represent the Cartesian components of F.

In a body-fitted coordinate system , n, Ç they can be written as

au
at
+
- (.mv
a

a
+ visc ) + - (.mv
an
+ visc J + aç- (iTi.
mv + H
visc )
= 0 (1.1.5)

with U = TJ/J

= k1 + G +
(1.1.6)
G = (n F + n G + n H)/J

= + G +
 l-7

where J = a(,n,ç)/a(x,y,z) is the Jacobian (matrix) of the transformation from

Cartesian to body-fitted coordinates.

It is noted that the energy equation for a viscous, heat conducting fluid
expresses that the local rate of change of total energy (kinetic plus internal)
is the result of the sum of transport of total energy through convection,
transport of heat through conduction and the mechanical energy delivered by
viscous forces.
Furthermore, the more detailed structure of the equations implies (ref. 1.2)
that the transformation of kinetic energy into heat through viscous forces
(dissipation) is irreversible. The same holds for the transformation of heat
into kinetic energy through heat conduction. As a consequence of these
irreversible processes the entropy of a (moving) particle of fluid cannot
decrease (2nd law of thermodynamics).

Numerical flow simulation based on the time-dependent N-S eqs. is often referred
to as direct simulation.

The practical and biggest problem of direct simulation is formed by the fact
that almost all flows of aerospace interest are turbulent. Turbulent flows are
characterized by the presence of time and space dependent fluctuations of all
the flow quantities. In many cases the level of' the turbulent fluctuations can
attain 10% or more of the mean values of the flow quantities. The biggest part
of the problem is, however, in the fact that space- and time-scales of the
turbulent fluctuations can vary enormously. The numerical resolution, with
sufficient accuracy, of the small-scale fluctuations in particular is a
formidable computational task. It requires such small space- and time steps that
for all practical purposes the computational effort is prohibatively large for a
long time to come, even for the biggest and fastest supercomputers.

Reynolds -averaging
The next (lower) level in the hierarchy of flow models therefore involves
application of a time averaging process to the turbulent fluctuations in order
to obtain laws for 'mean', time-averaged, turbulent quantities. The time-
averaging is to be done in such a way that the time-dependence of other
phenomena with time scales different from those of turbulence are not destroyed.
 1-8

For compressible flows density-weighted averages of quantities A are introduced

by setting

(1.1.7)
p

with

T/2
A = f A(x,y,z;t+T)d-t (1.1.8)
-T/2

and

A = A + A" (1.1.9)

pA" = 0 (1.1.10)

Hence, the time-average of the fluctuating part is set to zero.

Application of this time-averaging process to the (full) Navier-Stokes equations
leads to the Reynolds-Averaged Navier-Stokes equations (RANS).

The RANS eqs contain all the terms of the original time-dependent N-S eqs
applied to the mean flow plus a number of additional terms. The additional terms
arise as a result of the non-linear character of the N-S equations. They appear
where (vector) products of quantities are to be taken, as a consequence of the
fact that the average of a product of fluctuating quantities,

pA"xA" 0 (1.1.11)

even if

pA" = O

(The average of a product is not the saine as the product of averages).

In the momentum equations the additional terms are interpreted as Reynolds-

stresses (turbulent stresses). Reynolds stresses appear in the equations in a
 l-9

way similar to the viscous stresses. Hence, the RANS eqs can be written in a
form similar to eq. (1.1.4), (1.1.5). In the energy equation additional terms
appear that are related to (turbulent) heat conduction.

Unfortunately the HANS eqs no longer form a closed system of equations; the
number of unknowns is no longer equal to the number of equations because of the
appearance of additional terms of the type (1.1.11). Apparently, the time-
averaging process gives rise to 'loss of information'. This loss of information
must be compensated for by explicitly adding external information from other
(experimental) sources.
The 'art' of turbulence modelling is now to 'close' the system of equations by
relating the additional, turbulent terms (Reynolds stresses) in some way to the
mean flow quantities.

At this point it can be mentioned that there is an intermediate form of

numerical flow simulation based on the (RA)NS equations that is known as Large
Eddy Simulation (LES). In LES, the time-step and grid-scales are chosen such
that the large-scale turbulent phenomena are resolved implicitly by direct
simulation, while sub-grid scale turbulent phenomena are treated through a
turbulence model as in RANS.

The table below lists the flow phenomena that are modelled through the RANS
equations
 1-lo

Model eq(s) Phenomena that are modelled

Time-dependent As N-S,

Navier-Stokes with but only for phenomena with length and time

sub-grid scale scales of the order of computational steps

space/time averaging (grid) and larger
(tLarge Eddy Simulation, sub-grid scale phenomena through explicitly

LES) added turbulence model

Reynolds-Averaged As N-S,

(Time-averaged) but only for phenomena with time scales much

Navier-Stokes larger than time scales of turbulence phenomena

(RANS) transition and turbulence completely through

explicitly added models

It must also be mentioned that there exist two (simplified) subsets of the RANS
equations that are often used in (the emerging) practice. One involves a thin
shear layer approximation leading to the Thin layer N-S equations (TLNS). In
TLNS it is assumed that the dominating influence of the viscous and turbulent
terms come from the gradients transverse to the main flow direction, which would
be appropriate for thin shear layers. In terms of eq.(1.1.6) the TLNS eqs take
the form

aF. aG.
-+
au mv + + p- (H mv + FI . = O (1.1.12)
at a aç visc }

where Ç is in the body-normal direction.

The use of TLNS is justified when the streamwise viscous terms are small and/or
cannot be resolved on the computational grid. The penalty accompanying TLNS is
that strong (adverse) pressure gradients cannot be handled properly.
Another simplified sub-set of RANS is formed by the Parabolised Navier-Stokes
equations (PNS). PNS is based on considerations similar to TLNS but in addition
applies only to steady flow. The PNS equations can be written in the general
form

3F. 3G.
mv +
mv +
3
an inv1i visc )
. =0 (1.1.13)

They offer a reduction of computational effort in case of supersonic flow with

small cross-flow and no streamwise flow separation (slender configurations). In
such conditions the character of the equations is predominantly parabolic, which
opens the way for an efficient, marching-type of solution procedure.

Zonal Modelling

It appears that further simplifications of the (RA)NS equations are not possible
without distinguishing and separating viscous and inviscid parts in the flow
field. It was first recognized by Prandtl (1904) that at high Reynolds numbers
without significant flow separation the (direct) influence of the viscous and
turbulent shear stresses is limited to a thin layer close to the wall (boundary
layer) and that outside these layers the flow behaves as inviscid. In Prandtl's
theory a simplified boundary layer approximation (of the NS equations) suffices
for the determination of the viscous effects, while the (indirect) effects
thereof on the outer inviscid flow can be represented through the concept of
(boundary layer) displacement thickness. More recently Prandtl's theory has been
reformulated in terms of the theory of matched asymptotic expansions (see, e.g.
ref.l.3).
Clearly such zonal modelling requires some form of interaction between the
boundary layer computations and the computation of the outer inviscid flow; the
'inner' and 'outer' solutions must be matched at a common interface: the edge of
the boundary layer.
Prandtl's theory also teaches that for attached flow with thin boundary layers a
fully inviscid approximation is a valid and consistent one for many flow
properties.
 1-12

INVISCID FLOW MODELS

The Euler equations

The most complete inviscid flow model is obtained by setting the viscous and
heat conduction terms in the NS equations equal to zero, which gives (see eqs.

(1.1.4), (1.1.5))

aG. 8H.
-+
au mv +
mv +
mv -O (1 . 1 . 14)
at ax ay az

or

au
at
+ a

a
(.mv ) + - (,mv
an
+ !_
az
(T
) = (1.1.15)

The table below lists the phenomena that are modelled.

Model eqs Phenomena that are modelled

Euler eqs pressure and inertia forces

convection (advecticn)
rotation (vorticity)
compressibility (shcck waves)

It is important to note that by neglecting viscosity and heat conduction we have

lost the 2nd law of thermodynamics and the implicit modelling of separation
(vortex formation), as well as other viscous effects.
Hence we must anticipate that we may have to reintroduce such lost information
in some other, explicit manner to obtain a problem with a unique and physically
relevant solution. An example is the condition of Kutta-Youkowski which states
that the flow must 'separate smoothly' at the (sharp) trailing-edge of airfoils
and wings.
 1-13

The Full Potential Equation

The next lower level of approximation for inviscid flow is obtained by assuming
irrotationality, leading to the Full Potential (FP) equation:

+ .(pU) = o

p = p(U) (1.1.16)

U=
J
The most important aspect of the full potential equation is that it contains
only one dependent variable: the velocity potential . Recalling that the Euler
equations form a system of 5 equations with 5 dependent variables this is,
indeed, a considerable simplification. The price to be paid is a further
reduction of the number of phenomena that are (properly) modelled (see below).

Model eqs Phenomena that are modelled

Full Potential Eq as Euler eqs but

without rotation (vorticity)
shock waves only accurate as long as they
are weak

Linearised Potential equations and Laplace's equation

If, at the bottom of the hierarchy, in addition to the assumption of

irrotationality the flow field is assumed to consist only of a weakly perturbed
uniform flow we obtain the Linearized Potential or Prandtl-Glauert equation
 i-14

xx
+ +
zz
=0
yy

= i -

= - U.x (1.1.17)

Ihi « ü
This is a linear partial differential equation which offers a further
computational advantage as compared to the FP equation.

Model eqs Phenomena that are modelled

linearized as Full Potential Eq but with

Potential Eq compressibility effects modelled only in as

('Praridtl-Glauert Eq) far as they are linear (no shock waves)

For M= O eq.(1.1.17) as well as (1.1.16) reduces to Laplace's equation

,
xx
+$ +4
zz
=0
yy

or (1.1.18)

xx
+ +
zz
=0
yy

Hence Laplace's equation is valid for incompressible flow, irrespective of the

level of flow perturbation.
 1-15

Model eqs Phenomena that are modelled

Laplace Eq as Full Potential Eq but

without compressibility effects

Boundary Layer equations

As mentioned earlier the direct effects of viscous and turbulent stresses at

high Reynolds numbers are confined to a thin layer adjacent to the wall
(boundary layer). In such situation the velocity component normal to the body
can be argued to be much smaller than the components paralel to the wall.
Introducing this assumption into the TLNS equations (1.1.12) leads to the
conclusion that the equation for the normal component of momentum

a a
(puw) + - (puw) + - pw 2 + p + H
a a
(pw) + = O
aç visc
. ) (1.1.19)

reduces to

aç
(1 . 1 . 20)

or p(,iì,ç)

where PeRfl) represents the pressure at the edge between boundary layer arid
outer inviscid flow.
The implication of (1.1.20) is that the pressure no longer appears as a
dependent variable but as a known external driving force which is to be obtained
from an inviscid flow computation.
The set of equations obtained in this way is known as the Boundary Layer (BL)
Equations. In the symbolic notation of eq.(1.1.5), (1.1.6) they take the saine
form as the TLNS equations (1.1.12). However, the vector U of conserved
quantities and generalised flux vector components contain one unknown less (the
pressure p) and less terms.
 l-16

It is further noted that as with the N-S equations we may distinguish between
Time-dependent, LES and Reynolds Averaged Boundary Layer equations.

Model eqs Phenomena that are modelled

Time-Dependent as N-S,

BL Eqs but only on scale of b.l.

(i.e., e.g., no shockwaves) and only in as far
as compatible with pressure/velocity field
of' external inviscid flow

LES BL Eqs as N-S/LES

idem (as Time Dependent BL Eqs)

Re-Averaged as N-S/Reynolds Averaged

BL Eqs idem (as Time Dependent BL Eqs)

A further simplification of the boundary layer equations is obtained if certain

assumptions are made with respect to the velocity profile in the boundary layer.
In that case the BL Eqs may be integrated beforehand in the normal direction
across the boundary layer, leading to the Integral Boundary Layer Equations

Model eqs Phenomena that are modelled

Integral BL Eqs as Re-Ave BL Eqs

but only for selected class of velocity profiles
 1-17

Interaction between boundary layer flow and external inviscid flow

Solving the boundary layer equations is opportune only in conjunction with

solving the equations for the outer inviscid flow. In doing so the boundary
layer solution and outer inviscid flow solution should interact and match at the
edge of the boundary layer.
Matching requires that in the boundary layer computations separation takes place
at the same location as assumed for the outer inviscid flow computation. Broadly
speaking this requirement is more or less automatically met in the case of
attached flows about airfoils and wings in which both the boundary layer flow
and outer inviscid flow separate at the sharp trailing edge. (In the outer
inviscid flow 'separation' is effectuated through the 'Kutta condition'). In
this situation one speaks of 'weak interaction' between boundary layer and outer
inviscid flow.
In case, in the boundary layer computations, separation takes place upstream (or
downstream) of the point where it was assumed to take place in the outer
inviscid flow one speaks of 'strong interaction'. In this case viscous (finite
Reynolds number) effects have a large impact on the circulation and the outer
inviscid flow (with Kutta condition) no longer represents a valid first
approximation.
In the latter situation the boundary layer computations can no longer be excuted
in the classical way with prescribed pressure. However, it appears that they can
still be performed if the pressure is 'relaxed' in such a way that it can
'adjust' to the location of the separation point. (For details see part II.)

App1icabi1it

The various flow models mentioned above all have their own area of applicability
in aircraft aerodynamics. Qualitatively this can be indicated as in fig. 1.12.
The figure illustrates the areas in the a-Mach plane within the flight envelope
of a transport-type aircraft where the various flow models are valid.
The Navier-Stokes equations are valid in the whole of the a-Mach or CL_Mach
plane (provided the representation of turbulence is adequate). The other flow
models have more limited regions of applicability in the sense that the higher
the model ranks in the hierarchy of fig. 1.11 the larger is its region of
applicability.
 i-18

1.2.3. Numerical Models

In order to be able to 'numerically solve' a (system of) PDE('s) for given

boundary conditions we must first DISCRETIZE the problem. That is, the number of
unknowns (in principle infinitely large if the solution is to be known in any
point of the domain in which the solution is sought) must be reduced to a finite
number.
There exist several alternative ways of discretization, each with its own
advantages and drawbacks.

One popular class of discretization methods proceeds as follows

Define values of the function(s) that we wish to determine (e.g. velocity,
potential) on a regular spatial mesh or grid (fig. 1.13).
Approximate the derivatives of function(s) in the grid points by means of
finite differences
The requirements that the differential (or rather difference) equation(s) must
be satisfied in the grid points and that the boundary conditions must be
satisfied in the grid points on the boundary leads to a large (1O to 106)

system of generally non-linear algebraic equations. The matrices corresponding

with (linearisations of) this system of equations are sparse and exhibit a band
structure. Special solution techniques exist that take advantage of this
structure.
Methods of this type are called Finite Difference Methods (FDM).

Instead of the differential form we can also apply the integral form of the
conservation laws directly to a volume element or cell of the space in which the
solution is sought.
For example, in case of the mass conservation law, this leads to (fig. 1.114)

ff1 div p Ç'cIQ = 1f p '.n dS = O (1.2.1)

c e
Q S
e e

in which we recognize Gauss' divergence theorem.

 1-19

The surface integral in (1.2.1) can be approximated numerically. For example by

defining values of the (as yet unknown) flux vector pV at suitably chosen 'nodal
points' in the computational grid and expressing the variation of the integrand
over the cell faces by means of a Taylor series expansion around the nodes.
Satisfying eq (1.2.1) for all volume elements, while satisfying the boundary
conditions for the fluxes in the cell faces at the boundary, one obtains, again,
a system of algebraic equations with a bandstructured matrix.
Methods of this kind are called FINITE VOLUME METHODS(FVM).

The variation of the integrand along the cell faces can also be expressed in
terms of (suitably chosen) polynomials with the function values in the nodes as
parameters. This case is, sometimes, associated with the notion of FINITE
ELEMENT METHODS IFEMJ. However, the notion of finite element methods is
generally restricted to a class of methods that is based on so-called
variational principles.

As an example consider the Laplace equation

div. grad = O
(1.2.2)

It can be shown that solving (1.2.2) is equivalent with minimizing the volume
integral

f!! [grad ]z dQ
Q
(1.2.3)

The volume integral (1.2.3) can be approximated numerically by dividing

Q in
volume elements and choosing suitable (local) polynomials for with the nodal
values of as (unknown) parameters. Minimization of the discretized form of
(1.2.3) leads again to a system of algebraic equations with a band structured
matrix.
FEM's of the type just sketched have their origin in structural analysis.
They
are not very popular in numerical aerodynamics.

The main advantage of FEM's is that, unlike FDM's, they do not

require regular
or structured grids. This means that they are attractive for application
to
irregular geometries (as often occur in structures). Another and
associated
 l-20

advantage is that the manual as well as the computational effort required for
the generation of the unstructured spatial grids is significantly less than for
regular, structured grids. A disadvantage is that the structure of the resulting
matrices is also less regular than in case of FDM's. As a consequence the
computational effort associated with FEM's is generally larger.
FVM's can also be formulated such that they can make use of unstructured grids.
Such 'unstructured-grid' FVM's have gained in importance since the late 1980's.
For a more general introductory discussion on the differences and similarities
between FDM, FVM and FEM see e.g., ref. 1.2.

Certain PDE's, for which elementary solutions are known, can be transformed into
an integral equation (by means of Green's theorem, about which later, in section
2.1.6). The solution of a boundary value problem for such PDE's can then be
expressed in terms of integrals over boundary and volume distributions of
elementary solutions (sources/sinks, vortices) of, as yet, unknown strength;
e.g.

13
= JI oCx.) K(x.,x.)dS + III o(x k
S
)
13
K(x.,x.)dQ (1.2.4)

where o(x.) represents the elementary solution in a point x. and K(x.x) an

influence or distribution function associated with o(x.).
For homogeneous linear PDE's, such as Laplace's equation, the volume integral
disappears, so that a distribution over the bounding surfaces suffices.
The boundary integral can be discretized on a surface mesh; for example as in a
FEM (fig. 1.15).
Requiring that a boundary condition for is satisfied in each element of the

surface mesh leads to a system of (10 to 10') linear algebraic equations for
the unknown source/sink or vortex strengths. In this case the associated matrix
is full ( no zero entries).
Methods of this kind are called BOUNDARY ELEMENT or PANEL METHODS. In Numerical
(Aircraft) Aerodynamics they play an important role since about 1970.
 1-21

1.2.4 Relative computational efforts

Computational aerodynamics methods based on the (Reynolds-averaged) Navier-

Stokes equations, the Euler equations and the Full Potential equation are
usually of the finite difference or finite volume type. Those based on Laplace's
equation are of the 'Panel' type.
Numerical algorithms for the various classes of methods, in particular those for
the Euler and Navier-Stokes equations are (still) being improved continuously.
Fig. 1.16 provides a crude impression, in a relative sense, of the levels of
computational effort required for the different classes of methods. The figure
illustrates that the computational effort associated with panel methods and
RANS-methods differ by several orders of magnitude.

1.3 Aspects of application

Computational aerodynamics methods can also be classified according to their

role in the aerodynamic design process. On may distinguish:
DIRECT (or ANALYSIS)METHODS
developed for the purpose of computing the flow around a configuration of
given geometry (digital electronic windtunnel).
INVERSE METHODS
developed for the purpose of computing the (detailed) geometry required to
generate a given pressure distribution (design).
(AERODYNAMIC) OPTIMIZATION METHODS
('flow' solver plus numerical optimization algorithm)
developed for the purpose of computing the (detailed) geometry required to
obtain given aerodynamic characteristics (e.g. minimum drag).
For obvious reasons, cat. (1) methods are also called ANALYSIS methods while
cat. (2) and (3) are called DESIGN methods. Clearly, the
complementary pos-
sibilities of CFD (relative to those of windtunnel testing) are most
pronounced
for cat (2) and (3) methods (Ref. 1.3).
 l-22

The COMPUTER POWER required by the various classes of methods varies strongly
and depends on
- mathematical/physical model (type of PDE's)
- numerical model
- required numerical accuracy of solution
- turn-around time required

Fig. 1.17 gives an indication of the computer time required by various models
for one flow computation for a wing-fuselage configuration. A modest
'engineering' level of accuracy has been assumed and several types of existing
computers are considered.
Fig. 1.18 illustrates the computer capacity, in terms of processing speed and
memory required for a 'turn-around' time of about halve an hour, again assuming
a modest engineering accuracy.

Apart from a 'flow solver' the application of CFD requires extensive facilities
for PRE- and POSTPROCESSING such as
- GEOMETRY HANDLING
- GRID GENERATION
- GRAPHICS for VISUALIZATION of input (geometry, grids) and output (numerical
flow visualization)
The complexity of the complete CFD process requires an INFORMATION SYSTEMS
approach involving METHOD BASE and DATA BASE MANAGEMENT, EXECUTIVE and GRAPHICS
SOFTWARE, etc (fig. 1.19).
The COSTS associated with CFD development and application are substantial.
For instance, the investment associated with the development of a professional,
production oriented computer program is of the order of one to several million
US $.
The direct computation cost involved with a single flow computation (3D) may
vary between a few hundreds and several thousands US $.
 l-23

1 .4 State-of-the-art

The current (1992) status of numerical aerodynamics (in the Netherlands in

particular) can be summarized as follows. Routine applications (in industry) of
- Panel methods for subsonic flow (including complex 3D configurations)
- 'Full Potential' FD/FV methods for transonic flow (2D and simple 3D
configurations (wing-body))
- Boundary layer methods; on a basis of, both, 'weak' (2D and 3D) and 'strong'
(2D) interaction.
- 'Euler' FV methods for 2D and 3D flows with rotation.
In stage of development
- 3D boundary layer methods on a basis of strong interaction
- Reynolds-averaged Navier-Stokes methods.

A substantial portion of the current development effort is devoted to improved,

more efficient grid generation techniques.

FUTURE PROSPECTS
Because of continuing developments in computer technology, informatics and
numerical mathematics, Numerical Aerodynamics still has an enormous growth
potential.

Rapidly increasing possibilities for the numerical simulation of complex viscous

flows stress the need for improved turbulence models.

The need for both 'simple' (that is fast and cheap) and accurate (that is
computationally intensive and expensive) methods will remain, due to different
requirements with respect to processing speed and accuracy in the various stages
of the aerodynamic design process.
With the objective of' improving design integration computational aerodynamics
methods will, increasingly, be integrated with computational models from other
disciplines such as structures, flight mechanics, etc. (ref. 1.1).
 l-24

References

1.1 Holst, T.L. The NASA Computational Aerosciences Program -

Salas, M.D. Toward Teraf lop Computing AIAA-92-0558, Jan. 1992.
Claus, R.W.

1.2 Hirsch, C. Numerical computation of Internal and External Flows,

Vol. 1. John Wiley & Sons, 1988 (ISBN 0 47191762 1).

1.3 Van Dyke, M.D. Perturbation Methods in Fluid Mechanics, Parabolic

Press, Stanford, 1975.

1.4 Slooff, J.W. Windtunnel tests and aerodynamic computations,

thoughts on this use in aerodynamic design AGARD-CP-
210, paper 11, 1976.

1.5 Chapman, D.R. Computational Aerodynamics Development and Outlook.

AIAA J. Vol. 17, Dec. 1979, pp. 1293-1313.
 2-1

2. PANEL METHODS

2.1 Mathematical foundations

Panel methods are numerical methods for the solution of partial differential
equations governing potential flow, i.e. inviscid, irrotational flow. They are
based on surface distributions of sources and vortices or doublets. Within the
restriction of potential flows such methods can be applied to compute the flow
past complex geometries in two and three dimensions. These geometries are
approximated by a large number of surface elements or 'panels'. Depending on the
number of panels any degree of accuracy may be obtained in principle.

2.1.1 Basic equations

We restrict ourselves to steady incompressible flow and to steady compressible

(subsonic and supersonic) flow subject to small perturbations.
Steady, incompressible potential flow is governed by Laplace's equation,

xx
+
yy
+
zz
=0 (2.1.1)

where is the total velocity potential defined by the velocity vector V;

V = grad (2.1.2)

and the subscripts denote double differentiation with respect to x,y or z, the
coordinates of a right handed coordinate system. For steady compressible flow
the linearized potential equation may be written in terms of the perturbation
potential ' as

(1 - M2) + + = o (2.1.3)
xx yy zz

where is defined by the free stream velocity U and the total velocity
potential as

= Ux + (2.1.4)
 2-2

and M is the free stream Mach number. The linearized potential equation is also
called the Prandtl-Glauert equation. The potential equations (2.1.1) and (2.1.3)
follow from the law of mass conservation

div (pg) = 0 (2.1.5)

through the assumption of irrotationality. The latter is valid for homentropic

flow and is expressed by eq. (2.1.2). For incompressible flow in which the
density p is constant we find

div grad = 0 (2.1.6)

which is the same as Laplace's equation, eq. (2.1.1). In terms of the

perturbation potential given in eq. (2.i.14) we obtain

+ + = 0 (2.1.7)
xx yy zz

If the potential ' or is known, the pressure may be obtained from the momentum
equation. For incompressible flow this takes the form of the Bernoulli equation

p+ 1
p
(2 + + 2) = + ipU2
2
(2.1.8)
x y z

PP,
or, in terms of the pressure coefficient C -
pU2

+ 2 + 2 -
IvIz (2.1.9)
1 = 1
p U2 U_

In compressible flow with small perturbations where the perturbation velocity

components are written as

<< U (2.1.10)
u,v,w = , ,
X y ¿

the density in isentropic flow is given by (ref. 2.1)

1 (U+ )2 + +
X y z
= = [i M2 [ U2
(2.1.11)
 2-3

the components of p'? in the continuity equation (2.1.5) are then expressed in a
series expansion of 4) 4) , 4) as
X y z

4)

= p U [i + (1-M2) +
Uo

4)

p
y = p U [-y + (2.1.12)

4)
z
p
Z
= pU U[+
o

Substituting these expressions into eq. (2.1.5) we obtain

. [(i-w)
ax
4)

x]
+ay. [4)
j + 3z [ J =o (2.1.13)

which is the conservation form; in non-conservation form the Prandti-Giauert

equation, eq. (2.1.3) is obtained. Eq. (2.1.3) is valid for subsonic (M<1) as
well as supersonic (M>1) flow. For the present we consider only subsonic flow.
The pressure is obtained from eq. (2.1.11) as

)2 + 4)2 4)2
(U o + 4) +
= (i _iM2 I L
Y z]T-1
(2.1.14)
Po O

or from

- 1
(2.1.15)
P 1pU 1M2 P0
2 2°

Note: In the remaining of this chapter we will understand by and 4) the non-
dimensionalized quantities /U and 4)/U.
 2-4

Expanding eq. (2.1.14) in a power series

P=1_Mz
p
(2+ ..) 2 U (2.1.16)

we obtain the approximation for the pressure coeficient,

C = - 2 - + (2.1.17)
p U

Since Laplace's equation is the limit for M40 of the Prandtl-Glauert equation,
eq. (2.1.3), the latter has a larger range of applicability than the first.
Laplace's equation is valid for incompressible, irrotational flow, whereas the
Prandtl-Glauert equation may be applied in compressible, irrotational flow
subject to small perturbations and in incompressible, irrotational flow with no
restrictions regarding the perturbations. Therefore most panel methods are based
on the Prandtl-Glauert equation. In general the pressure is computed fom eqs.
(2.1.14) and (2.1.15), which reduce to the exact expressions (2.1.8) and
(2.1.9), respectively, for M+O. Note that this feature is not present if
instead of eqs. (2.1.14) and (2.1.15), eqs. (2.1.16) and (2.1.17) are used.

The Prandtl-Glauert equation (2.1.3) may be transformed into Laplace's equation

xx V V
+
yy V V
+
zz'
=0 (2.1.18)

by the transformation

X' =

=y (2.1.19)

= z

or by

X' = X
(2.1.20)
=

=
 2-5

where = 1_M2.

In the following we therefore confine ourselves first to the treatment of panel

methods for Laplace's equation; the description to be given may be transferred
directly to the Prandtl-Glauert equation through the (inverse) transformation.
We will return to the subject of compressibility effects in section 2.3.3.
Both Laplace's and the Prandtl-Glauert equation equations are linear partial
differential equations, which allow superposition of solutions to create new
solutions. The latter property is utilized extensively in potential flow theory
and as a consequence also in panel methods.

2.1.2 Boundary conditions and additional conditions; uniqueness of the boundary

value problem*)

Types of problems

In order to obtain a specific solution for the potential of Laplace's equation

one has to impose boundary conditions. The latter have to be chosen in such a
way that the solution represents the desired flow situation (in our case the
flow past an aircraft or part of it).
For Laplace's equation two types of problems may be distinguished in terms of
topology:
- problem for an interior domain (internal flow), fig. 2.1.1
- problem for an exterior domain (external flow), fig. 2.1.2
Also two types of boundary conditions exist:

- Neumann-type boundary conditions, expressing that the normal derivative of the

potential on the boundary, (which is equal to the normal velocity component
(=.n)), is prescribed, fig. 2.1.3. A problem involving this type of
boundary conditions only is called a Neumann-problem. The tangential velocity
component at the boundary and thus the pressure distribution may be calculated
when the solution of the boundary value problem is known. The Neumann-problem
arises if one wishes to compute ('analyse') the flow about a body of given
geometry. It is therefore called an analysis problem

*) General literature: ref. 2.2, 2.3

 2-6

- Dirichiet-type of boundary conditions, expressing that the potential on the

boundary is prescribed, fig. 2.l.k. Prescribing the potential implies
prescribing the tangential velocity component by the derivatives and

where s,t are coordinates along the boundary. In this case the normal velocity
component is obtained from the solution. The Dirichlet problem is related to
the problem of designing the shape of a body for a given velocity or pressure
distribution.

Well-posed problems

A condition for (numerically) solving a boundary value problem is that the

problem is well-posed or properly-posed (ref. 2.3), by which we understand that
a continuous and twice differentiable (C2) solution exists;
at any point in the domain where that solution is sought it is unique;
there exist neighbouring solutions, i.e. for slightly different boundary
conditions the solution changes only slightly.

The properties with respect to well-posedness of some typical Dirichlet- and

Neumann problems are summarized below. For a proof one is referred to text books
on partial differential equations or to ref. 2.2.
Well-posed is a
- Dirichlet-problem for an interior domain Q bounded by a closed surface S (not
intersecting itself), fig. 2.1.5.
- Dirichlet-problem for an exterior domain, fig. 2.1.6. This problem may be
obtained from the problem shown in fig. 2.l.6b by letting S*.
Not well-posed is a
- Neumann-problem for an interior domain bounded by a closed surface, fig.
2.l.7a. It can be shown that there exists a solution only if the additional
condition

(2.1.21)
1f an
dS = O
s

is satisfied, indicating that inside the closed surface S no mass is created

or destroyed.
The solution for the potential that exists inside S under the condition
(2.1.21) is not unique, because an arbitrary constant, say, satisfying
c
 2-7

Laplace's equation, may be superimposed without violating the boundary

condition = O (grad = O).
an c

- Neumann-problem for an exterior domain. This problem may be obtained from that
depicted in fig. 2.l.7b for S9. The not well-posedness of the Neumann-
problem for an exterior domain with = O on S. (fig. 2.1.7b) and S* seems,
at first sight, to constitute a severe handicap for applications to aircraft
aerodynamics where we are particularly interested in such problems. However,
we shall see shortly that this is not really the case.

The mixed or Poincar problem for an interior domain bounded by a closed surface
(fig. 2.1.8) may be shown to be well-posed. No additional condition for j needs
to be prescribed. The mass flow imposed on a part of the boundary by prescribing
is absorbed by the part where is prescribed. The latter guarantees
uniqueness of the potential and as a consequence the uniqueness of the problem.
Note that in order to satisfy Laplace's equation must be twice continuously
differentiable in points of the boundary (A and B in fig. 2.1.8) where the type
of the boundary conditions changes. Because of the relevance for the aerodynamic
computations encountered in the aeronautical design problems we will analyse the
mixed boundary value problem in some more detail.

Consider the earlier mentioned case of a multiply connected domain shown in

figs. 2.l.5b and 2.1.7b, but now as a mixed boundary value problem with a
Neumann boundary condition on S and a Dirichiet condition on S0, see fig.
2.1.9. Independent of the choice of on S, any solution satisfies the Neumann
condition (for example = O) on S..
i
The question now is how to select the potential on S in order to have the
desired model that simulates the flow around the closed surface S., i.e. how we
have to select = g(s,t) on S
where g(s,t) is a known function of the
coordinates s,t on S. To answer this question we consider the physics of the
flow past an aircraft configuration (fig. 1.10). From that we learn that at
large distances from S. the flow conditions approach those of the undisturbed
flow. In other words the normalized potential and perturbation potential at
S = S behave like
o
 2-8

= g(s,t) 4 X

or for S (S )
(2.1.22)
o
(4
s0 40

Whether this condition holds everywhere on S remains to be analysed. For

example nothing has been said about the trailing vortex wake of the
configuration.
The mixed boundary value problem with a Neumann boundary condition at the body

surface 5B = Si and a Dirichlet condition expressing zero perturbation potential

at infinity appears to be the appropriate problem to solve for our applications.
However, there is at least one important shortcoming, as we will discuss in the
next section.

2.1.3 Modeling of circulation/lift

In order to illustrate the shortcoming of the flow model mentioned in the

previous section we introduce a plane Q intersecting SB and S, fig. 2.1.10;
the intersections are aS and aS, respectively. We determine the line integral
along the curve as from A to B (fig. 2.1.10) as

B
!! ds = (2.1.23)
as
A

The double brackets ] denote a jump in the value of the quantity inside.
Realizing that a potential flow does not contain rotation, eq. (2.1.23) may be
recognized as the expression for the circulation r around a cross-section 35B of
(see text books on elementary aerodynamics).
5B
For a continuous potential along as (as well as 3SB) we have

r = hm 0 (2.1.24)
B9A

In other words there is no circulation and therefore no lift.

The only way to introduce circulation is to allow for a discontinuity in the
potential between A and B. However, this cannot be done without due considera-
tion, since it would be incompatible with the requirement that (s,t) is

continuous on 5. This inconvenience is circumvented by introducing a cut asc

 2-9

between the points A and B on as and the cross-section aSB of the body. Across
the cut the potential is allowed to jump from a level to according to
2'

= h(s,t) (on Sc) (2.1.25)

where we have formally assumed that the jump is a function h(s,t) of the local
coordinates. In fig. 2.1.11 a sketch is given of the model representation with a
cut. The circulation around the sectìon 35B of' fig. 2.1.11 is then

(F)as L (2.1.26)

Note that the introduction of the cut has reduced the original multiply-
connected domain into a simply-connected domain. Without proof we mention here
that a well-posed problem with circulation can only be formulated for a simply-
connected domain. Every cross-section of the type eQ needs a cut to produce a
well-posed problem with circulation; in three dimensions we thus obtain a

discontinuity surface Sc extending from the body SB to the surface at 'infinity'

S, see fig. 2.1.12. The discontinuity surface intersects the surface S at a
line s (fig. 2.1.12). In general it will not be possible for both the upper and
lower side of S to blend smoothly with the body SB. thus violating the
requirement that the bounding surface of Q must not intersect itself for a
properly posed boundary value problem (section 2.1.2). As a consequence we may
expect that special measures must be taken in order to cope with this situation,
about which later.

We may further notice, that since (s,t) should be continuous over the bounding
surfaces, we must require, at the edges of the discontinuity surface Sc, that

ir = o (2.1.27)

In addition it may be clear that because of the jump in the potential across Sc
we have to abandon the Dirichlet condition (2.1.22) on S in the vicinity of s;
we come back on this later.
 2-10

Boundary conditions on the discontinuity surface

The requirement to have a jump in the potential across S is insufficient to

determine the solution of the mixed boundary value problem uniquely. For a
complete description, either 4 should be prescribed on both sides of Sc, or the
normal derivative on Sc, should be given, or an equivalent of both. How can

we achieve this?
The physics of the flow require conservation of mass, momentum and energy. This
should also hold across S.
Conservation of mass across SC involves that no mass is created nor destroyed,
which can be stated as

=o

or, in case of incompressible flow

° (2. 1 .28)
=

(normal velocity component is continuous across Sc).

It can be shown further (see, e.g., chapter 3) that in an irrotational flow
which is uniform at infinity upstreai conservation of energy is satisfied
implicitly when conservation of mass is satisfied. Hence, it remains to consider
conservation of momentum.
Conservation of momentum across S may be formulated as

[p + (-) a a
]j
=o (normal momentum) (2.2.29a)

and

a
LE aíp) j:B1 = o (2. 2. 29b)

(tangential momentum)
2

LE (p) j Ill
a 3
0 (2.2.29c)

Using (2.1.28) this can be reduced to

(2.2. 30a)
= o
 2-il

and

a a
P =o (2.2.30b)

a a
p j =o (2.2.30c)

respectively.
Since

a a

a a
=

and because [] is not necessarily constant on S (if it were, it would have to

be zero, because of (2.1.27), it follows that (2.2.30 b,c) can be satisfied only
if both

(p Jl = (2.2.31a)

and

=
an2 (2 .2. 31b)

We note that eqs (2.2.31 a,b) express that S must be a stream surface.
Unfortunately we do not know the exact shape and position of a stream surface a
priori, without knowing the solution. What we do know, however, is that if the
flow separates somewhere on the body a stream surface eminates from the line at
which the flow separates into the flow field. In case of sharp-edged bodies,
like a wing the flow is known to separate at the sharp (trailing) edge. This
then suggest that a surface approximating the stream surface from the sharp
(trailing) edge represents a proper a priori choice for the discontinuity
surface S. In case of small perturbations streamlines in the flow will almost
be parallel to the undisturbed flow. In such conditions S may chosen to consist
of straight generators, parallel to the x-axis, eminating from the trailing
edge.
 2-12

With such an a priori chosen shape and position of S the conditions (2.2.30a)
and (2.2.31a and b) are one too many for a properly posed boundary value
problem. In addition we have the complication that the zero pressure jump
condition (2.2.30a) is a non-linear one. The general approach is then to replace
the two conditions (2.2.31 a plus b) by the single (weaker) condition (2.1.28).
This, of course, still guarantees conservation of mass.
Using the Bernoulli eq. (2.1.8) the non-linear zero pressure jump condition
(3.2.30a) can be written as

onS0 (2.1.32)

(which is the saine as ( (II)2 = O or O.

For small perturbations (see eqs. (2.1.16), (2.1.17) this reduces to the
approximation

cP = O on S (2.1.33)

which is equivalent with

= const. in x-direction on S (2.1.34)

The nonlinear character of condition (2.1.32), complicates matters considerably,

therefore, in practice, linear approximations are used of which eq. (2.1.33) or
eq. (2.1.34) is an example.
An alternative may be found in writing eq. (2.1.32) in terms of 'normal' and

tangential coordinates as

+ 2 + = O on Sc (2.1.35)
IL-n s t 1

or, using eq. (2.1.28),

(2.1.36)
lt + = O OflSc

These expressions may be used if we posess a more accurate estimate of the flow
direction on Sc than the free stream direction.
 2-13

The coordinates (s,t) may then be chosen such that f = 0 In this case the
s- direction, denoted by s, coincides with that of the component, in S, of the
average streamline of S (see fig. 2.1.13). Eq. (2.1.36) then reduces to the
linear expression

[L * :11=0 on Sc
s
or
*
it: Ji = const along s

The boundary conditions on S are summarized in fig. 2.1.14.

It is finally noted that matters simplify considerably in the case of two-

dimensional flow. In that case tIe one remaining tangential momentum condition,
as well as the normal momentum (zero pressure jump) condition, are automatically
satisfied when = 0; this irrespective of the shape arid position of S.

Kutta condition

In the preceding section we have shown that a jump in the potential across the
discontinuity plane Sc is necessary in order to model circulation. However, the
magnitude of the jump has not yet been stipulated. It appears to be sufficient
that the jump [ ] is given along SB the intersection of 5c with SB (see fig.
2.1.12); the distribution on Sc then follows from eqs. (2.1.34) or eq. (2.1.38).
considering again the physics of our problem (fig. 1.10) we observe that the
condition that the flow separates 'smoothly' from the (sharp) downstream end SB
of' the body has not been satisfied yet. Here, basing on heuristic, physical
arguments, 'smoothly' is to be interpreted as with finite velocity () that is
continuous when passing from the body onto S across the trailing edge.
This condition is known as the condition of Kutta-Youkowski or simply as the
Kutta condition. The situation suggests that JJ along SB be chosen such that
the Kutta-condition is satisfied.
At this point it should be noted that the more formal argument that Laplace's
equation (2.1.1) or (2.1.6) must be satisfied everywhere in the flowfield, the
line
5b on the bounding surface SB+SC included, leads to the same requirement,
namely that is finite and continuous when passing from onto Sc. Hence, we
 2-14

may conclude that if we make sure that Laplace's equation is satisfied at the
trailing edge we automatically will also satisfy the Kutta condition.
Recall also in this context, that, in section 2.1.3, we anticipated that special
measures were to be taken at the intersection of S with SB in order to restore
the conditions for a properly posed boundary value problem (that is that the
solution is continuous and differentiable twice ("C2") everywhere in the flow
field). Apparently, the Kutta condition serves as such.
Finite and continuous velocity requires in the first place a continuous
potential across the trailing edge, which can be expressed as (see fig. 2.1.15
for the position of the points P)

and . for P1, P, P2, P -* T

(Pj)

or

- (P1) ii: (2.l.39a)

and for P1, P, P2, P + T

+ (P1) (P) + (2.1 . 39b)

Because eq. (2.1.39a) directly involves the jump in potential across S we must
expect that this is the governing condition for determining the circulation.

Continuity of velocity implies continuity of across the trailing edge, thus

for P1, P, P2, P T

(P )
1
or
(P2) - (P) +
for P1, P, P2 P2 + T (2.1.40)

Subtraction of these expressions yields

- (P1) (P) - (P) = = vS- (2.1.41)

 2-15

where the latter equality results from the condition that ff = O across
Sc;V denotes the gradient in S. Eq. (2.1.41) puts a requirement on the jump of
the tangential velocity across S.
Addition of the expressions (2.1.40) results in

+ (P1) 9 (P) + (2.1.42)

which puts a requirement on the average velocity.

Multiplying the upper part of eq. (2.4.41) with + (P) and the lower
part with V(P1) + e(P) we obtain

- [(P1))2 +
(()}2 - (P))2 =
2 p,
2

On the surface of discontinuity we have, according to eq. (2.1.35),

2j
= O.

Hence, it follows that

((P2))2 - [(P1)J2 O (2.1.43)

The condition given by eq. (2.1.43) may also be obtained by requiring continuity
of pressure, because this means (according to the Bernoulli equation, eq.
(2.1.9)), see fig. 2.1.14

+ + 2 = (2 + Z
)
(2 + + 2)
(
x y )
zP2 s tP2 s t nP2
and
(
P1PjP2P2 T
(2 + + + z) (2 + +
x y zP ) (
s tP s t
2)
np ,

J
(2.1.44)

where on the body surface (p1, P2) there holds = O

Taking into account the jump condition, eq. (2.1.35) we find

(25 + 2) (2 + Z)
p' p 1,P24T (2.1.45)
tP

which expresses the same as eq. (2.1.43). These conditions are nonlinear; taking
for s the streamline direction (or a good approximation of it), eq. (2.1.45)
reduces to
 2-16

P1,P2 T (2. 1. 46)

( ;
s -

which, in small perturbation theory, may be approxinated as

4 P1,P2 T (2.1.47)
(o X ) Pl (0
X
)

Eqs. (2.1.143, 146, 47) are sometimes referred to as 'equal pressure' forms of the
Kutta condition. Note that they (pre-)assume that (2.1.39) is already satisfied.

In order to be able to appreciate the implications of (2.1.41) in some further

detail we will consider the individual components of in 'normal' and
'tangential' directions n, s and t respectively, where we choose the s direction
normal (Sn) and t parallel to the trailing edge. Introducing unit vectors n,
and respectively eqs. (2.1.41) may then be written as (see fig. 2.1.15)

a
(Pa) caS 02 (P2) sin ' (Pi)
as 2
n

a a (2.1.48)
j-;: (P2) sin 02 (P2) cas 02 (P
as
n

at
9(P2)a
j

and

cas 01 - ---- (P1) sin (Pi)

fl
n

--- (P1) sin 0

as
n
(P1) cos 01 9- as
n
(Pi) (2.1.49)

at

Because we have -- = O at P1, P2, eqs. (2.1.148/49) may be reduced to

 2-17

(P2) sin (Pa) (a)

2
-ç

---- (P2) cos 02 8s

9- (P) (2.1.50) (b)
n n

at
(P
2
9(P2) J
(c)

and

as (P1) sin 01 4--

an
(Pi) (a)
n

--- (P1) cos O

n
1
9
as n (Pi) (2.1.51) (b)

(E'1) (c)
at ' 1 -: -J

Because the cut surface must be a streamsurface (or at least a good

approximation to it) (p1, P2) must also be zero. It then follows from
(2.1.50a) and (2.l.51a) that

both

(P2) sin 82 0
as
' for p1, p2 T (2.1.52)
and (P1) sin 0 O
--;; J

This can be satisfied only under either one of the following three conditions

17
for Pl, P2 T (2.l.53a)
O
::n (P1)
 2-18

ii) sin 02 = 0
(2.1.54b)
sinO1zO ---(P)90 forP1-*T
as
n

iii) sin 02 z O
as
(P2) 40 for P2 - T
n
sin = O (2.1.55c)

In case i) the component of the flow normal to the trailing edge is stagnant on
both upper and lower surface and the flow may leave the trailing edge in a
direction anywhere between that of the lower and upper surface of the wing.
In case ii) it leaves the wing in a direction parallel to the lower surface and
in case iii) in a direction parallel to the upper surface.
It has been shown by Mangler & Smith (ref. 2.14) that which of the three
possibilities occurs is determined by the sign of the vorticity shed at the
trailing edge (or, in other words, by the sign of and
Unfortunately this is generally not known a priori so that an estimate must be
made of the direction in which the flow leaves the trailing edge. In case of
small trailing-edge angles a suitable estimate is the direction of the trailing
edge angle bisector. A requirement of this kind can be expressed as

(P') = O for P T (2.1.54)

an 2

and is sometimes referred to as the 'directional form of the Kutta condition'

(fig. 2.1.16).
We conclude the discussion on the Kutta condition by noting that, apparently, it
is not possible to, exactly satisfy all (sub-)conditions involved
simultaneously, unless the exact shape and position of the cut/stream surface is
known beforehand. In general, the latter is not the case.

Discontinuity surface, vortex sheet and vortical wake

At this stage we know in general terms how to formulate a well-posed potential

flow problem with circulation. The presence of a surface of discontinuity S
carrying a jump of potential was identified as an essential feature. In addition
 2-19

this surface of discontinuity should be a stream surface or a good approximation

to it.

Further insight into the significance of our mathematical model may be obtained
by considering once more the physics of the real problem, (fig. 1.10). In the
real flow we observe a vortical wake shed from the trailing edge of the body. In
our mathematical model a cross-section approximately perpendicular to the
undisturbed flow direction reveals the intersection of the discontinuity surface
Sc, see fig. 2.1.17. The circulation along a curve t around one of the edges of
S may be written as

B
dt (2.1.55)
= B - =

where A and B are points at the lower and upper side of S, respectively (fig.
2.1.17).
Since z 0, the circulation along t is unequal to zero and we conclude
that there exists vorticity inside the contour t. The only location where this
can be is at the discontinuity surface S, because everywhere else the flow is
irrotational. Hence, the discontiniuty surface is a vortex sheet, with vorticity
along lines of ¡t ji = constant. Apparently, the vortical wake (which is of
viscous origin) in the real flow is represented by the vortex sheet type of
surface of discontinuity in our potential flow model.

At this stage it is worth recalling (see Chapter I) that in stepping down from
the Navier-Stokes equations to more approximate flow models, ranking lower in
the hierarchy of fig. 1.11, essential information on the physics of the complete
flow problem is lost. For inviscid, potential flows such as described by the
Prandtl-Glauert or Laplace equation this concerns, a.o., loss of the ability to
implicitly model rotation (vorticity) and the generation of vorticity.
We may now conclude that we have, explicitly reintroduced the ability to model
rotation (and, through that circulation) through the introduction of the surface
of continuity (or vortex sheet) S and that the Kutta condition now, explicitly,
represents the machenism for generating the correct magnitude of vorticity in
the vortex sheet.
 2-20

2.1.4 Summary of the boundary value problem for wing-like bodies

We can now summarize the incompressible flow past a body as follows;

The flow is modelled by Laplace's equation for the perturbation velocity-
potential 0, defined by eq. (2.1.4), and the following boundary conditions (for
the nomenclature see fig. 2.1.18)

= f(s,t) on SB (2. 1. 56a)

In general the normal velocity component at the body surface is zero, this is
expressed by

=0
an
onS B
or
U -

an
onS B (2. 1. 56b)
Uo

Furthermore

= o on S
and

IL O = constant along an average streamline on S,

The jump on SB (fig. 2.1.18) must be determined in such a way that the
O

Kutta condition eq. 2.1.39/40 plus (eq. (2.1.45/46/147) or eq. (2.1.54) is

satisfied at the trailing edge.
Far from the body on S we have

0=0 on S o (2.1.59)

The surface S (vortex sheet) must be chosen such that it is a streainsurface to

a good approximation.
As has been mentioned earlier the Dirichiet condition O = O at infinity. (Sm)

cannot be satisfied in the vicinity of the intersection of S0 with S. Because

of this reason the part of S infinitely far downstream, denoted by ST in fig.

2.1.18, is treated in a special way, see also fig. 2.1.19.

 2-21

At x * we take
the cross-sectional shape of S invariant in x-direction; it then
becomes a cylindrical surface, usually a plane surface is chosen;
the vortex lines = constant on S parallel to the x-z plane,
resulting in perturbation velocities that have only components in a
plane normal to Sc;
the downstream surface ST normal to Sc;
instead of = 0 (on

=0 onS T (2. 1 . 60)

(Note: this is in agreement with 1, 2 and 3).

In most cases S is chosen as a plane surface parallel to the x-axis, so that

(2.1.61)
anhST - - ax

Bodies with sharp trailing edges (wings, tailplanes, pylons, winglets, stub-
wings, flaps etc.) offer a fixed position for flow separation and thus the
beginning of the vortex sheet is known a priori (of course as long as the
Reynolds number is high enough to guarantee attached flow upto the trailing
edge). Examples of such bodies and their vortex sheets are shown in fig. 2.1.20.

2.15 Models for other types of bodies

For bodies without sharp edges at their downstream end (fuselages, tiptanks,
external stores etc.) the location of flow separation is not known a priori; how
to deal with such cases?
In practice two different approaches are encountered. The first one is based on
the empirical observation that isolated smooth bodies with rounded fore and aft
ends do not generate much lift, at least at small angles of incidence. This
suggests that in such case it is not necessary to adopt in our mathematical
model a jump in the velocity potential at the downstream end of the body, or, in
other words, there is no necessity for a vortex sheet behind the body, see fig.
2.1.21. We then have to solve the boundary value problem with a boundary
 2-22

condition at infinity given by eq. (2.1.22); this problem is depicted in fig.

2.1.9.
The second approach is used when body-lift may not be neglected. Then the
location of flow separation (separation of the boundary layer) is estimated. Two
possibilities may be distinguished:
The flow separates from an open line segment SB (analogous to the separation
from the trailing edge of a wing) as is sketched in fig. 2.1.22.
On the separation line SB the Kutta condition (eqs. (2.1.39/40) plus
(2.1.45/46/47) or (2.1.53)) has to be satisfied.

Note that eqs. (2.1.53a, b, c) require that, in general the vortex sheet must
be tangent to the body at

The flow separation occurs along a closed line, fig. 2.1.23, creating a
closed vortex sheet that subdivides the domain behind the body: an interior
and an exterior domain. In this mathematical-topological model the interior
domain Q. has no physical significance because in reality the flow inside Q.

will be dominated by viscous effects. Hence, in a potential flow model we

have some freedom in the choice of boundary conditions inside Sc.
If we do not allow forces acting on Sc we have to maintain the no-jump
conditions. Also, in order to obtain a well-posed boundary value problem in
the interior domain Q. the potential (Dirichlet condition) has to be imposed
on at least part of the boundary. We have seen further (in section 2.1.4)
that prescribing O on 5T is not possible. Usually one imposes on the
(interior) boundary aQ:

° (2. 1 .62)

and

on ST(3Q.) (2.1.63)

A constant potential on SB (aQ.) implies zero tangential velocity component.

From the Kutta condition (2.1.44) it then follows that
 2-23

(2 + 2
+ (O + 2)
s t
+
nP np
for P1,Pj,P2,P + T (2. 1 .64)
(2 + + 2) (Z
s t np , s
+
tP

where P1,P2 are on SB at the outer and inner side of S, respectively, and Pj,P
are downstream of S8, at the outer and inner side of S, respectively (fig.
2.1.23). Note that here (P1) = O but P2, P) O
With the zero jump condition (2.1.35) it follows further that

(2.1.65)
(
+ tPl nP2

This mathematical flow model implies that the total pressure of the (fictitious)
potential flow in Q. is equal to the total pressure of the external flow.
A more general formulation is obtained if we assume a jump in total
Pt
pressure across S, according to

on S (2.1.66)
= pt

where, as before the derivative * is taken with respect to the average

s
streamline direction, see fig. 2.1.23. In this case the Kutta condition may be
written as

(2) +
(2. 1 .67)
s
+
tP n P'2 Pt i

The magnitude of [ p cannot be determined without knowledge of the real

(viscous) flow inside Q.. The usual assumption in practice that = O is
based on empirical and heuristic considerations (the results such as forces and
moments agree reasonably well with experiments) and also by knowing that the
previously mentioned case a) tends to case b) in the limit of closure of the
separation line.
 2 -2

Jet-engine nacelles, fig. 2.1.24.

In this case the normal velocity component at the entrance of the compressor
(fan face) is controlled by the engine adjustment. The normal velocity component
at the jet boundary is determined by the entrainment effect of the jet. An
alternative for the treatment of the jet boundary may be found in case b
described on p. 2.22.

Flow-through nacelle, fig. 2.1.25.

This problem can be treated as a ringwing.

Wing-body configurations, fig. 2.1.26.

Outside the fuselage the vortex sheet from the wing may be treated as described
earlier for a wing alone. The vortex sheet from the fuselage, if to be modelled,
has been considered in the paragraph on bodies without sharp edges. There is,
however, a complexity at the wing-fuselage junction. To model this flow it is
generally assumed that there is no separation (generation of vorticity) at the
junction. The intersection of the vortex sheet of the wing with the fuselage is
taken as a line = constant, where [ is determined by the Kutta
condition on the wing. This results in a jump in the potential behind the
fuselage (fig. 2.1.26) and therefore the fuselage carries circulation and lift,
eventhough the fuselage itself may not generate a vortex sheet. The
experimentally observed phenomenon of a fuselage that by itself does not
generates lift but does so in the presence of a wing (lift carry-over) has thus
been modelled.

2.1.6 Elementary solutions of Laplace's equation

Laplace's equation is a second-order linear partial differential equation. The

fact that it is linear implies that summations, integrations and
differentiations of particular solutions also form solutions of the equation
(superposition principle). This suggests that a complicated flow problem
number of
governed by Laplace's equation can be synthesized by adding together a
elementary flows which are also governed by Laplace's equation.
 2-25

Sources

One of such elementary solutions is the potential generated in a point P(x,y,z)

by a point source in a point Q(,n,Ç). The strength of the point source is
denoted by o(Q) and the distance from P to Q is given by

r(P,Q)= (x_)2 + (y-n)2 + (z - ç)2

Then, the perturbation potential induced in P is given as

a(Q)
(P)
- 4nr(P,Q)
(2. 1 .68)

The associated (perturbation) velocity field is obtained by taking the gradient

of' (P) with respect to the coordinates x,y,z

-1
= a(Q) VP[4nr(PQ)) (2. 1 .69)

Distribution of sources

For a distribution of sources on a surface S, with the source strength per unit
area given by a(Q), the potential induced at a point P (Q) is evaluated
integrating over the surface S

o(Q)
(P) = If dS (QS) (2.1.70)
4iir(P,Q)

and the velocity by taking the gradient in P of eq. (2.1.70)

-1
= JI o(Q) V dS (2.1.71)
4nr(P,Q)
S

In eqs. (2.1.70,71) the case P-*Q is to be excluded, since then r-0. This problem
is dealt with as follows.
Consider a circular region X with radius e around Q on the surface S (fig.
2.1.27a) and divide the integral (2.1.70) into two parts: one part over the
surface outside the circular region and the other over the circular region
itself. Now let P-*Q together with c-*0, then eq. (2.1.70) may be written as
 2-26

(P) = 11m [ JI d(S-Z) + If -o(Q') d] (2.1.72)

4nr(P,Q') nr(P,Q')
P-Q P.Q S-
c40

where Q' is the 'current' (variable) point on S,X.

The evaluation is done most conveniently if the point P is positioned on the
normal in the point Q on the surface S. The first term of (2.1.72) yields a
finite, in general non-zero, value, that we denote by

dS
14nr(F,Q')

(Cauchy principal value , obtained by excluding the immediate vicinity of Q+P

from the area of integration). The second term is evaluated conveniently by
using polar coordinates (p,0) with respect to Q in , see fig 2.l.27b. Then the
integral becomes

2n C

f j
-o(p,0)
pdpd0 (2.1.73)
0=0 p=Q L11 Jpz+nZ

where n is the distance along the normal on S from Q to P. For c0, n-*0 and a
continuous source strength o(p,O) this integral vanishes.
Thus the following conclusion may be drawn: the potential of a source
distribution on a surface S is continuous across S and takes the form

a(Q) QS (2.1.74)
'(PS) =
i4nr(P,Q)
dS ;

Associated with the continuity of the potential across S is continuity across

S of the derivative of ' along S and thus a continuity of the tangential
velocity on S.

In analogy with the preceding analysis the normal component of the velocity may
be derived as

()p4Q = (n.V)Q = hm [ f! o(Q') an 4r(P,Q' d(S-) +

P-,Q S-

a -1
+ 1f o(Q') an 4nr(P,Q')
dl] (2.1.75)
 2-27

where is the derivative in P along the normal in QS.

P
In this case the integral over I does not vanish. In case P approaches Q from
above (PQ), we denote P by P and we find

!_ (P) = o(Q') dS + o(P); P,Q' tE S (2.l.76a)

an
S
- 4nrQ')

and from below (PtQ), where P is denoted by P, we have

L
an
P) o(Q')
t'4iir(PQ')
dS - o(P ) ; P,Q' S (2.1.76b)
s

Subtracting (2.l.76a) and (2.l.76b) we find

o(P) = L
an
(P) - L
an
(P) =
+
(2.1.77)

Conclusion: The normal velocity component at a surface S induced by a source

distribution on S is discontinuous across S. The magnitude of the jump in a
point of S is given by the source strength per unit area in that point.

A summarizing picture of the velocity field induced by a source distribution is

given by fig. 2.1.28.

In the following we consider some specific examples of source distributions and

their characteristic properties. (The reader is encouraged to keep these in
mind!)

- Uniform source distribution on a sphere with radius R (fig. 2.1.29). Inside

the sphere the potential is found to be constant with value -oR, where o is the
source strength per unit area and R is the radius of the sphere. Outside the
sphere the potential is the same as for a point source with strength _4nR2o at
the centre of the sphere (total surface distribution concentrated in the centre
of the sphere). The perturbation velocity inside the sphere is zero and outside
= _4nR2 o, where r is the distance from a field point to the centre of' the
sphere.
 2-28

It can further be shown that for an arbitrary closed surface there exist a
nonuniform source distribution that, uniquely, induces a constant potential
inside the surface.

- Uniform distribution on an infinite plane surface (fig. 2.1.30) In this case

the perturbution velocity in the field as well as on the surface is directed
along the normal to the surface. The magnitude of the velocity is constant
throughout the whole field including the surface itself, but the direction of
the velocity vector at one side is opposite to that at the other side. The
perturbation velocity may be expressed as = ± o, where z is the direction
normal to the surface and the plus(minus) sign is valid for z >(<) 0.

- Uniform distribution on an infinitely long circular cylinder with radius R

(fig. 2.1.31)
Inside the cylinder the perturbation potential is constant (oR ln R) and there-
fore the perturbation velocity is zero. Outside the cylinder the perturbation
potential is oR ln r, where r is the distance from a field point to the axis of
the cylinder. The potential is the same as if the total source distribution on
the surface of the cylinder were concentrated along the axis.

Doublets
From the source (or sink) type of elementary solution for Laplace's equation
another elementary solution may be derived by differentiation in a given
direction. Such an elementary solution is known as a doublet or dipole (sources
and sinks are also known as monopoles). It follows, that unlike a source, a
doublet posesses directional properties.
The potential in a point P induced by a doublet in a point Q at a distance r is
given by

(P) = p(Q) (2.1.78)

3flQ t4nr(P,Q)1

where flQ denotes the direction of differentiation and ii(Q) is the doublet
strength.
From this the velocity induced by a doublet may be derived as

= P(Q) VP{anQ (2.1.79)

 2-29

where the subscript P in V, again denotes differentiation with respect to the

coordinates of the point P.
The flow field produced by a (three-dimensional) doublet situated in a point Q
may be characterized by fig 2.1.32, which shows the streamlines in an arbitrary
plane through the doublet axis flQ. The direction of iQ is considered to be
positive if it points from the negative to the positive pole. The complete
three-dimensional picture of the flow field induced by the doublet is a series
of streamsurfaces obtained by rotating the streamlines of fig. 2.l.32a about the
doublet axis.

Doublet distributions

Like in the case of monopoles we may consider here also a distribution of

doublets with density p(Q) on a surface S. In particular we consider the
situation in which the doublet axis points in the direction of the normal to the
surface, see fig. 2.1.32b.
The potential in a point P is obtained by integration of (2.1.78) over the
surface S

(P) = fi'p(Q)
a -1 dS
arìQ t4nr(P,Q) (2.1 80)
.

and the velocity from the gradient

= ffp(Q) v{ tITP,Q)11 dS (2.1.81)

As for the source distribution the potential may be evaluated for the limit P-*Q.
Again we have to distinguish two possibilities, PQ or P and PtQ or P. The
result is found to be

(2. 1 .82)

(2.1.83)
 2-30

Thus, we may conclude that the potential generated by a doublet distribution on

a surface S Is discontinous across S, the jump being given by

ENII = (2.1.84)

Note that the jump is given by the negative value of the doublet strength in the
point under consideration. In the case of a source distribution on a surface the
normal velocity component was found to be discontinuous and to jump by an amount
equal to c, see eq. (2.1.76). The minus sign in (2.1.814) is due to the
difference in the differential operator: -p--- =
8flQ

From eqs. (2.1.82) and (2.1.83) we may derive expressions for the tangential
velocity on S, which, consequently is also discontinuous

= i-'(Q) nr1 dS - 5p(P) (2.1.85)

and

a 1
= i(Q)
¡ (- 1) dS + (2.1.86)

Here denotes differentiation tangent to the surface S.

The normal velocity component on the surface is found to be

!__ a a 1
[(FeS)) = dS (2. 1 .87)
anP 8n anQ (
4nr

irrespective of the direction from which the surface is approached. In other

words: the normal velocity component induced b a doublet distribution on a

surface S is continuous across S. A summarizing picture of the perturbation

velocity induced by a doublet distribution on a surface is shown in fig. 2.1.33.

Some particular examples of doublet distributions are given in figs. 2.1.34 and
2.1.35. Fig. 2.1.34 depicts a uniform distribution on a closed body. It may be
shown that inside the body the perturbation potential is constant and equal to
minus the doublet strength per unit area (for an inward normal). Outside the
 2-31

body the perturbation potential is zero. The perturbation velocity is zero

everywhere.
A uniform distribution on an infinite plane surface S (fig. 2.1.35) produces a
perturbation potential that is equal to -i (thus constant) everywhere
i
above S and equal to p (constant) everywhere below S. Therefore, the
perturbation velocity is zero everywhere.

Vorticity distributions and vortices

Consider again the expression (2.1.81) for the perturbation velocity induced by
a doublet distribution. This may also be written as

-1
= if P(Q) VptIflQVQ (i(P Q))) (2. 1 .88)
S,

It turns out that eq. (2.1.88) can be transformed in a different but equivalent
expression that can be associated with yet another type of elementary solution
of Laplace's equation (see also ref. 2.5).
Making use of the following rules of the differential calculus of vectors

Va(ab) + Vb(a.b)

= (bV)a + bx (Vxa) + (aV)b + a x (Vxb) (2. 1 .89)

and

x(x) = V x (.x) x (x)

+ b

= (b7)a - (Va)b + (V.b)a - (a.V)b (2.1.90)

it follows that the essential part of the integrant of (2.1.88) can be written
as
=0
(t)) =
p Q+
=0
+ (1) flQ X
(1
Q P) Q r (2.1.91)
+
 2-32

and that =0 =0
'1
X flX v0()J = Q -
-
n +

_=0
;)JflQ
---1
- InQ.Vp)VQ(_] (2.1.92)

Hence, it follows that

(_l))
= - px [iQX (2.1.93)

so that (2.1.88) can be rewritten as

= - JI p(Q) [flUX VQ()) dS (2.1.94)

or, using iix = -(xii), as

== JI p(Q) (Q(--1 x dS (2.1.95)

The latter can also be written as (tbringing i under the gradient )

V1D = x E
JI
4
)X fl dS - II (Q(P) X flu) dSJ (2.1.96)

The first (surface) integral in (2.1.96) can be reduced to a line integral by

using the integral theorem

If (iixA) dS = f A ds (2.1.97)
S as

which is related to Stokes theorem. (aS representing the curve bounding the

surface S and ds an element of as, (see fig. 2.1.38).

It gives

=
X [ f
as
-+
ds
4nr
[QPxQ)ds] (2. 1.98)
 2-33

Finally, using the identity

x(f) = fx + fxa (2.1.99)

arid re-arranging the result gives

=
rxds 1
- p(Q) + Jf [V (p)xn ) x dS (2. 1. 100)
S r S r

In the first term of (2.1.100) we recognize Blot & Savart's law (ref. 2.1)
applied to a closed line vortex bounding the surface S with strength equal to
the local doublet strength at aS. In case of an unbounded surface this term
vanishes. In the second integral we recognize Blot & Savart's law applied to a
distribution of vorticity on the surface S with the vorticity 'r equal to

:;
= VQ(P) XflQ (2.1.101)

Taking the divergence of (2.1.101) it follows further that

=0
= (px) = cp) = O (2.1.102)

In (2.1.102) we recognize Kelvin-Helmholtz law (ref. 2.2) that in an inviscid

(barotropic) flow vorticity can neither be created nor destroyed.

We can now conclude that

- The velocity field induced by an unbounded doublet distribution with density p

is equal to the velocity field induced by an kunboundedl vorticity
distribution with strength [=J[, which satisfies the Kelvin-Helmholtz
theorem.

- The velocity field induced by a doublet distribution with density p on a

bounded surface S is equal to the velocity distribution induced by a vorticity
distribution on S with strength p plus a line vortex along the boundary of S
having a strength equal to the doublet strength at the boundary.
 2- 314

A more detailed discussion on equivalence of source, doublet and vorticity

distribution can be found in ref. 2.6.

For the velocity induced at a point P by an arbitrary vortex distribution on a

surface S (not necessarily satisfying Kelvin-Helmholtz's theorem) it follows
from Biot and Savart's law (and the preceding analysis) that

ff:;(Q)
(P) = x dS (2.1.103)

In the limit, for P-*Q, we find

= (Q) x dS + x i (QS) (2.1.104)

P 14nr
S
and

= (Q) x dS - x (QES) (2.1.105)

P 14îir
S

The latter two equations involve that the tangential velocity across a vortex
distribution on a surface is discontinuous, the jump in the velocity is equal to
the vortex strength T.

The normal velocity across the vortex distribution is continuous.

Examples of equivalent vortex- and doublet distributions

- A uniform doublet distribution on a closed surface (fig. 2.1.37) is equivalent

to a vortex distribution of zero strength (perturbation velocity zero
everywhere).

- A uniform doublet distribution with strength p per unit area on a bounded

surface S (fig. 2.1.38) is equivalent to a ring vortex of constant strength
along the boundary of S.

- A vortex distribution constant in strength and direction (fig. 2.1.39) is

equivalent to a linear varying doublet distribution.
 2-35

Doublet distributions were seen to satisfy automatically Kelvin-Helmholtz's law

of conservation of vorticity. Therefore, there is a preference for modeling
boundary value problems by means of doublet distributions rather than vortex
distributions.
On the other hand the evaluation and interpretation of the results (velocity
distributions) in terms of (equivalent) vortex distributions has certain
advantages that are related to the singular behaviour of doublet distributions
in the near field; this will be discussed later on.

Due to the inherent singular behaviour of sources, sinks, doublets and vortices,
they are also referred to as singularities and their distributions as
singularity distributions.

In the next section we will see how solutions of' the boundary value problems
discussed in paragraph 2.1.2 can be expressed in terms of singularity
distributions.

2.1.7 Greens's identity

In expressing solutions of Laplace's equation in terms of singularity

distributions an integral theorem due to Green plays a crucial role.
Green's theorem may be formulated as follows. If and are scalar functions
of the spatial coordinates (x,y,z) within a domain Q bounded by a closed surface
S (fig. 2.1.40), then the following integral identity holds

1ff
Q
s
-
s
¿) dQ = Jf (
S
(i') - s
)J dS (2.1.106)

where = V is the Laplace operator.

The proof of Green's theorem (which is related to Gauss's divergence theorem)
may be found in any text-book on partial differential equations (see also ref.
2.2, 2.3, 2.5).
The significance of the theorem becomes clear if we choose for and to be
solutions of Laplace's equation.
For we select a solution that is uniformly valid inside Q, whereas for the
potential of a single source in a point Q of Q is taken, i.e.
 2-36

i
(2.1.107)
s - 4ir'(P,Q)

This means that satisfies Laplace's equation everywhere inside Q, except for

the point Q, where D becomes singular. We now apply Green's theorem to the
domain Q while excluding a small spherical region w with surface and radius c

centered in P because of the singular behaviour of for Q-P (fig. 2.1.41).

Then we have

s
- s ] dQ = - ii r an
- a
an
1
dS +
Qw

+ L 1
r3n
-
anr1)) dz (2.1. 108)
ITTlz

Since and satisfy Laplace's equation in Q (w excluded for ), the volume

integral at the left hand side of eq. (2.1.108) vanishes. The first part of the
second surface integral vanishes also because it can be shown that

hm
c-o
If r an
dz = C (2.1. 109)

The second part of the surface integral evolves as

him
° _
L jj
a
(- ) 1
dz] =- (P) (2.1. 110)
z

(These results are obtained due to the fact that z c2.

Eq. (2.1.108) may now be reduced to

= L
4n
¡j
r an - an
(_
1))
r
dS (2. 1. 111)

This integral representation is known as Green's identity for the solution of

the equation of Laplace in a point P of a domain Q.
For a point P outside the domain Q, it follows directly from Green's theorem
that

-
ff((
S
- 1
ran an (- )) dS = O (2.1. 112)
 2-37

Additionally, if and are also defined outside Q (in particular in a

vanishingly small spherical region around P) it follows that

(PeQ) = 0 (2.1.113)

We conclude that the potential inside and outside Q is entirely determined if on

the boundary S of Q both the potential and its normal derivative are prescribed.

We already know that for a well-posed boundary value problem (either inside or
outside Q) we must prescribe on the boundary S of Q either the potential
(Dirichlet problem) or its normal derivative (Neumann problem). Apparently
Green's identity describes a wider class of problems than the boundary value
problems in which we are interested. The difference is that Green's identity,
describes the solution both inside and outside Q, whereas the boundary value
problem is formulated either for the interior region or for the exterior region.

2.1.8 Boundary value problems and equivalent integral equations

As we have pointed out earlier, the boundary value problems we are confronted
with are of the mixed type: we seek the solution of a boundary value problem in
a domain Q bounded by a closed surface S, where on a part of S the potential
is prescibed and on another part the normal derivative . In this section we
analyse how we can use Green's identity to find the solution of such boundary
value problems. Considering eq. (2.1.111) it is clear that prescribing only on
S (interior Dirichet problem) is insufficient to determine in a point P inside
Q, because is still unknown on S. Similarly a Neumann problem cannot be
solved because also has to be known on S.
To solve either two problems an additional and important step has to be made.
For this purpose we consider Green's identity for the case that the (internal)
point PEQ approaches the boundary S, see fig. 2.1i2. Then Green's identity
gives in the limit

i a
(EeS) = H- ) - r
dS, qs (2.1.1114)
 2-38

The factor 1/2 arises from the fact that now around P a half-sphere ci of radius
-o has to be applied; the intersection of the half-sphere with the boundary S

must be excluded from the region of integration.

For the Dirichlet problem ( prescribed on S) eq. (2.1.114) may now be

considered as an integral equation for the unknown -----, thus
Q

1 1 3
- dS = (P) + '(Q) -p--
3flQ r
dS (P,QES) (2.1.115)

The left hand side of eq. (2.1.115) may be recognized as the right hand side of
eq. (2.1.74) which was the potential on a surface S produced by a source distri-
bution with strength a on S. In the case of eq. (2.1.115) we have a =
Q

If the solution (on S) of the integral equation (2.1.115) is known we are

flQ
able to compute the potential and the velocity (components) in the entire region
Q.

Eq. (2.1.115) is known as a Fredhoim integral equation of the first kind. (The
(general) theory of integral equations may be found in ref. 2.7).
In general the numerical treatment of Fredholm integral equations of the first
kind causes difficulties and should be avoided (see ref. 2.8); we return to this
later.

Now consider the problem where is given on S (Neumann problem). This leads to

an integral equation for the unknown on S

1) dS, (P,QES) (2.1.116)

(P) + (Q) (- dS =
?n (
r
4ns
Eq. (2.1.116) represents a Fredholm integral equation of the second kind, which
in general lends itself well to numerical treatment.
Note however, that in the particular case we are considering here, the solution
of (2.1.116) can not be unique, since, as we have seen earlier, a Neumann
boundary value problem for an interior region bounded by a closed surface is not
unique.
The non-uniqueness may be explained as follows. The left hand side of eq.
(2.1.116) may be recognized as the expression (2.1.82) for the potential on a
surface S produced by a doublet distribution with a strength per unit are of i
 2-39

-' on that same surface S and with the dipole axis pointed inwards. We know that
a uniform doublet distribution on a closed surface induces a zero potential
inside that surface (see section 2.1.6), resulting in a zero (normal) velocity
on the inner side of the surface. This implies that we may superimpose any
constant on without changing the righthand side of eq. (2.1.116), thus
reflecting the non-uniqueness.

Prescribing on part of S and on the remaining part we obtain, of course, a

mixed system of integral equations.

The considerations just given show that a well-posed problem does not automati-
cally lead to the preferred type of integral equation (Fredholm integral
equation of the second kind) and also that a, at first sight, suitable integral
equation (of the second kind) does not necessarily have a unique solution.

The formulation resulting in eqs. (2.1.115/116) is known as the direct

formulation (ref. 2.9). It has the important feature that, everywhere outside
the domain Q, the potential is equal to zero. In the practice of panel methods
the direct formulation is hardly used. However, a variant of it, to be discussed
in the next section, is used frequently.

2.1.9 General formulation of boundary value problems in terms of source and

doublet distributions using Green's identity

Once more we consider Green's identity, but now with scalar functions and
defined for the interior region Q and scalar functions " and ' for the
exterior region Q' of a closed surface S, see fig. 2.1.43.
If satisfies Laplace's equation inside Q and ' satisfies Laplace's equation
outside Q, i.e. inside Q', the potential in a point P in Q may be expressed
as, see eq. (2.1.111)

(p)=Lff{
1- -r l)3a(
3n 3n
1))ds
r
(2.1. 117a)

Because of (2.1.112) there also holds, inside Q

o=fJ[(-1
S
r an'
a
3n' r
(2. 1. 117b)
 2-40

Here n' is the unit normal on S pointing into Q'. If ' should be a solution of
a well-posed boundary value problem in Q' we formally have to bound Q' by a
second closed surface S' and apply the integral (2.1.117b) also to this surface.
In doing so we will assume that S' is infinitely far from S and that ''-*o for
S'4w; in that case the integral over S' can be shown to vanish, see ref. 2.2.
Then, superposition of eqs. (2.1.117a) ai-id (2.1.117b) yields, after replacing
a a
j-- by - j-

i a a' a 1
(P) = N- ) (j-h - j---) - - 3n
;fl dS (2. 1. 118)

The first part of the integral represents the potential due to a source distri-
bution on S with strength o
a
j-- -
a' and the second part represents the
=
potential produced by a doublet distribution on S with strength p = -( -

Eq. (2.1.118) implies that a solution of Laplace's equation may be expressed in

terms of a source distribution plus a doublet distribution on the boundary S of
Q. We write eq. (2.1.118) as

(P) = L jj [(- -)
r
i
+ p an- (-iS
a
(2.1.119)
4TIS

In analogy with the preceding analysis we may express the potential ' in a

point P' of the exterior domain Q' as

-
1
i- J! ( 1 a'
(j--
a
- j-h-;-) - (' - aii'
dS (2.1. 120)
S

or

a
=
[G (_ r an
-fldS (2. 1. 121)

Eqs. (2.1.118/119) and (2.1.120/121) govern a broader class of boundary value

problems than those adressed by the direct formulation of eqs. (2.1.111) and
(2.1.112). The latter case is re-obtained if we make the specific choice ' = 0.

In the general formulation ' is the solution of an arbitrary, (but well-posed)

problem in the domain Q'.

 2-41

The integral representation (Green's identity) (2.1.118) for the potential (P)

in Q is valid independent of the choice of '. However, it is not unique as long

as ' is not specified. Specifying ' in Q' can be done by prescribing either '

or on S, because this is identical to formulating a Dirichiet or Neumann

an
boundary value problem for ' in the domain Q'.

In order to evaluate (P) from eq. (2.1.119) for a point P on S we repeat the
limiting proces that was also applied in obtaining eqs. (2.1.74) and (2.1.82).
The result is

1
(P) = (a(Q) (- -) + j.i(Q) dS - p(P), for PS = aQ
arlQ
S
(2. 1. 122)

where we have taken the + side (or Q-side) of the surface S. If (PS) is
prescribed (Dirichlet condition) eq. (2.1.122) can be interpreted as an integral
equation for o and p. However, the integral equation will not have a unique
solution as long as either o, or p or a relation between both is not prescribed.
The situation suggests that prescribing either a, or p or a relation between
both is equivalent to specifying ', about which later.
If (PS) is prescribed (Neumann condition) we may obtain an integral equation
for a and p by considering the integral representation for and by letting

P4S (see also eqs. (2.1.76> and (2.1.87)). We then find

N 1))
i
an
(P) =[o(Q) a

an (- -)
r
+ i-i(Q) -q---
j dS + o(P)

for PES = aQ (2.1.123)

For the saine reason as before the solution for o and p is not unique. In this
case it is not unique even for prescribed values of o or p, because of the ill-
posedness of the Neumann problem for an interior domain.

For the mixed boundary value problem a system of integral equations may be
obtained that has a unique solution provided that either o or p or a relation
between both is prescribed.
 2-42

Similar to the analysis that led to the expressions 2.1.122/123 for (P) and
expressions may be obtained for ' (P') and (P') for a point P' at
np
the Q' side of S (aQ'),

=
i
[o(Q') (- ) + i.i(Q'
anQ? (
1j
- r1
S
P'aQ' (2.1.124)

and
i
an, - 4n [(Qe) 8n
8 1
(-)+p(Q')
r
a
ann, 8nQ
a
))dS-'
r
G(P')

P'caQ' (2.1.125)

These expressions may be used for formulating a boundary value problem in the
exterior domain Q' (see fig. 2.1.43).

2.1.10 Uniqueness of the integral representation derived from the general

formulation

In the preceding paragraphs we concluded that the general formulation leads to

an integral representation for or that is not unique unless ' is specified
and to an integral equation for a and ji that does not have a unique solution
unless o, or ji, or a relation between the two is specified. We also suspected
that prescribing either a, or ji, or a relation between the two is equivalent to
fixing '.

In order to confirm this consider again the integral representations

(2.1.118/119) for the potential in the interior region Q and (2.1.120/121) for
' in the exterior region Q' . Let the source strength a be prescribed on S.
Then, the doublet strength p and the potential (P) on S can be determined by
solving a boundary value problem for in Q, i.e. by solving an integral
equation for p of the type (2.1.122) or (2.1.123).
Similarly, if p is prescribed, one can solve an integral equation for a. If a
and p are both known on S, also ' is known in Q' by means of eq. (2.1.121).

Conversely the doublet strength p and the potential ' (P) on S' can be
determined by solving a boundary value problem for ' in Q', i.e. by solving an
 2-43

integral equation of the type (2.1.124) or (2.1.125). Then in Q is known

through eq. (2.1.119).

In other words: for any boundary value problem in Q there exists, for any choice
of o or p, an equivalent boundary value problem for ' in Q'. Which equivalent
boundary value problem is concerned precisely, is determined by the jump
conditions associated with the particular choice of o and p.
In principle the number of possibilities for fixing o or p (and through that ')

is infinitely large. The following possibilities are the most convenient ones
and indeed used most often:

p=O, or =' on S
In case of a Neumann boundary value problem for eq. (2.1.123) then leads
to a Fredholm integral equation of the 2 kind for o.

o0, or
a'
an
=onSa
an
In case of a Dirichiet boundary value problem for ' eq. (2.1.122) leads to
2nd
a Fredholm integral equation of the kind for p.

A third possibility, as we have seen already, is to choose '=O everywhere

inside Q'.
This results in the direct formulation, discussed in section 2.1.8. There
we saw that this is equivalent to

0=-
an
(2.1. 126)

and

u=- (2.1.127)

In case of a Neumann boundary condition on aQ for one may use (2.1.126) to fix
Then (2.1.124), with '=O as boundary condition for the (equivalent) boundary
nd
value problem in Q , leads to a Fredholm integral equation of the 2 kind for
In case of a Dirichlet boundary condition on aQ for one may use (2.1.127).
 2-44

Then (2.1.125), with = O as boundary condition for the (equivalent) boundary

value problem in Q', leads to a Fredholm integral equation of the kind for

a.

The possibility to replace a boundary value problem for by an equivalent

boundary value problem for ' (as implied by iii)) is often used in practice.
Example: a Neumann problem in Q as shown in fig. 2.1.44. The source strength is
determined by eq. (2.1.126). This problem can be replaced by the Dirichiet
problem '=O in Q'.
It may be interesting to use this possibility in order to obtain a Fredholm
integral equation of the second kind and/or to reduce the computational effort
(see section 2.1.13).

2.1.11 Restrictions with respect to the choice of source and doublet

distributions

Although the number of possibilities for fixing o or j.i is unlimited, in

principle, there are restrictions.

As an example consider the Neumann problem for the exterior domain of a closed
surface S with a net mass flow through S, given by (fig. 2.1.45)

dS xO
an

It is not possible to model this problem by means of only a doublet distribution

on S, because such a distribution does not produce mass. Even for

dS = O

we run into problems, because we can add a constant to any solution i for the
doublet strength which does not affect the normal derivative (the potential
of a constant doublet distribution on a closed surface is zero outside the
surface *), section 2.1.6). If we consider this problem from a different
viewpoint we arrive at the same conclusion: with only a doublet distribution on

*) Note that because of this the external Dirichlet problem cannot be solved
either with only a doublet distribution!
 2-145

a
S the normal derivative - is continuous across S therefore the equivalent
an
interior problem is a Neumann problem of which we know that it does not have a
unique solution.

A second class of problems which imposes restrictions on the choice of

singularity distributions is formed by boundary value problems with circulation.
The latter is given by

B
r= f asdsxO
A

(s is the distance measured along the contour of integration, see fig. 2.1.46).
This problem cannot be modeled by only a source distribution on the surface S,
because the circulation around a single source or around a source distribution
is zero. In order to be able to do so a non-zero circulation a doublet
distribution or a vortex distribution is necessary on SB and S (fig. 2.1.46).
These examples illustrate that a sufficient wide class of boundary value
problems can only be modeled by singularity distributions consisting of both,
sources arid doublets.

2.1.12 Summary of well-posed and not-well-posed boundary value problems and

integral representations

Boundary value problems

- Dirichlet problem: Well-posed for both interior and exterior domains.

- Neumann problem: For an interior domain there exists only a solution if

- dS = O, if a solution exists it is not unique.
an
For an exterior domain there exists a unique solution if the circulation
is prescribed (this is established by a Dirichlet condition at infinity
and requires in general the introduction of a cut).

- Mixed problem: In general well-posed for both interior and exterior domains.
 2-'46

Integral representations

- A solution of Laplace's equation inside a closed surface S may be expressed in

terms of a source distribution plus a doublet distribution on S.

- The integral representation thus obtained is not unique unless we prescribe

either the source strength or the doublet strength (this is equivalent to
prescribing the solution outside S).

- For a given boundary value problem inside S and a given source strength or
doublet strength an equivalent boundary value problem outside S can be
formulated.

- In general it is not possible to model boundary value problems with Neumann

boundary conditions without using a source distribution. Also it is not
possible to model boundary value problems involving circulation without using
a doublet distribution.

Integral equations

For a well-posed boundary value problem in a domain Q it is possible to derive

(through a unique integral representation) an integral equation having a unique
solution satisfying the boundary conditions on the boundary S of Q.
The resulting integral equation is of the Fredholm type, first or second kind.
From a standpoint of' numerical treatment the Fredhoim integral equation should
preferrably be of the second kind.
In case of a problem with Neumann boundary conditions an integral equation of
2nd kind is obtained when the source strength is taken as unknown, whereas
the
in the case of Dirichiet conditions the doublet strength should be taken as
such. Note however, that it is not always possible to do so. For example in
lifting surface theory a doublet (or vortex) distribution is the only admissable
type of singularity while the boundary condition is of the Neumann type.
 2-47

2.1.13 The application of the integral representation to the boundary value

problem for the flow around an aircraft configuration

Consider the boundary value problem for the flow around a configuration as shown
in fig. 2.1.47 (depicted earlier in fig. 2.1.18).
In this section we will analyse how such a problem may be transferred into a
system of integral equations. As we have seen (in the preceding sections), the
solution of a boundary value problem can be represented by means of a source
distribution and a doublet distribution on the boundaries S of the domain Q
under consideration. Hence, the perturbation potential in a point P of Q can,
according to eq. (2.1.119), be expressed as

1 1 a
P) ((_ -) ' dS (2.1.128)
BCT
+ _
r r
S +S + +

The various parts of the boundary surface, SB$SC,S and ST, distinguished in
section 2.1.3, have been indicated in fig. 2.1.47.
Considering first the contribution of the vortex sheet we learn immediately
that the source strength on S is equal to zero because of the boundary condi-
tion = O. In addition the doublet strength is directly connected to the
jump of the potential across SC, because i=- 1.
Therefore the contribution of S to the integral is

a
p (-
1
;) dS=--ff]--(-1)dS (2. 1. 129)

Since is constant in s- (or x)direction it can be written as

B
the jump in potential at the trailing edge. For an orthogonal coordinate system
(s,t) on SC one may write

tR
1 a
ISB( f ) ds) dt (2.1.130)
SB

where tL and tR represent the left and right edges of S, respectively.

Next we consider the singularity distributions on S and
 2-48

Taking ' = O in the domain Q' outside S and ST we have, according to eqs.
(2.1.126) and (2.1.127),
p = - = O, and o = - on S
an
and
o = - = O, and p = - on ST

On S not only the perturbation potential vanishes, but the perturbation

velocity as well. So, on S the source strength as well as the doublet
strength is zero, leaving no contribution of S at all to the surface integral
in eq. (2.1.128).
At the surface at downstream infinity, ST, the problem is more complicated.
Although the source strength vanishes, this is not so for the doublet strength,
in particular in the vicinity of the intersection S of S with ST. The latter
may be understood by considering the perturbation potential ' at the outer side

of ST (according to eq. 2.1.1214)

' (FS,) =
TTS
pp
an (_ rJ
dS + p(P) + fi J dS

T
=0
,. if as + if pj-
a 1j dS (2.1.131)
- r
S S
C

In (2.1.131) the integral over SB vanishes, because v*, whereas, as noted

before, the integral over S is zero because of the boundary conditions.

However, the contribution to the potential on ST due to the doublet distribution
does not vanish. Therefore, the first integral is not zero, thus p5 z Q.
on S
T
On the other hand, since the perturbation potential in a point P of Q due to the
doublet distribution on ST may be written as

a lj
(P) dS (2.1. 132)
r
T

and since ST is at downstream infinity, the distance r between an arbitrary

point P in Q at a finite distance of 5B and tends to infinity, the
contribution of 5T to the perturbation potential in such a point P vanishes.
In summary, the expression (2.1.128) for the perturbation potential at a finite
distance from S may be replaced by
 249

1
((
1)
r
+ an
p a

r
dS - L jj ['t'i (- )
dS, pQ
4TTS

(2.1. 133)

If the point P is on the body SB itself, eq. (2.1.112) becomes

PSB)
= L [(_ r) +
i a

anQ
(_1))
r
as - 3i(P) - ¡f (- )
dS
4SB
(2.1. 134)

In that case the expression for the normal derivative of ' takes the form

= L [a
3n
3
i- ) +
i a

i
a

r )
) dS +
2
o(P) +
S

a
-
3flQ
1-)dS=-U
r
n
- (2.1. 135)

Here we have (re-)introduced the subscripts P and Q to discriminate between the

'fixedt and 'running' variables on S.

Eq. (2.1.135) represents an (incomplete) integral equation for the unknowns a, p

and i
. To complete it we have to construct an equation for the Kutta
B
condition and to make a proper choice for a or p on SB. If eq (2.1.135) is to be
2nd
a Fredholm integral equation of the kind one must adopt a as unknown and fix
p and [DJS in such a way that the Kutta condition is satisfied.
C

With respect to the prescription of a or p we noted already (section 2.1.10)

that the 'external' Neumann boundary condition on SB may be replaced by the
'internal' Dirichlet condition = O inside S
B
if wete a = -U n.
 2-50

In other words an alternative for eq. (2.1.135) is (similar to eq. (2.1.124))

= = (_ U ) (1) dS +
1'a1Q
(.1) dS +

i(F') -tj- a 1) dS (2.1. 136)

r

nd
Note that this is (also) a Fredhoim integral equation of the 2 kind (provided
is known).
C
A further important point to notice is, that in case of (2.1.136) the
perturbation potential on SB follows directly from the jump relation

Q - Q' = -i-'

Since G' = 0. it follows that

(2.1.137)
Q(PSB) = -

In other words the velocity distribution VG on SB can be determined immediately

from the derivative of the doublet distribution. Later on we will see that this
implies an advantage in terms of computing time.

With respect to the Kutta condition it has been pointed out in section 2.1.3
that, in the first place, this requires continuity of the (perturbation)
potential across the trailing edge of the body. This can now be stated as (see
fig. 2.1.48 and recall that [[]1 -L')

G(P)
for P1,F2,P,P + T

+ Q(P) = G(P) +

G(P1) + (2.1. 138)

or

From the jump conditions for a and p on SB it then follows

 2-51

i(P1) 9 (P2) = IL ; P1, P2, P, P T

- (2.1.139)

Hence (or ) can be expressed in terms of

This means that in case of .
ii
C C B
the 'internal Dirichlet' formulation (2.1.136) there is no need for an explicit
additional condition (other than (2.1.139)) to fix the circulation. This
situation is therefore commonly referred to as one with an implicit Kutta
condition. Note that (2.1.139) implies continuity of potential (or a finite
velocity) only. Continuity of velocity across the trailing-edge is not
necessarily satisfied. The discussion of section 2.1.3 suggests that this will
be the case only if the position and shape of the vortex sheet SC satisfies the
criteria formulated by Mangler and Smith.

The situation is different in case of the 'external Neumann' formulation

(2.1.135). Recall that in that case o (on SB) is the unknown and that p5 and
B
must be specified in advance. Eq. (2.1.139) implies that p5 and p cannot
c B C
be specified independently and that, as before, p5 can be expressed in terms of
C
or vice-versa.
B
(Apart from this we are completely free in the choice of
B
Unfortunately this is now not sufficient to fix the circulation ); the
C
reason being that now , rather than , is the unknown in the integral
B B
equation. Apparently an additional, explicit condition linking o, with jic, and

is required to complete the formulation of the integral equation problem. In

B
the light of the fact that it is not possible to excercise control over
circulation (jump in potential) by means of source distributions this is not
surprising. On the other hand the source distribution on SB (as well as p5 and
B
p5 ) will have an influence on the velocity at the trailing-edge, but not
C
necessarily in such a way that continuity of the velocity across the trailing-
edge is satisfied. This suggests that we must add some form of explicit Kutta-
condition expressing the requirement of continuity of velocity. As discussed in
section 2.1.3 either the 'directional form', (2.1.54), may be used, or the
'equal pressure (velocity) form', eqs. (2.1.43) to (2.1.47).
The directional version, obtained by differentiation of eq. (2.1.133), is
 2-52

1 111a 1 a
+ 1J: -
a
(-
1
) dS +
anPPT = an
SB

an
a
(
-j dS = - (U. "PP+T' lDC (2.1.140)

A sketch is shown in fig. 2.1.49.

The pressure version in its general form is given by eq. (2.1.43). Using the
perturbation potential we find

2 - 2

[ 5
(u x+')) - [V
5
(u x+)) -* O for P1,P2+T (2.1.141)
2

with eq. (2.1.113)

(P) = L [a(- r)
+
anQr2
a
li-
4SB
1 a

_fI3flQ C
1)dS
r
(2.1. 142)

and where denotes the gradient tangential to S.

s

Panel methods are based on the numerical threatment of integral equations of the
type given by eq. (2.1.135) or eq. (2.1.136); in both cases the system is
completed by the additional relations (2.1.139) plus, in case of (2.1.135), by
(2.1.140) or (2.1.141).
Since there are many possibilities for choosing the singularity distributions a
and i, together with different possibilities for the formulation of the Kutta
condition and different numerical treatments, a considerable number of different
panel methods have been developed over the years. The most important of these
are summarized in Table 2.1. The methods mentioned in the table are all
applicable to general three-dimensional configurations.
There also exists a large number of methods for the computation of potential
flows around more restrictive geometries (two-dimensional, axi-symmetric, thin
wing configurations etc.). These methods also differ not only in their basic
formulation but also in the numerical treatment, particular applicability etc.,
they will not be discussed here.
 2-53

The numerical aspects of general three-dimensional panel methods will be

discussed in the next chapter.
 2-54

Naine of panel method Type of singularity distribution

DOUGLAS-NEUMANN, 1962 (ref. 2.10) - o distribution (p=O, no circulation)

- external Neumann boundary condition:

BOEING (TE)A230, 1967 (ref. 2.11) - o and p distribution

NLRPANEL, 1969 (ref. 2.12) - p distribution (inside SB) fixed in
MBBPANEL, 1970 (ref. 2.13) terms of p5 = f'(t): = f(S)ps
''s
C B C

BAC, 1972 (ref. 2.14)

BAC, 1976 (ref. 2.15)

- external Neumann condition
- 'directional' Kutta condition

DOUGLAS-NEUMANN LIFTING, 1972 - o and p distribution

(ref. 2.16, 2.17) - p distribution fixed in terms of
p5 = f(t): = f(s)ps
'1s
-
I)
C B C

- external Neumann condition

- nonlinear pressure version of Kutta
condition

PANAIR, 1975 (Boeing, ref. 2.18, 2,19) - o and p distribution

VSAERO, 1981 (Anal. Meth, Inc. ref. 2.20) - o distribution fixed through o =

QUADPAN, 1983 (Lockheed, ref. 2.21) o, p

AEROPAN (NLR, in development)

MCAERO, 1980 (Mc. Donnell-Douglas, - internal Dirichiet condition

ref. 2.21), directional Kutta condition - implicit Kutta condition

Table 2.1. List of main 3-D panelmethods.

 2-55

2.2 Numerical aspects

2.2.1 General considerations on the discretisation of the integral equations

In the preceeding chapter we have seen that integral equations for boundary
value problems representing the flow about aircraft configurations can be
obtained by either applying an external Neumann boundary condition or by
applying an internal Dirichiet boundary condition.
The integral equation based on the external Neumann condition is given by eq.
(2.1.135). For a point P on the (external) surface SB of a configuration it
reads

a i a l
r(P,Q)1 dS +
i o(P) +
= o(0) -

a a i )dS+fJp(Q) a a i dS
(
an aflQ r(P,Q) n an anQ r(P,Q)
s

=-U ; (2.2.la)

The doublet distribution (p(Q))5 on the configuration SB is expressed in terms

B
of the doublet distribution on the vortex sheet S, with [p(Q))5 being
C C
constant in the streatnwise direction.
We may write

= f(s) [p(0)J5
p(Q)115 = f(s)p (t)
B B
[u(0)) =p (t) (2.2. lb)
C SB

where p (t) is the doublet strength at the trailing-edge (see fig. 2.2.1).
B
The system of equations is closed by adding as 'explicit' Kutta condition either
the 'directional' form (2.l.l1tO)
 2-56

30 1 a -1 a a -1
- f! o(0) Lr(P,Q)) dS +
r(P,Q )) dS
-

an 3flQ
SB

for PS -s
+ p
an0 r(P,Q)
dS = - U xi
P CB (2.2.lc)

or (an approximation to) the 'equal pressure' form (2.1.141/142)

The integral equation based on the internal Dirichiet boundary condition is

expressed by eq. (2.1.136). For a point P' on the (internal) surface 5B it takes
the form

1 a -1
= ° dS + 1_l(Q) dS +
r(P',Q) anQ r(P',Q)
SB SB

+ 1l(P') i(Q) dS = O (2.2.2a)

4ii r(P',Q)

where the source distribution [o(Q)Js on SB and the doublet distribution

B
are given by
C

[a(Q))5 = - (U
B Q)sB
(2.2.2b)

[p(Q))5 = 1_ls
(t)
C B

(t) is eliminated by the implicit Kutta condition (see fig. 2.2.1)

B

p
(2.2.2c)
= p2) -i(P1) + = sB(t) for P1,P2,P' + T OrI SB

Substituting eq. (2.2.lb) into eq. (2.2.la) and rearranging the result we obtain
for a point P on SB
 2-57

Lc a(0)
a

an (
r(P,Q)
dS + o(P)
SB

+ L jj
4T1
f(SQ)
a
(t0) -;:;;
an0
a -1
r(P,Q)
dS +
SB
SB

tR
a a 1
[ ds(Q)] dt0= - U (2.2.3)
an anQ r(P,Q)
tL

Analogously we may derive by substituting eq. (2.2.2b) into eq. (2.2.2a) for P
on SB

a -1
dS + p(P) +
an0 r(P,Q)
s
B

t
R
a i
)[ f anQ L r(P,Q)1 ds(Q)] dt0 =
B0 sß(t)

L (U J
r -1
dS (2.2.L)
p r(P,Q)
SB

For completeness we note that the functions

-1
potential of a source
r(P,Q)

a i
L
potential of a doublet
r(P,Q)

(2.2.5)
a i
normal velocity component due to
an '- r(P,Q)
a source

r -1
normal velocity component due to
an -r(P,Q)
a doublet

appearing in eqs (2.2.1/2) are often referred to as (singular) kernel functions.

In panel methods the integral equations (2.2.3) and (2.2.4) are discretised
through the following steps.
 2-58

First the surface of the configuration and its vortex sheet are subdivided in
a number of elements or panels, see fig. 2.2.2. In most cases it suffices to
approximate the geometry of the vortex sheet by a rather simple geometry on
which the doublet strength is assumed to be constant in streaniwise direction.
Therefore, the panel distribution on the vortex sheet can be limited to only
one or a few panels in the streaxnwise direction.

On each panel i a collocation point P. is selected in which the boundary con-

dition and hence the integral equation is satisfied. In addition, if
applicable, points PQ where the Kutta condition is to be satisfied also have
to be chosen.

The integrals in eqs. (2.2.3) and (2.2.14) are replaced by summations of the
integrals over the individual panels. For example
dS
r(P.,Q)
SB i

is replaced by

NB
'
Z f! a.(Q) (- r(P.,Q) J dS + o(Q) r(P.
dS (2.2.6)
n i
j=1 ÓS
jB i
iB i,Q)
ixi

and
dS
3flQ r(P.,Q)

is replaced by

NB
L
4n
jj (Q) aÇ (- r(P1Q)
dS +
j=l ô.S
ii

L a

aflQ r(P
-1
,Q)) dS (2.2.7)
4T1ÔS i
iB

Here ó.S and ä.S denote the surfaces of individual panels with index i
iB
and j, respectively. NB and N are the total number of panels on SB and S,
respectively. Note that because of the fact that we are dealing with
 2-59

principle value integrals the integrals over panels with index j=i must be
treated separately.

Because of the fact that, on Sc, 1.1 is constant in the streainwise direction
it is convenient to rewrite the integral over S as
tn
1
SQ [f anQ
1
tr(P.,Q)1 ds(Q)] dtQ
sB(t)

and replace it by

Nc
1
z f ds(Q)] dtQ (2.2.8)
k=l äkt
ij
530 aflQ r(P.,Q)
sB(t)

4. The integral over a panel is approximated in the following way.

a. The geometry of a panel is approximated by a polynomial representation,
e.g. of the type

= a00 '- a10 + a01ii + a202 + a11n + a022 + .. . .0(h3) (2.2.9)

Here are local orthogonal coordinates attached to the panel under

consideration (fig. 2.2.3). The local coordinate system is generally chosen
such that the origin coincides with the collocation point, with the plane c =

O tangent to the panel (a00= 0). The quantity h represents a characteristic

length of the panel (for example the length of a diagonal); the coefficients
a , etc. are constants.
00
The area of a panel is then given by

dS
dd 2 3P)2]1/2
= [i + +
an
ddn (2.2.10)

where n is the Ç-component of the normal in a point on the panel, see fig.
2.2.3.
 2-60

b. The distribution of the source strength o and/or the doublet strength p on

a panel (with index j) are also expressed in terms of polynomials

O = 000 + 0.10 + 0.01 n. + ....0(h2)

J J J J J J
(2.2.11)

Pi = + p.1° + p.° n. + ....0(h2)

.

Assuming continuity on the singularity strength distribution the polynomial

coefficients l 0, 00 1
, l 0, po i may be related to the values 00 0 and po o of
the adjacent panels. For example, interpreting (2.2.11) as a Taylor series
for o around ,n = 0 it follows that

00 00
- °j-1
= ( = 0, n = 0). + 0(h2)
j+1 -1

In general o.'° may be expressed as

0.10 =a.ci°.° +b. 00.0 +c. 00.0 + (2.2.12)

J J 31 3 3 j j+l

etc.

In most panel methods the geometry of a panel is simply taken as a plane

surface and the associated source/doublet distribution as constant.

c. Now the integral over a panel can be worked out in its local coordinate
system by substitution of eqs. (2.2.10) to (2.2.12) into eqs. (2.2.6) etc.
As an example we carry this out for the integral of eq. (2.2.6)

-1
JI o.(Q) ,Q)) dS
r(P
os
jB i

-1 dÇdn
II (o°° + 010 +
os j J {()2 + (n.-n)2 + (.)2)h/2 InI
j B

= 00.0 A°.°. + 10 A'.°. + 00.1 A°.' + . . . . (2.2.13)

J 13 j 13 J J

where
 2-61

If
- 1 ddn
6jBS (()o + (-n)2 + (çç)2)'/2
A°.°. =
lnI

ddn
()Z
-
= JI (2.2.14)
'J ô.S (R_2 i
+
iP )2)/
(Ç.Ç InI
ddi
A°.'.
'3
=
6 S [z + ()2 +
i
f! -
(ç -ç )Z)l/2 Inçi

As we shall see later the integrals can be evaluated analytically for

planar panels, whereas for curved panels the evaluation must be established
numerically or through series expansion.
Summing over all panels gives

NB
a(- 1) dS = (A°° 000 + A1° + A°' 0.1 + (2.2.15)
)
4ii r ii j ii j
SB i =1

where it is noted that, at least in principle, the case j=i has to be

treated separately because of the principal value involved.
Since o.°, 01 etc. can be expressed in terms of' values of 00
, in adjacent
panels, aY1, 0°.°l etc. (see (eq. (2.2.12)), eq. (2.2.15) may be
rearranged into

N
L0(
1-hi r
) dS
B
(A.. 00.0
13 3
+ ) (2.2.16)
S5

10
where A. . + 10
13
= A°.°.
13
+ A'.°. b.
13 3 i,j+1 j+l a. + A.
i,j-1 c.j-1
. +

The coefficients A. . are called the aerodynamic influence coefficients

13
(AIC) for the potential induced by the source distributions on the panels.

It is possible to derive similar AIC's for all the terms of eqs. (2.2.3)
and (2.2.4). Consider for example the perturbation velocity-components on
the i-th panel due to the source distribution on the j-th panel. The
 2-62

contributions of the j-th panel are given by the derivatives

They may be expressed as

(!_)
a.
-
f!

3B
a(Q)
if i
r(P., Q) ) dS = (A i.. o.
13 3

i
a(Q) dS = (A ).. o.
ar = r(P., Q) )
T1J 3
i
'

o.s i i

= o(Q) r(P.,Q» dS = (A)1 o

ôiS

Since
3n. j .
j
= n.
-

i
,-
(V
P.
j .

3
we may write
. ,a
j.
an. j
= n.
-
i
-
A.
13
. a., where A.
3
-
ij
=

i.
[(A
(a
j
ij.,
. (A
T)ij
) .,
.
(A j . .J.
Çij
Writing further that A
13
. =
1
A .
13
we arrive at

j . = A a..
an.j
.

13 3

As a result of analogous treatment the integrals in eqs. (2.2.3) and

(2.2.4) may be represented by

NB
a(Q) dS + a(P.) E A' a + (Truncation Error)
r(P.,Q) ii
i j=i
i
(2.2.17)

Nc
f(SQ)
aÇ k - r(P.,Q)
i
aS =
k=l
E B'
ik (sB)k + T.E.
i
(2.2.18)

tJ: Nc
1
a a
r(P ,Q)i dSQ] dt =
I C' T.E.
(tQ
sB(t)
an 3flQ
i k=1 ik'sBk
i
(2.2.19)

NB
dS + B p + T.E. (2.2.20)
p(Q) -p---
an _ r(P
Q)
ti(P) = I
Q i' i=1
 2-63

NB
flQ) H r(P.,Q) dS A.. (U + T.E. (2.2.21)
= 3-
B

tR w Nc
1
SQ sB(t)
3flQ r(P
i
Q)) dsQ dt
k=1
c
ik SßK )
+ T.E.

(2.2.22)

In eqs. (2.2.17) to (2.2.22) A. ., A ., B. ., B ., C. , C are aerodynamic

13 13 13 13 ik ik
influence coefficients. a. and p. represent the source strength and doublet
strength a°.° and p°.°, respectively, in the origin of the local coordinate system
1 1
('midpoint') of the j-th panel. (Note that terms p have been included in
the definitions of A. and B.. in eqs. (2.2.17) and (2.2.20), respectively)
The truncation errors in eqs. (2.2.17) to (2.2.22) are dependent on the degree
of the geometry and singularity representations in eqs. (2.2.9) to (2.2.12) and
on the accuracy of the integration procedure over a panel.

5. If the integrals in eqs. (2.2.3) and (2.2.4) are replaced by the summations
given in eqs. (2.2.17) to (2.2.22) we obtain for eq. (2.2.3)

NB N0
X A . a. + I (B:
ik
+ c:ik ) (p 1 =- U n.
1
,i = 1,2,3, . . . N
B
SB 1(
k=1
(2.2.23)

and for eq. (2.2.4)

N N N
B c B
X B..1.i.+ I C.(p
ik J
SB k
= I A..(U
w ).,i=l23.....NB
k=l j=1
(2.2.24)

In the case of zero circulation(p 0) satisfying eqs. (2.2.23) or (2.2.24)

B
in each collocation point P. of SB leads to a system of linear algebraic
equations that can be written as

or
'r a. = _(Uw.i:)j
[A: .1
j
(2.2.25)
 2-64

[B..] p. = -[A..]
13 3 13
(U ).J
m
(2.2.26)

where [A'..] and [B..] are square matrices of which the elements are the
aerodynamic influence coefficients, o. and p. are the vectors of unknown
source and doublet densities to be solved for.

In the general case with circulation a number of N equations representing

the Kutta condition must be added in order to obtain a closed system. We
first consider the case of the 'external Neumann' formulation (2.2.23).
With the (explicit) directional version of the Kutta condition the additional
equations are discretizations of eq. (2.2.lc). The discretized additional
equations look much the same as eq. (2.2.23), as will be clear from the close

resemblance of eq. (2.2.lc) and eq. (2.2.la). Then, in matrix notation the
complete system of algebraic equations becomes

A. Bk + o.
J
i,j = 1,2,3 . .
. NB

k,Q = 1,2,3 N (2.2.27)

A. Bk + 0Qk
-U ...

Here the aerodynamic influence coefficients A'13., A'QJ describe the influences
.

of the source distribution on the body in the collocation points P. and in

the 'Kutta points' P respectively, B!k + C!k, Bk + describe the

influences of the doublet distribution in the collocation points P. and the

Kutta points PQ. Note that the number of unknowns is NB + N.

In the case of the 'internal Dirichiet' formulation eq. (2.2.24) and the

'implicit' Kutta condition (2.2.2c) may directly be expressed in terms of

B
the p. adjacent to the trailing edge (say p*).
This means that we may incorporate the ok of eq (2.2.214) into the B..
belonging to these p.. If we denote the modified influence coefficients as
B* the system may be written as
 2-65

Ü.i;j\ i,j = 1,2,3 NB - 2Nc

'li
- (2.2.28)

pi * U ij = 1,2 . .
.
j

Note that the number of unknowns remains NB and that the right-hand side of
(2.2.28) is identical to that of (2.2.26).

The linear algebraic system (2.2.27) or (2.2.28) are solved using numerical
techniques that are either standard or are developed specifically for this
purpose (see section 2.2.9).

When the solutions o. and p. are known, the potential and its derivatives can
be computed throughout the entire flow field including the body surface SB.
In doing so one can, again, make use of aerodynamic influence coefficients.
For example for an arbitrary point P one may write

for the potential

= if o(Q) J dS + L jj a
dS
[
r(P ,Q) 4ii r(P,Q)
ni
SB SB+SC

= A a. + X B p (2.2.29)
j
.

1113 3
k
mkk

and for the velocity components

a(P m) 8A . aB
mk
u(P
ni
)
ax 3x ax
.
k

a(P m) aA . aB
mk
v(P ) = (2.2.30)
ni ay ay j ay
k

a(P ) aA . aB
mk
w(P)= in az az G +X az
- k
 2-66

For 'internal Dirichiet' methods the velocity vector on the bodysurface SB

follows directly from (see section 2.1.13)

=- 1.1

Since p can be approximated directly through numerical differentiation (finite

differences) it is, in this case, not necessary to make use of the aerodynamic

influence coefficient ---, etc.. In other words it saves us the matrix-

ax ax
vector multiplications implied by (2.2.30).
Once the velocity components are known the pressure can be calculated from
Bernoulli's equation and the forces and moments through numerical integration of
the pressure distribution over the body. For example the force component in x-
direction can be obtained from

NB
K
x
= - f! p n
S
x
dS = -
.

i=l
Xli
(p n ) .ó.S + (2.2.31)

From the velocity components and pressure distributions one can also construct
streamlines, isobars, etc.

2.2.2 Surface-grid generation ('paneling')

In order to execute the necesary computations efficiently it appears to be

necessary to make use of computer aided generation of the surface mesh on which
the integral equations are discretized. This requires the availability of a
continuous, analytical description of the surface of a configuration.
There exist a number of commercially available software packages for the
description of geometries, among them we mention
CATIA (Dassault)
AEROLIS, (formerly SIGMA), (Aerospatiale)
CALMA (General Elelctric)
CADAM (Lockheed)
PATRAN
 2-67

These geometry packages offer the possibility to describe surface parts (called
patches, fig. 2.2.4) analytically, either in terms of elementary geometrical
shapes such as plane surfaces, cylinders, spherical surfaces, ellipsoids, etc.
or in terms of more general surfaces described by polynomials (e.g. splines).
Once the analytical description of a patch is established, the coordinates of an
arbitrary point on it can be determined easily.

Panel shapes
In principle there are many options for panel shapes, e.g. triangular, qua-
drilateral (trapezoidal), pentagonal, etc. Triangular elements (fig. 2.2.5),
applied frequently in finite element methods, have the advantage that they are
easily patched together in irregular shapes. Also a continuous representation of
a surface is always possible (no gaps). The disadvantages are a rather
complicated administration of panels (i.e. indexing) and a (more) irregular
distribution of nodal or collocation points which is undesirable for aerodynamic
applications. In aerodynamic practice almost exclusively quadrilateral elements
or panels are used (fig. 2.2.6). The main reasons are: they generally fit
reasonably well into regular aerodynamic shapes such as wings, fuselages, etc.;
they can be administrated by a simple, cartesian indexing system and it is
always possible to construct a triangular panel as a special case of a
quadrilateral. A disadvantage is that not all surfaces can be represented
smoothly without gaps.

Organisation of the panel distribution

Panels are defined by their corner points; the coordinates of these points are
measured in a reference coordinate system fixed to the (aircraft) configuration
considered, see fig. 2.2.7. In general the x-axis coincides with the axis of the
fuselage and the x-z plane is a plane of symmetry. The corner points p. are
organised along contours C. as sketched in fig. 2.2.8. Then, corner points are
given as c. p. and a panel is defined by four neighbouring points connected two
by two along lines i = constant, j = constant.
Two successive contours having the same number of points p. define what is
called a strip; a number of successive strips from a segment, see fig. 2.2.8.
Segments may be joined in parts such as wings, fuselages etc., see fig. 2.2.9.
From the foregoing the following hierarchy is apparent: points, contours,
panels, strips, segments, parts.
 2-68

The direction of the normal to a panel

In computing the aerodynamic influence coefficients and setting-up the system of
algebraic equations to be solved the direction of the normal to a panel is an
important quantity. It is important to ensure that on each panel the normal
points into the flow. This may be assured by arranging the points on a contour
and the contours on a segment in a particular sense, illustrated in fig. 2.2.10.
The direction of the normal is then given by the vector product T. x T.1, where
and are vectors along the diagonals of a panel from the corner points
c.p. to C. 1p.1 and from cp.1 to c.1p., respectively. In fig. 2.2.11 an
example is given of normal directions on part of a wing.
In order to have a good procedure for checking or inspecting the correct input
of points on the contour with regard to their coordinates and sequence
(direction of normal), graphical inspection is indispensable. Special software
has been developed for this purpose; inspection may be made by
wire frames (with or without hidden lines), fig. 2.2.12
pin-cushion figures, fig. 2.2.13
solid modeling with shading, fig. 2.2.14

2.2.3 Local representation of the geometry of a panel

As we have seen in section 2.2.1 we need a local panel coordinate system in

which the geometry of' a panel is described and which enables us to work out the
aerodynamic influence coefficients. The local system must be obtained from a
mathematical description of the geometry of a patch or a segment in the
configuration reference system (x,y,z). For example a patch may be represented
by an equation

F(x ,y ,z ) = 0 (2.2.32)

of a certain degree in x ,y and z .

An alternative is to set-up a parametric

representation of the form

 2-69

X = x (s,t) i

= y (s,t) (2.2.33)

z = z (s,t)
p p j
or

= (s,t)

where X is the position vector of a point on the surface, see fig. 2.2.14,
(s,t) is a curvi-linear surface coordinate system and f(s,t) is a polynomial in

s and t.
The relations (2.2.32) or (2.2.33) describing a surface patch may be obtained as
follows:
directly from a general geometry software package of the type mentioned in
the preceeding section.
from a separate, local surface fit (polynomial) through the corner points
of a panel (and, if necessary, through those of adjacent panels).

The first option has the advantage that use is made of knowledge that is already
available. The second one may be benificial from an informatics point of view
since the part of the computercode describing the gemetry as it is required for
the aerodynamic computations is decoupled from the (complex) general geometry
software package.
Depending on the degree of the polynomial representation there are several
possibilities for the local representation of the geometry of a panel.
In figs. 2.2.15 and 2.2.16 representations are considered for planar and
quadratic panels. It is important to notice that planar representations do not
need information of neighbouring panels, whereas for the higher order repre-
sentations such information is needed, complicating matters significantly.

Once the local description of the geometry of a panel in the body reference
system is known the local coordinate system (,n,Ç) of the panel can be chosen.
In general this is done in such a way that the (,ri)-p1ane is tangent to the
panel in its 'centre' see fig. 2.2.17. According to the calculus of vectors the
 2-70

tangent plane is described by the direction vectors obtained from eq

(2.2.33). The Ç-axis is directed along the unit normal in the 'centre' point

a a?

n= (2.2.34)
a a

How, exactly, the 'centre' of a panel is defined is not of great importance. For
planar panels one usually takes the centre of gravity. For paraineterized panels
the centre of the rectangular panel in the parameter space is used.
The 'centre' of the panel acts also as origin of the local coordinate system
(,ri,Ç) and as the point where the boundary condition is satisfied.
The ,ri directions are usually chosen such that the ri-axis is normal to the x-
axis of the body reference system. The t-axis then follows automatically if a
righthanded coordinate system (,ri,Ç) is required. With these choices the
representation (2.2.9) of the surface in the local panel coordinate system
reduces to

= a20 + a11n + a02n2 + 0(h3 (2.2.35)

If the transformation (x,y,z) - (,n,Ç) is known, the coefficients a20,a11,a02

etc. can be expressed in terms of the coefficients of the polynomial given by
eq. (2.2.32) of the local patch considered in the reference system (x,y,z). An
alternative is to interpret eq. (2.2.35) as a Taylor-series expansion of
around the origin of the (,n.Ç) system, then

(P) i
c(n) -
- 2!
_
2
2
ta3no fl + n2 + 0(h3) (2.2.36)

The second derivatives are determined in the following way, see fig. 2.2.18. As

indicated are unit vectors along the x,y,z (reference) axes, respectively;
are the unit vectors along the ,n,ç (local) axes, respectively, i is the
unit normal in the centre of a panel. The ri-axis is taken normal to the x-axis,
thus I=ü. If the vector locates the panel centre (point of expansion in eq.

(2.2.36)), then, since

 2-71

çP = r1 (2.2.37)

TI = - o

ç = - x0).

it can, using differentiation rules and the calculus of vectors, be shown that

a2ç P 2 2 a2p
+ 2 n (2.2.38)
2 2 asat +
as at2
o

as 3t
The coefficients , - can be expressed as
a a

a
as 1 an
m (2.2.39a)
D at - D at

l3r _l ax P (2.2.39b)
a Das Das

where D is the Jacobian determinant of the transformation matrix s,t .+

ap.
as at as at
D= (2.2.40)

an
as at

The scalar products in eq. (2.2.40) can be written in terms of the components of
Q,m in the x,y an z directions as

ax1 az
£ Q
as = as + as + as

ax az
m
as
- -m
as x
+ as my + m as z
(2.2.41)

etc.
 2-72

ax1, ax a2x
The partial derivatives ... etc. can be approximated by
at ' as ' as2
finite differences, for example central differencing

-
(P) (2.2.42)
as i,j - 2s

etc., or from eq. (2.2.33).

a2ç
The other second derivatives in eq. (2.2.36), a"' (an2'o' etc. are ap-
proximated analogously.

It is important to notice that the boundaries of segments need special treatment

with regard to the representation of the geometry of segments in the reference
coordinate system and also with regard to the finite differencing.

Using eqs. (2.2.36) to (2.2.42) has the advantage that the only information
that needs to be stored in the computer memory are the coordinates of the panel
centre. If polynomials like eqs. (2.2.32) or (2.2.33) are used throughout, for
each panel also the polynomial coefficients have to be stored, which amounts at
least to twice as many numbers.

2.2.4 Representation of source and doublet distributions over a panel

In section 2.2.1 we have seen that the aerodynamic influence coefficients are
computed using polynomial-type representations in the panel coordinates ,n of

the source and doublet distributions.

+ + o°'n + 0(h2) (2.2.43)

G(,n) =

ji(,n) = + + p01 n + 0(h2

(2.2.44)

Eqs. (2.2.43) and (2.2.44) may be interpreted als Taylor-series expansions of o

and p about the panel center (0,0)

o(,n) = o(0,0) + -(0,0) + 2(0,0) n + 0(hZ) (2.2.45)

(2.2.46)
pR,n) = p(O,0) + -(0,0) + (0,0) n + 0(h2)
 2-73

The derivatives in these equations can also be expressed in terms of the

parametric variables s,t which are related to the panel coordinate system,
according to

ao as aa at ao
(2.2.147)
a a as a at

8o
an
= as
--ao + at--
an as
ao
an at
(2.2.48)

where, as before, (see eq. (2.2.39))

a
as
- D at
(2 .2 .49)

a
at
fi
a - - D as

and

as_ 1ap
anDat Q

(2.2.50)
a
at
anDas

D is again the Jacobian determinant of the transformation (s,t) -* (,n) given in

eq (2.2.40).

Similar to ase' etc. also , may be approximated by central finite

difference expressions, i.e.

Q. . - o.
1+1,3 i-1,j
- 2s (2.2.51)

etc.

The same procedure may be followed for the doublet distribution, ji.
 2-74

With the representation of the geometry and the singularity distribution for a
panel established, the next step is the evaluation of the integrals appearing in
the integral equations, in other words, the evaluation of the Aerodynamic
Influence Coefficients (AIC's) introduced in section 2.2.1.

2.2.5 Aerodynamic influence coefficients

The aerodynamic influence coefficients may be regarded as the result of the

numerical approximations of the integral expressions for the perturbation
velocity potential and perturbation velocity components induced by source and
doublet distributions (see section 2.2.1). In this section we will study the
derivation of these coefficients in some further detail.
We first consider the simplest case, i.e. that of planar panels with source or
doublet distributions of constant strength.

Perturbation potential induced by a constant-strength source distribution

The perturbation potential ô.. in a point P. induced by a constant source

distribution on a planar panel j is given by

= L -r-
o.jr.. dS = A. .o. (2.2.52)
ji ijj

Thus, if the point is the collocation (mid) point of the panel i, A..o. is

the contribution to the perturbation potential in that point induced by the

source distribution on panel j. From eqs. (2.2.13) to (2.2.16) it follows, for a
source distribution of constant strength, that

A.. =- L jj [()2 + +
ddn (2.2.53)
'j 411ôs
jB

This may be integrated analytically resulting in

¿7

I.
A.
13
.
= L F1
(k+1)
- F1
(k)
)
+
-
k+l)F2 - ki_k)F2(k) +

k=1

(k+1)
+B o (k+1) F - B' F31], k+1 : = 1 if k=4 (2.2.54)
 2-75

where

(k) ((_)2 + çZ)

F1
(k) ik ik -B1
i
+
1k +
ik
i
-fl k
F2' - sinh1
ikj2 + ç2Jh/2

1 ik + B (k) ( fi)
F3 =
[1 + (B1)2J*2 sinh
(k)2 (B1)2) ç2] 1/2
[(B +
o

B' = ik - B1
k
73

B=R
1
-
v.

k+1 -
-

The subscript arid superscripts k indicate the cornerpoints of a panel.

Perturbation velocity induced by a constant-strength source distribution

The perturbation velocity induced by a constant-strength source distribution on

a planar panel may be expressed by the (vector) relation

r
Ó.. = - f! o.v(---) dS
ll
dS = A. .o. (2.2.55)
j i 14n j r.. j
jB jB
where

= k.-) + (T1.-fiJrn +
(2.2.56)
rl = r = [k1-) +
()2 + ç]l/2
J'
 2-76

Hence

-- (A. .)
1J
= (A j
13 = ddn (a)
ôiS

n . -n
= ff ddn (b) (2.2.57)
(A
f13
) . .

4TTÓS
jB
ci
(A) JI -- ddn (c)
ç ij -

The integrals in eq. (2.2.57) can also be evaluated analytically. The result is

r +r -d
21
n -n
log
r +r -d
(1212) 22
-n
log
2 3 23)
- d r1i-r2+d + d r2+r3+d23
12 12 23

n4-n3
log
3434 log (2.2.58)
+ d34 r3+r4+d34 d41 r4+r1+d41

r+r-d r+r-d
4nA
12 log
(12 12) 23 log
2 23)
3
n d r1+r2+d + d r2+r3+d23
12 12 23

r4+r1-d41
log
3434 +
4l log (2.2.59)
+ d34 r3+r4+d34 d41 r4+r1+d41

m e -h 1 m e
12 2
-h 2)
-1 12 1
nA=tan tan
cr
il i2
m
23 2
e -h 2) m e -h 3)
23 3
+ tan1 - tan
çr
12 i3
m e -h 34e4 - h4
3 3) - tan1
+ tan' i3 Ç.r4

m41e4 - h4 m41e1 - h
(2.2. 60)
+ tan' j- tan' ç.r
ji4 il

with
 2-77

d12 = [k 2 - + (fl2 -
1

d23 = [k3 + (n3 -

-
(2.2.61)

d34 = [k4 - + (n4 -

d41 = [k1 - + (n1 -

n2 - n1 n3 - n2
m12 m23
= - - 3 -
(2.2.62)
n14 - n3 n1 - n4
ffl314 a23
- - - -

and

rk = [ki - k)2 + (n. - nk)2 + ç ii k = 1,2,3,4 (2.2.63)

e
k = Ç2.
i +k - k = 1,2,3,4 (2.2.64)

h
k
= (n.i - n k ) k.i
k = 1,2,3,4 (2.2.65)

With respect to eqs. (2.2.58)-(2.2.65) the following remarks are made:

the subscripts (k=) 1,2,3,4 again refer to the cornerpoints of a panel in
the saine sense as in eq. (2.2.54);
for tan defined from - to , eqs. (2.2.58), (2.2.59) and (2.2.60) also
represent the perturbation velocity in the collocation point of the
emitting
. . .

panel itself, including the jump of ±

. . o.in the normal velocity
component according to eq. (2.2.17); in other words

ç. o. o.
(A ).. o.
i
= ii JI o -
r
dS ± -
2
= ±
2
(2.2.66)
iB

(iii) eq. (2.2.60) expresses the spatial angle under which the panel j is seen
from point i.

In the case of Neumann boundary conditions we must know the normal velocity
component in the collocation point of panel i. This component is obtained from
 2-78

o.(
3 fi). =
i J i =
1
. (VA.
13
.)a.
J
(2.2.67)

where

VA..
13
= [(A )..,
ij
).., (AciJ
(Anij )..) (2.2.68)

Perturbation potential induced by a constant-strength doublet distribution

The perturbation potential induced by a constant-strength doublet distribution

on a planar panel is given by

(2. 2. 69)
ô.4. = -
3 1
ff
o.s aç r.. dS = B..
13 3
jB

where

ddn
B.. = - (2.2.70)
13
[k-J + (n.-
+

The expression (2.2.70) is, apart from the minus sign, identical to eq.
(2.2.57c), thus

B.. = - (A
13 j..
1J
(2.2.71)

In other words, the potential induced by a constant-doublet panel is equal to

minus the normal velocity induced by a constant source density on the same
panel.
Recall that the expression (2.2.60) represents the spatial angle under which the
panel j is observed from the point P.. This spatial angle jumps in magnitude by
2ii when crossing the panel surface. Hence eq. (2.2.69) also contains the term ±
p/2 present in eq. (2.2.20).

Perturbation velocity induced by a constant-strength doublet distribution

The perturbation velocity induced by a constant-strength doublet distribution on

a planar panel is, in vector notation, given by
 2-79

ó..
j i
= - 1f p.
j
pÇ (i--)
r..
dS = (
B..) p.
ij j
(2.2.72)
jB

where B.. is composed of

i 3(.-)Ç.ddrì
.---- (B = -1f (2.2.73a)
= .)z
óS i
+
(
i T))2 +

3(n
(B)
liii
=_L j r5
ddn (2.2.73b)
jB
3ç
(Bc).. = -
r
ddn (2.2.73c)
ô .s
jB

For reasons to be discussed in section 2.i6, a different approach, instead of

evaluating eq. (2.2.73), is used in general in which use is made of the
equivalence of doublet and vortex distributions (section 2.1.6). This means that
eq. (2.2.72) is replaced by

'r x XdS
j i
= 1f ' ddri -
-
j (2.2.74)
ôiSB kI aô.S Ii:.13

where see eqs. (2.1.100/101).

=
Upon introduction of a unit vector g

= g + gin

describing the direction of a line element ã along the circumference of the

panel c5.S3, so that

as = gds

and

gds = ds =

gds ari
ds = dri,
 2-80

the individual velocity components may be written as

(- c)
= L jj ddn + pg ds (2.2.75a)
,J
i LIuós
jB
nr3 aó.S
jB
r3

o.'
j n.
=- J!
os T!
'r
(- ç.)
r' ddn
i

ao.s
r3
(2.2.75b)

jB 3B
n.-n
o.) = - 'r ) ddn +
.
' o'
3B
-(n-n)
( 1 1
(2. 2. 75c)
pg + pg j ds,
4n r3 r3
ao s
jB
-

where'r =Qand1 =. T'

Note that in eq. (2.2.75) the surface integrals are the saine as those in eq.
(2.2.57). This means that the velocity components induced by a doublet
distribution on a panel may be expressed in terms of the velocity components
induced by a source distribution on the panel plus those induced by a line
vortex distribution along the circumference of the panel.

For a constant-strength doublet distribution we have = O and therefore

'r = O. Hence only the line integrals in eq. (2.2.75) are retained. In that case
the velocity components may be written as

= (s°).. p
(2. 2. 76a)

= (°).. T' ij .
j
(2.2.7 6b)

(2.2.76c)
3cl
= ijj
(B°)..
ç

where, from eqs. (2.2.7'4) and (2.2.75),

g ç
(B°)..
T' ds -
r3 - 4ii g; r3
(2.2.77a)
13
aó aó s
jB
jB
 2-81

r
ds= r g 13
(2.2.77b)
36 36 rL
jB JB
g -(n.-n) g
(B°).. -
ç jj - 4ii
L
r3
ds
+ -; L
r3
ds
36 S
jB ao .S
jB
_L
-
g
ÓS
L
n.-n
r3 LIn g
L
r
d (2.2.77c)
JB jB
The integrals in eq. (2.2.77) can be evaluated analytically, giving

LI r i-r +-
k k+1 k+1 'k
adS
L
r3
d =
k=l
X Qn
rr
k k+1
-
kk+1-k)
(2.2.78)
jB
rr + -
kl k +.1
k k+1 'k+1 'k
L -d
r3
X {r r in
rr (2.2.79)
k=l k k+1
jB

and a similar expression for L dn with replaced by n.

r

In eqs. (2.2.78) and (2.2.79)

+
rk = ik)2 + (2.2.80)

Multi-pole expansions

As we have seen in the preceeding paragraph the aerodynamic influence coef fi-
dents of constant-strength source and doublet distributions on planar panels
can be determined in an exact way by evaluating the integrals analytically.
However, their computation is costly due to the large amount of coefficients to
be determined (N*N, N=0(103)). In addition, these coefficients contain a large
number of logarithms and inverse tangents.
The computing time can be reduced considerably by using, instead of the exact
expressions, series-expansions of

1
r
= 1( _2 + (._)2
i + (ç1_çJ2]*2 (2.2.81)

in terms of
 2-82

L
r
= [2

i
+ ri2.
i.
+
ii (2.2.82)
o

The, truncated, series expansion may be written as

i
iL{'+-
r r
o
r2
o
+
112 +
o
r
o
(2.2.83)

where h is a measure for the panel dimension (ref. 2.10). The first term at the
righthand side of eq. (2.2.83) represents the potential of an elementary point
source (monopole), the second and third term that of a point-doublet with axis
in c-direction and Ti-direction, and with strengths and n, respectively. The
terms of o(P) represent quadrupoles, etc.

Note that, when integrating over the panel, the second and third term at the
righthand side of eq. (2.2.83) vanish when =n=o represents the panel centre.

Substitution of the series-expansion (2.2.83) into the integral expressions for

the aerodynamic influence coefficients yields integrals of the type

mri
= ddn ;
m,n,p = 0,1,2, (2.2. 8'fl

o
The evaluation of these integrals results in relatively simple expressions in
and n which are considerably less expensive to compute than the logarithms and
inverse tangents occurring in the exact aerodynamic influence coefficients. The
number of terms to be retained in eq (2.2.83) depends on the ratio -. In

practice it turns out that for first order methods the first term suffices if
< 0.25. This means that about 90% of all exact influence coefficients may be

replaced by monopole approximations. 0f the remaining 10% about one half can be
approximated by quadrupole expressions.
The multipole expansions are used in the following methods (see also table 2.1,
section 2.1.13)
DOUGLAS-NEUMANN (Hess & Smith)
BOEING TEA 230
NLR PANEL
MBB PANEL
BAC PANEL
PANAIR
Hess & Friedman
 2-83

Higher-order approximations; 'small-curvature' expansion

If the local geometry of the panels and the singularity distributions are
approximated by functions of a higher degree than piecewise constant, integrals
of the following type result (see eqs. (2.2.13) and (2.2.14))

2
m n r ]1/2
n I1 + (---) + (.---)
If +2 ddn (2.2.85)
[k i J
2

+ (.-)
i
2

+ (c.-c)
i
1p/2

where (see eq. 2.2.9)

= a202 + a11n + a02n2 + (2.2.86)

Integrals of the type given by eq. (2.2.85) can be evaluated only by either
numerical integration or by approximation through suitable forms of series-
expansion.
The first approach (numerical integration) is found in refs. 2.14 and 2.23. In
these methods numerical integration methods are applied (6-or more point
Gaussian quadrature (see also ref. 2.24). The computation of the aerodynamic
influence coefficients by numerical integration is relatively expensive.
The series-expansion approach (reIs. 2.16. 2.17 and 2.25) is, although less
expensive in terms of computing time, more expensive in terms of code develop-
ment than the direct numerical integration approach. Due to the more complex
derivations and computer coding involved also the possibility of errors
increases.

In the series expansion approach the denominator of the integrand in eq.

(2.2.85)

2 2

(n.-n) + (ç_ç)2]P/2 = r
[k-) + (2.2.87)

is expanded in a power series. We have seen in the preceeding paragraph (eq.

(2.2.83)) that for r>>h, where r is the distance between the field point
and the collocation point of the panel considered (eq. (2.2.82)), the
term 1/r can be expanded in terms of 1/r (multipole-expansion, eq. (2.2.83)).
 2-8h

Also for r,r = 0(h) an expansion in a power-series is possible. A convenient

expansion may be obtained by writing

r2 = k1_)2 + (_)2 + (ç_ç)2 = r. - 2Ç + Ç (2.2.88)

where rf is the distance between the field point (.,n.,Ç.) and the projection
(,n,0) of a point (,n,Ç) of the panel on its tangent plane. Then can be

expanded as

1 Ç2
(2.2.89)
r r rfrf 2 r1 -):-f::-

This series expansion is valid for Ç/rf « 1 and therefore called 'small
curvature expansion'. Note that for r = 0(h)

rf = 0(h)
= 0(1)

Ç/r1 = 0(h)

because

Ç = a202 + a11n + a02n + . = 0(h2)

For zero curvature (flat panels) Ç = O and rf = r the full expression

for hr is

recovered.

form
Substitution of eq. (2.2.89) into eq. (2.2.85) yields integrals of the

mn
I ddn (2.2.90)
mnp
r

Integrals of this type can be evaluated analytically (see ref. 2.17).

As the accuracy requirements become more stringent, the complexity of computa-
0(h)) one
tion increases rapidly: in a first-order method (truncation error
integral of the type (2.2.90) has to be evaluated for each near-field
aerodynamic influence coefficient; a second-order method (truncation error
0(h2)) requires six such integrals to be determined and a third-order method
,
twenty-three.
 2-85

However, the computing time increases considerably less rapidly with increasing
power of the truncation error as a result of the fact that only about 10% or
less of the aerodynamic influence coefficients are of the near-field type.

Survey of relations between aerodynamic influence coefficients

Recalling that the aerodynamic influence coefficients may be considered as the

potential (eqs. (2.2.52), (2.2.69)) or the velocity (eqs. (2.2.55), (2.2.72) and
(2.2.77)) in a (collocation) point i induced by a singularity-distribution of
th
unit strength on the panel, it appears, upon inspection and comparison of
the various expressions that

(i) the potential induced by a doublet distribution of given strength is equal

to the normal velocity component induced by a source distribution of the
same strength and also equal to the tangential velocity component induced
by a vortex distribution of the same strength
or, symbolically

(p) = an(o) - 3s
(2.2.91)

From this it follows also that

a2
(u)
as
- anas (G) - as
(2.2.92a)

and - (ti) = (2.2.92b)

= anas
an2

(ii) the normal velocity component induced by a vortex distribution of given

strength = the tangential velocity component induced by a source
distribution of the same strength;
or
= (2.2.93)

also
O(o) = J -(T) ds (2.2.94)
 2-86

The conclusions (i) an (ii) should not be confused with the equivalence princi-
ple of doublet- and vortex-distributions.

2.2.6 Truncation errors; consistent approximations

An important question to be answered is: what degree (of polynomial) must be

chosen for the local geometry representation (eq. (2.2.36)) and the singularity
distributions (eqs. (2.2.45) and (2.2.46)) in order to realise a given accuracy
of the solution of the discretized integral equation? In this context the notion
of consistent discretization plays an important role.
In numerical analysis, a (numerical) discretization is said to be consistent if
the local truncation error tends to zero when the step size or meshwidth h tends
to zero. Symbolically this can be stated as follows.
If

A =
(2.2.95)

represents the analytical or continuuum problem, where i is the vector field to

be solved for, ? a known vector field (forcing function) and A an operator

working on u, and

h-h -h (2.2.96)
Au= f

represents the discretized problem, then the discretization scheme is called

consistent if, for h-O

A_Ahi max (A_Ah) k ha, for h 0 (2.2.97)

I IÓAI i
I
o
i ui 1=1

Here ii H represents the L 2

or Eulerian norm and k and a are finite and
o
positive constants (ref. 2.24).
For panel methods the consistency aspect has been studied by, a.o., Hess (ref.
2.26), in particular with the purpose of finding out where the geometry
representation and singularity representations can be truncated if their
ha is to be of the sane power in h.
contribution to the local truncation error k
In this section we will address this aspect in some detail, but in a somewhat
more general sense than Hess did.
 2-87

The potential induced by a source distribution on a panel

The potential induced by the source distribution on a single panel ois is given
by (see preceding sections)

=
1 1ddn (2.2.98)
3 i - kn r InI
i

where
r = [(_)2 + ()2 + (2.2.99)

If
r = [2 + n2. + (2.2. 100)
0 1 ]

then
2n.n 2ç.ç
r = r
o
[i
r2 + r2 - r2 +
2

i r2
+ j] 1/2
2

(2.2.101)
O O O Oi O O

where F,n are 0(h) and ç = 0(h2) + higher order. The terms to the right of the
vertical dotted line are non-zero only in case of curved panels.
(hZ jh/2
Since 0(ä.S) = O J
= 0(h2) (i + 0(h2) (see eq. (2.2.10)), it
ç
follows that

0(ó) = O(h2 a (i +l0(h2)+ (2.2.102)

Now two situations may be distinguished

1. r(.,n.,ç.) = 0(h); 'near-fieldt (neighbouring panels).
In this case the terms to the left of the vertical dotted line in eq. (2.2.101)
are of 0(1) and those to the right of the vertical dotted line are of 0(h) and
0(h2), respectively (plus higher order).
Isolating the (panel) geometry dependent part of (2.2.102) it is found that

o(hz 1 1
= o[
h2
f10(h2)+ jl/2j
r ii h[0(1)+ 0(h) + 0(h2)+ ...0(hZ)...)'/2

(2.2.103)

For planar panels (ç = 0) we have

 2-88

1
o(hz = 0(h)
r inçi

with a truncation error of 0(h2).

For quadratic panels (ç = O(h2)) we find, through power series expansion of the
denominators,

i
o(h2 = 0(h) + 0(h2) + ... 0(h3)
r ini

Verify that the expression that is obtained contains some but not all terms of
2 2
3 3

...). For example, with ç=a2 ç +a11Çn+a 2

.0(h the
0(h ), (indicated by ) , .

term ç2/r gives rise to a term of 0(h3). However (truncated) cubic terms in the
geometry representation are not present.
Hence, the lowest order truncation error is 0(h3).
Multiplying (2.2.103) with the source density representation o = 0(1) + 0(h) +
the accuracy and truncation error (TE) levels of the approximation for ÓA'.
are found to be as follows:

For flat panels with a constant-strength source distribution:

0(&t) = 0(1) * 0(h) = 0(h)
0(TE) = 0(1) * 0(h2) + 0(h) * 0(h) = 0(h2)
Truncation errors due to geometry and singularity distributions are balanced
(both 0(h2)).

For flat panels with linear source distribution:

0(ó) = {0(l) + 0(h)} * {0(h)} = 0(h) + 0(h2)
0(TE) = 0(1) * 0(h2) + 0(h2) * 0(h) = 0(h2) +
Note that the truncation error is still 0(h ) . Hence, the terms 0(h2 ) in the

full expression for ô are superfluous (do not really help); contributions to
the truncation error of geometry and singularity representations are not in
balance. This combination does not offer any real advantage over the preceeding

one.

For quadratic panels with linear source distribution we find

0(64) = {0(1)+0(h)}{0(h)+0(h2)+ } = 0.(h)+0(h2)+

0(TE) = 0(1)*0(h3)+0(hz)*0(h)+ = 0(h3)

 2-89

superfluous terms 0(h3) and higher order in full expression for ó4; truncation
errors balanced.

For quadratic panels with quadratic source distribution

0(ó) = {0(1)+0(h)+0(h2)}*{0(h)+0(hz)+ }= 0(h)+0(h2)+

0(TE) = 0(l)*0(ha)+0(h3)*0(h)+ = 0(h3)+

Superfluous terms 0(h3) and higher order in full expression for ô; truncation
errors not balanced. This combination does also not offer a real advantage over
the preceeding one.

2. r
o
= 0(1), far field.
Here the magnitude of the subsequent terms in eq. (2.2.101) is found to be

r = r (1 + 0(h) + 0(h2) + 0(h) + 0(h2)+ 0(h2) + 0(h)...0(h)....)

so that
0(h2 )=o{ h2 *
r InI )1/2
r [i + 0(h) + 0(h2) + 0(h2) + 0(h )+. . .0(h3)...
o

* + 0(h2)+ J'} (2.2.104)

This is = 0(h2) * (1 + 0(h)+....) with truncation error 0(hz)*0(hz), for planar

panels;
* (1+ 0(h)
and = 0(h2) + 0(h2) + ....), truncation error 0(hz)*0(h3), for
quadratic panels.

In this 'far field' case where we have to do with many, 0(N), panels contribut-
ing to the total perturbation potential 4 = &D. Hence the summation over
0(N) = 0(j-) panels cancels the factor h2, and we find:

For planar panels with constant strength source distribution

0(e) 0(1){1+0(h)+.. . .} = 0(l)+0(h)+.

0(TE) 0(h)*0(l)+O(h2)*0(1) = 0(h);
 2-90

superfluous terms 0(h) and higher order in full expression for 0; geometry and
singularity related parts of truncation error not in balance;

For planar panels with linear source distributions

0(0) = {0(i)+0(h)}*{l+0(h)+....} = 0(1)+0(h)+0(h2)+....

0(TE) 0(hz)*0(l)+0(h2)*0(l) = 0(h2);

superfluous terms 0(h2) and higher order in full expression for 0; geometry and
singularity related parts of truncation error in balance.

For quadratic panels with quadratic source distributions

0(0) = {0(1)+0(h)+0(h2)}*{1+0(h)+0(h2)+....} = 0(1)+0(h)+0(h2)+...

0(TE) 0(h3)*0(1)+0(h3)*0(l) = 0(h3);

superfluous terms 0(h) and higher order in full expression for 0; geometry and
singularity related contributions to truncation error in balance;

From this 'error analysis' it appears that in order to have a local truncation
error of 0(h) in the potential it is (more than) sufficient to use planar panels
with a constant-strength source distribution. The analysis also suggests that
the computation of the near field influence coefficients as well as in the
computation of the far field influence coefficients can be approximated through
series expansions (such as multi-pole expansion).
In order to achieve a truncation error of 0(h2 ) planar panels with a linear

source distribution are (also more than) sufficient. In this case also
simplification of the integrals by means of power series expansions is possible
(in fact, planar panels with constant source density would be sufficient in the
near-field).
 2-91

The potential induced by a doublet distribution on a panel

In this case the elementary contribution due to a single panel is given by

= _L11 a 1
()
ddn ddn
(2.2. 105)
¡n r ÍnI
i J

Following a similar analysis as for a source distribution, we arrive at

ç.-ç
0(ô) = o [h2 [l+0(h2)+..j'/2]

Again two cases are considered:

1. r (.,n., ç.) = 0(h), near field. Then

o i i i

0(h2
1 0(h) + ;0(h2) + ... *
= o[h
r3
h{O(1) + 0(h) + 0(h2)+...0(h2)..j3/2

)1/2]
* (i + 0(h2) (2.2.106)

where the dotted lines have the same function as before. Then it is found that

ç.-ç
i
0[hz = 0(1) + T.E. of 0(h), for planar panels;
r3 InI
= 0(1) + 0(h) + + T.E. of 0(ha), for quadratic panels

and for planar panels with constant doublet density

0(d) = O(i)
0(TE) = 0(h)*0(l)+0(h)*0(l) = 0(h);

for quadratic panels with linear doublet density

0(ä) = 0(l)+0(h)}*{0(l)+0(h)+. . .} = 0(1)+0(h)+

0(TE) = 0(hz)*0(l)+0(h2)*0(l) = 0(h2)

for cubic panels with quadratic doublet density

 2-92

0(6e) = 0(1)+0(h)+0(h2)+....
0(TE) = 0(h3)

Note: A particular 'near field' case arise when the point (.,ri.,ç.) is located
on an (adjacent) panel. For a (continuous geometry we than have Ç= 0(h2),
rather than 0(h), and the degree of the truncation error changes
correspondingly.

2. ro = 0(1), far field

o(h2 i = o[h 0(l);+ 0(h2) + *

r3 ¿n 3/2
r3
o
[i + 0(h) + 0(h2)i+ (h2)+0(h4)+...0(h3)
jl/2
* [1 + : 0(h2) j (2.2.107)

= 0(h2) * {0(l) + 0(h) + . . + T.E. of 0(h2)), for planar

panels;

= 0(h2) * [0(1) + 0(h) + 0(h2) + . . + T.E. of 0(h3), for

quadratic panels.

Combination with doublet density representation and summation over N panels

gives

0(e) = 0(1) +.. . + T.E. of 0(h), for planar panels with constant-strength
doublet distribution;
= 0(1) + 0(h) +.. + T.E. of 0(h2), for planar panels with linear doublet
.

distribution;
= 0(1) + 0(h) + 0(h2) +.. + T.E. of 0(h3 ), for quadratic panels with
.

quadratic doublet distribution.

Conclusion: Planar panels with constant-strength doublet distribution are

sufficient for having a truncation error of 0(h). Quadratic panels with linear
doublet distribution are (more than) sufficient if a truncation error of 0(h2)
is pursued. It is possible to simplify near field as well as far field
integrals.
 2-93

A summary of the results for the integrals for the potential is given in Table
2.2.

Consistent representations

geometry source doublet local (summed) examples of

distribution distribution truncation error methods

MCAERO *)

planar 1)2) constant1)2) constant2) 0(h) VSAERO

QUADPAN

planar1 )2),source linear 1)2) linear )2) 0(h2) PANAIR *)

quadratic1 )2 ) ,doublet AEROPAN

quadratic 1 )2 )source quadraticl)2) quadratic' )2) 0(h3)

cubic')2) , doublet

simplified near-field integrals possible

idem far-field
*) quadratic doublet distribution

Table 2.2 Consistent representation of integrals appearing in the expression for

the potential
 2. 9

The velocity induced by a source distribution on a panel

The expression for the contribution of a single panel is:

=
1 ddn L O
r ddrt
lnI
(2. 2. 108)
r
os i
-
i

1
0(óø) = 0h2 o

1. r (.,n.,ç.) = 0(h), near field

o ] i i

o(h2 L
r2 Jn1
1
= o[hz h2{0(1)
+
1
0(h) + 0(h2)+ ..O(h2)...)
*

* 0(h2)+...J'/2)
[1+ (2.2.109)

= 0(1) + T.E. of 0(h), for planar panels;

= 0(1) + 0(h) + . . + T.E. of 0(h2), for quadratic panels.

Then
0(óV) = O(i) + T.E. of 0(h), for planar panels with constant-strength source
distribution;
= 0(1) + 0(h) + . . + T.E. of 0(h2), for quadratic panels with linear
source distribution;

= 0(1) + 0(h) + 0(h2) + . . + T.E. of 0(h3 ) , for cubic panels with

quadratic source distribution.

2. r = 0(1) far field

o

0(h2 L
r2
1
)
= o[h2
r2 [i+ 0(h) + 0(h2) +
1
0(h2) + 0(h )+. . .0(h)...)
*

n o
ç

* (1+ 0(h2)+...)'/2) (2.2.110)

 2.95

= 0(h2) * (0(1) + 0(h)+...) + T.E. of 0(h2)*0(hz), for planar

panels;

= 0(h2) * (0(1) + 0(h) + 0(h2) + ....)) + T.E. of

0(h2)*0(h3), for quadratic panels.

Combination with source density representation and summation over N =

panels gives

0(u) = 0(1) + . . + T.E. of 0(h), for planar panels with constant-strength

source distribution;
= 0(1) + 0(h) + . . + T.E. of 0(h2), for planar panels with linear
source distribution;
= 0(1) + 0(h) + 0(h2) + + T.E. of 0(h3), for quadratic panels
with quadratic source distribution.

Conclusion: Planar panels with constant-strength source distribution are

sufficient if a truncation error of 0(h) is allowed. Simplified far-field
integrals are possible. Quadratic panels with linear source distribution are
more than sufficient if a truncation error of 0(h2) is allowed. Simplified near-
field and far-field integrals are possible.

The velocity induced by a doublet distribution on a panel

Contribution due to a single panel

c.-c
ô..
j i-= - 4 li f1
i
ji
ji.V.
a

i
i dd
Jn1 -
i

i
j i r
i ddn
Inc1

- L i i ddn (2.2.111)
- 4ii
a.S j r5 + r3

i .c.
0(ó$) = 0(h2p + 0(hlJ
r3
 2.96

1. r = 0(h), near field

o i i i

ç.-ç
1 1 0(h) + 0(h2)+ ....
o(h2 - o[h * (2.2.112)
r - h0(1) + 0(h) + 0(h2)+...0(h2)..j2
1/2
* (i + 0(h)+...)

= o() + T.E. of 0(1), for planar panels;

= o(i) + 0(1) + . . + T.E. of 0(h), for quadratic panels.

= 0(i) + 0(1) + 0(h) + + T.E. of 0(h2), for cubic

panels.

1
The saine holds for the 0(h2i term.

Then
0(ô) = o() + T.E. of 0(1), for planar panels with constant-strength
doublet distribution;
= o(i) + 0(1) + .... + T.E. of 0(h), for quadratic panels with
linear doublet distribution;
= o[) + 0(1) + 0(h) + ... + T.E. of 0(h2), for cubic panels with
quadratic doublet distribution.

From this analysis we conclude that we do not have a consistent approximation of

our problem, i.e. the error does not vanish if' the meshwidth (panel dimension)
approaches zero. Moreover the integrals themselves lead to infinity for h + 0.
The problem can be circumvented, at least partially, by switching from the
doublet formulation to the equivalent vortex formulation introduced in section
2.1.6. In order to see this we rewrite eq. (2.1.100), for a single panel, as

=r. =ï.
3---
i
L
14n
rxds
r3
+-ff 1
L4
- r
{.pxn}x-
j r3
dçdn (2.2. 113)
j j
36 S
i i

and recall that the line integral represents the velocity induced by a line
vortex with strength p. around the circumference of the panel and the surface
integral represents the velocity induced by a distribution of vorticity with
strength p on the surface of the panel.
An order of magnitude analysis for the surface integral of (2.2.113) leads to:
 2.97

o{hz 1
i o {h
0(h) + 0(h2) + *
Il
Inl
h3[0(l)+0(h)+0(h2)+...0(h3)...)3/2

1/2
* (140(h2 )+. . .)

and

o{ ¡J ]
= 0(1) + T.E. of 0(h) for planar panels with constant
ô .s
J
vorticity (linear doublet) distribution

= 0(1) + 0(h) + ... + t.e. of 0(h2) for quadratic panels with

linear vorticity (quadratic doublet) distribution

= 0(1) + 0(h) + 0(h2) + ... + T.E. of 0(h3) for cubic panels

with quadratic vorticity (cubic doublet) distribution.

Performing a similar analysis for the line integral it is found that this is
still 0(1/h).
Hence, it seems that we have solved the problem only partially.
The remaining problem with the line integral can be solved by combining
contributions of adjacent panels. For this purpose we consider the contribution

rxds
(2.2. ll24a)
i ,j=1 j r3
J

of the line integral of all panels to V.. Noting that adjacent panels have
partially coinciding edges ad that the ã at the common edge are of opposite
direction it follows that the expression for can be replaced by

°) ;xa (2.2.114b)
i =j=1 SnaókS
r3

Note that this represents a summation over the (common) panel edges rather than
over the panels themselves.
The virtue of using (2.2.l14) is in the fact that each of the individual
integrals of (2.2.114) remains finite, even for h40. This is due to the fact
that for each (common) edge we have
 2.98

= 0(h) with TE = 0(h) for p = constant per panel

J
aa.sAaoks TE = 0(h2) for linear p-distribution
TE = 0(h3) for quadratic p-distribution.

It appears further, upon substitution of Taylor series representations for the

p. and at the common edge (with p expanded around the respective panel mid-
points) that the truncation error of the integrals (2.2.114b), is dominated by
that of the doublet distribution. This means that, in general, the
discretization is consistent only for linear or higher degree representations
u(o)
for p; for p = constant the truncation error of is still 0(1). From a more
detailed analysis for p = constant panels (which is left to the reader) it
appears that the normal component of Q (but not the tangential component!)
happens to acquire a truncation error of 0(h) provided the collocation point is
positioned at panel centres.

The property just discussed forms the basis of the so-called vortex-lattice-
(constant-strength doublet) method of thin-wing theory. The accuracy of the
vortex-lattice method can be increased further by shifting the position of the
vortex pair and the position of the collocation point to the 1/4 and 3/4 panel
chord point, respectively. Note, however, that the tangential velocity component
which should exhibit a discontinuous behaviour when passing through a panel, can
only be modelled properly through the local gradient of the doublet strength.
Hence, a linear doublet representation, at least locally, remains necessary for
a consistent evaluation of the near-field tangential velocity. For more details
on the vortex-lattice method the reader is referred to ref. 2.27.

2. r (.,rì.,ç.) = 0(1), far field

o i i i
In this case one finds

1 0(1) + 0(h2) + ... *

o(hz = o[h2
r (i + 0(h) + 0(h2) + 0(h2) + 0(h)+...0(h))
I 1/2
0(h2)+...} (2.2.115)

= 0(h2) * [o(i) + 0(h) + ...) + T.E. of 0(h2), for planar panels;

= 0(h2) * [O(i) + 0(h) + 0(h2) + ...) + T.E. o f O ( h 3

) , for
quadratic panels.
 2.99

The factor 0(h2) again disappears by summation over N panels. Hence

0(V) = 0(1) + + T.E. of 0(h), for planar panels with constant-strength

doublet distribution;

= 0(1) + 0(h) + ... + T.E. of 0(h2), for planar panels with linear
doublet distribution;
= 0(1) + 0(h) + 0(h2) + + T.E. of 0(h3), for quadratic panels
with quadratic doublet distribution.

Conclusion: In general planar panels with linear doublet distribution are neces-
sary to achieve a truncation error of 0(h). Special measures (equivalent
vorticity formulation) are needed in the near field representation to keep the
perturbation velocity bounded for h -+ 0. The far-field integrals may be
simplified by power-series-expansions.
In order to attain a truncation error of 0(h2) the panel geometry and the
doublet distribution should both be quadratic. Both near-field and far-field
integrals may be simplified.
A summary of the results for the integrals of the perturbation velocity is given
in Table 2.3.
 2.100

Consistent representation

geometry source doublet local (summed) example of

distribution distribution truncat. errors methods

DOUGLAS-NEUMANN
planar 1)2) constant2) linear2) 0(h) BOEING TEA 230+)
NLR+), MBB+), BAC+

quadratic1 )2) linear' )2) quadratic 0(h2) Hess & Friedman

1)2)

cubic 1 )2) quadratic cubic 0(h3) Robert & Rundle*)

1)2) 1 )2

(BAC)

simplified near-field integrals

simplified far-field integrals
+) constant-strength doublet panels (large errors for large panel dimension/wing
thickness ratio)
*) cubic source distribution

Table 2.3 Consistent representation of integrals appearing in the expression for

the perturbation velocity
 2.101

2.2.7 Convergence and stability

In the preceding section we have seen that the concept of consistency tells us
something about the accuracy of the discretized representation of the integrals
(local truncation error) for given distribution of the source or doublet
densities. However, nothing has been said about the accuracy of the solution of
the discretized integral equations. For a judgement of this aspect we need the
concept of convergence, (ref. 2.24).
A numerical method is called convergent (see also eq. (2.2.97)) if, for h O

Iläull = I
- u-h J <kh (2.2.116)

where is the solution of the continuous problem.

A = (2.2.95)

-h
and u that of the discretized problem

A
h-h
u =f-h (2.2.96)

and k1 and are finite, positive quantities.

h
The difference is called the global truncation error.

A second concept playing a role in the numerical solution of equations is the

concept of stability (ref. 2.24). A numerical method is called stable if

-h h-h
Ilu I
<k2 lIA u J
(2.2.117)

where k2 is also a finite positive quantity independent of h.

It may be shown that for a consistent discretization of an equation, a numerical
method is convergent if it is stable (Laxe equivalence principle, see ref.
2.24). Stability is a property of the algebraic (matrix) operator Ah and thus of
the integral equation as well as of the way of discritization.
From eqs. (2.2.96) and (2.2.117) follows

-h h h-h
H = lEA j-' hll < I[Ah]lH IA u (2.2. 118)
 2.102

Therefore, stability requires

(2.2. 119)
I

To understand that this is also a condition for convergence when we have a

consistent discretization, we make the following analysis. From eq. (2.2.116)

and eqs. (2.2.95), (2.2.96) it follows that

-h = [A'] h = [Ah]A for h 4 0

h
(for a consistent discretization for h 4 0)

Hence

- -h- hi-1 h -
óu=u-u 4LAJ r
IA-Aju for h 0 (2.2.120)

Taking the norm

IIoII = I[Ah] (AhA) II

{Ah] (AhA) II IIH (2.2.121)

o r

Ioii I
lull
{] 1AhA) H
[AhJ
I I
IAAhI I
(2.2.122)

A consistent discretization requires (see eq. (2.2.97))

I
lAAhl I
= kha, for h 0 , II =i

Hence

IloII Il{Ah]IIk0h0, for h O , HUH = i (2.2.123)

Comparison with eq. (2.2.116) yields that if

k
h°, for h 0 (2.2.124)
II{A']ll <
o
 2.103

the numerical method is also convergent. In particular there holds that the
global truncation error in the solution and the local truncation error in the
integrals are of the saine order in h if

h-1
[A ] H is bounded (arid non-zero) for h O (2 .2. 125)

Eq. (2.2.125) expresses the same as eq. (2.2.119) and we may conclude that
stability guarantees convergence to the same order in h as the local truncation
error.
It may be remarked also that the norm I[Ah]_lIj is related to the condition
number of the matrix Ah (ref. 2.24) by

cond (Ah) = j[Ah]l (2.2.126)

Hence, in order to have stability a well-conditioned matrix is required.

From the preceeding analysis it follows that a good estimation of the norm of
the inverse matrix [Ah] is an essential ingredient of the stability analysis
of a numerical method. Unfortunately, good estimates of' such norms are generally
not available in the case of panelmethods. If they exist it concerns idealized
situations in two-dimensions, in which the matrix Ah can be reduced to diagonal
from using Fourier transformations.
Because of this lack of information about the norm of the inverse matrix the
estimation of the global truncation error of panel methods is in general based
on empirical/heuristical arguments. For example, it has been established
empirically, that the stability problem of' panel methods corresponds, in certain
aspects, to the problem of approximating a function by means of' piecewise defined
polynomials. In both cases lack of stability manifests itself in high-frequency
point to point oscillations. In the following we shall consider a number of'
stable and unstable formulations. In that respect we will distinguish between
(see fig. 2.2.19)
- absolute stability, local errors are damped and do not propagate through the
entire solution
- neutral stability, a local error propagates but does not increase
- instability, local errors propagate and increase in magnitude.
 2.104

Analo between the solution of piecewise approximated integral equations and

the piecewise approximation of functions through polynomials ('splines')

Spline-fit problem (see fig. 2.2.20)

Approximate a function f(x) that is given by function and/or derivative values
in a finite number of points, by means of piecewise defined polynomials f.(x)
defined by

(x-x.) + ...; for x. x x. i=0,l,etc.(2.2.127)

f.(x) =
:1.
f°.
1
(x.)
1
+ f'.
1 1 i i+1
,

The problem leds to a system of linear algebraic equations for the coefficients
f'., etc.
f°.,
J_ i
Fredhoim integral equation of the second kind

Equations of this type have the form (see section 2.2.1)

pi ç'-c
Z pK2dS K2 - (2.2.128)
i os.
J

o. --
oK4dS ; K4 (2.2.129)
i os.
J ii3

Inspection of these expressions teaches that the terms p.12 and o/2 are
strongly dominant and, in fact, independent of the panel dimension h (K2 and K14
are zero for flat surfaces).
Then, if p,o are piecewise defined polynomials of the type f.(x), mentioned
above, prescribing is (almost) equivalent with prescribing p. and, hence,
.

corresponds with prescribing f. in the spline-fit problem. The same conclusion

is obtained in case of prescription of with o as unknown.

Fredholm integral equation of the first kind

Here we have equations like

oK dS K =1 (2.2.130)
' ;
1 1 r
jóS
i
 2.105

r.n
i
= p1K dS K3 (2.2.131)
ani jj5
J as.
J
- + Ii3
In these cases a clear, strongly dominant term is absent. In fact, inspection of
the kernel K1 teaches that (even for flat surfaces) all source panels contribute
significantly to ' with the largest contributions coming from the 'near-field'
panels. A more detailed analysis (see section 2.2.6) teaches that the local
contribution of o itself is weakly dominant (0(h); the contribution of o' being
0(h2)). The implication of this is that solving the source/Dirichlet problem
(2.2.130) is equivalent, but in a weak sense only, to solving the spline-fit
problem with f. given. The same conclusion is obtained with respect to the
doublet/Neumann problem (2.2.131). (Remember, from section 2.2.6, that the near-
field contribution of p is 0(1/h) and 0(1) for p').

In the equivalent vortex formulation we are dealing with an integral equation of

the type

= I f! X dS
II
( )
.j os

Analyses teaches that the local contribution of I to is zero (a constant

vorticity panel induces zero normal velocity at its center) and the local
contribution of I' is weakly dominant. Hence, solving the vorticity/Neumann
problem is weakly equivalent to solving the spline-fit problem with f'. given.

Summary of the most important stable and unstable formulations

Figs. 2.21-2.2.27 ilustrate a number of stable and unstable formulations. They

reflect the rule of the thumb that:
- formulations resulting in a finite non-zero local contribution of the
singularity distribution are usually stable
- formulations resulting in a zero local contribution of the singularity
distribution are usually unstable.
It is further noted that in the integral equation case it is not possible to
specify f. or f at panel edges because of the singular behaviour of the
integrants.
 2. 106

Examples of absolutely stable formulations:

source distribution - Neumann boundary condition

Iconstant
constant doublet distribution - Dirichlet b.c

L2 linear source dist. - Neumann b.c

L2 linear doublet dist. - Dirichlet b.c
Fredhoim
integral eqs. CL2 quadratic source dist. - Neumann b.c
of 2nd kind C L quadratic doublet dist. - Dirichlet b.c
o2

L2 quadratic source dist. - Neumann b.c

L, quadratic doublet dist. - Dirichiet b.c

Neutrally stable

constant source dist. - Dirichiet b.c

Fredholm doublet dist. (line vortex) - Neumann b.c

j constant
integral eqs. L2 linear source dist. - Dirichlet b.c
of ist kind L2 linear doublet dist. - Neumann b.c
C linear vortex dist. - Neumann b.c
2nd kind fC1 quadratic source dist. - Neumann b.c
C1 quadratic doublet dist. - Dirichiet b.c

ist kind C L quadratic vortex dist. - Neumann b.c

o2

Unstable

ist kind constant vortex dist. - Neumann b.c.

ist & 2nd kind C linear source dist. - Neumann and Dirichlet b.c
o
C linear doublet dist. - Neumann and Dirichiet b.c

ist kind C1 quadratic vortex dist. - Neumann b.c.

Note that ail the absolutely stable formulations are associated with Fredhoim
integral equations of the 2nd kind. This is the main reason for preferring 2nd
over ist kind equations.
 2.107

2.2.8 Factors affecting the magnitude of the truncation error

The actual absolute level of the global truncation error in a particular

application is not only determined by the degree of the representations of the
singularity distributions and geometry representation but also by
- the geometry and flow conditions
- the panel distribution.

A regular panel distribution with gradually changing mesh dimensions that is

adapted to local gradients of source and doublet distributions (and thus to the
surface curvature) usually gives the most accurate solutions.
Fig. 2.2.28 presents the L2 or Euclidean error norm of the tangential velocity
of ist and 2nd order panel methods applied to the flow around a Karman-
Trefftz wing section with curvature-adapted panelling. Fig. 2.2.29 illustrates
that a less favourable panel distribution used with a 2nd order method does not
necessarily yield better results than a ist order method (same flow problem as
in fig. 2.2.28). From both figures it appears that a ist order method requires
80 to 150 panels to achieve an accuracy of 1% in the velocity distribution. This
number is 20 to 80 panels for a 2nd order method. To achieve an accuracy of 0.1%
the number of panels is 200 to 1000 and 80 to 200, respectively. These numbers
illustrate that higher order methods are particularly interesting when a high
accuracy is required.

Figs. 2.2.30 and 2.2.31 give comparisons of results of ist generation, ist order
panel methods (NLR PANEL, BAC/H(unt)-S(emple)) with a high accuracy 3rd order
solution (BAC/Roberts) for a 3-D swept wing. For the thick (15%) wing case of
fig. 2.2.30 the agreement is quite acceptable for both 360 and 730 panels (per
half wing). For the thin (5%) wing (fig. 2.2.31) this is no longer the case for
the lower number of panels.
This is associated with the fact that the representation of the (near-field)
tangential velocity induced at the surface panels by doublet/vortex panels on
the wing camber surface (table 2.1) becomes inconsistent for vanishing thickness
(see also the discussion in section 2.2.6). Methods employing (surface)
vorticity distributions or internal Dirichlet boundary conditions do not suffer
from such behaviour.
Figs. 2.2.32 to 2.2.3L illustrate that particular accuracy problems arise in
case of internal (duct) flows modeled with source distributions. Part of the
mechanism is sketched in fig. 2.2.32 and is associated with the fact that the
 2.108

local truncation error in the normal velocity induced by flat source panels on
ducts with circular cross-sections is equal to /2u, where O is the subtended
angle of the panel. It is further understood that the global truncation error is
also affected by the fact that modelling with source distributions satisfying
discrete Neumann-type boundary conditions is not (numerically) mass-conserving.
This is due to the fact that sources produce mass and that the sum over all
sources is not necessarily exactly equal to zero. For external flows this is of
relatively little concern, because the total mass flow external to the body is
unbounded. For internal flows, in which the total mass flow is bounded, the
error may become relatively large, in particular for channels with rapidly
varying cross-section where the source strengths are locally large.
Numerical models utilising doublet distributions do not suffer from this effect,
due to the fact that doublets have zero net mass production.
Figs. 2.2.33 (flow-through nacelle or ring-wing) and 2.2.34 (duct with
contraction) illustrate the seriousness of the problem for ist order source-
based methods.

2.2.9 Solution methods for the resulting system of linear a1ebraic equations;
computational effort

We recall that after discretisation we are left with the problem of solving a
large (0(103)) system of linear algebraic equations for the unknown source or
doublet strengths. The system may be written as
A = b

where A is the (full) matrix of influence coefficients and the vector of

unknowns.
Two types of solution methods for such systems of equations may be
distinguished:
- direct methods
- iterative methods
In panel methods both types of solution procedures are encountered.
The two types will be briefly reviewed, details are to be found in textbooks on
numerical analysis (e.g. ref. 2.24).
 2.109

Direct methods

Gauss-elimination
In this method the original system A=6 is reduced to a more convenient one
Ux=b', that can be solved readily, by successive elimination of unknowns from
the equations.
Symbolically:

A =6 = 6' (2.2.127)

In this process the matrix A is reduced to the 'upper' matrix U by multiplying

and subtracting successive equations so that zeros appear in the lower triangle:

A stepl step2 U

/xxxx xxxx\ xxxx\

oIxxx
I
oxxx oxxx
9
oixxx oo;xx ooxx
\oxxx oo1xx/ 000x/
t sub matrix t diagonal: 'pivot points'

In this elimination process the righthand side also changes. The system Ux= 6'
is readily solved by proceeding from the bottom to the top of the system. A
necessary and sufficient condition for this method to work is that the pivot
points (elements on the diagonal of the matrices) are 0. A matrix that does
not satisfy this condition is called ill-conditioned.

Gauss-elimination with pivoting is essentially the saine as the basic Gauss-

elimination procedure with the addition that rows and/or columms are inter-
changed in such a way that the element in the first row of the submatrix (see
the scheme above) is the largest.
The advantage is that the method becomes more robust arid gives better results
for ill-conditioned matrices. However, the computational effort increases with
about 50% as compared to the basic Gauss-elimination method.

L-U decomposition (factorization)

Here the starting procedure is again the Gauss-elimination method
 2.110

A = . U =

Now, if a matrix L is defined as

fi o o
(2.2.128)
L= m21 i O

m31 m32 i O
m,43 i
\ml m142

a
where the m.. = are the multiplyers from the Gauss-elimination process
.11 ai*1)
define a vector
(the superscript denotes the step in the procedure), and we also

s such that

(2.2. 129)
Ls =

then we have

(2.2.130)
= L' = L'

hence
(2.2.131)
U =

and
(2.2.132)
= U1s = U'L1L

and therefore

(2.2.133)
(LU)x = b = Ax

(2.2.131-i)
or A = LU

and from
Once L is known, may be determined readily from eq. (2.2.129/130)

eq. (2.2.131/132).
required for
The LU decomposition procedure has advantages if solutions are
(several angles of attack,
several different sets of righthand side vectors
(once and for all)
several normal velocity distributions); with U and L known
* vector (L) multiplication for the
can be determined through simple matrix
various b.
 2.111

A disadvantage is the larger amount of storage required because of the necessity

that both matrices U and L have to be stored.
The method is frequently used, for example in: BOEING TEA 230, Douglas-Neumann
Lifting, PANAIR, MCAERO, QUADPAN, VSAERO (optional).

Computational effort for direct methods

The operations count (the number of multiplications) is as follows:
n3-n
determination of U,L: ops.
3

n(n-1)
s=L -1-
-
b: in
*
2
ops.

- -1- * n(n+1) OÇ)S.

x=U s: in
2

Total number op ops.: + mn2 - (2.2.135)

where m is the number of righthand side vectors for which the computation has to
be repeated.

Iterative methods

With iterative methods the matrix is not decomposed or factored (e.g. A = LU)
but it is split up, in some suitable way, as

A = N + p (2.2.136)

The system of linear algebraic equations

Ax = b

may then be written as

Nx = - (2 .2. 137)

This is solved by means of an iterative process symbolized by

= (2 .2. 138)

Here the superscript denotes the iteration number; is suitably chosen

starting solution (e.g. the null-vector).
 2.112

As an alternative, eg. (2.2.138) may be replaced by

(_1)j
N{x' = b - (2 .2. 139)

or

(v-1)
NC'

The method symbolized by eq. (2.2.139) is called the residual-correction

-(u-l) -(y) = b- - Ax-(v-l) the
formulation; C-(u) = x-(u)
.

- x is the correction and R

residual.
The matrix N, the iteration matrix, is chosen such that eq. (2.2.138) or eq.
(2.2.139) can be solved relatively simple (can be 'inverted' relatively simple).
That is to say N should be a diagonal matrix or a upper/lower matrix.
The requirement that the iterative proces for eqs. (2.2.138) or (2.2.139)
converges implies that the original matrix A must satisfy certain conditions
regarding its structure and also that there are certain restrictions with
respect to the choice of the matrix N.

The conditions for convergence can be analyzed by (re)formulating the iteration

process as follows.
Subtracting eq. (2.2.137) from eq. (2.2.138) gives

((l)
N(- ) = - p )
(2.2.140)

or

N'=
where is the error of the vth iterate.
Multiplying with N gives

(v-l)
- N1

= -
(2.2.1141)

The matrix M = N P is called the 'error' matrix.

 2.113

It can be shown (ref. 2.2k) that the iterative proces converges if the matrix M
is 'convergent', or

hm M' = O
'J4

The latter is the case if the spectral radius p(M), of M satisfies

p(M) max H11 < 1

where A. are the eigenvalues of M. Since p(M) < 11Ml I this may also be stated as

(2.2.142)
lIMlI <1

The convergence rate of the iteration process is defined by

I II (2.2.143)
(u-1)
ll

Considering eq. (2.2.141) and taking norms one finds

llIl = tlM'll IIMII 1l1)I1

or

I
lll
-1)
I
IIMH

Hence,

(2.2.145)

and convergence is assured if

p(M)IIMII < 1

In practice the determination of the spectral radius is not an easy task.

However, it follows from eq. (2.2.141), that the possibility of convergence
increases with increasing I IN and decreasing I II I . This requires that the
matrix N contains a sufficient number of large elements of the matrix A.
In the following some examples are given of iterative methods used for
panel methods.
 2.114

Jacobi-iteration
In the Jacobi-iteration, also called simultaneous iteration, the N matrix of the
preceding method is the diagonal matrix, such that

A = N + P

xxxx\ /x000\ /oxxx

xxxx1= 0x00 + x0xx
xxxxJ 00x0 xxox
xxxxl \000x/ XXO

Gauss-Seidel iteration
The Gauss-Seidel iteration procedure (also called successive iteration) uses as
N matrix the lower triangle of the A matrix, thus

A N P

xxxx\ x000' o X x x

xxxx = xxoo + ooxx

xxxx xxxo 000x
XXXXJ XXXXJ 000 0/
Both Jacobi-iteration as well as Gauss-Seidel iteration lead to a simple and
fast solution procedure for eq. (2.2.138) or (2.2.139).

Block-iteration
Panel methods for 3D flow problems generally result in large matrices, which,
possibly, do not fit in the central memory of the computer. In that case the
matrix is organized in blocks (fig. 2.2.35). The Jacobi-iteration or the Gauss-
Seidel iteration is then applied per block while the blocks on the diagonal are
inverted by means of Gaussian elimination. The blocks are chosen in such a way
that they represent the influence of strips or segments of panels. The blocks on
the diagonal represent the influence of a strip or segment on itself.
It is also possible to choose N such that in addition to the diagonal blocks it
contains blocks adjacent to the diagonal (block-tn-diagonal matrix).
Advantages and disadvantages of iterative methods
The main advantage of iterative methods is in the computational effort. The
operations count is given by k*n2 operations, where k is the number of
iterations and n is the number of panels. This figure must be compared to 0(n3)
for direct methods. For large n and sufficiently small k, the iterative methods
are superior to the direct methods.
However, in order to have a fast convergence (small number of iterations) it is
necessary that the diagonal (blocks) of N and A contain elements that are
'sufficiently large'. This is generally so for Fredhoirn integral equations of
the 2nd kind and generally not so for Fredholm integral equations of the ist
kind. The larger the blocks the greater the probability of convergence (in this
respect it may be noted that there even exist converging block-iterative methods
for ist kind Fredholm integral equations!)
The precise structure of the matrix and, hence, the rate of convergence is
determined by the geometry of the problem and by the way panels are organised in
strips, (which is reflected in the block structure). The Gauss-Seidel procedure,
in this respect, is somewhat more forgiving than the Jaconi procedure, although
both methods formally have identical conditions for convergence.
Panel methods based on Fredhoim integral equations of the 2nd kind require, for
not too complex geometries, about 10 to 20 iterations. In such a case the
iterative methods are faster than direct methods if the number of panels exceeds
200 - 300, as is shown in fig. 2.2.36. However, iterative methods are sensitive
to 'unfavourable' paneling, causing large elements outside the diagonal blocks
that have a negative effect on the rate of convergence. For complex
configurations this may be unavoidable. In other words iterative methods are
less robust than direct methods.
Examples of panel methods using iterative solution procedures are
DOUGLAS-NEUMANN (Hess & Smith, ref. 2.10), Gauss-Seidel
NLR Panel, block Gauss-Seidel
MBB Panel
BAC (Hunt & Semple, ref. 2.15)
VSAERO (optional)
DOUGLAS-NEUMANN Lifting 2nd order (Hess & Friedman, ref. 2.17)

Computing time and storage requirements

The computing time and storage requirements of panel methods will be shown for a
very considerably from one method to the other as a typical example. We consider
 2.116

a wing-body configuration with about 730 panels. The computation was performed
using the NLR panel method with an iterative solver on a Cyber 180-855 computer
(comparable with the IBM-3080 series):

CPU time: 165 sec (corresponding to n2 ops)

of which 80% for determination of' the aerodynamic influence coefficients,

15% for the solution of the system of equations,

5% for postprocessing purposes.

Central Memory: about los numbers

1(3*n2 velocity components, 1*na
Background Memory: 4*n2 2.1*106;
normal velocity components)
The last figure implies that I/O operations will take a considerable time
(several times CPU-time!). However, this is not necessarily the case when a
large(r) central memory is available.
A further comparison of the computational requirements of several methods is
shown in fig. 2.2.37.

2.2.10 Postprocessing; lift and drag

We have discussed in section 2.2.1 that, once the solution (in terms of the
source strength distribution o and/or the doublet strength distribution p) of
the discretized integral equation is known, we are able to compute the velocity
arid pressure everywhere in the flow field and thus also in the panel centres.
Forces and moments on the body may then be determined by numerical integration
of the pressure distribution. This should be done through an integration
rule with a truncation error of the same (or higher) order as the truncation
error of the calculated velocity or pressure distribution. Usually, the simple
mid-point rule suffices. For example one may use for the force in x-direction

NB
K
X
= - f! pn
S5
X
dS = - I
i=1
xi i
(pn ). ó.S (+ 0(h2)) (2.2.1146)

Similar expressions hold for K and K

y z
For the lift coefficient one can write

K NB
z 1 (2.2.1147)
1 (Cn) ô.S
CL=pU2S s i=l
i
2 ref ref
 2.117

and for the drag coefficient

K NB
X 1
CD=l (Cn) ó.S
i
(2.2. 148)
PUZ S S i=1
2 ref ref

In eqs. (2.2.147) and (2.2.148) Sref is a reference area, e.g. the wing area,
and C the local pressure coefficient. CL and CD may also be determined for a
(wing-) strip.
The pitching moment coefficient taken with respect to a given point (x,z) can
be written as

-J! {pn (x-x) +pn (z -z)] dS

z o x o
CM=
pU2S ref Q
ref
2

N
B

pzi (xO-x i
1
[ (C n J J
+ (C n
pxi )
(z -z J )o S
o i
(2.2. 149)
S Q i
ref refi=1

where is some reference length, such as the wing mean chord.

ref

The general experience is that the determination of CL and CM in the way

described above is sufficiently accurate even for moderate numbers of panels.
This is not so for the determination of CD. In the first place the accuracy
requirements put on CD are higher than for CL. For example an error of .01 is

usually acceptable for CL, but one would like to have an error of .00O1
(1"count") in CD. Secondly the strong curvature of the surface at the wing nose
results in strong pressure gradients that may produce large errors in CD unless
a very fine panel distribution is used, see the qualitative sketch of fig.
2.2.38.
According to the Kutta-Youkowski theorem the drag must be zero for a closed body
in a 2-D potential flow. Unfortunately, 2-D panelmethods do not give zero drag,
the error depends on the number distribution and of panels and on the angle of
attack; the error may be positive or negative. First order panelmethods using 60
to 80 panels around a wing profile may show errors ranging from a few counts to
lo counts depending on the airfoil shape and the angle of attack. For a typical
value of CD = 0.01 for a wing profile with turbulent boundary layer flow this
implies an error of few percent to 10 percent or more. For higher order methods
the errors are usually smaller, in general less than 5 counts. The error tends,
of course, to zero when the panel dimension h approaches zero.
 2.118

In 3-D potential flow the drag that is modeled by the panel methods is the
induced drag, in other words the drag due to the downwash of the vortex wake.
Hence, we may conclude that in 3-D potential flow the integrated pressure drag
is equal to the induced drag.
A better way of evaluating the (induced) drag (it may also be done for the lift)
is obtained by applying the momentum conservation law to a closed controle
volume the bounding surface of which is at large distance from the body, see
fig. 2.2.39.
For the component of momentum in x-direction we then find

' (p.) u] dS = 0 (2.2.150)

+
S

where S = SB + Sc + S + ST + S (see fig. 2.2.39), is the local velocity

vector and U the x-component of . From eq. (2.2.150) we find that the
(pressure) drag D can be expressed as

D = - 1f pn dS f! [pn + [p]U] dS (2.2. 151)

P
SB S_SB

where we have made use of the fact that

1f (pi)U dS = 0.

(Neumann boundary condition on

Eq. (2.2.151) may be simplified by letting S+; in that case p+0 and the
contribution of S vanishes. The contribution of the integral over S (the

vortex sheet) also vanishes if conservation of momentum is preserved over S.

For a 'frozen' shape of the vortex sheet, which is usually the case in
applications of panel methods, the S0-integral does not vanish necessarily,
see section 2.1.3. However, the experience is that the error is small if the
contribution of Sc is neglected. Thus we are left with the integrals over S -
(upstream infinity) and over ST (the 'Treff tz-plane').
These result in

D = 1f [p-p + (U-U)pU] dS (2.2. 152)

P
S+S
 2.119

It can be shown (see e.g. ref. 2.53) that the surface integral (2.2.152) can be
reduced to the line integral

D = - f
P
SQS,

or, since [ = r (circulation), see eq. 2.1.26,

D = - f r dt = D. (2.2.153)
P
SS i

where t is the direction along the intersection SCIÌST and the normal on S,
see fig. 2.2.39.
The integral (2.2.153) can be approximated numerically and the experience is
that evaluation of the (induced) drag along this way is more accurate than the
pressure integration procedure discussed previously. Fig. 2.2.40 provides an
impression of the differences in results for the two procedures for a simple 3-D
wing (the same wing as that of fig. 2.2.31) as computed with the NLR panel
method for 60 x 12 panels.
The Trefftz-plane procedure may also be applied to complex non-planar
configurations consisting of several lifting surfaces.

2.3 The Modeling of compressibility (subsonic)

As mentioned already in section 2.1.1 panel methods may also be used to obtain
approximate solutions of the linearized potential (Prandtl-Glauert) equation
governing compressible flow problems. It was also indicated that, because of the
fact that the Prandtl-Glauert equation can be transformed to Laplace's equation
by a simple coordinate transformation, everything said so far on panel methods
for Laplace's equation is also applicable to panel methods for the Prandtl-
Glauert equation. (In fact most panel methods solve the Prandtl-Glauert rather
than Laplace's equation.)
However, the precise way in which various panel methods handle (linear)
compressibility effects varies from one method to the other. It is useful,
therefor to discuss this matter in some detail.
 2.120

2.3.1. The classical Göthert rule

We begin with recalling that the Prandtl-Glauert can be written as

,2O
xx
+4> +4>
zz
=0
yy
or (2.3.1)

._(,24>J
ax X ay y
+ az-
a
(4>
zj
= 0

where 4> is the perturbation potential defined by eq. (2.1.4) and =

Using the transformation

x=x
(2.3.2)
=

we obtain

+ <p + = O (2.3.3)
xx yy zz

The boundary condition on the body is defined by

(2.3.4)
(Ux + = O

or

xx +4>n +4>n
zz =-Un
4>n (2.3.5)
yy
the
Recall also that eq. (2.3.1) stems from the mass conservation law under
assumption of irrotational flow and small perturbations, i.e.

4> ,4> 4> = 0(c), c << 1 (2.3.6)

X y z

Since 4=i, it folows from eq. (2.3.5) that there must hold

(2.3.7)
n = 0(c); n ,n
z
= 0(1)
y

Hence eq. (3.2.5) can be reduced to

zz =-Un 0(c2)
(2.3.8)
4>n +4>
yy
or
 2i21

+ n = - Un
oex
+ 0(c2) (2.3.9)
yy z

Since the (unit) normal vector transforms as

- - l( -n
1 1
-n (2.3.10)
n (n n- n ) = - nx
- F y z
X y z

where

F2 = ,2n2 + n2 + n2 (2.3.11)
X y z

we have
n = Fn
X -X
n = Fn (2.3.12)
y
= Fn

and eq. (2.3.9) may be written as

+ = - U= + 0(c) (2.3.13)

Introducing the transformed potential

= (2.3.14)

we find for eq. (2.3.13)

+ + (2.3.15)
= - Uco
yy
- -
zz X

where the derivatives are expressed as

X
X
= ct, (2.3.16)
y

z
=-
13
 2. 122

Eqs. (2.3.2), (2.3.114) to (2.3.16) represent the well known Göthert rule from
classical small perturbation theory. The Göthert rule enables us to obtain an
approximate solution of the compressible flow around a configuration from the
solution of an incompressible flow past a similar or 'analogous' configuration
that is determined by the transformation (2.3.2) with the boundary conditions
(2.3.15) (the term 0(c2) can be neglected).

2.3.2 Application of the Göthert rule in panel methods

For panel methods which utilize the 'exact' boundary condition (2.3.5) for the
solution of the incompressible flow problem, the question may be raised what to
do with the term n (= 0(c2)J in the compressible flow case when the Göthert
rule is applied.
One of the possibilities is to use, instead of eq. (2.3.15), the complete
boundary condition of the analogous incompressible flow, i.e.

+ + = - uÇ (2. 3 . 17)

xx yy zz X

Transforming back we then get

2n
xx +n
yy +nzz=-Un (2.3.18)

which is not the same as the real boundary condition eq. (2.3.5).
However, for small perturbations the error is only of 0(c2), thus, formally, in
the asymptotic sense (c-*0) the error is of the same order as that of the
linearized boundary condition given in eq. (2.3.8) or eq. (2.3.15). More
importantly, eq. (2.3.18) unlike (2.3.15) reduces to the exact boundary
condition (2.3.5) for M+0 (-l) . This suggests a preference for eq.
(2.3.17/18).
At positions where locally nx 0(1), such as at the wing nose, the small
perturbation theory on which the partial differential equation (2.3.1) is based
is violated anyway, and neither (2.3.17/18) nor (2.3.7/8) provides a valid
approximation for compressible flow.

The Göthert rule given in the form of eq. (2.3.2), (2.3.16), (2.3.17) is
sometimes referred to as Göthert rule I; it is applied for example in the NLR
panel method.
 2.123

A second procedure encountered in the literature (in particular in those methods

where Neumann boundary conditions are applied) follows the following steps.
1. First the perturbation potential and the perturbation velocity of the trans-
formed flow problem are formulated in terms of the (as yet unknown) source and
doublet distributions on the analogous configuration,

(P)
1
= - fi
ku
1
oH -) dS +
a
(-1)d (2.3.19)
r anQ r

(P) ku=ifa (-è) d (2.3.20)

r r

where the integrals are obtained from chapter 2.1, eq. (2.1.107) and its
gradient; the tildas denotes the transformed quantities according to eq.
(2.3.2).

Next the Göthert rule (eq. (2.3.16)) is applied to the components of the
perturbation velocity of the transformed problem. In vector notation this may be
written as

(2.3.21)

where
loo
[B] = O O (2.3.22)
oo
Finally the 'real' boundary conditions are imposed according to eq. (2.3.5),
thus

= - U n (2.3.23a)

or, with eqs. (2.3.10) - (2.3.12)

1---
-n --
+ n -- + n = - Un (2.3.23b)
xx yy zz X

The procedure just described is applied in the panel methods of MBB (Kraus &
Sacher, ref. 2.13) and of BAC (Hunt & Semple, ref. 2.15). It is known as the
Göthert rule II.
 2.124

The method has the formal disadvantage that the integral equation which is
eventually solved is not derived from the application of Green's theorem to the
governing partial differential equation, i.e. the Prandtl-Glauert equation. This
implies that the uniqueness of the solution cannot be proven in a formal way,
see also ref. 2.30; in practice this does not appear to be a problem.

From a viewpoint of thin-airfoil theory the Göthert rules I and II are

equivalent. For configurations with thickness their results differ, in
particular in the vicinity of blunt noses and leading edges. Since the Göthert
rule is based on the assumption of small perturbations, neither rule I nor rule
II gives valid results in the vicinity of blunt noses. However, Göthert rule II
appears to give a smaller absolute errors in such regions than Göthert rule I.

2.3.3 Alternative treatments, mass flux boundary condition

A third possibility for solving the compressible flow problem with panel methods
is based on a 'direct distribution' on the surface of a configuration with
elementary solutions of the Prandtl-Glauert equation in psysical space. In other
words, we assume

(P) = j- Jf oH )
dS + fi p 3flQ
1
)
dS (2.3.24)
s r 11s r
or
fi o(- dS + j-- if H dS (2.3.25)
P) = )
i --p-- )

S r 11S anQ r

where f ( ) as well as - (- are elementary doublet-like solutions of the

Prandtl-Glauert equation. From the integral representations, eq. (2.3.24) or eq.

(2.3.25), integral equations are deduced which satisfy the boundary condition of

eq. (2.3.5).
Like in the case of the Göthert II method also the integral representations
given by eqs. (2.3.24) and (2.3.25) are not obtained from the application of
Green's theorem to the Prandtl-Glauert equation. If we do this (see ref. 2.30)
we find

(P) = L ji
-
dS + L j QV3 C- dS (2.3.27)
4rT
S r S r

where
-
= -Bx + -ay+ -az=
a a a 2
[sW'] ([s])
 2.125

with the matrix [B] given by eq. (2.3.22). Then eq. (2.3.27) can be written as

= f! o[- )
dS + ff pQ.[B] ([B']) (- )
dS (2.3.28)
S r S r

Curiously enough, as far as known to the author, no use is made in practice of

this formal integral representation.

Instead of the velocity boundary condition eq. (2.3.5) a so-called mass-flux

boundary condition is applied in sorne methods. The latter condition is obtained
by setting to zero at the surface of the configuration the normal component of
the approximated mass-flux vector. According to eq. (2.1.12) we have

3- ( 2-Y ) M2
pe[u M22 M2 (2 + 2)
x = + (1-M2e)
x
-
2 ex - 2 e x z
+

y = e [ y - M2e xy + (2.3.29)

= pe [i - M2 +
z z xz

Then, the lowest order approximation to the normal component of the mass flux
vector becomes

((U + l32)i y: (2.3.30)

+

which leads to the boundary condition

n
xx + n
yy + n
zz ex
= - Un (2.3.31)
or
j2 [B] ([s]V) = - U n ex (2.3.32)

Note that this is the same as the boundary condition (2.3.18) of the Göthert I
model after back transformation to the physical space.
The mass-flux boundary condition is applied in the panel methods PANAIR and
QUADPAN.

2.3.4 Summary of linearized models for subsonic compressible flows

The integral representations and the imposed boundary conditions of all subsonic
compressible flow models mentioned above may be expressed in terms of quantities
 2.126

in physical space by using the relations (transformed quantities denoted by

tilda (-))

= [B] i (2.3.33)

[B] (2.3.34)
I[B'] $

= 2
I[B'} :;i dS (2.3.35)

= [B] (2.3.36)

a
=
-- [B
-1
] n
-
[B1] (2.3.37)
Bn [B ] n

where
i o o

[B1] = o o

o o
i

We then obtain:

Göthert I

(P) =
1 1ff(l)d+lff (-è) d]
r anQ r
s s

= L!
1
- If o
s
r
) [] I
dS + ff

s
[B1J .
{B1] VQ (-
r
)
dS

(2.3.38)

boundary condition

2
[B'] [81] = - U n (2.3.39)
 2.127

Göthert II
Integral representation as Göthert I;
boundary condition

= - U n
mx (2.3.40)

Direct I

'(P) =ff o (iJ +tffpQQ (1) dS (2.3.41)

4iï
r r
s s

boundary condition
'exact': eq. (2.3.23) or eq. (2.3.40)
or 'mass-flux': eq. (2.3.31) or (2.3.39)

Direct II

(P) = f! o -
1) dS + ([B] (1) dS
r I
r
s s
(2.3.42)

boundary condition as for 'direct I'.

Formal integral representation for the Prandtl-Glauert equation

1)
(P) = f! o - dS + lip [B] ([B]Q) (- )
dS
r r
s s
(2.3.43)

boundary condition as for 'direct I'.

The five models are equivalent in terms of subsonic thin airfoil theory. All
models provide (numerical approximations of the) solutions of the Prandtl-
Glauert equation.
For thick bodies those formulations are equivalent which work with the same
boundary condition. The elementary solutions forming the surface distributions
are different, but they become identical for M-0 as is also the case for the
velocity boundary condition and the mass-flux boundary condition.
 2.128

The expressions, in terms of a and p, for the jump of the potential and the
normal velocity component across the surface S are different for MO and
In particular the doublet-like distributions of eqs. (2.3.38/41/42/43) do not
only exhibit a jump in the potential but now also in the normal velocity. This
is due to the fact that the kernelfunction in eqs. (2.3.38/41/42/43) differ from
those in the incompressible case. This can be interpreted as that in the
compressible case the doublet axis does not coincide any more with the surface
normal. Note however, that the jump of the normal component of the linearized
mass-flux is zero for the doublet-like distribution of eq. (2.3.43). Therefore,
this integral representation is particularly attractive in combination with
mass-flux boundary conditions. The jump behaviour in the case of a source
distribution is, qualitatively, the saine as in the incompressible case.
The difference in jump behaviour has also consequences for the equivalent
internal Dirichlet boundary conditions.
Evidently all models fail in the vicinity of blunt leading edges. Moreover, also
in regions where the approximations are valid (nx small), the error for a given
geometry increases with increasing Mach number, which is due to neglecting the
higher order terms in eq. (2.3.29). However, the thinner or more slender the
configuration the higher the Mach number can be for reasonable results.

For some panel methods (NLR, BAC) semi-emperical correction methods for non-
linear compressibility effects have been developed, which decrease the errors at
blunt noses and at high subsonic Mach numbers. Such correction methods are based
on analytical solutions of higher order approximations of the potential equation
for the compressible flow past simple geometries (ellipses etc.), see ref. 2.31.
Fig. 2.3.1 shows some examples of the application of such a correction method.
(NLR-panel method, as applied to airfoils)

We finally note a problem that arises in the modeling of compressibility effects

for engine inlet flows (fig. 2.1.24). If the prescribed values of the normal
velocity component or the mass flux at the compressor fan face differs much from
U or pU (which is usually the case at take-off) the condition of small
perturbation with respect to the free stream is violated. In such a case the
internal flow problem may have to be treated separately in terms of
perturbations with respect to the conditions at the compressor fan face.
 2.129

2.4 Examples of applications of 'ordinary' 3-D panel methods

2.4.1 Status and role of panel methods in the process of aircraft aerodynamic
design

Panel methods are utilised in several phases of the process of aircraft

aerodynamic design. Fig. (2.4.1) depicts the major phases that can be
distinguished.

In the preliminary (or even conceptual) design phase (phase I in fig. 2.4.1)
panel methods are used at an increasing rate for the purpose of obtaining rapid
estimates of the forces and load distributions on (complete) aircraft
configurations. As such there is a tendency to replace (part of) the functions
of empirical data base type methods such as the well-known USAF 'DATCOM' and
ESDU 'Data Sheets'.
In phase II (Design Development) one is concerned with the detailed shaping of
(parts of) aircraft configurations. Fig. 2.4.2 depicts the design development
process in some detail.
In the process (A) detailed shapes are generated, which are analysed in the
processes (B) and (C) . Panel methods are used in process (A) as well as in
process (B).

In phase II (B) panel methods serve to analyse details of the flow pattern, load
and pressure distributions on parts of or complete configurations. In phase II
(A) they serve to solve the 'inverse' problem: given the pressure distribution,
what is the shape that produces that pressure distribution?

2.4.2 Low speed flow

Simple wing-body configuration (ref. 2.31)

An example of a simple wing-body configuration is shown in fig. 2.4.3. The

results of computations using the NLR panel method are compared with those of
experiments in figs. 2.4.4 - 2.4.9. Overall forces and moments are compared in
figs. 2.4.4 - 2.4.6, load and pressure distributions in figs. 2.4.7 - 2.4.9.
In considering these comparisons it should be realised that the main reasons for
discrepancies between theory and experiment are
 2.130

viscous (boundary layer) effects;

the assumption of a planar vortex sheet and the approximate Kutta condition
in the theoretical model (this effect is expected to be small in general);
numerical errors.

Note that:
- the viscous effects, at the relatively low Reynolds number of 0.65 x 10' cause
a decrease of the lift of about 20% and an increase of the moment.
- the increase of the viscous drag with the lift should be positive,
illustrating that the integrated pressure drag in fig. 2.3.6 is not correct.
- viscous effects increase rapidly when flow separation occurs (CL 0.65).

Wing-alone and wing-fuselage configuration of a small transport aircraft kref.

2.31)

Results of computations are shown in figs. 2.4.10 and 2.4.11. In the panel
method used, some 350 panels are applied on fuselage and fairing, and about 700
panels on the wing; for symmetry reasons only one half of the configuration is
considered.
It can be observed from the results that the effect of the presence of the
fuselage on the wing load and pressure distribution is considerable. Again the
main reason for the difference between theoretical and experimental results is
found in viscous effects.

Wing-fuselage-tailplane configuration of a small transport aircraft(ref. 2.31)

Ssee figs. 2.4.12 and 2.4.13. The following observations may be made.
Variation of the position of the vortex sheet from the wing does not notably
effect the solution, as may be observed from the lack of difference between the
computed results for an undisplaced vortex sheet (='wake') ('initial ist wake
estimate', i.e. straight downstream) and a displaced (corrected) vortex sheet,
see fig. 2.4.12b and 2.4.13a,b,d. The corrected vortex sheet was obtained by
integrating once the normal velocities on the vortex sheet to obtain a stream
surface.
However, the results reveal clearly the effect of a sting support, especially on

the stabilizer. This illustrates that one should be careful in the

interpretation of wind tunnel measurements in particular with respect to the
sting-induced flow angle at the horizontal tailplane.
 2.131

The remaining causes of the discrepancies between theory and experiment are
thought to be:
- viscous effects on the various aircraft parts (in particular separation on
the rear part of the fuselage);
- absence of the modeling of the vortex sheet (wake) leaving the rear part of
the fuselage.

Fuselage-wing-flap (research model) (ref. 2.28)

The configuration and paneling (1219 panels) are shown in fig. 2.4.14; the
schematized vortex sheet models and position of vortex sheet are depicted in
fig. 2.4.15. The results are shown in figs. 2.4.16 to 2.4.19.
Note that in this case the (correct) positioning of the vortex sheet of the main
wing relative to the flap has a significant effect. Note also that the absence
of the modeling of' flow separation at the side edges of wing tips and flap tips
is clearly reflected in the local pressure distributions (figs. 2.4.18/19a).

Fuselage-wing-flap-stubwing-nacelle of a small transport aircraft (Ref. 2.29)

See fig. 2.4.20. This example illustrates the possibility to model the effect of
nacelle and engine control on the lift distribution on wing and flaps. No
experimental results are available.

2.4.3 High speed, subsonic flow

Examples of application (NLR Panel Method) to high speed subsonic flow are given
in figs. 2.4.21 through 2.4.23 (from ref. 2.29). Fig. 2.4.21 provides a
comparison of calculated and measured pressure distributions for a wing-tiptank-
pylon(-stare) configuration. Note that for this, admittedly thin and slender
configuration the agreement is quite good, in spite of' the high Mach number
(0.8). Note also that the effect of the presence of the store on the wing
pressure and span load distribution is quite large and well predicted.
Fig. 2.4.22 gives results for a rear-mounted nacelle plus stabwing
configuration, for a civil transport at cruise condition M= .745). Purpose of
this application was to explore the risk of local supercritical flow and shock-
wave formation in the duet-like area between nacelle and rear-fuselage. As
indicated by both calculated and experimental results (which are in fairly good
agreement) there is a small area of local supersonic flow near the leading-edge
of the stub-wing.
 2.132

Fig. 2.4.23 gives results for a wing-fuselage configuration with underwing

pylon-nacelle at M= 0.75. Experimental as well as computational results have
been obtained for a so-called blown-nacelle configuration. With the blown-
nacelle technique the effect of the jet is simulated in the windtunnel by means
of high-pressure air that is fed into the nacelle through an auxiliary pylon-
type part. For this purpose the nacelle inlet must be faired-over. In this case
the agreement between measured and calculated pressures on the wing lower
surface and pylon is less good. Insufficient resolution of the complicated
geometry in the nacelle-pylon-wing intersection region is believed to be one of
the reasons for this (ref. 2.29). An interesting aspect of this kind of
application is the possibility to compare computational results for the real
(flight) configuration with those of the windtunnel configurations. The
difference between the two could be used to correct the windtunnel results for
the effect of the air-lead pylon and the inlet fairing.

2.5 Simulation of viscous effects

In the preceding chapter we have observed that viscous effects can lead to
significant differences between potential flow theory (panel methods) and
experiment. From classical boundary layer theory (Prandtl boundary layer
approximations) we know that for high Reynolds numbers and attached flow
viscous effects are limited to the boundary layer and its extension (wake)
behind the airfoil. In such conditions, viscosity has a relatively small
effect on the forces and pressure distribution acting on the body.
According to the Prandtl's boundary layer theory, the high Reynolds number flow
about a body can be modelled, approximately, by distinguishing two 'zones': an
iriviscid part governed by the potential eq. (or Euler eqs.) and a viscous part,
governed by the boundary layer equations. We call this zonal modeling, see also
section 1.3.2.
That such a subdivision is allowed follows also from more formal consideration
based on the theory of matched asymptotic expansions, ref. 2.33 of
characteristics of the Navier-Stokes equations at high Reynolds numbers.
In the theory of matched asymptotic expansions solutions of the Navier-Stokes
equations for high Reynolds number are expressed in terms of series expansions
with the Reynolds number (Re = UL/') as expansion parameter and the airfoil
chord (or body length) L and boundary layer thickness ô (- Re )
as

characteristic lengths.
 2.133

If x, y, z is a body-fitted curci-linear coordinate system (x, y tangential; z

normal), the solution away from the body or airfoil can be expanded as ('outer'
expansion for the velocity vector)

- 1) 1/2 (i)
u (x,y,z; Re) = u0 (x,y,z) + Re u1 (x,y,z) + ...; Re (2.5.1)

with similar expressions for the pressure and density. The leading (that is
first) term of (2.5.1) must satisfy the inviscid flow equations which are
obtained by substituting (2.5.1) into the Navier-Stokes equations and retaining
only the lowest order ternis. The boundary condition on the body is
(i)
u0(x,y,O) .n=0.

In the boundary layer the solution is expanded as ('inner' expansion)

(y) _i/z ('i)

u(x,y,z; Re) = u0 (x,y,z/ó) + Re u1 (x,y,z/ó) + (2.5.2)

again with similar expressions for the pressure and density.

The leading term of (2.5.2) must satisfy the boundary layer equations with, on
(y)
the body, the boundary condition u0 (x,y,O) = 0. Matching of the inner and outer
solutions requires that the outer limit of the inner expansion equals the inner
limit of the outer expansion, i.e.

(y) (i)
hm u0 (x,y,z/ó) = lu u0 (x,y,z) (2.5.3)
z/ó+ z+0

Extending the matching technique to the second order term in the 'outer
expansion' , it is found that this second order term represents the effect of the
displacement thickness of boundary layer and wake. The displacement thickness
influences the velocity and pressure distribution of the inviscid outerflow (the
outer flow has to be adapted to a new shape including the displacement
thickness). In turn, the modified pressure distribution has its influence on the
boundary layer development, etc. In this way an iterative proces may be thought
of, where, alternatingly, the inviscid outerflow and the boundary layer and wake
are computed until the pressure distribution and displacement thickness in
innerfiow and outerflow are matched together. In the block diagram of fig. 2.5.1
this process is illustrated. In the remaining part of this chapter we shall
confine ourselves to the modeling of the displacement effect in the outerflow of
 2.1314

boundary layer and wake. The computation of the boundary layer itself is not
considered.

The displacement effect on the outerfiow may be modeled in various ways, see
ref. 2.3k. Two of the most important methods are discussed below.
1. The displacement thickness of boundary layer and wake, in 2D flow given by

* o
0 f l U
pU
j dz (2.5.14)
O ee

is added to the body and the dividing streamline behind the trailing edge (see
fig. 2.5.2). In the expression for the displacement thickness, O is the boundary
layer or wake thickness; the subscript e refers to the edge of the boundary
layer or wake.
The potential flow is then computed for the thickened body and wake
(displacement surface). The boundary conditions are: zero normal velocity at the
surface of the displacement body and plus the additional condition of zero
pressure difference across the wake (j p = O). Also a modified Kutta condition
must be satisfied at the transition between displacement body and wake at the
trailing edge. The additional condition at the wake causes the problem to be
overdetermined if the position of the (thick) wake is kept fixed (see section
2.1.3); it requires that the wake position is determined iteratively (non-linear
problem).
2. An alternative is to model the displacement effect by prescribing a given
in/out flow at the original profile contour and at the (infinitesemally thin)
vortex sheet (transpiration model).
Through asymptotic theory it can be shown that the displacement effect is
simulated with the same accuracy as in the first method if the in/out flow is
prescribed as

2D:
an
= li
p
-
ds
U
ee
0*) (2.5.5)
e

or,

ee + ô[as (p ee
U
3D: p - () = p W J + (PeVe)] (2.5.6)

Here s,t are orthogonal coordinates on the surface and Ue ,Ve ,W e are velocity
components at the edge of the boundary layer (e) in the directions s,t and
normal to the surface, respectively. The 2D case is sketched in fig. 2.5.3.
 2.135

In linearized flow theory the lefthand side of eq. (2.5.6) may be approximated
by (see section 2.3.3)

-
3n
(p
n )
= (u +
x
)n
x
+ c
n
yy
+ n
zz (2.5.7)

Also in the case of the transpiration model we need an adapted Kutta condition.
The corresponding (jump) conditions on the vortex wake now also require a source
distribution on Sc.

The transpiration model is preferred over the first model in which the
displacement thickness is just added to the actual model and trailing vortex
sheet contours. The reason is that the (modified) boundary condition is imposed
at the original contour and therefore the aerodynamic influence coefficients
remain unaltered during the proces of iteration shown in fig. 2.5.1.

It is important to note that the asymptotic theory underlying the computational

process depicted in fig. 2.5.1 is valid only in case of attached flow. In that
case one speaks of weak (viscous-inviscid) interaction and weak interaction
models. In case of separated flow and at the trailing-edge of wings the theory
is no longer valid, at least not in the form sketched at the beginning of this
(i)
section. Considering eq. (2.5.1) and recalling that the first term u0 of the
outer expansion represents the solution of the inviscid flow equations, with
Kutta condition at the trailing edge, this is not surprising: if boundary layer
separation does not take place at the trailing edge the Kutta condition (at the
(i)
trailing edge) is not appropriate and u0 cannot represent a valid first
approximation. Even with attached flow up to the trailing edge there is a
problem: because of their parabolic nature the boundary layer equations do not
possess a mechanism that accounts for the upstream influence of the
discontinuity (no slip/slip) in the surface boundary condition at the trailing
edge. As a result the calculated displacement thickness at the trailing edge is
inaccurate.
For the reasons just mentioned the iteration process (fig. 2.5.1) of weak
interaction models is not necessarily convergent. Nevertheless, computational
methods have been developed based on the weak interaction model which are quite
useful. In these methods the trailing edge region is necessarily treated in an
emperical-heuristic manner (see for example ref. 2.35-2.37).
Recently considerable progress has been made with more proper modeling of
trailing edge flows and of separated boundary layers. These models belong to the
 2.136

category of so-called strong interaction models; a good insight of the theory

may be found in ref. 2.38-2.40.

Panel methods with simulation of viscous effects are commonly applied for the
analysis of the flow past wing profiles (2D flow) with flaps (airfoils in take-
off and landing configurations), see figs. 2.5.4 and 2.5.5. For 3D wings there,
to the author's knowledge, are only a few methods available in which viscous
effects are taken into account; results from ref. 2.16 are shown in fig. 2.5.6.

2.6 Special methods

In addition to the panel methods for subsonic flow past thick, lifting
configurations of given shape described uptil now there are also methods for
special applications. The most important of these will be discussed briefly in
the following sections.

2.6.1 Methods based on thin-wing theory

In thin-airfoil theory also one looks for solutions of the linearized potential
equation. However, instead of satisfying the full zero mass flux boundary
condition at the surface of the body an approximate boundary condition of the
form of eq. 2.3.8 is used, which is imposed at the chord of the profile (2D) or
at the wing reference plane S(3D), see fig. 2.6.1. In the latter case the
approximated boundary condition reads

at the upperside of S w of S w

+ + + + +

y y
n
z
n =-Un
z
(2.6.1)

and at the lower side, S w

n + n = - U oex
n (2.6.2)
y y z z

+ + +
Here n, n, n z are the components of the normal on the actual wing surface.
x y
Usually the coordinate system is chosen such that the x,y plane coincides with
or is parallel to the wing reference plane S. Then, for thin wings, in addition
to also yfly may be neglected (since then n =0(E)), and eqs. (2.6.1) and
(2.6.2) reduce to
 2.137

+
- +
+ n -
'
-
z
=-U =Uax x
+
aZ
onS w (2.6.3)
n
z

where represents the chordwise slope of the wing profile contour.

The boundary condition (2.6.3) can also be expressed in terms of the jump in
(jump in normal velocity component) across SW and the average of and
We then obtain

aZt
-ii =
ax
= 2U
ax
az - aZ
+
w
= +
z
+
=
U Caz
ax
+ 1
ax
= u ax

aZ
where is the chordwise slope of the thickness distribution of the wing and
ax
az
the slope of the camberline plus the angle of attack.
In order to model the lift we have, again, to introduce a discontinuity surface
extending downstream of S. The boundary conditions on S are the same as
those in the 'thick' wing problem.

From the application of Green's theorem it follows that the solution of the thin
wing problem may be represented in terms of a source plus a doublet distribution
on SW and a doublet distribution on S. The jump in normal velocity over S
given by eq. (2.6.4) can be related directly to the source distribution on S;
its strength follows from

+ aZ(x,y)
G(x,y) = FL(x,y)]J = 2U_ (2.6.6)

Hence, it is not necessary to solve an integral equation to determine the source

strength.
With the source strength known, the perturbation velocities due to the source
distribution are obtained by numerical evaluation of the integral

= J' a dS
(
r ) (2.6.7)
S
w
 2.138

With the source strength on S known directly from eq. (2.6.6) only the doublet
distribution on S + S remains as unknown. This unknown doublet distribution is
w C
determined from eq. (2.6.5) by expressing the mean normal velocity.

:; . = (2.6.8)
+
w 2 z z
as

w p
=
)4
If pFi
w p an -- (1) r
dS + fío w p (
1) dS
r (2.6.9)

wC w

Here the index p refers to the point on S + S where the gradient is taken and
ii is the normal on S + S see fig. 2.6.1. Note that if S is a planar surface
,

w w C w
the second integral (the source distribution) vanishes.
Substitution of eqs. (2.6.8) and (2.6.9) into eq. (2.6.5) results in the
following integral equation for the unknown p

w pan w
a
- r ) dS=U
aZ
W_jj0j. (-)dsr
S+S
wC S
w
(2.6.10)

Unfortunately but also unavoidable this happens to be a Fredholm integral

equation of the ist kind.
The possibilities for discretization and solution of eq. (2.6.10) are the saine
as those of the tordinary' panel methods, with the exception that for the
solution of the system of algebraic equations only direct methods can be
used. The latter is due to the fact that Fredhoim integral equations of the ist
kind do not lead to a sufficiently diagonally dominant matrix.
Examples of discretizations used in practice are;
- constant ti-distribution on each panel (vortex-lattice)
- linear/constant p-distribution (ref. 2.41); constant/line vortex distribution
(unstable!)
- quadratic/linear p-distribution (VSAERO, ref. 2.42)
- quadratic p-distribution (NLRAERO, ref. 2.43).

In practice, methods based on pure thin-wing theory are seldom used. On a large
scale, however, methods are used which combine thin-wing theory for lifting
parts (wings, tailplanes, winglets etc.) with distributions on the actual
surface of fuselage-like parts, see fig. 2.6.2. The integral equation that has
to be solved in that case is of a mixed type (ist kind with regard to the
lifting parts and 2nd kind with regard to the fuselage parts). The number of
 2. 139

unknowns is less than in the case of ordinary panel methods, because on the
lifting parts we have reduced the number of unknowns with one half.

A disadvantage of methods based on thin wing theory is the occurrence of

additional modeling errors, which are caused by the fact that not the real
boundary condition on the actual surface is satisfied but an approximation to it
on an average plane through the lifting surface. As a consequence the lift
computed on the basis of thin wing theory is usually less than the lift computed
on the basis of the real boundary condition; the error increases with increasing
wing thickness. Figs. 2.6.3/14 illustrate the magnitude of the error for a
configuration considered in ref. 2.43. The wing profile thichness is 9% (the

same configuration is considered in figs. 2.4.3 and 2.4.4).

By chance, the error in lift due to neglecting wing thickness is compensated

largely by the error occur due to not taking into account viscous effects. The
thicker the wing, the higher the lift in potential flow (at the saine incidence),
but also the larger the loss of lift due to viscous effects. Both effects are
about proportional to the thickness-chord ratio of the wing. Because of this,
the lift computed on the basis of thin-wing theory is often better in agreement
with experimental results than that obtained by applying the exact boundary
conditions.
Of course, this does not mean that methods based on thin-wing boundary
conditions are more accurate than methods based on the exact boundary
conditions. However, it does mean that thin-wing methods are particularly
suitable for obtaining a first impression of total forces and moments acting on
a configuration. Their attractiveness is found in the smaller computational
effort as compared to 'thick wing' panel methods and also in the more simple
treatment of the geometry (simpler paneling). As a result, the use of thin-wing
methods in the conceptual and preliminary design stage of aircraft is expanding
rapidly.
Early examples of this kind of methods are the Woodward codes (refs. 2.41 and
2.42) and the NLRAERO code (ref. 2.43). Also NPLS, a vortex lattice method
developed from NLR PANEL by Fokker, and PDAERO, a version of AEROPAN developed
by NLR belong to this category. Codes like Boeing TEA 230, PANAIR and QUADPAN
also contain thin-wing options.
 2.140

2.6.2 Panel methods for supersonic flow

In section 2.1.1 we have indicated that the linearized potential equation is

also valid for supersonic flows (M>l), provided only small perturbations of the
flow are involved. We have then

(M2 - i) - - = O , M >1 (2.6.11)

xx yy zz

which is a linear, hyperbolic, partial differential equation for the

perturbation potential 4.
In the discussion (section 2.1.2) of the linearized (elliptic) potential
equation the notion of well-posed boundary value problems was introduced. For
hyperbolic partial differential equations well-posed problems can be formulated,
only for initial value or mixed initial/boundary value problems (see e.q. ref.
2.54).
Obviously the problem of (linearized) supersonic flow about a wing or aircraft
represents a (mixed) initial/boundary value problem. The theory of initial/
boundary value problems for hyperbolic equations is more complex and less
developed than the theory of boundary value problems for elliptic equations.
However, this has not prevented the development of methods for solving the
problem of linearized supersonic flow about aircraft configurations. Without
going into the details of well-posed initial/boundary value problems for the
linearized equation for supersonic flow we note (see also refs. 2.44, 2.45) that
these methods respect the following important consequences of the hyperbolic
nature of eq. (2.6.11):
- disturbances are only propagated downstream, leading to the notion of domain
of dependence and domain of influence, see fig. 2.6.5
- disturbances are propagated undamped along characteristic surfaces (Mach
cones), F(x,y,z) = 0, defined by

3F2 3F2 = Q i
2
(;) - ,2 = M2 - (2.6.12)

- different boundary conditions at infinity.

Also in the case of supersonic flow the linearized potential equation allows the
solution, by means of application of Green's theorem, to be expressed in terms
of surface singularity distributions. In this case the integral relation takes
the form
 2.141

(P) =
1 1
If [o (- -)
r,
+p a

anQ
(- ir-)] dS (2.6.13)
where
r = [(x_xt)2 - 2 (yy,)2 - (z_z')
ji/2 (2.6.14)

- i
=
n = (n n, - n)
z
(2.6.15)

and x>x'

The integration has to be performed only over the part S' of S lying within the
forward Mach cone of P, resulting in a factor in front of the integral
instead of as in the subsonic case.
4ii
The asterix above the integral sign indicates that we have to take the finite
part of the integral (in the sense of Hadamard, see ref. 2.44 or 2.45) , because
on the surface of the forward Mach cone emanating from P the integrand of eq.
(2.6.13) is singular (this in addition to the singular behaviour in x=x' , y=y'
z=z').
Based on eq. (2.6.13) an integral equation may be formulated by satisfying the
boundary conditions on the configuration (normal velocity component zero). The
discretization and solution procedures are, in principal similar to those for
subsonic panel methods. (Note: due to the limited extent of the domains of
influence (see fig. 2.6.5) we no longer have full matrices!).

Almost all supersonic panel methods are based on the application of (supersonic)
thin-wing theory to lifting surfaces in combination with surface distributions
on, body-type components. The computer codes are usually organized in such a way
that subsonic as well as supersonic flows can be treated. Examples are the
earlier mentioned Woodward method (ref. 2.41), USSAFRO (ref. 2.42) and NLRAERO
(ref. 2.43). Also PANAIR (refs. 2.18 and 2.19) contains a supersonic option.
More recently a higher order subsonic/supersonic code has been developed by MBB,
(HISS-code, ref. 2.46).

In a numerical sense supersonic panel methods are less accurate than subsonic
methods. This is due to the property that 'physical' as well as numerical
discontuities are propagated undamped along the surfaces of Mach cones (fig.
2.6.6). It appears to be difficult to represent accurately the 'physical'
discontinuities, while the numerical discontinuities should, of course, be
absent. Because of the latter, higher order methods (ref. 2.46) are more
desirable in the supersonic case than in the subsonic case.
 2. 142

Examples of application of a (low order) supersonic panel method (NLRAERO, ref.

2.43) to a wing-body and a missile configuration are given in figs. 2.6.7 arid
2.6.8.

2.6.3 Methods for modeling vortex sheet roll-up

It was mentioned in chapter 2.1 that the conditions to be satisfied on the

surface(s) of discontinuity (vortex sheet(s)) S of lifting configurations,
are that they should be stream surfaces with zero pressure jump. This implies
that lines of constant trailing vorticity (lines along which ji = flJ = constant)

must be directed along (average) streamlines in Sc. The determination of the

correct shape and position of (the stream surface) Sc in this sense is a
nonlinear problem that must be solved iteratively.
As we have seen in section 2.4.2 it is sufficient for a large number of
applications to treat the vortex sheet in a linearized sense, i.e. to choose
shape and position in such a way that the boundary conditions on S are
satisfied approximately. However, as we have seen also, there are applications
(wing with flap, wing-fuselage-tailplanes) where a more accurate determination
of the shape and position of the vortex sheet is desirable or even necessary.
Examples of aerodynamic problems for which a correct determination of the shape
and position of the vortex sheet is essential are the following.

Slender wings, with separated flow from the side edges which influences
relatively large regions of the wing, fig. 2.6.9 (see also figs. 2.4.14-19).

Wings with large sweep kdelta wings with flow separation from the (sharp)
leading edges (fig. 2.9.10).

configurations with canards, where the trailing vortices from the canard
interfere strongly with the wing flow (fig. 2.6.11).

In all these cases the rolling-up vortex sheet has a dominating influence on the
pressure distribution and aerodynamic forces acting on the wing.

In the 1980's considerable advances have been made with the modeling by means of
panel methods of such type of flows, see for example ref. 2.47.
However, the advances have not yet led to applications in every day practice for
complicated practical configurations. Reasons for this are
 2.1143

- lack of robustness of the iteration process due to the strong nonlinearity of

the problem;
- relatively long computing times (partly because also on the vortex sheets a
large number of panels is needed);
- the general topology of the vortex sheet configuration must be known in
advance;
- the intersection of the vortex sheets with parts of the actual configuration
(wing-fuselage) is not known a priori; it is part of the solution and has to
be determined iteratively wich asks for readjustment of the panel
distribution.

A survey of panel methods having some form of capability for modeling of vortex
sheet roll-up can be found in ref. 2.47. They can be distinguished according to:

the dimension of the problem: conical flow, slender body approximation

( <' ,4 ), 2D time-dependent ('far-waket), 3D;
xx yy zz
the treatment of the edges of the vortex sheet (special measures are
necessary because of the singular behaviour of the vortex core that forms);
the degree of the representation of the geometry and doublet (or vortex-)
distribution on the configuration and vortex sheet;
the technique for iterating on the nonlinearity of the problem.

With respect to pt»4 a subdistinction may be made between

- relaxation methods ('wake relaxation'), and
- Newton-iteration.
In both cases a (good) initial estimate of the geometry of the vortex sheet and
of the doublet distribution (i.e. the direction of the vortex lines) on it is
needed. In the case of wake relaxation, at each iteration step, the doublet
strength is determined for a fixed estimate of the geometry of the vortex sheet.
Once this solution is known the velocity distribution on the vortex sheet is
computed (through multiplication of the matrix of aerodynamic influence
coefficients with the solution vector). Then, by integrating in flow direction,
a new approximation of the shape of the vortex sheet and the direction of the
vortex lines is made. This process is continued until the desired accuracy is
obtained. The procedure is depicted schematically in fig. 2.6.12.
The wake relaxation method is not a very robust technique and because of this it
is mainly applied to problems that are only 'mildly' nonlinear (vortex sheets
from the trailing edge of conventional wings). Fig. 2.6.13, as an example, gives
 2.144

results obtained with a special version of the NLR PANEL method of the computed
shape of the vortex sheet behind a swept wing after 4 iterations.

Newton iteration represents a more general method for solving systems of non-
linear equations that may be found in text boóks on numerical analysis, see for
example ref. 2.24. In this technique the boundary conditions on the solid part
of the configuration (wing, fuselage) and those on the vortex sheet, i.e.
- zero normal velocity and
- vortex lines are (average) streamlines,
lead to a system of nonlinear algebraic equations with as unknowns the doublet
(vortex) distribution and the parameters describing the geometry of the vortex
sheet.
This nonlinear system can, symbolically, be written as

Ñ() = 0 (2.6.16)

Here Ñ is a discrete nonlinear operator and a vector representing the unknown

doublet strengths and wake geometry parameters.

(n)
Let be an estimate of , satisfying

(n) (-(n)'
x J
=R-(n) (2.6.17)

(')
n1) of may be
where is the residual vector. Then, an improved estimate
obtained as follows. In the neighbourhood of the solution the residual vector R
can be expanded in a Taylor series. For each component R. of , we write

i
aR.
dx (n)
) + o[(o x (2.6.18)
n) j
j ix=x

or, with eq. (2.6.17) and neglecting second and higher order terms

aN.

--(n) dx.j (n) )=-R.i (n) (2.6.19)

i "ix=x
.n that drives R. to zero.
where ô x. is the correction to
j i
eq. (2.6.19)
Since the partial derivatives can be determined for

represents a system of linear algebraic equations in

I'
ô x.''. The matrix with
aN.
elements i is called the gradient matrix or Jacobian.
ax.
j
 2.145

The system (2.6.19) may be solved with one of the methods discussed in section
2.2.9. If ô is known a new estimate of follows from

-(n+1)
X
(n)
=X-(n) +ôx (2.6.20)

(n+1)
and a new residual vector may be computed from eq. (2.6.17).

Note that according to eq. (2.6.18) a condition for the convergence of the
Newton procedure is that the initial estimate should be sufficiently close to
the solution, a condition which may be difficult to fulfil in practice.

We further note that an iterative solution of the linear system (2.6.19) allows
a blending of iteration techniques. For example, if the Gauss-Seidel iteration
technique is applied to eq. (2.6.19), it is possible to perform already a Newton
step after one or a few Gauss-Seidel sweeps; this is called Quasi-Newton. Even
during the Gauss-Seidel process which is carried out row by row, one can
aN.
determine, for each individual row, ---i values based on the most recent
ax.
'J

information with respect to x. This procedure is called Successive (Over)

Relaxation-Newton. The wake relaxation process discussed earlier can be
considered as a special case of SOR-Newton. The application of Quasi-Newton or
SOR-Newton has computational advantages (less memory requirements than full
Newton), but the rate of convergence may be affected unfavourably.

Fig. 2.6.14 shows an example of aplication of a panel method to curling vortex

sheets in which the Quasi-Newton iteration technique is used. The example is
concerned with the flow over a delta wing with leading edge separation (ref.
2.48, NLR VORSEP code based on thin airfoil theory and quadratic doublet
distributions).
Methods for treating vortex sheets roll-up of the type described above are not
(yet) suitable for situations in which the vortex sheet intersects parts of the
configuration (see fig. 2.6.11). An additional problem then is that during each
iteration sweep parts of the configuration have to be repanelled also.
 2.146

2.6.4 Other methods

For completeness we mention that panel-type methods have also been developed for
- unsteady (periodic) flows described by Helmholtz' equation; in particular for
the purpose of prediction of unsteady airloads for flutter analysis
- evolutionary time-dependent flows
- the simulation of massively separated flows
- the prediction of flows in a rotational frame of reference (propellers,
helicopter rotars)
- the prediction of' propulsion-airframe interaction effects
- compressible potential flows (see part II)
The scope of this course does not permit to go into any detail of such
developments. The interested reader may find entry to the literature through
reference 2.0.

2.7 Inverse methods

Sofar we have discussed methods to compute the incompressible or subsonic

potential flow past configurations of' given shape. However, panel methods may
also be used to find (approximate) solutions to the inverse problem, i.e. for a
given pressure distribution (= given lift distribution) on a wing, find the
detailed geometry of the wing sections (including camber, twist, angle of
attack) that generate the given pressure distribution.
The nature of' this problem may be sketched as follows for the simplified case of
a 2D incompressible potential flow past a profile, fig. 2.7.la. We ask for the
solution of Laplace's equation for the potential in a region around a closed

contour f(x,z) = 0, satisfying the (boundary) conditions:

on f(x,z) = O (2.7.1)
= O
an

and
C(s) = g(s) on f(x,z) = O (2.7.2)

With the Bernoulli equation, eq. (2.1.9) we have

= i (2.7.3)
U2
 2.147

and because of (2.6.22) eq. (2.7.2) can be replaced by

= h(s) on f(x,z) = O (2.7.14)

8s

where

h(s) = U (± Ji-c (2.7.5)

(+ sign for s>s, - sign for s<s, where s represents the position of the
stagnation point (C =1) at the leading edge).

Note that the boundary condition (2.7.4) may be transformed by integration into
the Dirichiet condition

s
(s) = f h(s') ds' + (2.7.6)
o
s, =0

where is an integration constant.

o

The problem as sketched above is nonlinear in its boundary conditions and,

consequently, must be solved through iteration. It has a strong resemblance with
the problem of vortex sheet displacement and roll-up discussed earlier. In both
cases we have Neumann + Dirichiet-like boundary conditions on a surface the
geometry of which is part of the solution. The main difference is that the
inverse problem is concerned with a closed surface, whereas in the vortex sheet
problem an open, bounded surface is considered.
As in the case of linear (Neumann or Dirichiet) boundary value problems also the
nonlinear inverse problem raises the question if and/or what additional
conditions are necessary for a well-posed problem with a unique solution. It
appears that in the case of the non-linear inverse problem this is considerably
more complicated than in the case of the linear direct (given geometry) problem.
This is in particular true for the 3D inverse problem.
For a fairly detailed discussion of the problems associated with well-posedness
and uniqueness we refer to ref. 2.29 and ref. 2.49. Here we confine ourselves to
mentioning that at least the following conditions must be satisfied:
 2.148

- The total arclength or chordlength of a wing section has to be prescribed.

- In the lifting case a discontinuity surface behind the wing is necessary, over
which the jump in potential is equal to the local circulation. The latter is
known if the velocity distribution or the pressure distribution is prescribed
along the contour of the section.
- The prescribed pressure distribution should at least satisfy the conditions
(2D; similar conditions for 3D wings)

C <1
p-
C = 1 in one or two (stagnation) points (2.7.7)

C (s=O) - C (s=l) = O Kutta condition

p p

The situation with only one stagnation point occurs for a wing section with a
cusped tail were we have a stagnation point at the leading adge and no
stagnation point at the trailing edge (fig. 2.7.16).

The conditions given above are necessary but usually not sufficient to guarantee
a well-posed inverse problem with a unique and acceptable solution. Foir the
case of two-dimensional flow it is known (ref. 2.49, 2.55) that additional
conditions must be satisfied, explicitly or implicitly to obtain a valid and
unique solution, i.e.
- a 'closure' condition (or conditions)
- a 'regularity' condition.
A closure condition is required because the Neumann condition eq. (2.7.1) by
itself does not guarantee a closed contour or a given trailing edge thickness.
For example shapes like that of fig. 2.7.lc may satisfy (2.7.1) for any width of

the trailing edge gap. If there is a gap, there is also a resulting nett mass
flow from the airfoil to infinity. The level of this mass flow and, hence, the
width of the trailing edge gap, can be controlled through the level of the
potential at the airfoil which is determined by the constant of integration in
eq. (2.7.6). Alternative ways for imposing closure (through free parameters in
the pressure distribution) are discussed in ref. 2.55.

A regularity condition must be imposed at the leading edge to avoid singular

behaviour of the solution at the stagnation point at any stage of the iterative
solution process. This can be illustrated as follows.
 2.149

(n)
Let f (x,z) = O, where x = x(s), z = z(s), and j-- (s) = V(s), j (s) =
Vt(s) represent a uniformly valid approximation to the solution of the inverse
problem. The slope e of the airfoil contour is then given by

dx
= cos 8
ds
(2.7.8)
dz
= sin O
ds

and the local flow angle a by

V
n
a = 8 + arctan (2.7.9)
t

From the requirement that the airfoil is a stream surface it follows that

V
arctan - O
t

Since we have Vt = O at the stagnation point this requires that, at the

stagnation point

Vn sO

at any stage of the iterative solution process.

Without imposing the regularity condition the solution may tend to produce
shapes like that of fig. 2.7.ld. Imposing the regularity conditions requires a
free parameter in the specified pressure or velocity distribution (ref. 2.55).

Moreover, if a solution exists that satisfies all the conditions discussed

above, this solution is not necessarily acceptable from the point of view of
structures design or because of other aerodynamic requirements, etc.; for
example the wing profile may have unsuitable values for thickness, tail angle,
nose radius etc.
In 3D flow the situation is even more complex and, to the author's knowledge,
not (yet) fully understood.
Because of reasons discussed above the inverse problem is often formulated as a
minimization problem in which the prescribed pressure distribution and
requirements concerning the geometry are satisfied in an approximate sense, for
example according to the method of (weighted) least squares. A 2D example of
 2.150

such a formulation is a method in which a discretized form of a functional of

the type

i
F = J [w (c -c )2 + w (z-z )2 ds (2.7.10)
p g o
s=0 Pt

is minimized. Here C is the actual pressure distribution and C the prescribed

p Pt
or target pressure distribution; z and z represent the actual and required
geometry, respectively, and w, W are weighting functions to be chosen by the
designer.
Like the direct (Neumann) Laplace problem, the inverse problem allows the
solution to be represented by source- and doublet distributions on the (as yet
unknown) surface of the wing and a doublet distribution on the discontinuity
surface behind it. The possibilities for discretization are the saine as for
'ordinary' panel methods.
Eventually a nonlinear system of algebraic equations for the unknwon source-
and/or doublet distributions and the unknown geometry parameters is obtained.
This system has to be solved iteratively. In the literature two possibilities
are encoutered:
- Newton iteration (see section 2.6.3)
- linearization of the boundary conditions ('shape relaxation').
Ref. 2.50 describes a method utilising Newton iteration.
A simple example of application is given in fig. 2.7.2: only three iterations
are required for the (re)construction of a circle.

In the second method the problem is linearized by imposing the Dirichiet

boundary condition (2.7.6) on an estimate of the shape of the wing profile or
wing. The solution of this Dirichlet problem yields a normal-velocity
distribution over the contour or surface. Integration of the flow inclination
angle along the contour or surface results in a correction to the shape. In 2D
this becomes

s
dÇ(s) J (-/-1 ds (2.7.11)
s
o

where dÇ(s) is the correction normal to the profile contour as a function of the
arc length.
 2.151

The iteration procedure is sketched in the diagram of fig. 2.7.3. In general

only a few (3 to 4) iteration sweeps are necessary in order to obtain

convergence within the range of engineering accuracy.
The PANAIR code of Boeing (refs. 2.18, 2.19) contains an inverse module of the

type just described

A variation of the procedure depicted in fig. (2.7.3) has been developed at NLR
(code INSYST, ref. 2.51). In ref. 2.51 the geometry correction at each iteration
cycle is determined by means of inverse thin-wing theory type panel methods
(source-and doublet/vortex distributions on the wing reference plane). The
driving force in the inverse thin-wing module is the difference (residual)
between the required (target) pressure distribution and the pressure
distribution on a current estimate of the geometry. The latter is determined by
means of an 'ordinary1 panel method. The approach just described is called a
residual-correction method. A flow chart is shown in fig. 2.7.4.
The method allows for geometry constraints in a weighted least square sense.
Examples of' application of the NLR INSYST code are given in figs. 2.7.4/5. Fig.

2.7.4 illustrates a wing-alone case, without (a) and with (b) geometry
constraints. Fig. 2.7.5 deals with the same wing, but now attached to a body.
These examples illustrate, a.o.:
- that convergence to engineering accuracy is obtained in about 3 to 4
iterations
- the necessity to be able to apply geometry constraints.
 2.152

References
2.0 Katz, J. Low Speed Aerodynamics; From Wing Theory to Panel
Plotkin, A. Methods.
McGraw-Hill, Inc., 1991
ISBN 0_07_0S0L46_6.

2.1 Moran, J. An Introduction to Theoretical and Computational

Aerodynamics, John Wiley & Sons, New York, 1984.

2.2 Lamb, H. Hydrodynamics, Cambridge University Press, 7th ed.,

Cambridge/New York, 1975.

2.3 Kellog, 0.D. Foundations of Potential Theory, New York, Dover, 1953.

2.4 Mangler, K.W., Behaviour of the Vortex Sheet at the Trailing Edge of

Smith, J.H.B. Lifting Wing, RAE Rep. TR69049, Royal Aircraft

Establishment, Farnborough, 1969.

2.5 Batchelor, G.K. An Introduction to Fluid Dynamics, Cambridge University

Press, Carnbride/New York, 1970.

2.6 Hunt, B. Relationships Between Volume, Surface and Line

Distributions of Vorticity, Source and Doublicity, BAC
(MAD) Rep. Ae/3814, 1977.

2.7 Mikhlin, 5.0. Integral Equations, Fergamon Press, Oxford/London, 196'4.

2.8 Delves, L.M., Numerical Solution of Integral Equations, Clarendon

Walsh, J. Press, Oxford 1974.

2.9 Morino, L., Subsonic Potential Aerodynamics for Complex Configura

Kuo, C.C., tions: A General Theory, AIAA J.12, p. 191-197, 1974.

2.10 Hess, J.L., Calculation of Nonlifting Potential Flow About Arbitrary

Smith, A.M.O. Three-Dimensional Bodies, J. of Ship Research, September

196, p.22-1414.
 2.153

2.11 Rubbert, P.E. A General Three-Dimensional Potential Flow Method Applied

Saaris, G.R. to V/STOL Aerodynamics, SAE J.. Vol. 77, Sept 1969.

2.12 Labrujère, Th.E. An Approximate Method for the Calculation of the Pressure
Loeve, W., Distribution on Wing-Body Combinations at Subcritical
Slooff, J.W. Speeds, AGARD CP, No 71, 1970.

2.13 Kraus, W., Das Panelverfahren zur Berechnung der Druckverteilung von
Sacher, P. Flugkörpern in Unterschallbereich, Z. Flugw., 21 Jahrg.,
Heft 9, September 1973.

2.14 Roberts, A. Computation of Incompressible Flow About Bodies and Thick

Rundle, K. Wings Using the Spline-Mode System, BAC (CAD) Rep. Aero
Ma 19, 1972.

2.15 Hunt, B., Economic Improvements to the Mathematical Model in a

Semple, W.G. Plane/Constant Strength Panel Method, Paper presented at
Euromech Collog. No. 75, 1976. Also: VKI Lecture Series
1978-4.

2.16 Hess, J.L. Calculation of Potential Flow About Arbitrary 3-D Lifting
Bodies; Final Technical Report, Mc Donnell-Douglas Rep.
No. MDC J 5679-01, 1972. Also: Comp. Meth in Appl. Mech.
Eng., Vol. ¿4,
1974.

2.17 Hess, J.L., An Improved Higher-Order Panel Method for Three-Dimen

Friedman, D.M. sional Lifting Flow, Mc Donnell-Douglas Rep. No. NADC-
79277-60, 1981.

2.18 Johnson, F.T., Advanced Panel-Type Influence Coefficient Methods Applied

Rubbert, P.E. to Subsonic Flows, AIAA Paper 75-50, 1975.

2.19 Johnson, F.T. A General Panel Method for the Analysis and Design of
Arbitrary Configurations in Incompressible Flow, NASA CR-
3079, 1980.
 2.154

2.20 Maskew, B. Prediction of' Subsonic Aerodynamic Characteristics - A

Case for Low Order Panel Methods, AIAA Paper 81-0252,
1981.

2.21 Youngren, H.H. Comparison of Panel Method Formulations and its Influence
Bouchard, E.E. on the Development of QUADPAN, an Advanced Low Order
Coopersmith, R.M. Method, AIAA Paper 83-1827, 1983.
Miranda, L.R.

2.22 Bristow, D.R. Development of Panel Methods for Subsonic Analysis and
Design, NASA CR-3234, 1980

2.23 Botta, E.E.F. Calculation of' Potential Flow Around Bodies, Ph.D.
Thesis, Groningen, 1978.

2.24 Isaacson, E. Analysis of Numerical Methods, John Wiley and Sons, New
Keller, H.B. York, 1966.

2.25 Lötstedt, P. A Three-Dimensional Higher-Order Panel Method for

Subsonic Flow Problems, Description and Applications,
SAAB-Scania Rep. L-O-1 R100, 1984.

2.26 Hess, J.L. Consistent Velocity and Potential Expansions for Higher
Order Surface Singularity Methods, Mc Donnell Douglas
Rep. MDCJ 6911, 1975.

2.27 James, R.M. On the Remarkable Accuracy of the Vortex Lattice

Discretization in Thin Wing Theory, Mc Donnell Douglas
Rep. DAC 67211, 1969.

2.28 Joosen, C.J.J. Application of' the NLR Panel Method to a Swept-Wing/Body

Sytsma, H.A. Combination with Part-Span Flaps at Low Speed and

Comparison with Experimental Results, NLR TR 80030 U,
1979.

2.29 Slooff, J.W. Application of Computational Procedures in Aerodynamic

Design, AGARD Rep. 712, Paper 7, 1983.
 2.155

2.30 Hunt, B. Von Karman Institute Lecture Series, Rhôde-Saint-Genèse,

Belgium, Series 5, 1980.

2.31 Labrujère, Th.E. An Approximate Method for the Calculation of the Pressure

Loeve, W. Distribution on Wing-Body Combinations at Subcritical

Slooff, J.W. Speeds, NLR MP 70014 U, 1970.

2.32 Labrujère, Th.E. Aerodynamic Interference Between Aircraft Components;

Sytsma, H.A. Illustration of the Possibility of Predictions, NLR MP
72020 U, 1972.

2.33 Van Dyke, M.D. Perturbation Methods in Fluid Mechanics, Parabolic Press,
Stanford, 1975.

2.34 Lighthill, M.J. On Displacement Thickness, J. Fluid Mech. Vol. 4, p. 383,

1958.

2.35 Piers, W.J. Calculation of the Displacement Effect in Two-Dimensional

Slooff, J.W. Subsonic Attached Flow Around Aerofoils; Examples of
Calculations Using Measured Displacement Thicknesses, NLR
TR 72116 U, 1972.

2.36 Piers, W.J. Calculation of the Flow Around a Swept Wing, Taking into
Schipholt, G.J. Account the Effect of the Three-Dimensional Boundary
Vai-i den Berg, B. Layer; Part 1: Wing with Turbulent Boundary Layer, NLR FR

75076 U, 1975.

2.37 Oskam, B. A Calculation Method for the Viscous Flow Around Multi-
Component Airfoils, NLR TR 79097 U, 1980.

2.38 Oskam, B. Cotnputation of Viscous/Inviscid Interactions, AGARD CP

No. 291, 1981.

2.39 Oskam, B. Applications of Computational Fluid Dynamics in

Aeronautics, AGARD CF No. 412, 1986.
 2.156

2.40 Oskam, B. Recent Advances in Computational Methods to Solve the

Laan, D.J. High Lift Multi-Component Airfoil Problem, NLR MP 84042U
Volkers, D.F. 1984.

2.41 Woodward, F.A. Analysis and Design of Wing-Body Combinations at Subsonic

and Supersonic Speeds, J. Aircraft, Vol. 15, No. 16,
1968.

2.42 Woodward, F.A. An Improved Method for the Aerodynamic Analysis of Wing-
Body-Tail Configurations in Subsonic and Supersonic Flow,
NASA CR-2228, 1973.

2.43 Hoeijmakers, A Panel Method for the Determination of the Aerodynamic

H.W.M. Characteristics of' Complex Configurations in Linearized
Subsonic or Supersonic Flow, NLR TR 80124, 1980.

2.44 Ward, G.N. Linearized Theory of Steady High-Speed Flow, Cambridge

University Press, 1955.

2.45 Robinson, A. Wing Theory, Cambridge University Press, 1958.

Laurmann, J.A.

2.46 Fornasier, L. 'HISS-A Higher-Order Subsonic/Supersonic Singularity

Method for Calculating Linearized Potential Flow, AIAA
Paper 8141646, 1984.

2)47 Hoeijmakers, Numerical Simulation of Vortical Flow, Von Karman

H.W.M. Institute, Rhôde-Saint-Genèse, Belgium, Lecture Series
LS-1986-1O, 1986.

2.48 Hoeijmakers, Methods for Numerical Simulation of Leading-Edge Vortex

H.W.M. Flow, NLR MP 85052 U, 1985.

2.49 Slooff, J.W. Computational Methods for Subsonic and Transonic

Aerodynamic Design, AGARD Rep. No. 712, Paper No. 3,
1983.
 2.157

2.50 Labrujère, Th.E. MAD-A System for Computer Aided Analysis and Design of
Multi-Element Airfoils, NLR TR 83136 U, 1983.

2.51 Fray, J.M.J. A Constraint Inverse Method for the Aerodynamic Design of
Slooff, J.W. Thick Wings with Given Pressure Distribution in Subsonic
Flow, AGARD CP No. 285, Paper 16, 1980.

2.52 Hoeymakers, Computational Aerodynamics of Ordered Vortex Flow.

H.W.M. Ph.D.Thesis, Delft University of Technology, 1989.

2.53 Bos, H.J. LR-521, Delf t

Induced drag for non-planar wings. Report
University of Technology, faculty of Aerospace
Engineering, 1987.

2.54 Courant, R. & Supersonic Flow and Shock Waves. Interscience Publishers,
Friedrichs, K.O. Inc., N.Y., 1967.

2.55 Volpe, G. Inverse Airfoil Design. Applied Computational

Aerodynamics, Chapter 11, edited by P.A. Henne. Progress
in Astronautics and Aeronautics AIAA, 1990.
operating on a graphics termina(

Fig. 1.1
Exmp1e of d surface (computationat) grid.

Fig. 1.2
 DATE uç/II/2g. TUlE
X X X s EPRUG s 05/I /29.
X
ItRPIES USERSBYTSP5A

PART SEOPtNT COIITOUR PART SEßt)4T co4ToL

X X I X K 12

Examp'e of a surface (computationaL) grid.

Fig. 1.3
Examp'e of a spatia[ grid.

Fig. 1.4
 $IVI$FA A

Ì 'q $J?44d_____
I

tZ
I

pi,:,.

Example of a spatial grid.

Fig. 1.5
 UPPER SURFACE ISOBARS
A1I CATE TDE
177'
J! ALPHA- .ß.4 . CQTaR DirERVN.- . TICAL LT)

Example of graphical output.

Fig. 1.6
Fig. 1.7
Fig. 1.8
 Hardware developments continu fo be the driver for CFO
-persistent performance trend factor 10 er 9 years
106
Q conventional computers
Performance peak
(M flops)
O massively parallel
sustained

tera lo'

peak
/ xl000
speed up
lo'

giga
sustained

102

peak scalar
10°
peak scatar

10_2 I. I

195O 1960 1970 1980 1990 2000

Year introduced

-we are currently seeing the advent of massively parallel

computerswhich are expected to lead fo a quantum jump
of a factor 1000 at the turn of the century.
-(work stations)

Growth of computer speed and cost

Fig. 1.9
 z

Fig. 1.10

NAVIER - STOKES EQS.

(TIME DEPENDENT)
SPACE - TIME AVERAGING

SUB GRID SCALES ON LY

LARGE EDDY SUB - GRID SCALE

ALL TURBULENCE SIMULATION TURBULENCE MODEL
SCALES I

I
REYNOLDS AVERAGED
NAVIER STOKES EQS.
IN N N /

SYMPTOTIC THEORY
ZONAL MODELLING COMPLETE
INVISCID TURBULENCE
MODEL

EULER EQS.

IRROTATÇONALITY

FULL BOUNDARY SIMPLIFIED

INTERACTION
POTENTIAL LAYER TURBULENCE
LAW
EQ. EQS. MODEL
LINEAR COMPRESSIBILITY

LINEARIZED
POTENTIAL
EQ. HIERARCHY OF MATH./PHYSJCAL MODELS

Fig. 1.11
 ONSET OF B.L. SEPARATION

t LARGE - SCALE SEPARATION

ANGLE
OF
ATTACK BOUNDARY FOR TRANSONIC EFFECTS
OR CRITICAL T'tACH NUMBER
LI FT

LINEAR INVISCID FLIGHT ENVELOPE

OR EULER FULLÀ
POTENTIAL
OR
IR

(WEAK VISC. / INV. INTERACTIOÌ

MACH NUMBER

Areas of applicability of various f(owmode[s (transport aircraft)

Fig. 1.12
 Fig. 1.13

¡JI div, pd .fld 5e°

vol. e
surf. Se
Fig. 1.14

Fig. 1.15
 RE-AVERAGED N-S

RELAT IVE \
COMPTL. \
EFFORT - EULER

'O FULL POTENTIAL

-
i
LINEAR INVISCID
(PANEL)
I I I I
.1
1970 1975 1980 1985 1990 YEAR
(NL SITUATION)

Reduction of conputaton1 effort due to more effident Igorithms.

Fig. 1.16
 X
L1
LJ
o'rn.
L) I
'- NX
QE
Eo zL)o
>-LJ
>
ro
>_" >-, VI

.-
LJ L) L) ro r.j

10'
- TR-
io7 rn
C ompu t a Hon
time per flow -MO-
z
(seconds) 3-D wing-body Reynolds average N-S
10'
-WK-

-DAY- full potential

linearized potential

io'
-FIR
1/2

10-
min
- Euler
102

min
I i
10
0.1 1 10 102 io3 io' 10'

Computation speed Mflops

Computation time versus computation speed for a fuselage,
wing configuration ("engineering accuracy)

Fig. 1.17
 10

memory
(words)
1015

lo

Fut) Pot.
W Wing
Euler AC Aircraft
io

Reynolds averaged N-S

.. Cray YMP/832
l0
Nec SX-3 at NLR
Cray i
/<PaneL '
i'
J J

102
1018 1021
io3 io' lO9 1012

speed (flops)

Computercapacity required for a computation of 1/2 hr approximately

Fig. 1.18
 HOD BASE

PREPROCESSORS
(GEOMETRY/MESH FLOW SOLVER(S) POSTP ROC ESSO RS
GENERATION)

APPLICATIONAL ASPE[T Necessity for (informatics) system approach.

Fig. 1.19
 Fig. 2.1.1

Fig. 2.1.2

V.

f (s,t)

Fig. 2.1.3

Fig. 2.1.4
 or

G) b)

t=g(s,t)CONTINUOUS on S and Zx DIFFERENTIABLE

Fig. 2.1.5

or

Fig. 2.1.6

=;7. =f )s,t)
òn
onS
(CONTINUOUS I

Fig. 2.1.7

Fig. 2 .1.8
 n g(s,t)

so
Q

òn
(st)

Fig. 2.1.9

Fig. 2.1.10

g(st) òS
(st)
ds

ÒQ
(f:gffì,,j:,4 B'

A
òSc

= f(s,t)
òn

Fig. 2.1.11
 Fig. 2.1.12

Fig. 2.1.13

Fig. 2.1.14

Fig. 2.1.15
 òn
P, T

Fig. 2.1.16

t
Sc

Fig. 2.1.17

Fig. 2.1.18
 [P]1const. r
I s3

'j Sc

X st

Fig. 2.1.19

Fig. 2 .1. 20
 Fig. 2.1.21

if
o .ònJi

Sc

Fig. 2 .1. 22

IT ò
òn òs

T Sc

Pl

Fig. 2.1.23
 con st
òn

Fig. 2.1.24

Fig. 2.1.25

òn - n

ILine []const.

lines []= corst.

Fig. 2.1.26
 n

/
alo, i
/
source strength per unit area

Fig. 2.1.27a

Fig. 2.1.27b

Fig. 2.1.28 0= const.

Fig. 2.1. 29

0= const.

Fig. 2.1.30
*- -
Fig. 2.1.31
 Fig. 2.1. 32

tangential normal

Fig. 2.1.33

Fig. 2.1.34
 Fig. 2.1.36
PzO

no Fig. 2.1.37

S iconst.

C-
Fig. 2.1.38

I z const.

Fig. 2.1.39

Fig. 2.1.40
 Fig. 2.1.41

Fig. 2 .1. 42

-1
n
.Pl

i:i

/ s
Ql

(2
I
Q

\
//
Fig. 2.1.43

Fig. 2.1.44
 dSO
s
t net rncisspraduction)

Fig.,2 .1. 45

B
'
òs

Fig. 2.1.46

R) determincted by kuttQ condition

SB

Fig. 2.1.47

's
Fig. 2.1.48

Fig. 2 .1. 49
Fig. 2.2.1

Fig. 2.2.2

Fig. 2.2.3
 Patch
Grid Point Network

Conti g u ration

Fig. 2.2.4

Fig. 2.2.5

Fig. 2.2.6

Z ref

7 'T'ref

Fig. 2.2.7
 Fig. 2.2.8

Upper W-B Fairing

Lower W-B Fairing
Lower
Body
Inbd Wing

Upper
IT I. 3

Body
'u Lower
I" I"
I
4 3 AiLeron
III
'u 'ulu
In Ti
2
Aft- Body

2
Mid Win
Upper
Aileron

Fig. 2.2.9
n

vector (13)
vector (2. 4) }
is the positive normal
vector defined in the
panel collocation point
pI
/4
1,
/ ..--
c 4/-

Fig. 2.2.10
 Z ref
A

' ref

Fig. 2.2.11

Fig. 2.2.13 Example of "pin-cushion" figure.

 a) ithcut "hidden tines'.

b) with 'hidden Eines'.

Fig.2.2.12; Example of "wire" figures.

 Fig. 2.2.14

PLANAR representations

Flat surface

Straight
Line

Pure planar Parameterized Segmented planar

(least squares fit through (hyp erb o t o idal)
center of gravity)
Resultant surface nearly continuous Resultant surface continuous Resultant surface continuous
Slope discontinuous Slope discontinuous e Slope discontinuous
Totally local i.e shape depends Totally ocal i.e shape depends Totally local i.e. shape depends
on four corner points only on four corner points only on four corner points only
(Douglas Boeing TAE 230; NLR, PANAIR, MCAERO)
MSS VSAERO OUADPAN)

Fig. 2.2.15

O.UADRAT(C representations

Quadratic Surface

Plane
Quadratic
Curve Edge

Pure paraboloidal P ara me te rized Explicit (segmented)

Resultant surface nearly continuous Resultant surface continuous Resultant surface continuous
Slope nearly continuous Stooe nearly continuous Slope nearly continuous
Localbut shape depends on Local,but shape depends on Local,but shape depends on
neighbouring points neighbouring points neighbouring points

Fig 2.2.16
 X

- -k - unit vectors

Fig. 2.2.17

Fig. 2.2.18

Fig. 2.2.19

)
 fIx)

X2 X3 X4 X5

Fig. 2.2.20

SINI3ULARITY SPLINE CONSTRUCTION

s Constant sp(ine

f(x)

f7
f6 -a f8
s u
f 5-_
u- \ slope f expressed as slope of
I

f2
tine through neighbouring values
f4
f3

X0 xl X2 X3 X' X5 X6 X7 X

TotaLLy LocaL

Discontinuous

Stable for f specified at midpoints

Absolutely stable t source/Neumann doubtet/Dirichiet

(Douglas-Nn (VSAERO,UAOPAN,
Boeing TEA 230 AEROPAN)
NLRMBB
Hunt-Semple (BAC)

NeutraLly stable: doublet/Neuman (vortex lattice) source/Dirichtet

Unstable for f' specified at midpoints: vorticity/Neumann

C linear spline )derivative is constant spline)
fo

fix)

I I

X0 X2 x3 x X5 X6 X7 Xe

Local

Continuous

Absolutely stable for f specified at corner points : Not usable

Neutrally stable for f0 specified and f' specified at midpoints: vortex/Neumann

Generally unstable for f specified at midpoints with any other conditions:

source/Neumann; source/Dirichlet doublet/Neumann; doublet/Dirichiet

Fig. 2.2.22
Least square linear spline

fix)

t N X-
fo

slope obtained from slope of

line through neighbouring values

X0 xl \ X2
f3

X3 X4 X5 X6 x7 X0

Good and simple

Nearly continuous

Absolutely stable for f specified al midpoints: source/Neumann (Hess)

double t/Dirichlet

Neutrally stable source/Neumann; Doublet/Neumann

Fig. 2.2.23
 C q uadratic spline (derivative is C0 linear sptine)

f(x)

X2 X3 X' X5 X6 Xl Xe

Not recommendabLe
Continuous with continuous slope

Neutrally stable for f specified at midpoints and x0 and xe source/Neumann;

doublet/Dirichiet

Unstable for f0 specified along with f' at midpoint and X8: vortex/Neumann
Fig. 2.2.24
s C0 corner point least square quadratic spline (derivative is Least sq. Linear sp(ine)

f(x)

Parabola which interpolates f6 and f7 exactly

and f5 and f8 ri least square sense

X0 Xl X2 X3 X4 X5 X6 X7 X8

Continuous with nearly continuous slope

Absolutely stable for f specified at corner points: Not usable

Neutrally stable for f0 specified and f' specified at midpoints: vortex/Neumann

e Oenerally unstable for f specified at midpoints along with any other condition:
source/Neumann ; doubtet/Oirichlet : source/Dirichlet Doublet/Neumann

Fig. 2.2.25
 define f at corner points by least square
C0 least square quadratic spline
fitting parabola to 4 neighbouring values
(heavy weights on 2 closest values)
then fit parabola to f7 and resultant
fIx) corner values.

Good but rather complex

Continuous with nearly continuous slope

Absolutely stable for f specified at midpoints and x0 and x1

source/Neumann; doubtet/Dirichlet (PANAIR, MCAERO).

e By forward weighting least square procedure, can be made stable for

specifying f0 along with f' at midpoints and x8: source/derivative Dirichtet;
vortex/Neumann

Fig. 2.2.26
[east square quadratic spline

f(x)

fi

Parabola which interpolates f6,f7 ,f1

f3

X0 X1 X2 X3 X x5 X6 X7 X0

Good and simple

'4early continuous with nearly continuous slope

Absolutely stable for f specified at midpoints and x0 and x0

source/Neumann; doubtet/Dirichlet

Fig. 2.2.27
 100

-9 L2-E
C
(Vt)
-8
1st ORDER METHOD

-7 I 21 ORDER METHOD

-6

-5

10
-4

-3

-2

-1 AREA 0F INTEREST

o X/C
near-optimal (curvature adapted) spaci

10-'
100 iO1 102
2
- -
NUMBER 0F PANELS

100
Fig. 2.2.28
L2-E
(Vï
st ORDER METHOD

1 A 2nd ORDER METHOD

10.2

AREA OF
INTEREST

non-optimal spacing

1°-4
100
J

101
I ill 102 N io3
NUMBER OF PANELS

Fig. 2.2.29
cP

- 0.8

-0.4

x/C
0.0

-s-. ROBERTS (OATUM) error

0.4
O NLR (60x12) 1%
X NLR (30x12) 1.5%

0.8
CONVERGENCE STUDY
CI-IOROWISE PRESSURE DISTRIBUTION
RAE WING TIC = .15. a= 5.0 .q= .549
1.2

Fig. 2.2.30
CI,

-0.4

-0.2

x/C
0.1 0.2 03 0.4 05 O6 0.7
0.0

0.2 -- ROBERTS (DATUM) error

f u H-S(SHEETS) (60x12) 1%
BAC
'L x H-S(SHEETS) (30x12) 3%
o NLR (60x12) 1.5%
0.6 (30x12) 4%
+ NLR

CONVERGENCE STUDY

0.6 CHORDW!SE PRESSURE DISTRIBUTION

REA W)NG.T/C .05. u5.0.t.549

Fig. 2.2.31
 influence of da(0)
Curved panel w=2.[°i
2 remote sources dx
subtending an
angle o
This source induces
normal and
tangential velocity Error analysis for linearly
components at varying o strengths
boundary
point

i9) [i+_P_l_error 0(8)

2
L 2TTJ

curvature
effect
Error in w is 0(0) and uncouptes from
dx
Error in u is 0(82) and uncouptes from a
-!2iS?2 [ i - error =0(82)
2n dx
NoterRoles of u & w reversed for vortex panels L

curvature
Remark Flat source panels present a serious effect
problem for circular duct simulation.

Fig. 2.2.32

[p
3

-- ROBERTS (DATUM) .)
-0.8 o H-S(SHEETS (60x10) Ì
0 NLR (60x10)J

-0.4
outside
0.8
00

0.4

3rd
) order source-based method
0.8 1St
order source-based methods

1.0 FIRST ORDER METHOD COMPARISON

CHOROWISE PRESSURE DISTRIBUTION
NACELLE :C/DE=3.3333 . s =0.0 . 8= 0.000

Fig. 2.2.33
 CALCULATED SURFACE VELOCITY DISTRIBUTIONS IN A CIRCULAR DUCT

HIGH ACCURACY SOLUTION (AXISYMMETRIC)

3-DIMENSIONAL HU3HER-ORDER SOLUTIO4 756 EFFECTIVE SOURCE PANELS (42x18)

3-DIMENSIONAL iST ORDER METHOD SOLUTION 756 EFFECTIVE SOURCE PANELS (42x18)

3DIMENSIONAL iST ORDER METHOD SOLUTION '+200

EFFECTIVE SOURCE PANELS (62x100 OR 84x50x1

2 3 4
X

PLANE OF
S Y MME TRY

vo

AXIS OF SYMMETRY

Fig. 2.2.34
 A N P

EED 000
[]flo
0EDD
rioDD
E1D r' r' r'
I I
i.._i L_J
+ r-i II
I
e---'
L....J LJ
I

r--iIir--, r-, o
EDDD I I I
L_J L_J L_J
I I I
L_J L._J t_J i I

Fig. 2.2.35

loo

C
E Boeing TEA23O
ni
60
.4-
C-
o
I, DIRECT METHOD
In
w
e 40
C- ITERATIVE METHOD
o.
X X
X

20
X NLR Panel

o
o 500 1000 1500 2000
Number of singularities
Comparison central processor times on CDC 6600 computer of vortex source method
for a direct and for an indirect method of solution.

Fig. 2.2.36

2000 c' HESS (istorder direct) 1.0 0 HESS(15t order direct)

A VSAERO VSAERO
'J In
ni t UAOPAN o
L)
0.8 UADPAN
1500
w MCAERO MCAERO
E PAN AIR
I.- . 0.6 O PAN AIR
1000 E
0.4

500
0.2

200 400 600 200 400 600

Total panels Total panels
comparison of CYBER-175 CPU execution time Comparison of CYBER-175 job costs

Fig. 2.2.37
 z

s--

Fig. 2.2.38

0.0 08

0.006
60x12 pane's

¿ Trefftz-pLne integr. C0.

o Pressure integr. C

(COpO.00O2 at CLO)

2
0.1 LL

Fig. 2.2.39
 - - - - Corrected slon. of lin, pot. eqn.
-- -- Corrected soin, of tin, pot. eqn.
Exact( SELLS) full potential equation
Exact (NIEUWLANIJ) full potential equation
- - Solution linearised potential equation (64 contour
-1.0 -1.0 points) GUTHERT II
- Solution iinearised potential equation (68 contour
Cp p
points ) GOTHERT if local1
-0.8 -0.8
Cp*(M(ocaL:: 1)

r
-0.6 -0.6

\
-0.4 -0.4

-0.2 -0.2

0.0 00
0.2 0.4 0.6 0.8 1.0

0.2 0.2

0.6 0.4

0.6 0.6
//
0.8
I Symmetrical non-lifting aerofoil 0.8 ¡ Lifting aerofoil

1.0 1,0

auasi-etliptical aerofoit section t/c=0.17; M=0.714; cx=0° NPL 3111 aerofoit section t/crO.34 M0.667;a'1.2°

Fig. 2.3.1
 DESIGN
/1OBJECTWES

PHASE I PRELIMINARY
DESIGN

/ MAJOR
( DIMENSIONS
IPLANFORM
ETC.)
L.
i
DESIGN
PHASE U DEVELOPMENT CURRENT FOCUS

/ DETAILED
/ GEOMETRY
T

/
Major phases of aerodynamk design process.

Fig. 2.4.1
 // /PERFORMANCE GOALS,
CONSTRAINTS

MAN - IN - THE - LOOP

.ThI
J I

f/ AERODYNAMIC
/ OBJECTIVES,
CONSTRAINTS

AERODYNAMIC
DESIGN CONTROL
SYSTEM

CONTROL
AERODYNAMIC
ANALYSIS
SYSTEM

MAN - BESIDE - THE - LOOP

WIN DTU N N E L PROCESS CONTROL

SYSTEM

-J

/GEOMETRY,
AERODYNAMICS OBTAINED

TypicaL flow chart for aerodynamic design process.

Fig. 2.4.2
 r)

A
C

3
PRESSURE
WII 413 SECTION
PLOTTING 2
RA E 101,90/ THICK
SECTIONS
BOIJY PRESSURE ON WING
PLOTTING SECTION
-V

Fig. 2.4.3

Ø experiments

NLR panel method 200(body) +240(wing) paneLs

half-configuration

/
CL

0.6
/4 onset of boundary

// Layer separation
in experiment

0.4-
/o

0.2 0

-8° -4° 8° 12°

-
16°

-0.4

o
-06
o
,

Fig. 2.4.4
 Ø experiments

NLR pane method (200+240 p3nes)

C
25

0.4

02 0

of separation

CL

:I:

Fig. 2.4.5
r2
LL
Trefftz integrated

00
pressures
plane
o
¶
0.6
I

0.5
4

i
/ viscous drag

X
0.4

X
0.3 I I

'I I

02
X

0.1 X

o 0.1 0.2 0.3

C0- C00 (ex p.

C01(comp.)

O experiments

CDi(p) NLR panel method (200240 panels)

--x--x-- COitTr)NLR pane method (200+240 panets)

Fig. 2.4.6
 PRESSURE
PLOT TI N G WING SECTION
SECTIONS RAE 1019I. THICK
BODY PRESSURE ON WING
PLOTTING SECT IONS 31
5
6
i

X NLR METHOD
body '6 o EXPERIMENT DFL
OEwing a6 Rec s 0.65 x iO6

0.6
C

0.4

0.2

0.2 0.4 0.6 0.8 1.0

- 780 paneLs on wing

MEASURED AND CALCULATED LIFT DISTRIBUTION ON A STRAIGHT

WING-BODY COMBINATION IN INCOMPRESSIBLE FLOW

Fig. 2.4.7
 M0
3
PRESSURE
PLOTTING OEbodyr 6G
SECTIONS 2
NLR METHOD
ON WING o EXPERIMENT DFL
Rec z 065 x106
-x
Cp
-2.0 SECTION 1 -2.0 SECTION 2

X IC

SECTION 3

+1.0

COMPARISON OF MEASURED ANO CALCULATED PRESSURE DtSTR8UTIONS ON WiNG

,i o
 M'O
bodyl6 awjnga6°

PRESSURE -. - NLR METHOD

BODY PLOTTING
PRESSURE PLOTTING SECTIONS o EXPERIMENT DFL
SECTIONS 31 ON WING
RecO.65X1O6

BODY SECTION BODY SECTION

Cp
-0.4
i lo

0.4

COMPARiSON OF CALCULATED ANO MEASURED BODY PRESSURES

Fig. 2.4.9
 PRESSURE
PLOTTING
SECT N S
4 ON WING

BODY PRESSE
PLOTTING SECTKDNS
Y

II__i,i,IIIIHIIUUIIIIUI____

Ma O4 WIflQ.body I CALCULATKDN
a .637 -wIflg alon. J
EXPERIMENT
Re.l 3x1O6

C{ 8

t6
4

o +
ßOOY SIDE TIP

Panet arrangement for a typical short range,subsonic Spanwise Lift distribution for a typical short rangesubsonic
transport configuration. transport configuration.

Fig. 2.4.10
 PPSSuRE
PWTTI3
SEC
ON WING

X
Ma .04 wriQ.body
U :837 - - wig alone ) CALCULATON
EXPEPIMENT Iz
Peel 53x106

No.04 - CAtC1LAt4(wç .by)

CD a .637 ttRS.4P4T
-06 (9.body *.4.,J

-04
/S(Cl t (TOP)
.02

.02

.04

CD
06

.Ø4

-02

.02

Wing pressure distributions for a typical short range Body pressures for a typical short range.
subsonic transport configuration. subsonic transport configuration.
Fig. 2.4.11
 -

Generai arrangement and panel plot of

wing-body-tail configuration.

a)

WAKE POSTIOP4( --
f--- U637
t - 1ITIAL EST IMAÎE
eo
z
ST AB Ij E R
40

o
o 200 300 400 5Q0
'Y

-40

-60

Wake positions at horizontal stabilizer

b)

Fig. 2.4.12
 - .T_.... - 200 panels on sting

WITHOUT STING }WAKE DISPL.

WITH STING 1bCALC
5t WAKE ESTIMATE (WITHOUT STING)J
o EXPERIMENT

015 C
0.12
CM U US(LAOC Sa1ER
0.12

0.06
,rnuu
.. u...
....
oe
CL

02 0.4
0.04
0.04

0.4
L

......rn
-- - - ein STPG
o - W1Th40Ut STP4O } 5pt.ACCO
0.2 0.4\ z
O o8 1.0
E
lt WA.E ES ATE IW1T4JT STING)
}CMCULATI0.
CL

-0.04 b)
Coritriution to the pitching oorent
of parts of the corifigu.ratior..
-006

-0.12 0.4
WITH STING
a) WITHOUT STINGS
DISPLACED WAKE
1

CALC
Total pitching coer.t versus tota..]. -. - 1st WAKE ESTIMATE(WITHOUT STING) J
lift. e L EXPERIMENT

-o
r 0.4
WITHOUT STING }CALCULATION
--WITH STING
EXPERIMENT
02
Ct
0.1

o
r63

200
AAA
120 160 2
08 OY

-0.1
04

O
-02
Y
c)
CI.an.ise liftdistributior. or wiri. -0.3
u u
d)
Spanvis .]]ftdistritutlon on
j ter.

Fig. 2.4.13
 TOP VIEW

1500

0 200

ELLIPSOID
R= oo\yH ._(.__)
11%C 0.5%

CROSS - SECTION C-C

(ENLARGED SCALE)

k\oo
MOMENT REFERENCE
POINT

300

PARABOLOID
R 100 [1

INCLINATION OF WING CHORD TO

50
BODY AXIS SECTION A - A
SECTION B - B 1°

BOX-TYPE ADDITION TO
FUSELAGE

GeometricaL data of the configuration. Dimensions in mittimeters.

Fig. 2.4.14a
 (I

'h

'A
'lIIwaMffD'DDMM.DMW 1I111

_M'li4 ////,,,,,,,, ,,,,, 'u//hIll,

Ríii,..u.u.urn...u,.a...
------------.-.-.-.-.-u.---.--
1IIItIIIIIIIIIIIIIIIIUIIIILIIUiII

Example of panel distribution (configuration W8 * F0)

Fig. 2.4.14b
 Body centerline

Section A-A

a. Calculation model (with gap)

Body centerline

Section A-A

b. Calculation model (trailing vortices aligned with flap)

Configuration WBF25

MODEL I

4/

WING VORTEX SHEET LOCATION OBSERVED WITH

VORTEX INDICATOR AT 1S, BEHIND FLAP
TRAILING EDGE

Vortex sheet Locations used in the calculations for 6F 250

Fig. 2.4.15
 CALCULATION MODEL

CONFIGURAT!ON WBF25

WING + FLAP
1.5 /
7)<
CLW + CLF,
7
7
CLw 7
7 X WING
CLF 7
.7
7X7
1.25
7
7 -,
7
7
.7 7/
7 ,-
7
X 7
7
7- /4<
X
0.75 7
7 X
7
7
X

0.5

0.25
FLAP

O
-2.5 0 25 5 75 10
a

Wing and flap lift versus angle of attack (calculation models

IB and II)

Fig. 2.4.16
 CALCULATION MODEL

Is CALCULATION MODEL

U IB
CONFIGURATION W8F25
EXPERIMENT X n
EXPERIMENT X
CONFIGURATION WBF25

SEPARATION AT FLAP TIP EDGES

1.25 -- s..
X
-- X

[a a7i°
0.4

1.5

0.2
--- X

L25
0.25 05 0,75
1

X 0.4 La =
X

SEPARATION AT
Xi TIP EDGE
0.2
0.75
a \
II
X
\
\
I
I

8 o
O 0.25 0.5 0.75
\ r)
(a = 87501
0.5
0480
0.4 La 1
CD
\0
0.25 La r 0480
0.2

O
O
o 0.25 05 0.75
o 0.25 05 0.75
7) '1

Spanwine wing lilt distribution

Spanwise flap lift distribution
 CALCULA TI ON kIOOE L

TE
PRESSURE SECTION
-4 u
EXPFRIMENT 7
CD

-3

CONFIGURATION WBF20

77 = 0.208 a = 4.6201
-2

X'X'
-t

o 0.25 0.5 0.75 I O 0.2 0.4

-4 X/C X/C

-3

CONFIGURATION W8F25

77 = 0.500
-2

cP

I I I I

-3 C 0.25 0.5 0.75 1 0 0.2 0.4

X / C_ X/C
CONFIGURATION W8F25

-2 77 = a = 4.620

SEPARATION / VORTEX
FORMATION AT TIP EDGE

o 0.25 0.5 0.75

X/

Fig. 2.4.18
 CONFIGURATION WBF?S
CALCULATION MODEL

L PRESXURE SECTIoNS

CR EXPERIMENT CU j
= 4.62°

-2

= 0.407 q u 0.743

SEPARATION/VORTEX FORMATION AT SIDE EDGES

I I I I i t I I I
0.2 0.4 0 0.2 0.4
o 0.2 0.4 0
X/C X/CW X / C,
Finp prennuro diWtributiona

a)

2.8 I 308
j I T
CONFIGURATION WBF25
(CI.W CLF)2 CONFIGURATION WB

2.4
0.8 (j. 0A EXPERIMENT i
2 CALO. MODEL T
LW1

0.4

1.6
044 0 0.02 0.04 0.06
CO,

f, CALCULATION MOOSL

is

-N
1.2

CONFIGURATION WBF0
EXPERIMFNT

0.8
4 4

0.4

o
0
44- 0.02 0.04 0.06 0.08 0.10 0.12 0.16 0.18
C0
0.2

Iflducod dra0 aa a function of the lift

b)

Fig. 2.4.19
 .-j...uIU..UU ulîi;; -

I//I"
I/I/fill
'I//fl,', a) CONFIGURATION PANELING
11111111

/1/lili!
ill"
/1111111
.V///I/llu,lIf

,,,/I,I,,lI,u
ii.
IIIIuIa___
liuLIilLil_Il -a..

NACELLE OFF
CZ
rn.f.r. 2.3 NACELLE
m.f.r. 0.4 ON

2.5
/ \\
I
m.f.r. = MASS FLOW RATIO

4/
// \°-. .

2.0

b ) SPANWtSE LIFTDISTRIBUTION .

1.5

CENTRE-LINE NACELLE

1
o .2 4 .6 .8 10v

Effect of naceL'e flow ori lift (take-off configuration a = 18°) (Ref. 2.29).

Fig. 2.4.20
 M..O.e cL.3
wrTH PYLOI I STO4E } EXR (Rt.14..40A1
r WITHOUT PYLOII I STORE
wrTH PYLON ¡STORE
WITHOUT PVWN/ } CALCULATION

Panel arrangement of wing-tiptank-

pylon-store configuration

Spanwise liftdisfribution on wing

cp
- 0.4
PYLON
-0.2

o
A
o X

4_
o X d UPPERSOC
SS

Chordwise pressure distributions ori

wing (Ma0.8,x:3°)
o X

0.2
(from: Ref. 2.32)
0.4
Store pressure distributions
(Ma 0.8 30)

Fig. 2.4.21
 a)GEOMETRICS OF CONFIGURATION

t etC...I.tIe., s. S a.
4..'. 3
n, .a..,e -.

(meet VIO. ret.. .,t*

cp*

- 1-

(from: Ref. 2.29)

b) COMPARISON OF MEASURED AND CALCUIATED PRESSURES

Rear-mounted high-bypass ratio naceLLes a locaL interference problem

Fig. 2.4.22
 F LIGHT CONF.

WINDTUNNEL CONF

LOCATION OF TUNNELWALL

Cp

a ) CONFIGURATION PANELING

o
C-

Cp

0.2

1.0 X/C

b) WING LOWER SURFACE PRESSURES

Under-wing pylon/naceLle interference study (Ref. 2.29)

Fig. 2.4.23a
 C.

RO'

ROS

75 i
C EXP
ROS RO? - CALC

O NACELLE
STAR N R O)

NACELLE PRESSURES

Under-wing pylon/nacelle interference study (Ref. 2.29)

 POTENTIAL FLOW

r
PRESSURE DISTRIBUTION

BOUNDARY LAYER COMPUTATION

DISPLACEMENT THICKNESS
I

POTENTIAL FLOW WITH. DISPLACEMENT EFFECT

Fig. 2.5.1

new body
IpI=O

oriqinat body

ADDING ON DISPLACEMENT THICKNESS

Fig. 2.5.2

- dividing streamLine

SURFACE BLOWING

Fig. 2.5.3
cp

NLR 7301 PLUS 32 PERCENT FLAP

a 60 , Re 2.5x io6
o EXP. UPPER SURFACE REF. 33
EXP. LOWER SURFACE McOl9
- INVISCID / VISCOUS THEORY

I I i I I I I

0 0.2 0.4 0.6 0.8 0.9 1 1.1 1.2

X X

FLAP DEFL. 20°

GAP 2.6
OVERLAP 5.3 %
0.1
Y
o

0.1-
I I I I I I.
0 0.2 0.4 0.6 0.8 1

NLR 7301ftap configuration; comparison of pressure distribution (Ref.2.37).

Fin
 3.8
NLM 7301 + FLAP

3.4
INVISCID THEORY
C,
riREDICTEO
TURBULENT
f3 o BOUNDARY 3.6
'I LAYER NIR 7301 + FLAP
SEPARATION
C,
WING
UPPER
3.0 - t3.4
SURF ACE THEORY
g0-.----- F LAP WAKE
o
2.8 - 3.2 WING WAKE

2.6 -
/ 3.0

2.4 - 28
/ 0.-o

2.2 - 2.8 /
/

2.0 2.4 -o
R. 2.5 MILLION Re 2.5 MILLION

.. VISCOUS THEORY (N_ -01 VISCOUS THEORY IM0I

1.8 2.2 4-
-O-O- EXP DATA (M_. 0.1851
WITH BLOWING EXP. DATA M -01851
1.8 2.0
-V-V- EXP.DATAIM..- 0.1851 -O--O- SPANWISE AVERAGE
WITHOUT BLOWING
SPANWISE VARIATION

1.4e 1.8
4 B 12 18 lB 2(8) 300 400 500 500
o 100
(D EGH E E S) . 1OC

Corinparlion of lIft-drag curveS between viscous theory end enperlesni

Coopaelson of lift coeffIcIent versus angle of attack between VISCOUS for NLRI3OI plUs flap
theory and copen i.ent for IILR73OI plan flap

Fig. 2.54b
 CONFIGURATION 2

t . . .
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6
rX

Wing with double-slotted flap; denoted by configuratIon 2

4.5

CONFIGURATION 2

Re = 2.6 MILLION

INVISCID THEORY
M= O

4.0-

O
t
3.5-

3.0-

..- VISCOUS THEORY (M0, = O)

O---0 EXP. DATA (M,= 0.19)

2.5
0 2 4 6 8 10 12
- a(DEGREES)
Comparison of Lift versus angle of attack between viscous theory and
experimental data for confiquration 2 (Ret.2.0).

Fin. 2.5 S
 4.0
CONFIGURATION 2
C4 Ri - 2.6 MILLION

5.

8
7. i1
3.5

VISCOUS THEORY
=0

EXP. DATA
/ M = 0.19

3.0
/6 OSPANWISE AVERAGE

/ I ISPANWISE VARIATION

I I
2.5
o 200 400 600 800
- 1O Cd
Comparison of lift-versus-drag curve between viscous theory and
experimental data for configuration 2
Ci

CONFIGURATION 2

-4 ReS 2.6 io: a=g

- vISCOUS THEORY (M= 0)

-3 O O EXP. DATA lM0 = 0.19)

-2
SEPARATION

Pressure distribution on vane and flap of configuration 2

Fig. 2.5.5b
A SUPERCRITICAL TRANSPORT CONFIGURATON

THE COMPLETE CONFIGURATION

AIRFOIL SECTION OF THE WING

Fig. 2.5.6a
 SPAN WISE DISTRIBUTIONS OF SECTION LIFT
COEFFICIENT ON THE WING FOR A
SUPERCRITICAL TRANSPORT CONFIGURATION
0.7

0.6

0.5

0.4
, /p
CQ

CL
0.582 - INVISCID CALCULATION
0.484 --- VISCOUS CALCULATION. DI'LACEMENT

0.1 1°" 0.469 0 EXPERIMENTAL DATA. R. - 4.07 x 106

20 40 60 80 100

PERCENT SEMISPAN

PRESSURE DISTRIBUTIONS AT 50% SEMI-SPAN ON

A SUPERCRITICAL TRANSPORT CONFIGURATION

CL

0.625 INVISCID CALCULATION

0.518 VISCOUS CALCULATION. DISPLACEMENT
0.515 0 EXPERIMENTAL DATA, Re - 4.07 z 106

cp
DISPLACEMENT
THICKNESS WAKE
NOT MODELED

Ref. 2.16)
 z wìng-reference pLane SW

Fig. 2.6.1

WING PANEL (B)

BOOY PANEL (A)
WING SEGMENTS WAKE PANEL (C)

BODY SEGMENTS

RING

tr.f

INTERNAL
LIFT -CARRY -OVER
STRIP
STRIP

Fig. 2.6.2
 O experiments

NLR paneL method 200(body) +240(wing) panels

haLf-configuration
WLR AERO

/
CL "thick" wing theory
'thin" wing theory

10.8
¡0
/0
0.6

/ / onset of boundary
Layer separation
in experiment

_ L __L0
/
- 02 o
/
40 40

/
-8° 8° 12° 16°

0.2
o

-0.4

/0
-0.6

or

Fig. 2.6.3
 O exp ?riments

NL R panel method 200 (body) +240 (wing) panels

- - - NLR AERO

CM25

0.8

'thin' wing theory

0.2 "
"thick' wing theory
onset of sepaafion

-0.6 -0.4 -02 .- 0.2 0.4 0.6 0.8

0.2

0.4

Fig. 2.6.4
 flow direction

sin i= 1/M_
domain of L domain of influence
dependence

characteristic cone

Fig. 2.6.5

Nonanatytic Character of Supersonic Flows Renders Goal Numerics Much More

Difficult to Achieve. Methods have Lower Numerical Reliability.

x/C

Mach Wave
Panel Corners
Internal Wave Create Spurious
Reflection s Mach Waves Unless
Singularity Strengths
are Continuous

Fig. 2.6.6
 AIRFOIL : NACA 65A004
L=36.5, Dmax=333
s = 12.0 SUP. L.E.
Yref c=8.89 SUP. T.E.
Ct = 2.0
15.0

-® y/s=.9428
10.0
.7448
-ay', w
ç,,"
- .5469_,.-

5.0
15. 5925 49.40° ..349
-N,
- ------I---- ref
5.0 10.0 15.0 25.0 30.0 35.0 40.0

ALL DIMENSIONS IN INCHES

Swept wing-fuselage configuration at M = 2.01

Fig. 2.6.7a
 PRESENT METHOD

O UPPER SIDEÌ
EXP. REF. 1451
D LOWER SIDEI

Pressure distribution on the wing

Fig. 2.6.7b
 AIRFOILS FLAT PLATE
L -54.86 Dm = 3.048 /168 PANELS ON BODY
WING s = 597 c cr7.42 c=l .45 25 PANELS PER WING SURFACE
TAIL s = 4.62 Cr =7.09 Ct=O.O 20 PANELS PER TAIL SURFACE
Az rei
lo - lo

«:p= 00 (= 450

10 rei

YreI

WING
L_HINGE LINE WING
i/ill
&iia TAIL

6.86 F221
"la
lo
23.21
20
'II 30 40
MOMENT REFERENCE POINTx
re (=28.87
60

ALL DIMENSIONS IN INCHES

Body-wing-tail configuration. Panel scheme

Fig. 2.6.8a
 6pitch = angle of wing/fuse'age

0
-- pitch
=100
PRESENT
METHOD
O = Ø0 EXP. REF.1471
o =100

LO

CA.CA1' 001
A CACA«) = 0°'

.2 -
/ D

.2

'D
00
o 2 4
( .'l
D
_,zj
2
I

4
o
M 1,5

Normal force, pitching moment and axial force coefficients

for 4 = 00 with and without pitch control

Fig. 2.6.8b
Fig. 2.6.9

Fig. 2.6.10

Fig. 2.6.11
 "WAKE RELAXATION"

CONFIGURATION +1st WAKE ESTIMATE

/
COMPUTE SOLUTION

/ SINGULARITY STRENGTH (p)

I

COMPUTE WAKE STREAMLINES

I,
IMPROVED WAKE ESTIMATE
I

Fig. 2.6.12
 i

30
1'te
CONSTANT CHORD

Computed vortex wake behind AR = 5.0 swept wing,

C = 0.82

CL :0501
o CL :082
..________tALCuLA0OM

.0
-025
000000000 o
0 025 050 075 10 0
2 _-025 "

025
4 o
O
0000 00 oc o

Fig. 2.6.13
a) THREE_DIMENSIONAI VIEW

56
D cx' xijv*'cL 57
- vO'Sc

Pressure distribution for AR = 1.0

delta wing at = 20 deg. Comparison
with turbulent flow experiment

Fig. 2.6.14
 = O, Cp f(s)

a)

stagnation point
cusp
ç- stagnation points
/
b)

C open tait

c) d)

Fig. 2.7.1
 1,400

1.200

1.000

0.800

0.600
0.4 00

0.200

0.000

- 0.200
-0.400

-0.600
-0.800
-1.000 ITERATION STEP
- 1 .200

- 1.400
-1400 -1.000 -0.600 -0.200 0.200 0.600 1.000 1.400
-1.200 -0.800 -0.400 0.000 0.400 0.800 1.200

Reconstruction of circular cylinder (contour)

-3 -
ITERATION STEP
cp

-2 - Q TARGET

0-

-0.5 0.0 0.5 1.0

-1.0
Reconstruction of a circular cylinder (pressure distribution)

Fig. 2.7.2
 / Cp(s)-..èj
/ os
/ STARTING
/ GEOMETRY

SOLUTION OF
WIRICHLET)
BOUNDARY VALUE
PROBLEM

CORRECTION
ON
GEOMETRY

NEW GEOMETR'Y

Fig. 2.7.3
 I START GEOMETRY
/ C TARGET
/ WEIGHT FACTORS

NON PLANAR
PANEL METHOD

/ Cp PANEL

DETERMINATION
RES IDUAL

I
fi Cp TARGETCp PANEL
I

(LINEAR)
INVERSE METHOD

NEW GEOMETRY

FLOW DIAGRAM OF ITERATION PROCESS

Fig. 2.7.4
 lop ìsw and p.n,I .rrang.msnt of .x.mpl. win;

STARTING GeOu(TRY/PRESSURE DISTRIBUTION CL rn STARTING GEORETRY/PRESSURE DISTR.

- TARG'T P f FINAL RESULT OF FIG. S. I CL
55LIRE DISIIIBUTION CL
- - - TARGET PRESSURE DISTRIBUTION CL . .5
FINAL RESULT (4 ITERATIONS) CL
- FINAL RESULT I 3 ITERATION I CL

C, - C,
If. .347

q .075
q. 57S

q. .274
Ip. .274

q. .012
q. .012

.) HIGH WEIGHT QN TRAILING - EDGE ANGLE I0) b) EFFECT OF CONSTRAINTS ON THICKNESS AND TWIST
LOW WEIGHT 034 LEADING-EDGE RADIUS I2.5)
ZERO WEIGHT 034 THICKNESS AND TWIST

ReuIt Tor wnq alone (M = Q.?)

.1 WITHOUT CONSTRAINTS 034 THICKNESS AND TWIST

I WITH CONSTRAINTS OW THICEI4ESS AND TWIST

Front view of wing configuration Cz-;c;). enlarged Sx)

Fig. 2.7.5
 Top end front view and panel arrangement of wingbody configuration (body at
zero angle of attack)

STARTING GEOMETRY/PRESSURE DISTR. STARTING GEOMETRY IPRESSURE DISTa.

SAIE WING GEOMETRT AS RESULT 0F FIG. Sb) C .191 FINAL RESULT OF FIG. I.
LUE TI .51 I
TARGET PRESSURE DISTRIBUTION C
1-NETT
.5 ---- TARGET PRESSURE DISTRIBUTION
CLNETT
CLUE11 - .5
FINAL RESULT 3 ITERATIONS) C .011
I4ETT FINAL RESULT 2 ITERATIONS)
LNETT -

.847
'wET1 - .847
'NETT

NETT -
'tNETT

.271 '1HETT .274

1NETT

.012 11IETI .0)2

11HETT

b) EFFECT OF CONSTRAINTS AN THICKNESS AND TWIST

HIGH WEIGHT ON TRAILING - EDGE ANGLE 10I
HIGH WEIGHT OH TBAILIHG - EDGE RADIUS 2.ST
ZERO WEIGHT Ola THICKNESS ANO TWIST

Rrul(s for wnç-body confryraraIwn IM O.7f

WITHOu CONSTRAINTS RN THICKNESS ANO TWIST

bI WITh CONSTRAINTS Ow TIIICENESS AND TWIST

Front view of wings of wing-body configurafions.

Fig. 2.7.6