An Unsteady, Accelerated, High Order Panel Method with

Vortex Particle Wakes

by
David Joe Willis
Submitted to the Department of Aeronautics and Astronautics
in partial fulfillment of the requirements for the degree of
Doctor of Philosophy, Department of Aeronautics and Astronautics
at the
MASSACHUSETTS INSTITUTE OF TECHNOLOGY
June 2006
© Massachusetts Institute of Technology 2006. All rights r PcMASSACHUSETTS
prTm
ta VI
INSTnITE
OF TECHNOLOGY

JUL 10 2006

Authr... .. ...- ....

a... ..... ............. LIBRA PRIES
X Departen of Aeronautics and Astronautics ARCHIVES

'A'
i',-~Jaime
Pera
I fryai 1Afl AnK
IVAlLy , LVVU

Certified by

Professo, Department of Aeronautics and Astronautics

/7 - Thesis Supervisor

by..
Certified ..... .....................
Jacob K. White
Professor, Department of Electrical Engineering and Computer Science
,4 17 4n Thesis Supervisor
Certified by ............. ,.........................................
Mark Drela
Professt Depatmerit Rf Aeronautics and Astronautics
Committee Member
Acceptedby.....
Jaime Peraire
Chairman, Department Committee on Graduate Students
2
 An Unsteady, Accelerated, High Order Panel Method with Vortex
Particle Wakes
by
David Joe Willis

Submitted to the Department of Aeronautics and Astronautics

on May 30, 2006, in partial fulfillment of the
requirements for the degree of
Doctor of Philosophy, Department of Aeronautics and Astronautics

Abstract
Potential flow solvers for three dimensional aerodynamic analysis are commonly used in
industrial applications. The limitation on the number of discretization elements and the
user expertise and effort required to specify the wake location are two significant drawbacks
preventing the even more widespread use of these codes. These drawbacks are addressed
by the hands off, accelerated, unsteady, panel method with vortex particle wakes which is
described.
In the thesis, an unsteady vortex particle representation of the domain vorticity is cou-
pled to several boundary element method potential flow formulations. Source-doublet,
doublet-Neumann membrane (doublet lattice), and source-Neumann boundary integral equa-
tion formulations are implemented. A precorrected-FFT accelerated Krylov subspace iter-
ative solution technique is implemented to efficiently solve the boundary element method
linear system of equations. Similarly, a Fast Multipole Tree algorithm is used to acceler-
ate the vortex particle interactions. Additional simplification of the panel method setup is
achieved through the introduction of a body piercing wake discretization for lifting bodies
with thickness.
Linear basis functions on flat panel surface triangulations are implemented in the accel-
erated potential flow framework. The advantages of linear order basis functions outweigh
the increased complexity of the implementation when compared with traditional constant
collocation approaches. Panel integration approaches for the curved panel, double layer self
term are presented. A quadratic curved panel, quadratic basis function, Green's theorem
direct potential flow solver is presented.

Thesis Supervisor: Jaime Peraire

Title: Professor, Department of Aeronautics and Astronautics

Thesis Supervisor: Jacob K. White

Title: Professor, Department of Electrical Engineering and Computer Science

3
4
Acknowledgments

I would like to express my sincere gratitude to my thesis supervisors Professor Jaime

Peraire and Professor Jacob White. Their suggestions, guidance, feedback, discussions
and support throughout my time at MIT is deeply appreciated. I would like to thank Pro-
fessor Mark Drela, who provided significant insight and technical support. I would like to
thank the readers of the thesis, Professor David Darmofal, and Professor Kenneth Breuer
for their abundance of useful suggestions and comments. In addition, I would like to thank
Jean Sofronas and Chadwick Collins without whom nothing would be organized. I would
like to thank the National Sciences Foundation, The Singapore-MIT Alliance, and the Nat-
ural Sciences and Engineering Research Council of Canada for their support of this work.

My time at MIT has been influenced and enhanced by the many friends I have met. I
would like to thank them for the friendships we have formed. Ivan Skopovi for his encour-
agement and positive outlook. Jay Bardhan for the long hours spent discussing BEM in the
bowling alley office. Carlos Freire da Silva Pinto Coelho for his magic with C++ and claims
of bovine world dominance. Junghoon Lee for his matlab skills and for the valuable techni-
cal discussions. Sudeep K. Lahiri for keeping my one foot planted in the FDRL/ACDL lab
and for the many discussions both technical and otherwise. Christopher Sequiera for his
hard work in coupling his 6 D.O.F. dynamics engine to the flow solver. Zhenhai Zhu and
Ben Song for the first pFFT++ implementation which was used in early versions of the flow
solver. Xin Wang for helping me settle into the corner office. To Luke Theogarajan for his
discussions of Feynman and science in general. A sincere thank you to (Professor) Luca
Daniel, Nader Farzaneh, Xin Hu, Tom Klemas, Shihhsien Kuo, Dimitry Vasilyev, Michal
Riewenski, Homer Reid, Bradley Bond, Bo Kim, Per Olof Persson, Krzysztof Fidkowski,
Garrett Barter, Brian Taff, Anne Vithayathil, John Rockway, Lily Kim, Kin Sou, Du Ke,
Lei Zhang, Tarek Moselhy, Laura Proctor, Joe Kanapka, the triathlon club, and all the other
members of the MIT community I have interacted with.

Finally, I would like to thank those closest to me, without whom this thesis would
not be possible. To my wife, Kayla, I thank you for your patience, encouragement, love,
and support through all the ups and downs of the research. Your ability to make the hard

5
times bearable and the happy times even happier will be forever cherished. To my parents,
Owen and Tina and my two sisters, Natasha and Suzie(QT), who have taught me more
than education can teach me. For their love, and support throughout my life and academic
pursuits. Finally to the Riccio family, John, Maria and Gina who have become my second
family. Thank you all.

6
Contents

1 Introduction 17
1.1 Panel Methods: A Brief History ................ . . . 18
1.1.1 BIEs in the Pre-1960s. . . . 18

1.1.2 1960s - 1980s ...................... . . . . . . 18

1.1.3 1980s - 1990s . . . . . . . . . . . . . . . . . ... . . . 19

1.1.4 1990s - present ..................... . . . 20

1.2 Vortex Particle Methods: Background and Previous Work . . . . . . 21

1.2.1 Combined Panel Method - Vortex Method Approaches . . . 22

1.3 Challenges with Panel Methods ................. . . . 22

1.4 Thesis Outline .......................... . . . 23

2 The Governing Equations 25

2.0.1 The Domain. 25
2.1 The Governing Flow Equations ........................ 26
2.1.1 The Boundary Conditions. 26
2.1.2 Velocity Definition .......................... 27
2.1.3 The Scalar Potential Relationships ................ 28
2.1.4 The Vector Potential Relationships ................ 28
2.1.5 The Integral Equation Relationships for the Vorticity in the Domain 29
2.2 Boundary Integral Equations for the Potential Flow Equation . . . . . . . . 30
2.2.1 Derivation of the Green's theorem BIE ............... 30
2.2.2 BIE Formulations ........................... 32
2.3 The Pressure-Velocity Relationship ................ ... 40

7
 2.4 The Wing Trailing Edge Kutta Condition ................... 43
2.5 Conclusions .................................. 45

3 High Order, Curved Panel Integration 47

3.1 Boundary Element Discretization ....................... 47
3.1.1 Discussion of Orders of Approximation . . . . . . . . . . . . . . . 48
3.1.2 Orders of Approximation ....................... 50
3.1.3 Brief Review of Flat Panel Integration Strategies . . . . . . . . . . 52
3.2 Integration Approaches for Quadratic Basis Functions on Quadratic Curved
Panels. 52
3.3 Panel Integration Approaches for Non-Self Term Integrals ......... 53

3.4 Panel Integration: The Self Term Integrals .................. 55

3.4.1 Single Layer Integrals: The Self Term Integral ........... 55
3.4.2 Double Layer Integrals: The Self Term Integral ........... 61

3.4.3 Potential Flow Around the Unit Sphere ............... 67

3.5 Conclusion. 70

4 Wake Details 71

4.1 Body Piercing Wakes ............ .......... .72

4.1.1 Comments on Body-Piercing Wakes . . . . . . . . . . 77

4.2 The Vortex Particle Method ......... .......... .78

4.3 Conclusions ................. .......... .80

5 Implementation Details 81
5.1 Boundary Element Method Implementation . . . . . . . . . . ..... . 81
5.2 Wake Implementation details .................. ..... 82
5.2.1 Kutta Condition Implementation . . . . . . . . . . . . ..... 83
5.2.2 Doublet Sheet, Vortex Sheet and Vortex Particle Wakes ..... 85
5.2.3 The Farfield Approximation Model ......... ..... 85
5.3 Putting it all together ...................... ..... 87
5.3.1 Summary of Solution Steps .............. ..... 87
8
 5.4 Acceleration Algorithms . . . . . . . . . . . . . . . . . ........ 88

5.4.1 precorrected-FFT . . . . . . . . . . . . . . . . . ........ 89

5.4.2 Fast Multipole Tree Method ................... .. 90

5.4.3 A Brief Discussionof the AccelerationRoutines . . . . . . . ... 94
5.5 Implementation Details of the BIE Formulations . . . . . . . . . . ..... 95

5.6 Conclusions .................................. 96

6 Simulations and Results 99

6.1 Validation Computational Experiments . . . . . . . . . . . . . . . . . 99

6.1.1 NACA-0012 Non-Lifting Wing ................... 99

6.2 Steady NACA 0012 Lifting Wing Simulations ................ 100
6.3 Unsteady Computations . ........................... 101
6.3.1 Unsteady Startup Flow. ....................... 101
6.3.2
Unsteady
Finite
Aspect
Ratio
Heaving
Wing
............ 103
6.4 Wing-Body Simulation Example . ...................... 108
6.4.1 Rigid Body Piercing Wakes ..................... 108
6.4.2 Combined Body Piercing Wakes and Vortex Particles ........ 109
6.5 Example: Flapping Flight . . . . . . . . . . . . . . . . . ...... . 111

6.5.1 Flapping Wing Example ....................... 112

6.6Some
Comments
onSolution
Time
and
Accuracy
.............. 115
6.6.1 Doublet Neumann Formulation . . . . . . . . . . . . . . ..... 116

6.6.2 Source Doublet Formulation. . . . . . . . . . . . . 116

6.6.3 The Body Piercing wake Formulation ................ 117

6.6.4 Rigid Vs. morphing Bodies .................. 118
6.6.5 Comments on the use of Galerkin linear basis approaches . . . . .. 118
6.7 Conclusions
................... ............... 119

7 Conclusions 121
7.1 Contributions. . 121

7.1.1 The general panel method framework ............... . 121

7.1.2 Body Piercing Wake Formulation ................. . 122

9
 7.1.3 Higher Order Integration Approaches .............. 122
7.2 Future Work .................................. 122

10
List of Figures

2-1 The domain of interest includes all fluid external to the aircraft ....... 26
2-2 The fluid domain considered in the derivation of the Green's Theorem BIE . 32
2-3 A plot of the potential due to a unit strength wake . . . . . . . . . . ... 42

3-1 A comparison of the pressure around a sphere ......... . . . .. 50

3-2 An illustration of the inconsistency of using constant doublet representa-
tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...... 51

3-3 The relationship between the quadratically curved panel and the flat para-
metric
triangle
.. ......... . .............. 53
3-4 The 6-quadratic basis functions on the reference triangle. .......... 54
3-5 The orthonormal projection shown for a large selection of points on the
curved
panel
tothe
tangent
plane
............ ...... .... 56
3-6 The panel after it has been appropriately transformed ............ 58
3-7 The points which are used to compute the polynomial approximation of the
smooth integrand ............................... 59
3-8 The orthonormal projection shown for a large selection of points on the
curved panel . . . . . . . . . . . . ........ ....... 60

3-9 The radial projection of a flat panel onto the unit sphere centered at the
evaluation point, x. ............... .. . ........ . 63
3-10 The 6-subtriangles illustrated on the curved panel .............. 65
3-11 The radial projection of the points on the panel and the edges onto a unit
sphere centered at the evaluation point. . . . . . . . . . . . . . . . . . 66

3-12 An illustration of the solution X computed on a 32-panel sphere ....... 68

11
3-13 The error convergence of the potential .................... 69

4-1 A pictorial representation of the Green's Theorem formulation of the Bound-

ary Integral Equation ............................. 72
4-2 The effect of dipole wakes on the potential at a wake-body interface..... 73
4-3 Representing the potential wake as vortices. ................. 74
4-4 A demonstration of the body piercing wakes. . . . . . . . . . . . ... 75

4-5 A schematic of the sphere example ...................... 76

4-6 The total potential computed due to a body conforming wake, (o), and body
piercing wake, (*). . . . . . . . . . . . . . . . . ... . . . . . . . . . . 77

5-1 A schematic illustrating the conversion of dipole sheet panels to vortex

particles .................................... 86
5-2 The pFFT algorithm presented graphically . . . . . . . . . . . . . ... 91

5-3 The domain is divided into several concentric spherical domains ...... 93
5-4 A schematic illustrating the body piercing wakes .............. 96

6-1 The pressure distribution around the midsection of an Aspect Ratio 12 finite
wing with NACA 0012 profile at a zero angle of attack .......... 100

6-2 The pressure distribution around the midsection of an Aspect Ratio 12 finite
wing with NACA 0012 profile at a zero angle of attack .......... 101
6-3 The pressure distribution at the midspan of a Aspect Ratio 12 membrane
wing . ........................... 102
6-4 The time evolution of the Z-component coefficient of the force due to a
sudden startup of a series of finite wings ................... 103
6-5The
startup
flow
around
anaspect
ratio
8wing.
...... .. ....... 104
6-6 The Theodorsen lift deficiency function as a function of the reduced fre-
quency. The function is presented in terms of the magnitude (C(k)) and
the phase shift in radians. .......................... 105
6-7 The z-Component Coefficient of Force resulting from heaving oscillations
at a reduced frequency of . 5 .
.............
. . . . . . . . . . . . . . . . . . . . . . . .
107

12
6-8 The z-Component Coefficient of Force resulting from heaving oscillations
at a reduced frequency of . .... 108
6-9 The z-Component Coefficient of Force resulting from heaving oscillations
at a reduced frequency of . ................... 109
6-10 Plots of the vortex wake structures behind oscillating wings. . . . . . ... 110
6-11 The geometry considered for the wing-body example. . . . . . . . . .... 111
6-12 Results for the surface potential and surface pressure considering a rigid
body-piercing wake representation. . . . . . . . . . . . . . . . . . . 112

6-13 The geometryconsideredfor the wing-bodyexample. . . . . . . . . .... 113

6-14 A demonstration of the use of a wake only optimal vorticity distribution
result to construct an efficient three-dimensional flapping wing geometry
forsimulation
inanunsteady
panel
method
solver.
. . . . . . . . . ... 115

13
14
List of Tables

15
16
Chapter 1

Introduction

The solution of unsteady fluid dynamic flow around morphing bodies remains a challenge
for computational scientists. Several algorithms for handling this class of flows have been
proposed [1, 2, 3, 4, 5, 6, 7]; however, excessive computation time and power required
for moderate to high Reynolds Number flows as well as long setup times for complex
geometries limits the widespread use of the methods for practical design and analysis.
In this thesis a potential flow solver for simulating unsteady potential flows is presented.
Although unsteady potential flow approaches have been proposed in the past [35, 19, 21, 13,
14], the current method has a collection of distinct advantages for the solution of unsteady
flow around morphing bodies:
* Wake Body Intersections: A formulation which allows for body-piercing wakes is
presented, thereby relaxing the strict requirements of body conforming wakes.
* Particle Wake Representation: The approach makes use of a time dependent Vortex
Particle Method (VPM)[54, 57, 88] to represent domain vorticity rather than a dipole wake
sheet [21].

* Hands-off simulation of complex geometries: Due to the combined use of a vortex

particle wake representation and the implementation of body piercing wakes, hands off
simulation of potential flow around complex geometries is possible.
* Fast Acceleration Approaches: The precorrected-FFT algorithm (p-FFT)[84], and the
Fast Multipole Tree Method (FMT) [110, 70] are used to reduce the solution computational
complexity to O(N log(N)) (where N represents the number of unknowns in the system).

17
 * Steady State Wake Location: The method has advantages for steady applications. In
particular, it does not require a wake geometry to be prescribed by the user.
* High order basisfunctions: New integration methods for the curved panel dipole are
presented. In addition, a method for computing the single layer self term integral over
curved surfaces similar to previous work in the field is implemented.

1.1 Panel Methods: A Brief History

Aerodynamics Panel Methods were first investigated in the late 1950s. Since their initial
development, they have been instrumental in the design, optimization and analysis of air-
craft and aerodynamic bodies [15, 16, 107, 18, 22, 23]. A brief outline of some salient
history of panel methods is presented:

1.1.1 BIEs in the Pre-1960s

Prior to the use of digital computers, basic analytical solutions to the potential flow Bound-
ary Integral Equations were employed [24, 25, 26]. The principle of linear superposition
of fundamental solutions such as point sources, and point doublets was used regularly to
solve potential problems[10]. The field of panel methods was born in 1958, when Smith
and Pierce from Douglas Aircraft Company used a discrete form of the boundary integral
equations to solve for the potential flow around bodies of revolution [29].

1.1.2 1960s - 1980s

With the success of the initial panel methods, the Smith group received support to con-
tinue development of panel methods for both two and three dimensional flow [78]. They
pioneered the panel method solution to the lifting body problem in 2-Dimensions [13] and
in 3-Dimensions [12]. The development continued to include higher order discretizations
of the BEM approach in 2-Dimensions[ll]. The Douglas group panel methods were al-
most exclusively of Neumann type, using either source or vorticity distributions over the
surface [10]. In the 1970s, the Green's Theorem perturbation potential based Dirichlet

18
problem was introduced by Morino [31]. There were also several variations of different
complexity of the surface singularity boundary element method/membrane lattice approach
[30, 121, 123].

The early panel methods were limited by computer memory and processing power.
Some alleviation of computational complexity was achieved by using multipole expansions
in place of analytical expressions for panel integral expressions for farfield evaluations;
however, the methods still required the solution of a dense linear system.

1.1.3 1980s - 1990s

During the 1980s, several low order three dimensional panel methods were developed
[34, 35, 107]. In addition to the low order methods (low order here referring to the constant
basis function approximation of the solution), several high order implementations were also
developed. These high order methods were developed for the benefits of increased solution
accuracy as well as for satisfying the solution continuity requirements imposed by super-
sonic flow applications. A combined Boeing and NASA effort resulted in PANAIR/A502,
a quadratic basis, flat-sub-element high order panel method [116, 36]. Additionally, HISSS
[37] a panel method based on PANAIR was developed. In the late 1980s PMARC [106]
was developed at NASA-Ames Research Center and was later released as a controlled
access computer program. Although the 1980s brought with them great advances in com-
putational power, limitations on computational time and memory still prevented large-scale
panel method solutions. Solutions with several thousand panels were routinely performed
on large computers; however, due to the coarseness of surface discretizations, limitations
on the practical use of panel methods existed. In addition to developments in three dimen-
sional solvers, two dimensional panel methods were being developed and used heavily for
inverse airfoil design [38, 27, 39, 28]. Furthermore, the use of boundary layer coupling was
investigated for incorporating viscous effects[38, 39].
In the 1980s several algorithmic developments were also made which have had a signif-
icant impact on the development of panel methods. These developments included iterative
solution methods, most notably for this thesis the development of Krylov subspace iterative

19
solvers [100, 42, 43]. In addition to iterative solvers, several sparsification and acceleration
routines were also developed in order to facilitate the rapid computation of matrix vector
products of the dense BEM linear systems. The first category of fast methods involved
multipole expansion approximations of the farfield influences [110, 70, 86]. The second
category of methods relied on rapidly approximating farfield interactions using a Cartesian
background mesh[102, 50, 68].

1.1.4 1990s- present

By the 1990s panel methods had largely given way to higher fidelity Navier-Stokes and
Euler solvers [105, 51]. Although Eulerian reference frame domain solvers were being
heavily investigated, several Lagrangian based approaches were developed. The Vortex
Particle Method was refined and further investigated for the simulation of largely vortical
flow [54, 59, 56, 57, 60, 67]. The section which follows describes some of the history and
development of vortex particle methods. Despite the promise of Navier-Stokes solvers,
accurate viscous drag prediction remained an elusive task. In the 1990s several researchers
started to consider the problem of 3-Dimensional Integral Boundary Layer Methods[40, 41]
with some success.

The Fast Multipole Method [110] was used and further developed in practical boundary
element method solvers for many diverse disciplines [44, 49, 46]. In the early 1990s, the
precorrected-FFT algorithm [84] was developed. The precorrected-FFT approach provided
a kernel independent framework for the acceleration of BEM and N-body problems[45, 83].
In the 1990s and 2000s several panel method codes continued the advancement of higher
order approximations to the boundary integral equations[1 15, 117, 118, 119, 120, 123, 124,
125, 111], however, due to the complexity involved with higher order methods and the lack
of robust and efficient integration techniques for higher order approaches, their adoption in
the BEM community is limited in comparison with the popular constant collocation type
approaches.

20
1.2 Vortex Particle Methods: Background and Previous
Work

Like the panel method, the vortex particle method (VPM) has a rich and long history which
is presented in several references [54, 57, 88]. Only some of the brief points of interest are
outlined in this section. The framework for Vortex Particle Methods originated well before
the digital computer age. In 1930, Rosenhead performed a dynamical vortex calculation
using singular point vortices [53]. The vortex particle method in 3-Dimensions was largely
derived from the 2-Dimensional method. Initial investigations in 3-Dimensions used vortex
filament approximations to account for the domain vorticity [58, 55]. The difficulty with the
vortex filament and sheet methods was the need to store the connectivity information. Beale
and Majda [57] proposed the first vortex blob method in which the vortex blob positions and
strengths were updated in a Lagrangian manner with no connectivity between the blobs.
Several advancements have been made to the vortex particle method since its inception
which can be found in the following works [54, 57, 59].

The use of vortex particles, vortex filaments as well as vortex sheets has been explored
for representation of domain vorticity in many diverse applications [60, 54, 59, 65, 67, 68].
Many current vortex particle methods are used for the simulation of turbulent flows. The
method has been found to be quite useful for the prediction of jet flows [66], turbulent
bluff-body wake flows [59], internal flows [67], etc.

Much like panel methods, the vortex particle approaches require the evaluation of an
N-body interaction problem at each step of the simulation. Since the particles advect in a
Lagrangian manner, their positions change at each timestep, requiring the re-computation
of inter-particle influences. With the development of O(Nlog(N)) algorithms such as the
Fast Multipole Tree [110], and particle in cell [102] methods, the particle based approaches
became viable for the large number of particles required to adequately resolve the flow
physics.

21
1.2.1 Combined Panel Method - Vortex Method Approaches

Vorticity representation in panel method aerodynamics has been traditionally limited to

vorticity sheet and vortex filament approaches [113, 107, 10]. Although vortex particles
are used commonly for representing vorticity in unsteady 2-Dimensional aerodynamics
computations [61, 113, 65], the use of vortex particle methods to model lifting body vortex
wake sheets in 3-Dimensional aerodynamic panel methods has been limited to a small
number of researchers [62, 63, 64]. The previous use of the combined 3-Dimensional
panel method and vortex particle approach has been limited to applications in which a
lifting surface or series of lifting surfaces is considered isolated from the body (sails [64]
and wind turbines [62]). Combined panel method and vortex particle approaches require a
sufficiently high density of vortex particles to practically and accurately model the vortex
sheet; therefore, when acceleration methods are not considered, these approaches suffer
from large computational times and low accuracy. The vortex particle approach has also
been coupled to several boundary element potential flow methods in disciplines other than
aerodynamics; however, this coupling often involves approaches for modeling the near wall
viscous effects which are beyond the scope of this thesis [67, 69, 59].

1.3 Challenges with Panel Methods

Despite the comprehensive development of panel methods over the past several decades,
there are several drawbacks to existing methods which hinder their more common use,
namely:

* Stringent requirements for wake discretizations: There are two primary drawbacks
with current wake discretizations.

1. Due to the necessity to impose a potential jump in the wake region, there is a need to
have a body conforming wake surface mesh. This poses significant challenges in the
problem setup and meshing.

2. For unsteady problems, the development and accurate advection of the wake is of-
ten compromised when wing-body simulations are considered. Common problems

22
 include: intersection of the wake with downstream surfaces, the difficulty in wake-
body intersection meshing, and the time and effort required to compute full unsteady
flows.

By providing a panel method framework which incorporates an automatic wake generation

and advection approach, the difficulties associated with user expertise and user interference
will be reduced.
* Memory and Processor Imposed Limitations: Physical computer memory constraints
as well as unrealistic solution times hinder the applicability of panel methods to practical
problems. This limit on the number of elements prevents the interfacing of panel methods to
traditional CAD and CFD grid generation tools. In addition to the panel count limitations,
specific surface panel types (such as quadrilateral) often prevent the use of CAD and CFD
compatible discretizations. By ensuring a panel method can use the same surface grid
as a CFD tool will reduce the significant preprocessing workload for aircraft designers.
Often, preprocessing such as surface grid generation and geometry description is more
time consuming than solutions.
* Discrete Approximation to the Continuous Problem : Most panel methods in use
consider low order approximations for both the geometric discretization as well as the
solution basis representation. In particular, a common approach is the constant collocation
approximation on flat polygonal elements. Although the approach works well for many
simulations, it is shown in Chapter 3 that the method has several drawbacks. Increased
accuracy and faster convergence rates are a direct consequence of appropriately designed
higher order methods.

1.4 Thesis Outline

In the following chapters of this thesis, solutions to each of the drawbacks of the tradi-
tional panel method are presented. The resulting solution framework which is implemented
provides rapid simulations of potential flow simulations which are accurate, and easy to
compute. Chapter 2 outlines the applicable theory and Boundary Integral Equation formu-
lations considered. In chapter 3, a novel integration approach for high order curved panel

23
integrals is presented for the double layer potential (the single layer potential approach is
also demonstrated and closely resembles traditional approaches). In chapter 4 approaches
for handling the wake vorticity are presented. In chapter 5, the details of the implementa-
tion of the panel method is presented. Finally in chapter 6 and chapter 7 validation, results
and conclusions are presented.

24
Chapter 2

The Governing Equations

In this chapter the governing fluid dynamics equations are presented.

2.0.1 The Domain

Consider the domain illustrated in Fig. 2-1. A point position, R(X, Y, Z, t), in space at a
given time, defined in the fixed-in-space global reference frame is:

Rp =RG + Gp, (2.1)

where RG is the position of the body frame origin, and rp is the position of point p relative
to the local body frame. A given point on the body will have a velocity with respect to the
global reference frame given by:

VP = VG + VG + (Q X rGp), (2.2)

where VGrepresents the velocity of the body frame origin in global coordinates, Q is the
angular velocity, and VGPrepresents the relative motion of the surface due to deformation
of the body (eg. deflection of a control surface). For clarity, body velocities are denoted by
V and fluid velocities by U.

25
 z

Figure 2-1: The domain of interest includes all fluid external to the aircraft surface. The
particles trailing the wing section, the vertical tail and horizontal stabilizer represent regions
in which vorticity exists due to the lifting surface trailing shear layer.

2.1 The Governing Flow Equations

In the paragraphs which follow, the governing equations and assumptions are presented.

2.1.1 The Boundary Conditions

At any point on a solid surface in the domain, a no penetrating flux boundary condition

(flow tangency condition) is given by:

(2.3)

where, n is the outward unit normal vector on the body at a given point R on the body
surface. The rate at which the perturbations in velocity decay with distance from a non-

lifting body is:

R~oo
.......... (1IIRII3 )
lim U(R, t) ::; 0 - .....
- . (2.4)

26
Similarly, for a lifting body, the velocity radiation condition is:

lim (R, t) < ( 2 (2.5)

The slower decay of the velocity in a lifting body simulation is due to the presence of the
trailing vorticity in the wakes.

2.1.2 Velocity Definition

The flow is assumed to be inviscid, incompressible and have constant density. Any vorticity
in the domain is localized on the thin wake regions trailing the lifting surfaces. Otherwise,
the flow is assumed to be irrotational. These assumptions greatly simplify the form of
the governing equations. The fluid velocity, U(R, t) satisfies the following two farfield
conditions:
V. U(oo,t) = 0, (2.6)

by the continuity requirement, and:

V x U(oo, t) = 0, (2.7)

by the farfield velocity decay boundary condition. As a consequence of the above fluid
velocity field properties, the fluid velocity, U(R, t) at a given point in the domain can be
expressed as the superposition of a scalar potential component, U(R, t), and a solenoidal
vector potential component, U (, t), using a Helmholtz decomposition [77]:

U(R, t) = O(, t) + V(F, t) = + V x P. (2.8)

The scalar potential component of the velocity is irrotational and any rotational effects
are captured in the vector potential component. The decomposition of the velocity field
into a scalar potential and a vector potential component is not commonly considered in
most panel method implementations. By performing the decomposition, the traditional
scalar potential boundary element method formulations can still be used; however, with

27
the addition of the vector potential, wake vorticity can be modeled directly using vortex
distributions (volumes, sheets or points) rather than indirectly through the use of dipole
sheets.

In terms of the vector and scalar potentials, the boundary condition of equation 2.3 is:

ii (V + V x T) = (VG+ VGp+ X rap). (2.9)

Note that due to the inviscid flow assumption, only the normal velocity boundary condition
is applied at the body surface.

2.1.3 The Scalar Potential Relationships

The governing continuity equation for a constant density fluid is expressed in differential
form as:

V U=0,

Substituting the velocity as defined in equation 2.8 into the continuity equation the resulting
mass conservation equation is:

V. (VO + V x ) = V. (V) = V2 = 0. (2.10)

Which is the Laplace's equation for the scalar potential.

2.1.4 The Vector Potential Relationships

The Vector Potential - Vorticity Relationship

The vorticity in the domain, wc(R, t), is defined as the curl of the velocity [8]:

V x U= (2.11)

The velocity component due to the Vector Potential, I, is:

Vx = UT (2.12)

28
Substituting the vector potential relationship (equation 2.12) into the definition of vorticity
(equation 2.11), and choosing the vector potential to be a solenoidal vector field (V -· = 0),
results in:

V2~ = -a, (2.13)

which is a vector Poisson equation relating the vector potential to the vorticity.
The Vorticity Evolution Equation
The vorticity evolution equation is derived from the incompressible Euler equations [8],

au- Vp
+ U-VU= -- (2.14)
at p

where, p, is the fluid density, and p is the pressure. Taking the curl of eqn. 2.14, the
resulting equation for the vorticity evolution in the domain is [8],

Dt4
Dt = at + U V = VU (2.15)

where the term ac. VU on the right hand side represents the vorticity stretching (or how the
magnitude and direction of the vorticity changes as it is exposed to velocity gradients in the
flow field). The left hand side of the equation is simply the total derivative of the vorticity
with respect to time.

2.1.5 The Integral Equation Relationships for the Vorticity in the Do-
main

The vector Poisson equation governs the velocity vector potential (see equation 2.13). In
integral form, the vector potential due to the vorticity in the domain is:

I(i4,t) = 4j _ dV', (2.16)

29
where, V is the volume of the domain in which the vorticity, Ac,exists. The vorticity induced
velocity is determined by taking the curl of equation 2.16:

u = V x (r,t)= V x I 1 .- , 11 dV'. (2.17)

The resulting expression in equation 2.17 for vortex filaments is the familiar Biot-Savart
law [113, 8]. Similarly, the associated component of the gradient of the velocity term used
for the vorticity stretching in the vorticity evolution equation is determined by taking the
gradient of equation 2.17.

2.2 Boundary Integral Equations for the Potential Flow

Equation

The derivation of the potential flow boundary integral equations is briefly presented in this
section. First, a general derivation is presented, following which is a closer examination of
the integral equation formulations.

2.2.1 Derivation of the Green's theorem BIE

The following derivation is similar to that presented in [94, 127, 95, 96, 113]. Consider the
potential governed by Laplace's equation:

V2 0(rl = o. (2.18)

A particular fundamental solution (Green's function) to Laplace's equation is:

G(r-, i)= i 1=11" (2.19)

30
The following statements are a result of integrating by parts (for r # rf') [25]:

l]/ dV
|j[ 1V~ (V ir]
j=I .)v'+
I[|| _ ihln. V dS, (2.20)

and, similarly:

-lII l i dS, (2.21)

where, ' is the integration variable position on the boundary surface, S represents the
domain boundary surfaces under consideration, and n, represents the boundary surface
normal (directed into the fluid domain) at the integration variable position. Combining
equations 2.20 and 2.21, results in:

[()v ] dV= [ 1i dV
V -V ( _ )]n dS' (2.22)
W+SB+S+SSphere
[(l- l' - (ir ) dS'. (2.22)
Figure 2-2 illustrates the boundaries of the fluid domain as well as the spherical exclusion
surrounding the point f' = . The volume integrals in equation 2.22 are identically zero.
Hence the integral equation reduces to a boundary integral equation:

Sw +SB+So [ ,.
_+SSphere - ) V i, dS' = 0. (2.23)

Integration over the farfield boundary reduces to zero due to the radiation boundary con-
dition. The integration over the surface of the spherical exclusion region simplifies to
[113, 95]:

[V (f,)
LsSphere -(ji1r ,)V ] *.dS' =-4r(r. (2.24)
31
 Sinjinity

@S'iphere

SilljillilY Sill/illil)'

SB

Figure 2-2: The fluid domain considered in the derivation of the Green's Theorem BIE

The result is the Greens Theorem boundary integral equation for describing the potential at
an evaluation point r inside of the fluid domain:

2.2.2 BIE Formulations

In this work, the following Boundary Integral Equation formulations are considered:

1. The Greens Theorem Formulation: The traditional Greens theorem formulation [113,
95, 127] involves solving for the surface potential on the body boundary by specify-
ing the normal derivative according to the boundary conditions.

2. The Source-Doublet Formulation: The source doublet formulation [31, 113] is a com-
mon panel method formulation consisting of an interior and an exterior domain over
which potentials are defined. The jump in the potential is prescribed using doublets
and the jump in the normal derivative of the potential is prescribed using surface
sources.

32
 3. The Source-Neumann Formulation: The source-Neumann formulation [13] satisfies
the Neumann boundary condition directly. The potential and velocity in the domain
are computed using source singularities. In traditional source-Neumann formula-
tions, lift cannot be modeled; however, in this thesis a Source-Neumann method is
used with a wake piercing formulation to solve for the flow around lifting bodies.

4. The Doublet-Neumann Formulation: The doublet-Neumann formulation [113] con-

sidered is equivalent to a doublet lattice method. Doublets are placed onto an in-
finitely thin surface and their strengths are adjusted in order to satisfy the Neumann
boundary condition.

In the following paragraphs the above methods are presented.

The Direct Green's Theorem BIE

The Green's Theorem boundary integral equation for computing the potential at a point,
(rf, in the fluid domain due to a non-lifting body is [95, 127, 113]:

(r = 4- ()I '- dSB 4 rSB - a (2.26)

The lifting-body potential flow problem requires the addition of a potential jump in the
trailing wake region to account for the physical vorticity[1 13]. In equation 2.25 the fol-
lowing Green's Theorem boundary integral equation expression resulted when a wake was
included in the domain:

( 47 ISB n [r - fi] dS- 47r +S- [- i'l dSB+w. (2.27)

In this case the potential, 0, includes a component due to the body, 5B, and a component
due to the wake, bw.

A slightly different expression for the lifting body potential arises from a superposition
principle, in which the potential due to a dipole sheet trailing all lifting surfaces is added to

33
the Green's Theorem representation of the non-lifting potential flow (equation 2.26):

- 4
l q ,
aB )
I
1
dS - -ll
B
OB W)
_s 1
l dSBl 4s
(2.28)
1
Where w represents the added domain potential jump influence due to the wake. The
potential q(r) represents a contribution from both the body and the wake:

(r =qB(f)
+ w(r-- (2.29)
Notice, in equation 2.28, that the unknown is the body surface potential (B). The wake
potential, 1 w can be expressed as a function of the prescribed wake potential jump, uw as:

+s
1= [1 dSwv (2.30)
Iw(r')()
In order to uniquely solve equation 2.28 for the potential, O(f'), the following steps are
taken:

1. Prescribe a wing trailing edge Kutta condition which relates the unknown wake
strength(s) w to the surface potential or its derivatives.

2. Specify the normal derivative of the potential at the boundary. Equation 2.3 is used:

i VB(R, t) = - (VG+ VGp+ X rGp-V X (R,t)-V w() ), (2.31)

In order to solve for the body potential B(i, the problem must be reduced to one in
which the influence of all of the other potentials comprising the solution become boundary
conditions.

Depending on the wake representation, either of the expressions (equation 2.25 or equa-
tion 2.28) for a lifting body problem may be used. In this work, the wake is often considered
as a potential influence close to the body (in order to satisfy the Kutta condition) and as a
velocity influence further away from the body (using a vorticity representation).

34
The Indirect Source-Doublet BIE

The source-doublet formulation is commonly presented in place of the Green's theorem

formulation in many aerodynamics texts [113, 20]. The source-doublet approach considers
the two domains which are separated by the surface of the body. The first region is the
domain interior to the aerodynamic body under consideration. The second region which is
considered is the domain exterior to the body boundary surface (which includes the fluid of
interest). Inside of each of these domains a potential can be defined.

Setting up a potential flow inside of the body yields the following interior potential
Oi(r) expression for an evaluation point inside of the body:

i1 isB
and() i dS-
B

4i | sBir a i[ ] dSb (2.32)

Similarly, for an evaluation point in the region exterior to the body (inside the fluid of
interest), the influence of the internal potential is:

47rs
anf,(v) Il f dS-

47r#SB()i [2ldSB,'
one, I-r ll
(2.33)

For a point exterior to the body, the potential O(f) is described by:

W_ a 1
47r an llr-r, I-1
Adding the interior potential expression (equation 2.33) to the expression in equation 2.34
(while appropriately adjusting the surface normal definition to point into the fluid domain

35
(exterior to the body) ), the following integral equation results:

r =4 VJsI() dS-

1rSB
47 ( ) ( 5n9)) [l2r fil] dSB +4w, (2.35)

By compactly representing the surface singularity distributions, the following equation re-
sults:

W%J
if L
=
d-
p~ [~r
¢11
1
Iu(J -
1
'SB
s fa(S)B
If -
dSB +,Dw. (2.36)

The doublet, represents the interface potential jump:

(r = (r - -e i(ir, (2.37)

where the superscripts e and i represent the limit as the point r is approached from the exte-

7(-) (2.
rior domain, and the interior domain respectively. The source, a represents the magnitude
of the interface discontinuity in the normal derivative:

( )

The application of boundary conditions to the source-doublet formulation is similar to that

in the Green's Theorem formulation; however, care must be taken to ensure that the inte-
rior and exterior problems are appropriately considered. A benefit of the source-doublet
method is that the interior potential can be carefully chosen in order to simplify the par-
ticular problem being solved. In the following paragraphs a brief discussion of two of the
many possible interior potential choices is presented. The two formulations represent a
subtle difference between the traditional source-doublet method and a useful variation of
the source-doublet method which is investigated in this thesis:
Internal Potential Case 1: In the first example of the source-doublet approach, a zero

36
internal potential is examined:

OV)47T f 1 1
AWBdS'
ani jIIII
u(ri') dSB + w. (2.39)
47 JsB

In this case the source strength has a value:

(
5(F)an ) = nf
((G+ G,+ x rp - V X (R, t)) , (2.40)

and, similarly, the doublet strength can be written as:

A(i = ¢~(r) - oi(f)

= 0(), (2.41)

By inserting these values into the source-doublet BIE, the result is:

/( ) ISB () l dS -
T1--
On, WI I
1JSJSB
4?r
I i~. (" +-VGp+ X tCp - V x (, t)) _ Pl dS (2.42)

Dw.

hence,

B() =147/ (r) _ dSB

- -

47FJSB
hp,.VG+ VGPQ Cp-VX (R,t)) 1 dSB
+ X ,11 (2.43)

Solving equation 2.43 gives the potential (r) (which is also the doublet strength) at the
surface of the body. This expression is similar to equation 2.25. The potential and velocity
at points in the domain is a simple evaluation of the BIE in equation 2.39.
Internal Potential Case 2: In the second example of the source-doublet approach, an
interior potential corresponding to the wake induced potential is examined:

1f
4I c(i') dSB + w (2.44)
1SS)4/3 14'r )nrB
+s On I lf 1 ;fl dSB 471 JSB II- ?II

37
where, the source strength has a value:

r ane)
= . (VG + VG + X rpG - V x TR, t)-Vw(r1), (2.45)

and, similarly, the doublet strength can be written as:

u(r = i = e(r-' w()

(r- - i(r wB(r,= (2.46)

By inserting these source and doublet values into the source-doublet BIE, the result is:

-)I
(~=47 B(ranf llr- ilB + __w-
_dS
47r , (VG.+ VG + Q X rGP-V X (R,t)-V () 1 fll dS
(2.47)

hence,

(r =47r (
1/~B~r
an llr--14FJSB
l dS$-
IB-

+P + xrc-x (R,t)- (fi I *'d-

(2.48)

Which is similar to the expression presented in equation 2.28.

In the first internal potential case the boundary integral equation absorbed the wake
potential into the surface singularity strength representations; however, in the second ex-
ample, the wake was treated as a superimposed potential field and as such the potential
on the body. In the second case, one effectively sets up an internal potential which cor-
responds to the wake potential. The difference between the two formulations is subtle.
For the work presented in this thesis, the ability to treat the wake using a velocity and/or
potential representation is useful.

38
The Indirect Neumann Source BIE

The remaining two BIE formulations which are considered impose Neumann boundary
conditions. The Neumann-Source formulation potential is [113]:

4Ss() T1 dS, (2.49)

where, o-(') is the boundary surface distribution of fluid source strength. To add generality,
an arbitrarily defined potential 'Tw can be added:

q(r)= OB(r + = J (/) lr- d

11 W. (2.50)

The gradient of the equation is taken to yield a boundary integral equation for the velocity
in the domain:

=V 4r|
V+(z) C(d)
1' dS+ (2.51)
When the boundary conditions are applied to the resulting equation is:

f*(VG
+VG~p+Qx
rGp-VXV =i *( ( f|1u(') fi dS + w) . (2.52)
Solving for the strength of the source distribution a(x'), one can back substitute into equa-
tion 2.49 or 2.50 to obtain a relationship for the potential in the domain.

The Indirect Neumann-Doublet Membrane BIE

The Indirect Neumann-Doublet BIE is an identical formulation to that used in a doublet

lattice type code[113]. The potential is defined as:

(M = 47rI )n [r-'I dSB, (2.53)

where, u(') is the boundary surface distribution of dipole or doublet strength. In order
to once again add generality, one can explicitly incorporate the superposition of any other

39
solutions of Laplace's equation:

X(/-) = 4B(r)
B )11
+W = 47
1
I (0)n
f
I(-
_I

dSB + w (2.54)

The boundary conditions given in equation 2.3 are applied to give:

n, *(VG+ VGp+
Q X rip -V X ~) = L, V (4|4.. (- (_ l) dS+ w)-
(2.55)
The indirect Neumann-Doublet approach is particularly useful for lifting membrane appli-
cations. For example, thin fabric or membrane structures are difficult to model well with a
source-doublet type approach due to the need for panel elements which are of comparable
size to the thickness of the membrane. Instead, the Neumann doublet membrane formu-
lation approximates the thin surface using a single sheet of doublet singularities which do
not have strict limiting constraints as the Dirichlet approach.

2.3 The Pressure-Velocity Relationship

The Bernoulli Equation is used to determine the forces and pressures on the body. Since a
potential-vorticity approach is used, the applicable unsteady Bernoulli equation is derived[ 113].
The incompressible Euler equations are:

ac - Vp
+ U. V =U
- (2.56)
ot p

If those regions of the flow which have zero vorticity (all of space excluding the trailing
vortex wake region) are considered, the resulting equation is:

+
1VIdl2
_0 VP (2.57)
at 2 P

The definition of the velocity, given by equation 2.8, can be substituted into equation 2.57
resulting in:
o(V + V x ) 2
2 Vp
atV+2) VIV +l+ V (2.58)
dt 2 ~~P

40
Collecting like terms, and re-arranging:

a(Vxt) + (2.59)
At At + VIV+V x +V
12 (P) =0.

Integrating eqn. 2.59 along the streamline from surface point xil, to a farfield reference
point at oo (where at R = oo, Ut=o to = 0, and p = po) results in:

a(V x I) 1V +V X x ) poo - px (2.60)

at dC+ at 2 xl P

Furthermore, note that the term 90in equation 2.60 is defined in an Eulerian reference
frame. The change in potential with respect to time for a point on the body surface can be
computed by converting to a body Lagrangian reference frame:

Eulerian body -
I= (VG + VGP+ ( X rGP)) V (2.61)
at

The overall Unsteady Bernoulli equation is therefore:

Poo- Px = fPX o( x )
dC + JIbody-
POO-Pa
l (o at at
(VG + VGp+(Qx rGp)) V + 22IVq + V X j12.
.
(2.62)

The unsteady term due to the domain vorticity:

V .d,
1xa(t (2.63)

is difficult to handle in the form written above. By considering the contribution of the
vortex wake as an analogous contribution due to a dipole sheet, one can write:

0l a(v x ) dC =
IPX1 PXl 1ao
Iwake '
d
t. (2.64)
f0 at

41
Integrating the expression for the wake potential cp is simple:

(2.65)

where, \7 cp is simply the velocity due to the wake. The overall Unsteady Bernoulli equation
is therefore:

The above expression for the change in potential due to the wake can be examined from an
order of magnitude argument. For finite wing-only simulations the term is small (see figure
2-3 for a pictorial argument). Although the contribution to the pressure from the vortex

positive
potential,
wake
infuence

r
1
negative
potential,
wake
influence

Figure 2-3: A plot of the potential due to a unit strength wake trailing behind an airfoil.
The wake extends 50 chord lengths behind the airfoil and is not shown in entirety. Notice
that the potential contribution to the airfoil is nearly zero due to the airfoil lying nearly
in plane with the wake sheet. For airfoils, the change in potential due to the wake in an
unsteady computation is a perturbation of an already sInall potential influence. In a scaling
argument, one can argue that the change in potential due to the wake is negligible for many
simulations (especially for thin airfoils).

wakes in many cases is small, for certain cases the order of magnitude of the unsteady

42
wake potential contribution is similar to the other contributions in the pressure calculation.
This tends to occur when the wake passes above or below a downstream surface. In these
cases, one can compute the unsteady pressure contribution of the domain vorticity on the
body by solving the linear BIE system for the interior to the body wake potential:

01i 1 - 1 _ _ 1
(P(r-=4 j a( _|1 1- p( _1 adSb, (2.67)

where, the potential p represents the potential due to the wakes at the body surface. The
normal derivative of the wake potential is known at the surface by the velocity expression
for the vorticity in the domain. By solving the boundary integral equation, the potential
contribution due to the wakes on the surface of the body is determined. As a result, the
unsteady Bernoulli's equation for the pressure at the surface of the body is:

Poo- PX1 _' a

p t
= Ibody
+ ActIbody-

(VG + VGp + ( X rTGp)) [V + V X (P]+ 2IV

2
+ V X i2. (2.68)

From this pressure-velocity relationship, the forces can be computed by integration:

F(t) = j (pOO
- p(t))fdS'(t). (2.69)

Similarly, the moments about the body reference frame origin can also be computed by
integration:

M(t) = i r/GpX [(Po - px,(t)) hi]dS'(t). (2.70)

2.4 The Wing Trailing Edge Kutta Condition

Most surfaces with sharp geometric cusps will produce forces when the body moves relative
to a fluid. This net force is due to the viscous nature of fluids which prevents flows from
traveling around these cusps. Rather, the flow will tend to separate from a sharp corner
and shed a trailing shear wake. Since the flow is assumed to be inviscid, the production of
forces requires an additional condition on the governing equations. In this thesis a Kutta

43
condition [9] is applied at all wing trailing edges. The Kutta condition prescribes a finite
surface potential jump discontinuity across the geometric cusp, thereby resulting in smooth
and finite trailing edge velocities. In order to assure the potential jump at the trailing edge,
a fictitious wake surface is introduced into the domain. The wake surface is considered as
a cut in the fluid domain capable of sustaining a potential discontinuity. Both a linear and
non-linear Kutta condition have been considered in different formulations presented in this
thesis.

The non-linear Kutta condition requires that there be no pressure jump across the trail-
ing edge [72, 73, 74]:

Pupper - Plower = 0, (2.71)

Which, for an unsteady flow is:

t body od dy- (VG + V Gp+ ( X p)) [Vb + V x ] + 2 Vq - V X 2]u

77dt~~~~o~dt~~~ +2~~ 7~IT~otupepper
~
r ibody + tIbody - (VG
ad~
+ VGp + (
~~~~~
~~It
X TGp))' [Vo + V X ] + 2 IV + V X f
2lower a~ ~~~~~~~t
0,

(2.72)

and for a steady flow, the above expression simplifies to:

p|Utupper- 2-uer
P = 01 (2.73)

where, U in equation 2.73 is the total velocity of the fluid relative to the trailing edge. A
linearized version of the steady pressure continuity condition at the trailing edge is also
considered for use in certain formulations. The linearized, steady, Kutta condition provides
a rapid means to compute the potential jump [123],

Oupper- lower -= A¢cake. (2.74)

Here, the subscripts upper and lower refer to points on the upper and lower surfaces of the
trailing edge of the wing. This linearized Kutta condition assumes the fluid flows in a direc-
tion normal to the trailing edge cusp. For most applications this is an adequate assumption,

44
however, for highly loaded, low aspect ratio, or highly swept wings this linearization will
be inaccurate. Furthermore, although this Kutta condition performs well as a first approxi-
mation for unsteady lifting flows, it should be cautioned that highly unsteady flows will not
be accurately modeled.
For membrane flows, a condition which imposes no trailing edge vorticity has been
implemented. The condition requires that the lifting surface potential jump at the trailing
edge be equal to the potential jump imposed by the wake at the trailing edge:

/lmLwing= Allwake- (2.75)

Equation 2.75 requires that there be no net vorticity capable of turning the flow around the
trailing edge line.
For unsteady flows, any increase in bound vorticity on the wing must be balanced by
an equivalent increase in vorticity in the wake. The formal statement for this condition
(attributed to Kelvin [9]) is:

[dspan rspan (2.76)

dt wing dt wake

where F represents the circulation strength of the wing and body. In this relationship,
rspan on the wing represents the integral of the bound vorticity. Therefore, in order to
satisfy the above equation, the rate at which body bound vorticity increases must be equal
in magnitude (but opposite in direction) to the rate of of vorticity shed into the wake.
Satisfying the trailing edge Kutta condition of choice, and representing the domain vorticity
using a vortex or doublet wake representation will result in the bound vorticity increase on
the wing being appropriately accounted for in the wake.

2.5 Conclusions

This chapter presents the theoretical background for the potential flow solution framework.
Considering the wake as either a potential influence (due to a dipole sheet) or a veloc-
ity influence (due to domain vorticity) proves valuable in the implementation of the panel

45
method. The Helmholtz decomposition of the fluid in the domain provides a means to rep-
resent the scalar and vector potential components of the flow independently, thus providing
flexibility in the discretization.

46
Chapter 3

High Order, Curved Panel Integration

The justification for the use of curved panel, high order integration approaches is examined
in the first section of the chapter; following that integration techniques for the self term,
higher order, single and double layer potentials are presented.

3.1 Boundary Element Discretization

Discretizing the boundary integral equation involves approximations to both the boundary
surface and the solution. In most boundary element methods the body is discretized into a
series of triangular or rectangular panels. Over each of the panels, a basis function is used
to represent the unknown quantity. Due to the discrete approximations, panel integrals
such as the ones presented in equations 3.1-3.4 need to be evaluated for the discretization
of interest:
(JpPanel 7 77) 1
(3.1)
Palnel I 1 1dpanel

(3.2)

1
/,Panel =
( ( t
uP7anel
'
4 9'i
/)aPae
[HI) - x1dS e) (3.3)

Panel I (3.4)

47
Due to the singularity of the integral equation kernel under consideration, expressions for
computing the above integrals when the evaluation point and the panel are coincident are
not simple by nature. This section discusses the limitations of lower order basis functions
as well as lower order geometry representations. In the sections which follow, several tech-
niques for improving the discrete approximations using higher order techniques to remedy
these drawbacks are presented.

3.1.1 Discussion of Orders of Approximation

In order to justify the investigation of higher order approaches a brief presentation of the
limitations of some lower order approaches are presented.

1. Convergence Rates: One would expect that more accurate approximations of the
boundary integral equation will produce more accurate computations. Care however
must be taken to ensure that the surface discretization is consistent with the order of
approximation of the solution in order to achieve optimal convergence for a given
order of representation. Figure 3-13 illustrates the various convergence properties
for a Green's Theorem solution of the potential around a unit sphere with succes-
sive grid refinements. A flat panel discretization of curved surfaces has an area
which converges to the true area as O(N-1), where N is the number of panels in
the discrete approximation. As such, flat panel discretizations of curved surfaces
at best will converge with an O(N-1) rate due to the area's role in the discrete ap-
proximation. To achieve higher rates of convergence, a coupled approach in which
higher order surface representations must be used in conjunction with higher order
solution approximations. Several researchers have investigated higher order BEM
approaches[115, 116, 117, 118, 119, 120, 121, 123, 124, 125, 111]; however, diffi-
culties in evaluating the panel integrals prevent the more common use of these meth-
ods. In this chapter a set of high order curved panel integration approaches on curved
panels will be presented.

2. Velocity Calculations: The accurate computation of velocity is important for comput-

ing the forces and moments. This is especially true for the potential based Green's

48
 Theorem and Source-Doublet formulations. There are two approaches for computing
the velocity:

(a) Numericaldifferencingapproachesfor computingsurfacederivativesof thepo-

tential. For the potential based computations (Green's Theorem and Source-
Doublet formulations), the unknown is a potential. The surface gradient of the
potential is the tangential velocity due to that potential. Computing the surface
derivatives of the potential requires that either:

i. The potential is represented using linear or higher order basis functions

such that gradients in surface potential amount to differentiating the basis
function representation of the solution,

ii. The discrete elements are ordered in such a manner that finite difference
approximations can be made.

(b) Directforward evaluationof the integralequations.The forwardevaluationof

the boundary integral equation for the velocity is another option for velocity
post processing. Figure 3-1 illustrates the evaluation of the velocity integral
equation post-solution for a constant collocation solution. The velocity com-
puted using the Green's Theorem BIE is visibly wrong. This arises due to the
nature of the dipole representation. The constant dipole is equivalent to a vortex
ring. The vortex ring represents the normal velocity accurately; however, close
to the surface of the body the computation of the tangential derivative is increas-
ingly different than physical reality. Figure 3-2 demonstrates the inconsistency
with the constant basis dipole used for a tangential velocity calculation.

As a result of the desire for more accurate solutions with fewer panels as well as accuracy
in the computation of velocities, pressures, forces and moments, higher order BEM approx-
imations to the BIE's were investigated. This chapter presents some methodology which is
a result of the investigations.

49
 -1

&.0.5
=3
4)
o
0
4) 0
0n

0.5

·1
-1.5 -1 -0.5 0 0.5 1 1.5
x-position

Figure 3-1: A comparison of the pressure around a sphere, computed using a source Neu-
mann velocity formulation and a Green's theorem formulation. The velocity was computed
using a forward evaluation of the boundary intergral equations. The result deomstrates the
diminished accuracy of the Green's theorem approach.

3.1.2 Orders of Approximation

These approximations are made as follows:

1. Surface Discretization: The surface discretization is performed using a watertight

surface triangulation. The majority of the results presented in this thesis were per-
formed with a flat panel discretization of the surface, although, initial work presented
in this chapter demonstrates several advantages of using curved element discretiza-
tions.

2. Solution Discretization: The solution discretization is performed using piecewise

linear basis functions with solution continuity enforced. These CO continuous ele-

50
 Dipole Strength

/vortices\
~UiV~

(a) (b)

~
+
I
I
~
I
V(x) r------I V(x)
I I
r-----.J
I
I I
I
I
+I
I
+
I
I
I I I I
I I
t (c)
t (d)

Figure 3-2: An illustration of the inconsistency of using constant doublet representations

when modeling the tangential velocity close to a body when a forward evaluation of the
BIE is used. (a) shows the linear dipole and the equivalent constant vorticity distribution
on a surface, (b) shows the same for a constant dipole and the equivalent point vortices,
(c) shows the tangential velocity computed a small E above the surface for the linear dipole
distribution, while (d) shows the tangential velocity due to the constant distribution of dou-
blets. Notice, in (d) the tangential velocity at the wall of the body is zero at all locations
except for locations where there is a discontinuous jump in doublet strength. There the
velocity is infinite. As the evaluation point tends away from the surface, the velocity in
(d) will approach the velocity in (c), however, on average this occurs after several panel
lengths away from the surface. The difficulty with computing the tangential velocity arises
due to the poor approximation of tangential velocity self term for constant dipole panels.
As such, linear or higher order dipole distributions are recommended

ments are similar in nature to those used in the finite element community [126]. In
addition quadratic basis function approaches are presented in this chapter.

51
By discretising the surface geometry into a set of flat panels with associated compact sup-
port basis functions, the Boundary Integral Equations presented in Chapter 2 form a Bound-
ary Element Method linear system of equations.

3.1.3 Brief Review of Flat Panel Integration Strategies

In their paper on panel methods, Hess and Smith [13] proposed several analytical formulas
for computing panel integrals for sources and doublets. Later, Newman [112] proposed a
more general framework for the computation of analytical formulas over flat panel poly-
gons. The framework proposed by Newman included a proposed recursive algorithm for
the integration of higher order polynomials over flat polygons. The recursive formulas are
presented in more detail by Wang [111]. Although Hess, Smith [13], Newman [112], and
others have proposed formulas for the analytical integration, numerical approaches do ex-
ist. Based on adaptive quadrature schemes, these methods are typically more expensive for
computing the self term integrals. Some examples of these approaches for both flat and
curved panels can be found in [128, 129, 130, 131, 132, 133].

3.2 Integration Approaches for Quadratic Basis Functions

on Quadratic Curved Panels

The geometry is described using parametrized triangular quadratic patches [126]:

X(, r) = ax 2 +b 2 + cqr d+ + e + fx (3.5)

Y(S,:r) avy2 +br 2 + cyrt]+dd + er + fy (3.6)

Z(, r) = a, 2 + br 2 + cZrl + dz + erl + f (3.7)

The coefficients ai through fi are determined from the nodes defining the curved panel
patch. The representation of the parametrization is shown in figure 3-3. Parametric quadratic
basis functions are used such that the patch geometry definition is consistent with the basis

52
 ",,
oa I
(0.0,1.0)
0('1
1
i
1
o~.!
1
o 2~
I
D,I

iI J'j
04 i
1 (0.0,0.5)
I
.utJI

.u~.i
1
1
1 ~:(---

\
.o~'\
\
U\

Q~\.~ '~e\.o"~~~~-;-~-\
'08
(0.0,0.0) (0.5,0.0) (1.0,0.0)

X-}(.Z: Coordinate System Quadratic Patch s-t: Parametric Space Reference Triangle

Figure 3-3: The relationship between the quadratically curved panel and the flat parametric
triangle. The curved panel in this example is a ~ segment of a sphere.

function definition [126]. These basis functions for a given panel are presented pictorially
in figure 3-4.

3.3 Panel Integration Approaches for Non-Self Term In-

tegrals

High order methods currently consider adaptive quadrature schelnes[115], expansions[ 119]
and semi-analytical approaches [111, 121] for the evaluation of curved panel integrals.
Nearby non-se~f term panel integrals can easily be computed using the lnethod developed
by Wang et a1.[111]. The far field interactions can be computed using a quadrature approx-
imation. Although the Wang et a1. method is used for nearby panels, it should be noted that
the methods presented in this chapter for the self term integrals can be extended to compute
the non-self term integrals; however, they are typically less efficient than the Wang et a1.
approach. The Wang et a1. method does not work well for the self term evaluation in this
quadratic panel method due to:

53
 ,",
'-....:~¥

:'~ ;~
. " .... ~ .~
- ..,~ ..."-:1
8(u'tJ'# Z B;;-;;/lJ

BastJ#4

Figure 3-4: The 6-quadratic basis functions on the reference triangle.

1. The inability to form a ratio between the double layer kernel for a flat panel and the

double layer kernel for a curved panel. This is shown below in the equivalent flat

panel representation of the curved panel integral:

<p(x) = is
r I
~l(X)f}n
f} ( 1
Ilx-x'll ) (~ (~)
2..-(_l_),JI )
dB!,
I
(3.8)
f ! on Ilx-x~II

but:

]( a
_
(1) - ----
o~ (lIx-1xf II) (3.9)
an;' ~1 (IIX~X~ II )

Is undefined at x = x' and zero elsewhere. Hence, it is not possible to fit a polynomial

accurately through this representation.

2. The difficulties associated with smoothly mapping non-constant radii of curvature to

a flat, straight edged reference triangle. It is challenging to determine a mapping in

which the self term curved-to-flat panel transformation is easily fitted with a polyno-

mial in these cases. Although the Wang et al. method can be applied, unless more

54
 complicated mapping approaches are considered, the accuracy of the approach will
typically be diminished.

3.4 Panel Integration: The Self Term Integrals

Methods for quadratic basis function single and double layer self term integration have
been developed for quadratically curved panels. These methods can easily be extended to
higher orders.

3.4.1 Single Layer Integrals: The Self Term Integral

Although the self term integration for the single layer which is used is similar to the Wang
et al. approach and those traditionally used[l 15], the differences introduced in this compu-
tation, albeit small, are important to improve the accuracy of the quadratically curved panel
integrals.

Background Theory

The single layer potential self term integration is performed numerically after appropri-
ately transforming the curved panel into a flat panel (with curved edges). Following the
transformation, the panel integral is evaluated using quadrature in cylindrical coordinates.
The transformation to cylindrical coordinates to desingularize the integrand was presented
by Hess and Smith[122] in their analytical computation approaches for flat panel integrals,
and the ideas have been further used in numerical implementations for higher order[l 15].
The background theory of the method used in this approach is described below:

1. Consider the integral governing the self term single layer on a curved panel:

@(z) = j (1, i
II· 'H
idS', (3.10)

where a(J', I') is the basis function representation of the single layer strength on the
panel.

55
 PotnLv are pnjeciedfrnm the cuzrvedpanelto thex-yplane,
such that theprojectionis orthogonalto thexy plane.

_0 .I * ~C .rO f·

Figure 3-5: The orthonormal projection shown for a large selection of points on the curved
panel to the tangent plane. Notice that the curved panel is projected onto the tangent plane
in a manner that is normal to the tangent plane.

2. Project the curved panel geometry using an orthonormal projection onto the tangent
plane at the evaluation point as shown in figure 3-5. As a result the integral over the
flat tangent plane is:

T(x) /
,fo
,7~C(c
'Cc) 1~
I,
-x IIIf ~d~
IJldSf, (3.11)

where IJI is the determinant of the Jacobian due to the orthonormal mapping. The
value of lJI evaluated across the flat panel will be smooth assuming large distortions
are not encountered in the mapping.

3. As in the Wang et al. method [ 111], multiply and divide the integrand by 11xx,11:
f1-;I

(x)= /7(/, ?77)l J I) dS. (3.12)

s , I~x - Xf Hix- X/ 1

4. Combine the basis representation polynomial (ua(', r'))

Cl with
CrlrlII VLI the polynomial repre-
CC ~l

56
 sentation of the mapping, into a single polynomial representation Q (r, Of), and the
integral becomes:

· (x)
= I Q(r 0') dS (3.13)

5. Equation 3.13 is then re-written in cylindrical coordinates as:

(z) j
0=27r

Jf=O
r=R(0)

Jr=O
Q(rW,[O)'
1

11
Ijr Id8'dr' (3.14)

6. Which, when simplified becomes:

T(x)- f j Q(r, Of)dO'dr' (3.15)

7. The resulting integral is simple to evaluate by first analytically integrating the poly-
nomial Q (r, 0.) in r, and then using one dimensional quadrature, for the remaining
0-integral. In cases where curved panels are significantly distorted or stretched, care
should be taken as the polynomial fit will likely require significantly higher order
polynomial representations for an accurate solution than when the panels are regu-
larly shaped.

Implementation Approach

In order to implement the above approach, the following steps are used:

1. Transform the curved panel such that the tangent plane at the evaluation point lies
in the XT - YT plane, and the panel normal at the evaluation point is in the +ZT
direction, as shown in figure 3-6.

2. Determine a series of points on the curved panel which will be used to fit polynomial
approximations to distribution on the panel. One can use any appropriate set of
quadrature points; however, the polynomial fit will occur in R - e coordinates, and

as such, an appropriate quadrature scheme might be similar to that used over 4-sided
regions(with the sides of the rectangle being the R and 0 limits). In order to achieve
a more accurate approximation to the integral, subdivision of the original panel into

57
_ I----
PanelNormal Iectorat evaluationpoint.

EvaluationPoini

_1-·-)---^--I-^- ---XI--_-----_l--- --I-^I_--_-PI··P-L·--i- --^_IIIX·llll-_-----_-111111111111111 · 1I_1IIII

Figure 3-6: The panel after it has been appropriately transformed such that the tangent
plane at the evaluation point is the XT - YT plane, and the panel normal at the evaluation
point is in the +ZT direction.

three panels should be performed, where the subdivision occurs by connecting lines
between the vertices and the evaluation point. Figure 3-7 demonstrates the points
which will be used for the polynomial fit, as well as the division of the panel into
three subpanels.

3. Orthogonally project the polynomial fit points onto the XT - YT plane. One can also
project the edges of the panel onto the XT - YT plane. The projection of the curved
panel edges and polynomial fit points onto the XT - YT plane forms a flat reference
panel for the integration. This projection is demonstrated in figure 3-8 and figure 3-5.

4. Compute the integrand, Q(r), O ) at the polynomial fit points, which is product of
the following values:

(a) The determinant of the Jacobian of the mapping between the curved panel and
the orthonormal projection of the panel onto the XT - YT plane. The expression

58
 Points on the curved panel which are lued
to npproximare the integrand as a polynomial.

02

I
.1

Figure 3-7: The points which are used to compute the polynomial approximation of the
smooth integrand. Notice the polynomial fitting points are arranged in a manner that re-
sembles a quadrature scheme for a degenerate quadrilateral. This apparent degeneracy will
disappear once R - 8 coordinates are used. In addition, one should note the division of
the panel into three subpanels. The layout of quadrature points and division of the panel is
performed on the ~ - TJ parametric triangle.

of the determinant can be found analytically for quadratic panels.

(b) The ratio between the kernel evaluated on the flat panel to the evaluation point
and the kernel evaluated on the curved panel to the evaluation point.

(c) The Basis function value, o-(rj, 8j).

5. Compute a polynomial approximation to represent the integrand in the R - 8 co-

ordinate system. This polynomial approximation is computed by a simple linear
system solve (The right hand side vector contains the known values of the function,
the system matrix is a Vandermonde matrix where the entries are evaluated at the
polynomial fit points (in R - 8 coordinates), while the vector of unknowns contains
the unknown polynomial coefficients).

59
 Points on the curved panel and their
/ orthonormal projection onto the x-y
/ plane.

IU
...... ,

.. -'<:;~'E'1~~ J

-t2

'( 1

'.0

Tire x-y plane onto which the curved panel

is projected

Figure 3-8: The orthonormal projection shown for a large selection of points on the curved
panel. Note that the panel edges will not be straight lines in the general case. This is
the fundamental difference between the Wang et al. method and the current method for
computing the self term integration of the single layer.

6. Integrate the resulting polynomial in R analytically. The integrand of:

w(x) =
r
J8=O
8=27r lr=O
r=R(O)

Q(rj,8j)drd8' (3.16)

is a polynomial representation in R-8. As such, the integral is simply an appropriate

coefficient multiplication and exponent augmentation. The integral becomes:

w(x) = L:: 2n

U(Rj(Bj), Bj )dB', (3.17)

where U(Rj(8j), 8j) is the polynomial resulting from integration. Rj(8) is the dis-
tance from the evaluation point to the flat panel edge for a given angle 8. This integral
is merely an integration of U (Rj ( 8), 8j) in 8 around the edges of the panel.

7. Detennine a quadratic parametric relation for the edges of the curved panel (a one-
dimensional quadratic parametrization of the projected edge), such that in the param-

60
 eter t the 0-integral becomes:

r0=27 rt=1
T(x)= U(R (0), ')d'= U(R'('(t)) Os(t))IJotIdt. (3.18)
JO=0 t=O

Since the above integral is computed numerically using quadrature, place quadrature
points on the panel edges. These quadrature points are placed by mapping the Gauss
quadrature points from the parametric representation to the curved edges.

8. Integrate the line integral numerically using the panel edge based quadrature routine.

The scheme presented for the single layer self term can be applied to arbitrary curved panels
with high order basis functions. In addition, the method does not have the same limitations
as the Wang et al. method due to the lack of restrictions on the reference panel edge shape.
In the current method, the reference flat panel edges can be straight or curved lines. The
curved panel self term single layer integration using the above process is simple from both
an application and conceptual point of view.

3.4.2 Double Layer Integrals: The Self Term Integral

As with the single layer integration, the double layer self term integration approach we
present is a conceptually simple numerical integration. Similar to the single layer, the
approach we present for the double layer is not restricted simply to self term panels.

Background Theory

The double layer self term integral is performed through the implementation of some funda-
mental ideas presented in Kellogg [114] and later exploited in Newman [112]. The double
layer expression is considered from a slightly different perspective as shown below:

@(x 1
W)=On jx - xl ldS', (3.19)
ch~ xl:- xc' dS',

61
where b(x') is the basis function representation. If we re-write the expression as:

<>()=
-V [t(x)h
(i X)1
' dS' (3.20)

the expression in the square brackets can be regarded as the velocity at the panel surface
integration point x' due to a point charge at x. Furthermore, if we assume that for now
bt(x') is a unit constant distribution over the panel:

(X) f=t. (V( 11 - x') dS' (3.21)

which is exactly the flux through the panel S' due to a point charge at x.

Since the point charge at x emits a divergence free velocity in a radial manner, we can
equivalently say that the flux due to the point source which passes through the panel, is
identical to the flux which passes through any radial projection of the panel (where the
radial projection is centered at the point source which is also the evaluation point). If the
panel is radially projected onto a unit sphere centered at the evaluation point, the integration
is equivalent to determining the flux through a portion of a sphere defined by the radial
projection of the panel. This is shown in figure 3-9.

The integral for the flux through the spherical patch formed by the radial projection of
the panel (where the radial projection is centered at the evaluation point which lies at the
origin) is expressed as:

@(x)= (V( dS I dS' (3.22)

Which is identically the area of the sphere onto which the projection of the panel acts. In

order to compute higher order distributions a similar radial projection can be used. In the
higher order approach, the basis function is considered as a weighting function on the flux
through the panel due to the source point charge. In other words, both the panel and the
basis function are projected onto the unit sphere in a radial projection. This is shown by
considering equation 3.20.

The interpretation of the dipole integral here is a flux through the panel surface S, due

62
 ):1 ~ ",~;1f('!J:.~I';Y
Panel over which double layer integral
is being computed,
0.5.1 1'" .. ;:

Evaluation point: x
l).1
"
I
(l JJ

_--<:'-
-.----...;-----'
U~

,~ -;:- () Radial Projection of pa1U!/ onto the unit sphere

centered at the evaluation point

Figure 3-9: The radial projection of a flat panel onto the unit sphere centered at the evalua-
tion point, x.

a point charge located at x, weighted with the basis function value p,(X') at the point x' at
which the flux is being evaluated. Again, since the flux due to the source point charge is
divergence free, this weighed flux for each point x' on the panel can be projected onto the
unit sphere in a similar radial projection as before. The resulting integral becomes:

(3.23)

Which is merely a calculation of the basis function weighted area of the radially projected
panel onto the unit sphere.

A Note on the Self Term Double Layer Potential Integral

When the self term of the double layer potential integral over a surface is considered, care
must be taken. Consider the double layer integral:

(3.24)

63
At x = x', the integral is undefined. As such the point is excluded from the domain and the
Cauchy Principle Value is evaluated, such that:

(x) = ±2t(x)+ (x')On -1 xdS', (3.25)

where the ±2-u(x) term is the Cauchy Principle Value of the integral, the sign of which
depends on the direction in which one approaches the panel (from above or below).The
remaining integral is evaluated over all x' :7 x points on the panel. In the section which
follows, the evaluation of the integral is of interest and will be described.

Implementation Approach

The implementation of the integration of the double layer over curved panels is presented
for quadratic basis functions; however, it should be noted that arbitrary order basis func-
tions can be used. Furthermore, it should be noted that the constant basis computation is
much less computationally complex than the higher order computations, since it amounts
to computing the area of the sphere between the equator and the projection of the panel.
This is a direct result of the fact that there is no varying basis distribution on the panel.
The following steps are performed to compute the double layer integral.

1. Transform the curved panel such that the tangent plane at the evaluation point lies
in the XT - YT plane, and the panel normal at the evaluation point is in the +ZT

direction, as shown in figure 3-6. Furthermore, translate the transformed panel such
that the evaluation point lies at (Xe, Ye,: e) = (0, 0, 0).

2. Split the curved triangle into 6-sub triangles. These sub triangles are formed by:

(a) Dividing the panel along lines joining the evaluation point and the vertices.

(b) Dividing the panel along the lines joining the evaluation point and the nearest
point on each of the panel edges (in the current implementation this is measured
on the reference parameter based triangle).

3. Determine quadrature points on each of the six sub triangles. In addition to the
points themselves, determine the basis function value at each of the points. These

64
 quadrature points are determined based on quadrature rules applied to the ~ - T}

parameter triangles and then mapped to the surface of the curved panel. As in the
single layer self term computation, these points are used for the determination of
a polynomial approximation, hence any appropriate distribution of points will be
adequate. Again, we choose to use the quadrature rule for four-sided domains. Figure
3-10 demonstrates the splitting of the curved triangle and the placement of quadrature
points.

['
I
tT
J
.:i
l
I
,L

(a) Split the cun'f!d triangle inlo 6-sub Iriangles (b) Determine Function Evaluation Points on each triangle

Figure 3-10: The 6-subtriangles illustrated on the curved panel. The illustration on the
right shows the points which will be used in the polynomial approximation of the integrand
after it is projected to the unit sphere.

4. Radially project the points on the panel onto the unit sphere. The points can be
projected onto the unit sphere by:

P(x,y,z)
Ps ( X, y, z ) = ----- (3.26)
IIP(x, y, z) II

Note, since the integral in question excludes the point x = x', the projection of the
portion of the panel directly encircling the origin (which is the evaluation point), will
lie on the equator of the sphere.

65
 By radially projecting the panel onto a unit sphere centered at the evaluation point,
the integral can be represented as:

cp(x) =
i
S
J-l(X')-a
a
n x
1
II _
x
'lidS' -7
1
8=0
8
=27r
1'1>=4>(8)

'1>=0
jl(cjJ,e)-a
a
n
(1)
-
l'
r2sin(cjJ)dcjJcLe
(3.27)
Which, when we consider the unit sphere (1' = 1), becolnes:

rO=27r r4>=4>(O)
cp(x) = }o=o } 4>=7f Jl( cjJ,e) sin( cjJ)dcjJde (3.28)

If, the sin( cjJ) term in the integrand is combined with the basis representation to give
P( cjJ, B), the resulting integral expression is:

(3.29)

5. Determine the cP - 8 coordinates of each of the points.

Figure 3-11: The radial projection of the points on the panel and the edges onto a unit
sphere centered at the evaluation point.

66
 6. Sort the projected points such that they lie within the O angles of the original curved
panel subdivision.

7. For each of the 6 sub-triangles determine the polynomial approximation, P(, 0),
which best represents the basis function value as it is radially projected from the
panel points to the points on the sphere.

8. Integrate the inner integral of the polynomial, P(q, 0), analytically in cg:

(x): = ] P(, =)dqdO

W(x) - Q(q0(o),O)dO (3.30)

9. Integrate the remaining line integral using quadrature. The edge of the panel over
which the line integral is being performed must be expressed as a quadratic line
based on a parameter t. The integral to be computed then becomes:

t=l
4 (x) = j t=O
Q(((t)), (t))lJIotdt, (3.31)

which is easily computed using moderate orders of Gaussian quadrature in one di-
mension.

10. Since we are concerned only with evaluating the Cauchy Principle Value integral in
this computation, in order to get the full influence of the self term, we add to this
result the appropriate +27rp Cauchy Principle Value (CPV).

3.4.3 Potential Flow Around the Unit Sphere

In order to demonstrate the quadratic panel method arising from application of the integra-
tion techniques described, a series of solutions for the flow around increasingly refined unit
spheres is presented. The potential flow around the unit sphere has an analytical solution
for the perturbation potential which is given by[ 113]:

¢=Uoocos0 22 (3.32)

67
The simulations considered were run at U = 1, over unit radius (R = 1) spheres. The
panel method approximations which are considered in the convergence comparison are:

1. Constant collocation on planar panels

2. Linear basis Galerkin on planar panels

3. Quadratic basis Galerkin on planar panels

4. Quadratic basis Galerkin on quadratic curved patches

The results of the quadratic curved panels are visibly superior to the constant collocation as
presented in figure 3-12. The convergence of each of the methods to the analytical solution

n,
U.

TheoreticalPotential
0. 4 , Quadratic-BEM
O Constant-Collocation-BEM
0.3 - ConstantCollocationBasisFunctions

0.2 . . . .. . .. . .... ..........

! _... .... . ... ......... ...............

..
0..1 ....... .... .. . . .... .. ....... .. . .. ....... ... .......... ..... ................

C
I
0S C i : Z m~~

-0.1

-0.2

-0.3 s....

-0.4

71 t
..........' .'-0.8
.. .. . . i i
I
i
I
. .
i
i
i
i
g- I -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8
_11 ": X-Coordinate

Figure 3-12: An illustration of the solution ¢ computed on a 32-panel sphere comparing

constant basis function representations with quadratic basis functions on quadratic patches.
The theoretical result is also plotted and coincides with the result from the high order panel
method. For clarity in the depiction of the constant basis function results, both the potential
value at the centroid and the basis function itself are represented on the plot. This further
demonstrates that the constant basis approximation will have poor error convergence.

is demonstrated in figure 3-13. The error was computed for each panel and integrated
using high order quadrature schemes. From figure 3-13, it can be seen that the highest

68
 10

l100 -- . ...........-- ............... .... .. ........ . . ..... ................................ ..... ..........-.- ....................................
..............
. .. ......

'z
10

lo-'
_-
l LI Error:QuadraticBasis,Curved Panels .............. ......... .......... ...........
--- L2 Error: QuadraticBasis,Cured Panels
L I Error: QuadraticBasis,FlatPanels
-- L2 Error: QuadraticBasis,Fla Panels
10-3 -- L1 Error: Linear Basis,FlatPanels
- L2 Error :Linear Basis,FlatPanels
-. -LI Error: ConstantBasis,Flat Panels
-- L2 Error: ConstantBasis,Flat Panels
1 0° i I i i1 2 i 3
10 101 10 10
Nmnberof Panelsin DiscreteApproximation

Figure 3-13: The error convergence of the potential for constant basis representation on flat
panels, linear basis representation on flat panels, and quadratic representation on curved
panels. The plots show L1 and L2 error integrated over the sphere using the underlying
basis representation for the integration. Notice that the constant basis functions on flat
panels converge at a rate of O(N-2 ), the linear basis on flat panels at O(N- 1 ), the quadratic
basis on flat panels at O(N- 1 ), and the quadratic basis on quadratic patches at O(N-2 ).

convergence rate for flat panel discretizations is limited by the geometry convergence. This
is realized when comparing the quadratic, linear and constant basis solutions when flat
panel discretizations are used. The flat panel discretization surface area converges at a
rate inversely proportional to the number of panels (O(N -1 ) = O(h2 )), hence becomes
the limiting convergence rate for the case. It is interesting to note that the results for the
linear basis function representation are the most accurate of the set of computations done
on flat panels. Although this may appear strange, it is consistent with the approximations
being made. Quadratic basis functions will provide a more accurate solution to the discrete
geometry (which is based on a flat panel representation), which in the case of a sphere or
curved surface, is typically a worse approximation than the linear basis on flat panels.

The convergence results for the curved panel high order basis functions demonstrate
that increased convergence rates can only be achieved if approximation order increases are

69
realized in both basis function approximation and surface discretization. If we consider the
panel methods described, it would take a discretization on the order of 105 - 106 constant
basis flat panels to achieve the error levels displayed by a 128 curved panel high order
method. Similarly, it would take on the order of 104panels to reach that level of accuracy
using linear basis functions on flat panels. From these results it can be seen that low panel
counts using accurate geometry representations and high order basis functions can often
lead to accuracy levels hard to achieve with flat panel lower order approaches.

3.5 Conclusion

The work presented in this chapter highlights the methods involved in computing panel
integrals for the high order BEM representation and proposes approaches for integrating the
self term integrals. It was found that the high order method, when compared on a number
of panels basis was always significantly more accurate than the low order approaches.

70
Chapter 4

Wake Details

The explicit modeling of wake shear layers in panel methods often presents challenges.
More specifically, the following challenges exist:

1. Ensuring the wake position is representative of a physically realistic vorticity dis-

tribution: The wake should evolve with the local freestream velocity, otherwise it is
a force producing surface. Many panel methods require user specified wakes[106,
116], while several others [35, 107] incorporate wake roll-up and wake evolution
routines. The wake roll-up routines often suffer from limitations associated with the
use of panels to discretize the wake potential jump or vorticity.

2. Satisfying a potential jump along the downstream portion of the body: The wake
represents a potential jump in the domain. This potential jump must also be present
at the body-wake intersection. Since the potential jump in a domain discontinuity,
meshing challenges will typically exist at wake-body intersections.

3. Using appropriate models to deal with wake intersections and downstream lifting
surfaces: During unsteady motions, wakes may readily intersect with downstream
surfaces. When a doublet wake intersects with a downstream surface, the solution of
the potential flow will lose physical meaning. Although wake-body intersections are
non-physical, they can often occur due to the nature of the discrete representation.

In this chapter, several techniques for more advanced and automatic wake modeling and
evolution are presented.

71
4.1 Body Piercing Wakes

In traditional panel method approaches, the wing trailing potential jump is represented us-
ing a dipole sheet extending behind the wings and conforming to the downstreatn portions
of any bodies, as demonstrated in figure 4-1. The dipole wake sheets induce a potential

Figure 4-1: A pictorial representation of the Green's Theorem fonnulation of the Boundary
Integral Equation.

jump in the domain across the thin trailing wake. For potential based solutions, such as the
source-doublet formulation, wake-body intersections must be carefully treated to preserve
and accurately model the potential jump. In order to do this, the wake sheet must be abutted
to the body and form an impermeable intersection (The wake should align exactly with the
fuselage). The result is a large potential jump along the wake-body intersection. This is
illustrated in figure 4-2. When the velocity is computed, the effect of the potential jump
along the side of the body cancels with the effect of the potential jump due to the wake.
As such, the velocity at all points on the body is smooth. When determining the velocity
near the wake-body intersection, a limiting process must be used such that the velocity
evaluation point never coincides with the triple point intersection of the wake and fuselage.

Significant constraints on wake-body meshing are imposed through the use of dipole

72
 Figure 4-2: The effect of dipole wakes on the potential at a wake-body interface.

potential jump sheet wakes. Since the potential jump in the wakes is usually large, one
must ensure that the dipole wakes align with the body discretization exactly. If the wake
does not align with the body discretization, the edge of the wake sheet next to the body
will act as a strong vortex. In turn, this will cause erroneous values of the potential along
the body. For the wake piercing formulation presented in this thesis, the wakes are treated
as equivalent vorticity distributions which produce a velocity influence. By examining the
traditional dipole wake sheet implementation as analogous vorticity influences, one finds
that there are three strong line vortices at all wake-body intersections (see figure 4-3). This
set of strong vortices is not physical, rather it is a mathematical consequence of preserving
the wake potential jump right up to the fuselage. The combined-body and wake line vortices
are of equal magnitude, in the same positions, but of opposite signs, thus, in the velocity
evaluation, they completely cancel each other (to preserve the physical reality of Ininin1al
vorticity at the wake-body intersection). The net velocity at wake-body intersections is
therefore smooth.

The concept of the body piercing vortex wakes is a natural extension of using an equiv-
alent vorticity wake to produce a velocity influence on the body, rather than a potential

73
 Figure 4-3: Representing the potential wake as vortices.

Jump. Additionally, consider that the wake is merely a means by which the body bound
vorticity (or equivalent circulation) is shed into the domain. The body bound vorticity can
be mathematically represented in many diverse ways in a potential flow. As such, it is
convenient to extend the vortex wakes into the body and represent the circulation of the
body and wing using an embedded or body piercing vorticity distribution. By enforcing
the influence of this body piercing wake as a velocity influence on the wing-body system,
the potential jump in the solution for the body surface potential is not explicitly lTIodeled,
rather it is implicitly modeled via the velocity influence vorticity body piercing wakes. By
accounting for all contributions to the total or overall potential, the potential jump at the
wake-body intersection is re-introduced via the potential influence of the vorticity wake
(which is equivalent to a dipole sheet). This body piercing wake formulation is illustrated
in figure 4-4. Although the body piercing wake seems highly un-natural, the portion of the
wake inside of the body is merely a means of representing the body bound vorticity. The
body bound vorticity is required to produce the circulation around the lifting surface.

The body piercing wake formulation is implemented as follows:

1. Setup a solution to the wing-body component of the potential flow solution using any

74
 Wake inside wing &
body (circulation).

Body-wake piercing
intersection

Figure 4-4: A demonstration of the body piercing wakes.

potential flow formulation.

2. Prescribe a wake which pierces the body. The wake should have a gradual variation
of doublet strength from a value of zero at internal wake boundary points to the value
prescribed by the Kutta condition at points where the wake pierces the boundary from
interior to the exterior. For wings this piercing wake must intersect the wing at the
trailing edge at which a Kutta condition will be applied.

3. Express the wake influence in terms of a contribution to the velocity boundary con-
dition. Solve the system of equations, ensuring that the boundary conditions and
trailing edge Kutta condition are satisfied. Be careful to recognize that the quantity
being solved for has a different representation than it would if one were to consider
a traditional panel method formulation (ie. if one is solving for the doublet strength,
realize that the doublet strength no longer represents the perturbation potential, rather
it represents the perturbation potential minus the wake potential !).

4. If the potential is desired, compute the potential influence due to the freestream,
wakes and body and add them together.

75
Although, it is believed that the body piercing wake formulation is novel for three din1en-

sional potential flow applications, initial proposals for body piercing wake concepts can

be seen in SOlne of the two dimensional panel method investigations by Hess, Pierce and

Smith [13]. In two dimensions a single vortex point inside 2-Dimensional airfoils (repre-

sented using Neumann source panels) is used to provide the circulation necessary for lift

generation. The idea was not pursued for long, since Hess, Pierce and Sn1ith soon found

that the vorticity panel Neumann approach was Inore convenient for their applications.

A brief example of body piercing wakes is presented to illuminate the application.

(a) Body Conforming Wake.

Unit sphere test ca<;e.

(b) Body Piercing Wake.

Unit sphere test case.

Figure 4-5: A schematic of the sphere example with (a) Body conforming wakes, and (b)
body piercing wakes. The body conforming wake is represented using a potential jump
inducing potential jump, while the body piercing wake is represented using an equivalent
vortex ring (shown using thick lines).

A dipole sheet enclosing the trailing half hemisphere of a sphere: The exmnple of
a body piercing wake is the flow around a unit sphere centered at the origin. A trailing

wake sheet of unit strength is prescribed. In the traditional formulation, the wake sheet

is described using a dipole sheet of constant strength which conforms to the sphere as

shown in figure 4-5a. In the body-piercing-wakes formulation, the wake is represented

using a quadrilateral vortex ring which pierces into the sphere as shown in figure 4-5b. The

76
 0.6
• Body Piercing Wakes Formulation
,) Body Conforming Wakes Formulation

•••
-0.2
"
-0.4 .o:r
I
-6.1:
-0.6 .• L. ~_._~._.~_.__
., -2$ .: -1$ -~ -46 0

-0.8
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8
Streamwise Position

Figure 4-6: The total potential computed due to a body conforming wake, (0), and body
piercing wake, (*).

resulting total potential for the two separate simulations is presented in figure 4-6. The
total potential is the superposition of the body surface potential, and the potential due to
the wake sheet.

4.1.1 Comments on Body-Piercing Wakes

For unsteady panel methods, a simple, and general approach to handle the wakes is a ne-
cessity. The body piercing wakes presented in this chapter have the potential to simplify
panel methods greatly. Without the need for complicated wake-body intersection griding,
it is possible to reduce user interaction times with panel methods by orders of magnitude.
The body piercing wake concept can be applied to all non-membrane panel method formu-
lations by appropriately incorporating the correct boundary conditions. Due to the conver-
sion from a potential jump influence wake sheet to an equivalent velocity influcence, the
discontinuities at wake-body intersections associated with the potential jump requirement
are eliminated.

The body piercing wake formulation is effective for higher order distributions with

77
large panel counts. Although the method will work for low order, coarse discretizations,
the effect of the internal wake will tend to reduce the accuracy of the solution (since the
wake influence becomes a velocity boundary condition). The use of at least linear strength
basis functions for the solution of body piercing wake problems is recommended.

4.2 The Vortex Particle Method

The vorticity in the wakes can be represented in a number of ways. In this thesis, a vortex
particle method is investigated for the representation of the domain vorticity. Although
a vortex particle method is chosen in this work, one could easily implement many other
forms of vortex representation including sheet or filament vorticity.

The domain vorticity is represented using discrete vortex particles which capture the
strength and magnitude of the vorticity. The overall domain vorticity is represented as the
summation over all of the discrete vortex particles in the domain, which is written as[88]:

j(R, t) = A, 3(t)volp6(R
- Rp(t)) = Ao5(t)6(R- (t)),
P P

where, cZp(t)volp is written as dpv(t). Here, volp refers to the volume of the fluid domain
which is represented by the vortex particle.

For a vortex particle representation of the vorticity, the discrete vector potential is:

1 1
qp(R,t)= 4r A a(R, t) -

The vorticity induced velocity is:

1 1 x~(1~,
V x Tp(R, t) = 4-E V (t) (R)t). (4.1)

Similarly, the gradient of the velocity term used for the vorticity stretching in the vorticity

78
evolution equation is:

V (VX(/k ) =-4EV Vt x (Rt). (4.2)

A Lagrangian reference frame is used for the evolution of the vorticity, such that the posi-
tion Rp(t) of a discrete vortex particle at any given time is:

d
dt Rp(t) = Up(R(t), t) (4.3)

The evolution of the vortex particle strength as it travels through the domain is:

D
Dt - p Vp()() t),t). (4.
Each of the vortex particles, or vortons, has an associated core in order to mimic the phys-
ical vortex core as well as to reduce the numerical instability of the vortex interactions. A
collection of vortex particle core expressions are presented in [88]. The evaluation of the
vorticity induced vector potential can be represented as a matrix vector product:

[C] [a]=

Here the vorticity is known and a single matrix vector product results in the vector potential.
Additional information about vortex particle core functions and vortex methods in general
can be found in [88].

Discrete Form of the Vorticity Evolution Equation

The evolution of vorticity is computed by discretizing the ODEs governing the vorticity
evolution and computing the vortex particle position at time, t + 1. A simple approach to
this would be to use a forward Euler equation:

R(t + 1) = R(t) + Up(R(t),t)At (4.5)

79
and then update the strength of the vortex as:

ap(t + 1) = Gp(t)+ dp(t). VUp((t), t)At. (4.6)

4.3 Conclusions

In this chapter, methods for simplifying and automating wake modeling have been pre-
sented. The approaches presented in this chapter help minimize the user involvement in the
solution of lifting body flows by:

1. Allowing CAD and CFD surface grids to be used for lifting body simulations. The
development of body piercing wakes eliminates the task of prescribing wake-body
intersections.

2. Using vortex particle methods to evolve the domain vorticity in time. The vortex
particles are not constrained by connectivity requirements. This permits automatic,
hands of wake evolution in complicated environments.

With the flexibility of vortex particle methods, coupled with the generality of body-piercing
wakes, the lifting body problem is significantly simplified.

80
Chapter 5

Implementation Details

In this chapter the implementation details of the panel method formulations are described.
The boundary integral equation is discretized using a boundary element method [89,
127, 90], while the vortex volume integral equation is represented using a vortex particle
method [59, 88, 57].

5.1 Boundary Element Method Implementation

The Galerkin Boundary Element Method is used to discretize the various potential flow for-
mulations. Solution basis functions of linear order are considered for all examples shown.
The Boundary integral equation is written in terms of a residual R:

(BIE - RHS) = R. (5.1)

As an example, the Neumann-Source problem is written as:

[N
E
_j=l
ub
1
dS~~"·:~
+ 4)1V - ( (T + GP ++ Q x Gp - x )) - R.
(5. 2)
The weighted residual is minimized using a Galerkin formulation [127]. For each panel i
on the discrete surface, a minimization of the residual is achieved by first multiplying by a

81
target basis function bi and integrating the product over the panel:

J bi (BIE - RHS) dST = J biRdST (5.3)

for the example of the Neumann-Source formulation (for all target basis i):

EJSTi[b fIr jb
(4rJSJi l dS + (D dST -

JSbi [(n( + P+ X rp-V x ;))] dST, (5.4)

=i| biR dST.

In effect, this formulation minimizes the residual by ensuring it lies in a space orthogonal
to the solution basis. In the panel method presented, the integration over the test basis
functions is performed using numerical quadrature schemes[134, 135, 136]. In addition, the
testing basis functions bi are of identical form to those used for the solution representation
bj.

As a result of the discretization and formation of a residual minimization statement, the

boundary integral equations become linear systems of equations which are readily solved
using any number of techniques. Although the Galerkin formulation statement is robust,
the computational effort is increased over the simpler collocation based approaches [127]
due to the necessity of evaluating the integral over the target panel numerically. This outer
integral evaluation results in an increase in complexity by a factor proportional to the num-
ber of quadrature points used in the integral evaluation.

5.2 Wake Implementation details

The vortex particle method wake representation was presented in chapter 4. Some ad-
ditional details pertaining to the implementation of the wake model are presented in this
section.

82
5.2.1 Kutta Condition Implementation

The addition of wakes involves adding new unknowns to the system. These unknowns can
be determined by adding a Kutta condition to the linear system. The Kutta condition can
be linear or non-linear. The linearized potential Kutta condition [123] is implemented by
adding appropriate rows and columns to the the influence coefficients matrix. The nonlinear
Kutta condition [137, 117, 142] is implemented using a fixed, approximate Jacobian in a
damped Newton Method [139] solver. The use of Quasi-Newton methods such as the BFGS
[140] algorithm will be investigated in future implementations. The advantage of a Quasi
Newton method implementation is that an approximate Jacobian matrix is constructed with
minimal time investment, thus reducing the cost involved in constructing the full Jacobian.
The linear system including a wake and Kutta condition can be written as:

[B2B] [W2B] 1 {B} BCbdy l

L [B2W] [W2W] {W} T. EKutta} J

where, B2B represents the body to body potential or velocity influence matrix, W2B rep-
resents the wake to body potential or velocity influence, B2W represents the body to wake
influence at the Kutta condition forcing points, and W2W represents the interaction be-
tween the wake and the Kutta condition points. For a simple linearized Source-Doublet
formulation the linear system is:

[B2Bm] [W2B,] 1 {abody} 1 [B] {Ubode}

[(Du- qL] [I] t {Pwfake} JI {0} J

This is a compact manner to write the linear system including the Kutta condition. For the
non-linear Kutta condition, an iterative procedure is used to solve for the surface singularity
strength:

{
1. For the current iteration wake strength, Wn}, solve the following intermediate lin-
ear system for the body singularity strength iB)n }:

[B2B] {}, + [W2B]{W1,} = {B.Cbody} (5.5)

83
 Since, {W} is prescribed for the current step, the linear system to be solved is:

[B2B] {Bn} = {B.Cbody}- [W2B] {Wn} . (5.6)

2. Determine how well the Kutta condition is satisfied using the current wake strength,
but the updated body surface singularity distribution:

[B2W] {Bn} + [W2W] {Wn} - {T.EKutta}= {R}. (5.7)

3. If the residual R from equation 5.7 is less than some given tolerance, stop the iteration
loop, else go to 4.

4. Update the wake strength based on the Kutta condition, and the current solution. For
the non-linear pressure Kutta condition, this amounts to computing the change in
potential based on the current trailing edge pressure jump:

[ (pu- P-P)i
]{A=Wn} {Z(PU P)n}. (5.8)

Solving this linear system, the resulting wake strength 1Wis:

Wn+l= fn + IWn (5.9)

5. Go to step 1. Repeat until converged.

The above weak iteration can be used for both linear and non-linear equations. For linear
equations, the iteration will take a single step to converge, however, for non-linear Kutta
conditions the approach may take several steps to converge, depending on how well the
Jacobian matrix in equation 5.8 is represented.

84
5.2.2 Doublet Sheet, Vortex Sheet and Vortex Particle Wakes

Consider a discrete dipole sheet representing a wing trailing wake. The discrete doublet
wake can be equivalently represented using equivalent vortex sheets and vortex rings. For
piecewise continuous doublet representations, the equivalent vortex ring surrounding the
doublet panels is nullified by the neighboring panels. As such, the use of higher order,
piecewise continuous doublet wake panels is preferred. In order to convert the dipole panels
to equivalent vortex particles the following steps are taken (see figure 5-1 ):

1. Determine the equivalent vortex representation for each dipole sheet to be converted
to vortex particles. The equivalent vortex sheet is computed by examining the tan-
gential derivative of the potential distribution.

2. Determine the number of vortex particles to be emitted from each panel.

3. Divide the panel into equal area segments and generate a vortex particle at the cen-
troid of each of the segments.

(a) The particle strengths are computed by multiplying the area represented by the
panel segment by the strength of the vortex panel.

(b) The vortex particle radius is determined using the maximum radius of the panel
segment and adjust it by a user specified factor.

5.2.3 The Farfield Approximation Model

Two levels of farfield approximation exist:

1. Fixed in space particles: When vortex particles travel sufficiently far away from the
bodies in the domain their exact position becomes less important. As such, once the
vortex particles have reached a sufficient distance from the body, they can be frozen
in space. This effectively reduces the inter-particle computations since no induced
velocities are required for the farfield-fixed-in-space vortex particles. The vortex
particles which are fixed in space do however continue to move with the freestream
velocity relative to the body in question.

85
 (a) Linear Variation of (b) Equivalent vortex panel
Doublet Strength and ring representation

~,~ .~

(d) Equivalent vortex particle

representation
(c) Equivalent v0l1ex particle and
vortex ring representation

Figure 5-1: A schematic illustrating the conversion of dipole sheet panels to vortex par-
ticles. In (a) the linear variation dipole distribution is shown for a particular wake panel.
In (b) and equivalent vortex representation is illustrated. Notice, there are strong vortices
around the circumference of the panel, and there is a constant strength vorticity distribution
on the panel. In (c) the constant strength vorticity distribution on the panel is converted into
a series of vortons or vortex particles. Finally, in (d) the cancellation of the circumferential
panel vortices with the circumferential vortices of neighboring panels is illustrated. This
cancellation occurs due to the piecewise continuous linear basis representation of the dou-
blet strength. If the representation were not continuous, the circumferential vortices would
not cancel with their local neighbors.

2. Multipole Expansion Lumping Once particles reach extrel11edistances frol11the body,

their influence can be treated as a lUI11psum influence. The vortex particles are
lumped into multipole expansions and each l11ultipoleexpansion is treated as a simple
multipole particle. These multipole representations have no inter-particle influence
in the farfield; however, they do have influence on the near field vortices and the
body. Since these multipole representations have no inter-particle influence, they
simply convect with the freestream flow. This farfield simplification is made in order
to reduce the number of particles in the domain in order to l11aintaina fast sinlulation.

The farfield models are used to reduce computational complexity.

86
5.3 Putting it all together

The following section describes the solution approach implemented in the panel method
framework presented in this thesis.

5.3.1 Summary of Solution Steps

In order to simulate a body shedding wakes, it is necessary to compute velocities and use
those velocities to advect the wake. In addition to convecting the wake, the velocities due
to the wake and any other domain vorticity must also be computed as an influence on the
body boundary condition. In the current implementation this is accomplished by using a
time stepping procedure. The following general steps occur during each time step in the
solution process regardless of the potential flow formulation being used:

Step 1) If the simulation is a morphing body simulation, load the current geometry def-
inition. Solve the potential flow boundary integral equation, to determine the unknown
singularity strengths on the surface of the body. Also determine the unknown potential
jump in the new wake region for the current timestep. Any of the potential flow formula-
tions from chapter 2 can be used for determining the potential flow solution. The velocity
boundary condition should be computed based on the current body translational and rota-
tional velocities, as well as any morphing of the body surface. The current domain vorticity
influence, U,, is also included in the boundary conditions. Included in this flow solution
is the determination of the potential jump across unknown buffer wake, which enforces the
potential Kutta condition.

Step 2) Determine the strength of the vorticity that needs to be released into to domain.
This new vorticity is the converted vorticity from the previous timestep known buffer wake
layer.

Step 3) Determine the velocity and gradient of the velocity influence from the body onto
the vortex particles in the wake using the boundary integral equations. This involves eval-

87
uating the gradient of the potential flow integral equation at each of the vortex particle
positions. This is a single matrix vector product.

Step 4) Determine the velocity and gradient of the velocity influence of each of the wake
particles on each of the other wake particles using the integral equation for the vorticity-
velocity relationship. This is a matrix vector product.

Step 4 a) If necessary, compute the pressures and forces acting on the body, prior to updat-
ing the position and strength of the wake. This is an application of the unsteady Bernoulli
equation, to determine the pressure from the velocity.

Step 5) For each vortex particle in the wake, update the particle position and vortex parti-
cle strength using the Lagrangian vortex evolution equations. If the vortex particles are in
the farfield region, simply advect the particles with the local freestream velocity relative to
the body. For the nearfield particles, this step involves determining the new position and
strength by solving the ODEs governing the evolution of the vorticity.

Step 6) Compute the wake to body influence based on the new wake voterx particle po-
sitions and strengths. This is done by evaluating the matrix vector product represented by
the Poisson equation governing the vorticity-velocity relationship.

Step 7) Start over at (1) unless the iteration stopping condition has been reached.

5.4 Acceleration Algorithms

Two matrix vector product(MVP) acceleration techniques are implemented for the solution
of the potential flow. The precorrected-FFT [84] is used for accelerating the GMRES [100]
iterative potential flow solution, while the Fast Multipole Tree approach is used for eval-
uating the vorticity stretching and velocity influences. Following the introduction of the
acceleration methods a brief discussion of the use of the methods is presented. Although

88
the acceleration methods are briefly outlined in the following two sections, more detailed
information can be found in the original literature [110, 109, 70, 86].

5.4.1 precorrected-FFT

The precorrected-FFT method used is based on the original p-FFT algorithm [84]. Through
a collaborative effort, a generic version of the pFFT algorithm was implemented[47]. In
the pFFT algorithm, the matrix vector product:

[A]5
-= (5.10)

is accelerated using the following steps:

1. Projection: Through a polynomial interpolation procedure the charge distribution r

is projected onto a 3-Dimensional Fast Fourier Transform Grid. In order to make the
process generic, the projection operator takes the moments and centers of discrete
charge representation and projects those moments onto the FFT grid. By taking
moments and centers, the use of arbitrary charge sources is possible. The process is
summarized as:

[P]~ = 59 (5.11)

where, P is the projection matrix and ag is the FFT grid equivalent charge. The
projection operator is setup once per given geometry per FFT grid.

2. FFT, Convolution, FFT - 1: An FFT is performed on the equivalent grid source

charges. Also, an FFT is performed on the kernel evaluated on the grid. The FFT is
performed using the FFTw code[48]. Once in Fourier space, a multiplication of the
transformed charge and transformed kernel is performed (the original physical space
operation was a convolution). Finally, an inverse FFT is applied to the result of the
multiplication to give the grid potentials, b9.

[H]dg = bg (5.12)

89
 3. Interpolation: The grid charges are interpolated back onto the evaluation points
through a polynomial interpolation. The interpolation is achieved using the same op-
erator as in the projection with different moments and centers possible (if the testing
basis functions are different or the inner/outer operators are different). The interpola-
tion matrix is determined by taking the transpose of the result (in reality a transpose
multiply operation is performed to reduce the number of memory accesses):

~[I]qg = X-~ ~(5.13)

4. Precorrection: Finally, the approximation of the nearfield interactions are subtracted

and replaced by directly computed interactions (using an appropriate integration
scheme).
([D]- [I][H][P]) -= PC (5.14)

The end result is a matrix vector product which is represented as:

[I][H][P]
+ ([D]- [ili[HP]) = 0 (5.15)

The matrices [I],[P] and [D] as well as the Kernel transform a need only be computed
once per given fixed geometry. The matrix vector product is rapidly computed during the
iterative solution process by multiplying the stored matrices. The generic framework which
was developed enabled the separation of the operators.

5.4.2 Fast Multipole Tree Method

The fast multipole tree method [110, 70] is implemented for accelerating the evaluation of
the vorticity stretching and velocity influences. The method uses a multipole expansions
but does not go into the complexity of using local expansions [46]. The following steps
highlight the FMT algorithm which is implemented:

1. Octree Domain Decomposition: The following steps are used to setup the octree
and multipole representation:

90
 (4)

Figure 5-2: The pFFT algorithm presented graphically. First a uniform 3D-FFT grid is
generated. The domain charges are projected onto the grid. Following that the convolution
is performed in Fourier space to yield the grid potential. Following that the grid potentials
are interpolated back to the evaluation points, and lastly the nearby interactions computed
directly and subtracted from the grid based computation.

(a) Create a root cell enclosing the entire geometry under consideration

(b) Get NumElements, the number of elements in the current cell

(c) If NumElements > J11axNu1nElements then split current cell into 8-children

cells and cycle over each of the 8-children cells by going to step b.

(d) If NumElements < J11axN'umElernents then construct leaf cell multipole

representation by clustering elements in current cell into a multipole mOlnent

expansion at the centroid of the current cell. Go to b. and continue decomposi-
tion with the next cell in the list.

(e) Translate and add the current cell multipole expansion to the parent cell multi-
pole expansion. The parent cell is the one level higher cell to which the current

91
 child cell is associated.

2. Potential Evaluation: The potential evaluation is the second step in the multipole
process. The following steps are taken to evaluate the potential at a point R =

(x, y, z):

(a) Start at the root cell level:

(b) Evaluate the distance, Rpr between the cell centroid and the evaluation point.
Compare this with the cell sidelength S. Compare the separation to sidelength
ratio with a user specified minimum:

K = Rp/S (5.16)

(c) If K < Kser: The current cell multipole representation is too close to the
evaluation point. If the current cell is not a leaf cell, go to the children cells
of the current cell, cycling through each starting at step (b). Nested loops can
continue multiple levels until each interaction is accounted for. If the current
cell is a leaf cell, directly evaluate the interactions between the elements in the
current cell and the evaluation point.

(d) If K > Kser: Compute the potential due to the current cell to the evaluation
point:

i. If the current cell is a leaf cell, the potential is computed by cycling through
the elements in the cell.

ii. If the current cell is not a leaf cell, compute the potential at the evaluation
point using the multipole expansion expressions in the cell.

In the current implementation of the Fast Multipole Tree approach, a heirarchical distance
specification is used. One specifies a center of the domain and one can specify multipole
tree evaluation properties in different zones of the domain defined using increasingly larger
spherical zones. For example, near the body, one might wish to ensure that the majority
of interactions are directly computed, while far from the body, the more frequent use of

92
 Farticld Region
Farticld Region

Farlield Region
Farticld Region

Figure 5-3: The domain is divided into several concentric spherical domains. In each do-
main, the minimum separation distance between evaluation point and cell center can be
specified. This allows for high resolution of the interactions in certain regions of the do-
main, while other regions can have a more approximate, more rapidly computed interaction.

multipole expansions can be considered. This is shown in figure 5-3 The various controls
on the Fast Multipole Tree acceleration in the current code are as follows:

1. Octree Decomposition Setup Controls

(a) Buffer for the bounding box. The bounding box can be enlarged by a sInall
aInount to encolnpass additional space if necessary.

(b) Multipole Order. Second order is maximmn expansion order.

(c) Maximum number of elements per octree cell. The number of elelnents per cell
enables the octree decOInposition. It has been found through trial and error that
10-20 elelnents are the best cOInpromise between accuracy and cOInputation
time.

2. Evaluation Controls

(a) Minimum cell separation for farfield multipole approximation. This is based
on the ratio of cell-sidelength to the distance between the cell-centroid and

93
 element-centroid. Typically, separation distances of 2.5 cell sidelengths are
a good compromise between accuracy and evaluation time.

(b) The farfield radius. The farfield radius dictates where the motion of particles
is based only on the freestream velocity relative to a fixed in space body. The
magnitude and direction of the vorticity associated with each of the particles is
not adjusted.

(c) Nearfield Subdomian Radii. The nearfield subdomain radii are used to deter-
mine the minimum cell separation value based on the distance away from the
body.

5.4.3 A Brief Discussion of the Acceleration Routines

The pFFT and FMT approaches have different strengths for the problem under consider-
ation. For this reason, both of the approaches are used for their appropriate application.
The pFFT is used for computing the matrix vector product inside of the GMRES iterative
solution method, while the FMT is used to compute the matrix vector products for one off
computations (such as the velocity and stretching computations). This is discussed briefly
here.

1. The pFFT is used for the solution routines, since, once the matrices are set up for
a given geometry the matrix vector product is extremely efficient. Since the Fast
Multipole Tree method which is implemented only considers an upwards sweep of
the Octree structure (no local expansions), the matrix vector product is less efficient.

2. The pFFT is not used for the wake interactions, since the wake may extend a long way
behind the body. In order to accurately approximate the wake, a large equi-spaced
grid must be used. The FMT is more appropriate due to the hierarchical nature of the
octree grid. This permits rapid evaluation of the influence.

3. The pFFT is not used for the wake particle to particle interactions due to the expense
of precorrection. The precorrection cost becomes significant component of the setup
time when the direct computations are simple point to point interactions.

94
 4. The FMT is used for computing the majority of the velocity influences, since the
pFFT would typically require a larger direct interaction stencil and/or a larger pro-
jection and interpolation stencil and order. This increases the direct and precorrection
components of the matrix vector product setup significantly. Since the FMT can be
tuned, it is ideal for computing the velocity and stretching influences.

5.5 Implementation Details of the BIE Formulations

In this section a brief discussion of the three boundary integral equation formulations which
were presented in Chapter 2 are discussed. The application of these formulations for solv-
ing the body surface singularity distributions is examined more closely in the following
chapter; however, for completeness they are presented here.

1. Source-Doublet Formulation with Potential Doublet Buffer Wakes: The source-

doublet formulation used for solving the potential flow involves the use of sources
and doublets on the surface of the body, with a dipole buffer wake sheet used for
satisfying the linear Kutta condition. The no flux boundary condition is satisfied
through the use of the source distribution. In addition to the body motion and morph-
ing velocities, the vortex particle wake velocity is also used as a boundary condition
on the potential solution. The time stepping procedure is similar to that presented
earlier.

2. Doublet Lattice Neumann Formulation with Velocity Inducing Doublet Buffer

Wakes: The doublet lattice formulation is similar to a traditional doublet lattice
solver except for the linear order doublet strengths and the use of vortex particles
in the wake rather than a traditional wake sheet. The Kutta condition is the simple
linearized no trailing edge vorticity Kutta condition.

3. Source Neumann Formulation with Body Piercing Doublet Wakes: The imple-
mentation of the source Neumann formulation with body piercing wakes is presented
in a little greater detail that the previous two formulations. The surface singularity in

95
 this formulation is a source distribution of piecewise linear strength. The source dis-

tribution is treated as a velocity influence. The wakes are represented using doublet

sheets which pierce the body at the trailing edge. For a single, silllple wing, Figure

5-4 illustrates the piercing doublet wake. The doublet wake sheet is treated as a ve-

locity influence. The combined velocity influence must satisfy the surface boundary

condition (which is used to determine the source strengths) and the Kutta condition

(which is used to determine the doublet strengths). The non-linear Pressure Kutta

condition described earlier is implemented.

TE-Wake
Intersection

,
_'I(
Linear Variation To Internal I External
Zero Strength ButTer BuOer

Figure 5-4: A schematic illustrating the body piercing wakes demonstrates the linear vari-
ation of the internal to the body wakes. For effective computations, it has been found that
a linear distribution of doublet strength inside of the body results in sufficiently SIllooth
source distributions on the surface of the body. For efficiency the source-Neumann formu-
lation is the approach considered in this work (note: the body piercing wake concept can
also be used for source-doublet type approaches).

5.6 Conclusions

The various implementation details for the potential flow solution framework have been

presented in this chapter. In the chapter which follows simulation results and examples are

96
presented for the various potential flow formulations investigated.

97
98
Chapter 6

Simulations and Results

This chapter demonstrates the solution of a selection of relevant sample potential flows
using the formulations presented in this thesis.

6.1 Validation Computational Experiments

6.1.1 NACA-0012 Non-Lifting Wing

A rectangular wing with a symmetric NACA 0012 airfoil at zero degrees angle of incidence
is examined in this section. The discretized wing comprises 4096 triangular panels and
2145 vertices. The wing is a lattice of ordered right angle triangles (32 panel widths in
the spanwise direction by 64 panel widths in the chordwise direction) with leading edge
refinement. Although an analytical solution is not known, the NACA 0012 airfoil is a
common test case and as such, the comparison to computational data from the 2-D panel
method XFOIL [38] is made. The surface pressure distribution for the source-Neumann
formulation and the Source-Doublet Potential formulation is illustrated in figure 6-1. Since
the flow is at a zero incidence angle, the wake sheet strength is identically zero.

99
 c

xlc

Figure 6-1: The results are illustrated for both the Neumann-Source approach and the
Source-Doublet approach. Comparison is made to the 2-D results from the panel method
code XFOIL [38] and show a good match.

6.2 Steady NACA 0012 Lifting Wing Simulations

In order to determine the steady state solution to a given lifting problem, there are two
options. First, the solution of a potential flow problem with a very long trailing wake. This
solution is determined in one step and does not require iteration in time. The second ap-
proach is more costly yet more accurate. It involves solving for the steady state conditions
by iterating some startup flow up to steady state. The second option, though more costly is
chosen for this example. The lifting body wing centerline pressure distribution is presented
in figures 6-2 and 6-3. Figure 6-2 illustrates the pressure distribution around the midsection
of a NACA 0012 with Aspect Ratio of 4. Since bodies of thickness are presented, only the
source-doublet and source Neumann-formulations are considered. The results demonstrate
good agreement with each other, as well as a good agreement with a reduced angle of at-
tack 2-D simulation. In addition to the thick body simulation results, a plot of the wing
centerline pressure difference due to a flat plate with an aspect ratio 4 wing is presented in
6-3. The result demonstrates the large leading edge contribution to the pressure due to the

100
strong leading edge vorticity in this case.

-0.8 .. ... ... .. .. ......... ...... ource eumann,AoA= 2.5

: ' ~b;+
SourceDoublet,
AoA= 2.5
-0.6 :_ ..................
........................ \ - XFOIL,AoA= 1.650 .........
0
..... XFOIL2D,AoA
= 2.5

0.2

0..
0.2 . .... .. . ... ..

0. . ..

-1 -0.5 0 0.5 1 1.5 2

Figure 6-2: The pressure distribution around the midsection of an Aspect Ratio 12 finite
wing with NACA 0012 profile at a zero angle of attack

6.3 Unsteady Computations

In this section results for wings under unsteady motions are examined.

6.3.1 Unsteady Startup Flow

The simulation of an unsteady startup flow is presented in this section. The wing in this
example is subject to a unit step velocity increase, such that the acceleration of the wing
at time t = 0 is an impulse function. As a result of the sudden start up, several interesting
flow features can be seen. First, an infinite force will act on the body due to the sudden
acceleration of the body (in the discrete simulations, this force is large, not infinite). A
starting vortex is formed due to the satisfaction of the Kutta condition at the initial startup.
As the distance between the starting vortex and wing increases, the effect of the downwash

101
 a

x/c

Figure 6-3: The pressure distribution at the midspan of a Aspect Ratio 12 membrane wing.

from the starting vortex is reduced. As a result the lift will gradually rise and after a
long time, settle at a steady state value. Finally, since the wing is three-dimensional, one
can readily see the wing tip vortices and the wing tip vortex roll-up. The startup flow was
examined for three different aspect ratio wings (Aspect Ratio 4, Aspect Ratio 8, and Aspect
Ratio 12), each with a NACA 0012 profile, except for the membrane which is an infinitely
thin plate with no camber. The non-membrane wings are discretized using 4096 triangular
panels and 2145 vertices, while the membrane wing is discretized using 2048 panels and
1089 vertices. In these simulations, there is significant wing leading edge refinement in
order to accurately capture the large gradients of singularity strength at the leading edge.

Figure 6-4 shows the time evolution of the Z-component of the resultant force for the
wing at an incidence angle of 2.5 degrees. The results are compared with data extracted
from a vortex lattice startup flow computation presented in Katz and Plotkin [113]. Figure
6-5 shows the resulting vortex particle trail behind the wing. The unsteady startup flow
demonstrates that the methods developed in the solution framework accurately model time
dependent flows. The unsteady lift increase with time compares well with the results pre-
sented in [113]. The only noticeable difference between the results is the early transient

102
 45
0. r- .. ....... .. ...... .......... ...... .... ....... ...... ..... ..........

AR12-M
0.4, . - --.AR12-G.......
* AR12-K.P.
-e-AR8- M
0 .35 .................. ... ......... . ... .............. AR - G
AR8-N
* AR12-K.P.
0.3H·.......
0.3 ··- ·· · ·-- · ··· ; ··· ....
··- ··. . -: · ·· · · ~·· · · s AR4-M
· -.........-.
AR4 - G
AR4-N
AR12-K.P.
0.25

0.2 ......
. .. .....

0.15

0.1. . ... .. .. ........... ............

.....

0.05 . . .........
........
.......
-- i

0 1 2 3 4 5 6 7 8
0 1 2 3 4 5 6 7 8
UT/c

Figure 6-4: The figure illustrates the time development of the Z-Component coefficient
after a sudden step function velocity change. The wing is instantaneously set to 2.5 degrees
angle of attack at time t=O. The results are plotted for an aspect ratio of 4,8, and 12. In
addition results are plotted for the membrane solution approach (M), the source doublet
approach (G), and the Source-Neumann approach (N). The plot also illustrates the results
presented in Katz and Plotkin [113].

time response of the Neumann-Body-Piercing Wake formulation result. These deviations

between the Neumann Source formulation code and the other formulations are attributed to
the satisfaction of a nonlinear pressure Kutta condition rather than a linear Kutta condition.

6.3.2 Unsteady Finite Aspect Ratio Heaving Wing

The unsteady heaving wing problem provides an opportunity to validate the unsteady time
dependent flow solution capabilities of the panel methods examined. The heaving wing
demonstrates the ability to compute unsteady flows in which wake-body interactions play
an important role in the overall production of forces. The heaving wing simulations are
compared to analytical expressions for an oscillating 2-Dimensional plate which are pre-
sented by Theodorsen [104] and Katz and Plotkin [113].

103
 -001

(a) The vortex particle representation of the vor- (b) The vortex particles trailing an aspect ratio 8
ticity behind an aspect ratio 8 wing after a sudden wing after a sudden unit step velocity startup. The
startup. The figure illustrates the vorticity distri- figure clearly demonstrates the startup vortex roll
bution in the wake. As expected, there is a strong up as well as the commencement of the wing tip
concentration of vorticity trailing the wingtips, vortex roll up.
which eventually forms two vortices. In addition,
the figure clearly illustrates the startup vortex, as
well as the gradual increase in vorticity due to the
sudden startup of the wing.

Figure 6-5: The startup flow around an aspect ratio 8 wing.

Theodorsen Theory For 2-Dimensional Flat Plate Heaving Motions [104]

In the Theodorsen theory, oscillations are assumed to have a small amplitude. The Z-

Component of force produced during heaving oscillations of a 2-Dilnensional flat plate is

expressed as [104]:
2
dh c d h)
Fz = -npc ( UC(k)-dt + -- (6.1)
4 dt2 '

The first part of equation 6.1 represents the generation of traditional steady state lift mul-

tiplied by a lift reduction factor C (k), while the second part of the equation represents the
unsteady acceleration component of the lift. The lift reduction factor C (k) can be written
as a function of the reduced frequency k as[1 04]:

C(k) = F(k) + iG(k). (6.2)

104
The expression for the real component, F(k) is:

F(k)- J(k)(J(k) + Yo(k))+ Y(k)(Y(k)- Jo(k)) (6.3)

2 +(Y1- Jo(k))
(J(k)+ yo) 2

while, the expression for G(k) is:

(k)
- Y (k)Yo(k) + J1 (k) Jo(k)
(6.4)
(Jl(k) + Yo)2 + (Y1 - Jo(k))2'

A plot of the circulatory lift reduction factor and phase shift due to Theodorsen's lift defi-
ciency factor is presented in figure 6-6. The reduced frequency k is defined as:

Figure 6-6: The Theodorsen lift deficiency function as a function of the reduced frequency.
The function is presented in terms of the magnitude (C(k)) and the phase shift in radians.

k= (6.5)
2U'

The Theodorsen lift deficiency factor C(k) is a direct result of the downwash on the main
wing produced by the unsteady wake development. The unsteady release of vorticity causes
a reduction in the circulatory lift while also causing a lag in the development of the lift.

105
Katz-Plotkin Theory [113]

Katz-Plotkin derived the following similar expression; however, considered different as-
sumptions [113]:

L = -rpc
(U dth + 3cd
4 dt
h\
2
2 (6.6)

To arrive at equation 6.6, Katz and Plotkin [ 113] neglected the effect of the wake downwash

on the flat plate. Although this simplifies the derivation, this results in significant deviations
from the Theodorsen [104] result and the computational results.

Simulations of Heaving Wings

Simulations are performed for several different frequencies. Each simulation has a heaving
velocity defined by:
d(t= -0.05sin(wt). (6.7)
dt
The computational results are presented and compared with the Theodorsen Theory [104]
and the Katz and Plotkin theory [113] in figures 6-7 - 6-9. The results for the heave oscil-
lations compare more favorably with the theory presented by Theodorsen [104] than Katz
and Plotkin [113] as would be expected. This is due to the panel method accurately repre-
senting the unsteady vorticity in the wake. As the heave oscillations increase in frequency,
the importance of more accurately capturing the wake downwash effect on the heaving
wing is apparent.
The computational results demonstrate slight differences when compared to the theo-
retical predictions. There are several reasons for this:

1. The computational results consider a finite aspect ratio wing. As such, the com-
putational results expectedly demonstrate a slightly lower Z-force prediction due to
the 3-Dimensional nature of the flow. This effect is clearly captured in the lower
frequency oscillations where tip vortices play a significant role in Z-force reduction.

2. The computational results consider several different Kutta conditions, whereas, the
Theodorsen theory [104] considers a finite-trailing-edge-velocity Kutta condition.
As such, the computational results would be expected to differ from the theoretical

106
 0.4

0.1

-0.1

-0.3

-0.4
o 10 12 14 16 18 20
Time

Figure 6-7: The plot illustrates the z-Component Coefficient of Force resulting from heav-
ing oscillations at a reduced frequency of ~. Notice, at this higher reduced frequency, the
unsteady component of the Forces begins to have an effect on the overall forces. This
can be seen from the phase shift which appears between the velocity and the resulting lift.
Computational results are presented for the source-doublet, the doublet Neumann mem-
brane lattice, the source-Neumann wake piercing formulation with body piercing wakes
(for both a steady and unsteady pressure Kutta condition). Comparison is made to theories
presented by Theodorsen [104] and Katz and Plotkin [113].

predictions, especially at higher frequencies. For example, in figure 6-7 results have
been presented for steady and unsteady, linear and non-linear Kutta conditions. It can
be seen that differences exist between the various simulations. The Inost notable dif-
ferences occur when the Kutta condition considered differs from the one considered
by Theodorsen.

The simulation results for the heave oscillations demonstrate the versatility and ease
of use of the cOlnbined panel method-vortex particle method approach. In figures 6-10
vorticity wake structure plots are presented for both fast and slow heave oscillations. The
plots illustrate both the vorticity structure as well as the vorticity evolution with time.

107
 0.4
-<>- Source-Doublet
-+- Membrane
-.0.- Neumann
0.3 - Theodorsen
-.- .. Katz-Plotkin
--- Velocity

0.2

" ,"
'..-- ........ \.
i
-0.1
., .
i,
"

i
.,
-0.2 ...,
.,.,
-.. I
-0.3 .- ,-;,"

-0.4
o 10 15 20 25
Time

Figure 6-8: The plot illustrates the z-Component Coefficient of Force resulting from heav-
ing oscillations at a reduced frequency of ;0'
At this frequency, the Z- force is of opposite
phase to the velocity as would be expected in lower frequency motions. Computational
results are presented for the source-doublet, the doublet Neumann membrane lattice, the
source-Neumann wake piercing formulation with body piercing wakes (for a steady pres-
sure Kutta condition). Comparison is made to theories presented by Theodorsen [104] and
Katz and Plotkin [113].

6.4 Wing-Body Simulation Example

In this section, a practical simulation of a wing-body configuration is demonstrated. The

example demonstrates the ability to reduce the problem setup time through the use of the
novel approaches presented in this thesis. The use of a body piercing wake allows the
wing-body mesh to be imported from a standard mesh generation tool.

6.4.1 Rigid Body Piercing Wakes

The first example presented demonstrates the use of rigid body piercing wakes. The ge-
ometry under consideration is presented in figure 6-11. The wing trailing edge vertices are
prescribed, and from that the rigid wake geolnetry is constructed. The wake is extended
into the wings through the trailing edge, as well as through the fuselage. The results of the
wing body rigid-body-piercing-wake-only simulations are presented in figure 6-12. The re-
sults presented in figure 6-12 demonstrate the body-piercing wake approach for a practical

108
 0.4
-<>- Source-Double!
-+- Membrane
-....- Neumann
0.3 - Theodorsen
--" Katz-Plotkin
--- Velocity

0.2

o~ 0

-0.1

-0.2

-0.3

-0.4
o 10 15 20 25 30 35 40 45 50
Time

Figure 6-9: The plot illustrates the z-Component Coefficient of Force resulting frOln heav-
ing oscillations at a reduced frequency of ;0'
At this frequency, the Z-force is of opposite
phase to the velocity as would be expected in lower frequency motions. Computational
results are presented for the source-doublet, the doublet Neumann membrane lattice, the
source-Neumann wake piercing formulation with body piercing wakes (for a steady pres-
sure Kutta condition). Comparison is made to theories presented by Theodorsen [104] and
Katz and Plotkin [113].

geometry.

6.4.2 Combined Body Piercing Wakes and Vortex Particles

A second example is considered in which both body-piercing wakes and vortex particles are
used in the same simulation. The body piercing wakes are used to represent the vorticity
at the trailing edge of the wing for the timestep under consideration. The wake is then
converted to vortex particles in preparation for the next timestep. Results are presented in
figure 6-13. The results demonstrate the pressure distributions as well as the vortex particle
wake evolution. Tip vortices are seen forming downstremn of the lifting surface. The
effective lift of the wing-body configuration at 3 degrees angle of attack was LeI I = 0.618,
while an aspect ratio 6 wing alone had an effective lift of LeI I = 0.629 and an aspect
ratio 5.2 wing was found to have an effective lift of LeI I = 0.516. The lift coefficient
for the wing-body was determined to be CL = 0.168, while for an aspect ratio 6 wing

109
 10 25

(a) An illustration of the vortex particle structure (b) The vortex particle structure behind a rapid Iy
behind a slowly heaving wing. The appearance heaving AR8 wing. The figure illustrates the vor-
of distinct vortex rings is clear. These rings are a tex particle structure behind a rapidly heaving
direct result of the heaving motion during which wing. The appearance of distinct vortex rings is
the Z-component of the forces oscillates between clear. These rings are a direct result of the heav-
a positive and negative value. ing motion during which the Z-component of the
forces oscillates between a positive and negative
value.

(c) An illustration of the vortex particles in the (d) An illustration of the vortex particles in the wake
wake trailing a slowly heaving wing. The wake trailing a rapidly oscillating wing. The vortex wake
roll-up during the oscillation can be seen. self influence is seen to further deform the wake as the
wake propagates downstream.

Figure 6-10: Plots of the vortex wake structures behind oscillating wings.

the lift coefficient was CL = 0.171 and for an aspect ratio 5.2 wing CL = 0.162. The
results illustrate successful fuselage lift carryover is achieved with the body-piercing wake

110
 0.' o.
02 02
o o
-02 -02
-0 •
-0 •
-3 -2

(a) Side View. (b) Front View.

(c) Top View.

Figure 6-11: The geometry considered for the wing-body example.

formulation.

6.5 Example: Flapping Flight

In this section an example of the application of the panel method framework to a flapping
flight application is presented.

111
 O't
02

:~~ _.- :.wu .._n ~

7 8
-005

(a) Surface potential 3D view. Notice the body (b) Surface potential side view. Notice the body
piercing wake approach captures the discontinu- piercing wake approach captures the discontinu-
ous potential along the fuselage when the overall ous potential along the fuselage when the overall
potential cp is determined. Since the body pierc- potential cp is determined. Since the body pierc-
ing wakes are represented as a velocity influence ing wakes are represented as a velocity influence
the solution of the body surface potential cp B is the solution of the body surface potential cp B is
continuous. Note that the surface potential repre- continuous. Note that the surface potential repre-
sentation is based on the potential at the centroid sentation is based on the potential at the centroid
of each of the panels, therefore, the wake poten- of each of the panels, therefore, the wake poten-
tial jump does not appear along the wake shed- tial jump does not appear along the wake shed-
ding line. ding line.

02

01

(c) Surface pressure plot. Notice here that the (d) Surface pressure plot. Notice here that the
pressure is continuous. pressure is continuous.

Figure 6-12: Results for the surface potential and surface pressure considering a rigid body-
piercing wake representation.

6.5.1 Flapping Wing Example

The method presented by Hall et al. [143] is used to determine the minimum power wake
vorticity distribution for a particular flapping wing lifting line. In this case a simple hinged
112
 (a) The surface pressure distribution as viewed (b) The surface pressure distribution as viewed
from above. from below.

-2 -3 -3

(c) Surface pressure plot including vortex particle (d) Surface pressure plot including vortex particle
wakes. wakes.

Figure 6-13: The geometry considered for the wing-body example.

flapping motion is considered. Considering the lifting line flapping motions, the wake im-
posed downwash and the resulting minimum power vorticity distribution, a thin melnbrane
flapping wing geometry representation replicating the minimum power wake vorticity dis-
tribution is constructed as follows:

1. Construct a reference wing platform which flaps with the correct parameter depen-
dence while shedding the least amount of vorticity possible into the domain. The
wing can be designed to have a zero vorticity shedding throughout the wingbeat by
instantaneous local modifications of section camber and angle of attack. In this ex-

113
 periment, however, the wing is assumed to have a zero camber and the wing section
is aligned with the flow at the leading edge. This assumption will produce a reference
wing which has non-zero vorticity, especially in regions near the tip where incidence
angles are important; however, as a first order approximation for demonstrating the
translation process this is adequate.

2. Compute the downwash at lifting line positions along the wake considering only the
vorticity downstream of the lifting line.

3. Compute the local section angle of attack for the entire wing beat at each spanwise
section of the wing. This is computed considering the optimal vorticity distribution
in the wake from the Hall et. al. simulations.

4. Adjust the reference wing platform for the local downwash effects due to the wake
as well as the local sectional angle of attack. The resulting geometry will have sig-
nificant twist as the wing span is traversed during the flapping motions.

5. Mesh the flapping wing as either a thick body representation (using an airfoil profile)
or as a thin membrane. In this paper the thin membrane is considered due to the ease
of meshing the geometry (especially near flapping joints).

The membrane doublet lattice code is used for these simulations due to the relative
ease of meshing the morphing geometry with a thin surface approximation. The minimum
power vorticity wake result from the Hall et al.[143] simulation approach is presented in
figure 6-14(a). The panel method simulation results based on this result are presented in
figures 6-14(b) - 6-14(d). Although the wake resulting from the panel method differs in
fine details from the desired vorticity distribution, the result of the panel method simulation
demonstrates a good overall agreement between the target vorticity wake structure and the
resulting simulated vortex particle wake structure. In addition, the panel method simulation
demonstrates the effects of wake roll-up during wingstroke segments where large vorticity
release is expected (such as during the downstroke portion of the flapping cycle).

114
 0.6

.,.........-- - - ' '\ ,~" nm

I
r

-~""'"--_.---~ .
T--"-'
-1_1-_
3 2.5
.- " J
'5 1 05
a 05

(a) The reference wake from the wake only mini- (b) A top view of the resulting FastAero simula-
mum power vorticity solution. Flight path is from tion of a flapping geometry. Flight path is from
right to left. right to left.

02
0._
-.j
-0.2 -.• ._

1.5
--r.--._ . .......-<:.
......
< -0.4
1 r- _ >-- -< -02
05 05 " a
. 02
O'

(c) An illustration of the vortex particle wake (d) An illustration of the vorticity in the wake
trailing a flapping wing, simulated using Fas- trailing a flapping wing, simulated using Fas-
tAero. tAero.

Figure 6-14: A demonstration of the use of a wake only optinlal vorticity distribution result
to construct an efficient three-dimensional flapping wing geometry for simulation in an
unsteady panel method solver.

6.6 Some Comments on Solution Time and Accuracy

The methods presented in this thesis have all been accelerated to scale nearly linearly with
the nUlnber of discrete elements. Despite the acceleration of the methods, the computa-
tional burden can still be extensive. In this section the computational cost is discussed
briefly. Although the current code is not completely optimized in terms of computational

115
structure and implementation, the discussion of the computational times will provide in-
sight into the relative cost of each of the three panel method formulations.

6.6.1 Doublet Neumann Formulation

The doublet Neumann formulation is consistently the fastest approach for simulating un-
steady lifting surface flows. This is due to the following reasons:

1. The number of surface elements needed for a particular simulation is typically one
half of that needed for a similar thick body lifting surface simulation.

2. The problem involves the setup of a single linear system matrix. This results in a
single accelerated system being setup per geometry.

3. The Aerodynamic Influence Coefficient Matrix is typically very well conditioned due
to the strong diagonal dominance of the matrix.

4. The implementation of a linear Kutta condition means that the solution for any given
timestep is known after a single solve. This does not involve additional solves to
refine iterations.

5. The forces, moments and surface pressures can all be computed based on the dou-
blet strength solution and surface gradients of the doublet strength. As a result, no
forward evaluation of the integral equations are necessary.

As a result, the Doublet-Neumann method yields rapid solution times. For example, the
unsteady rigid wing simulations in this chapter were all simulated in less than one hour on
a 3.2 ghz Intel Pentium, 2GB RAM memory personal workstation.

6.6.2 Source Doublet Formulation

The source-doublet formulation is the next fastest approach for simulating the unsteady po-
tential flow. The advantages of this approach are similar to the membrane doublet Neumann
code with a few differences:

116
 1. The rigid body problem involves a single setup of two linear system matrices. This
results in a set of accelerated systems being setup per geometry.

2. The implementation of a linear Kutta condition means that the solution for any given
timestep is known after a single solve.

3. The forces, moments and surface pressures can all be computed based on the dou-
blet strength solution and surface gradients of the doublet strength. As a result, no
forward evaluation of the integral equations are necessary.

The source-doublet formulation performs solutions in a comparable time as the Membrane

Doublet Neumann code, however, is slower due to the comparably reduced conditioning of
the linear system matrix and the increased complexity of thicker bodies. The thick wing,
unsteady simulations were computed in under an hour.

6.6.3 The Body Piercing wake Formulation

The body piercing wake formulation is by far the most costly of the implementations. This
is due to several reasons, but most notably:

1. The body piercing wake formulation involves and iterative pressure Kutta condition.
Although a time savings is realized in using an approximate Jacobian, there is still
significant computational cost in solving the linear system multiple times to deter-
mine the wake strength and surface singularity distribution.

2. The post-processing of the solution involves the setup and evaluation of an integral
equation for both the body induced velocity and the body induced potential. In a
future version of the method these matrices should be setup on a one-off basis for
rigid body simulations. The setup and forward evaluation of the integral equation for
post-processing purposes imposes a significant time penalty on the method.

Although the source-Neumann-body-piercing-wakes formulation presents increased solu-

tion times, it should be noted that the reduction of the solution setup time and grid gener-
ation complexity should easily offset the cost of the solution. In addition, when compared

117
on an even level with the other methods (rigid body simulations implementing non-linear
iterative pressure Kutta conditions) the source-Neumann approach will not be significantly
unfavorable.

6.6.4 Rigid Vs. morphing Bodies

The computational framework has the capacity to simulate both rigid bodies as well as
morphing geometry bodies. The two options prove to have significantly different solution
times. For the rigid body, the body system matrices for the are setup once. At each iteration
the solution involves solving for the particular boundary conditions on the given geometry.
In the case of the morphing geometry problem, the solution is performed on a different
geometry at each step. In the current implementation this means that the system matrix
setup time for each timestep takes a significant portion of that particular step's time. As
a result, morphing body geometries are typically significantly more costly that their rigid
body counterparts. As a result, future work will involve reducing the linear system setup
time in the pFFT algorithm through thorough optimization of the computational routines.

6.6.5 Comments on the use of Galerkin linear basis approaches

In each of the formulations implemented in this thesis, a Galerkin, linear basis approach
was considered. The additional computational complexity due to the high order galerkin
method compared with a lower order constant collocation approaches is seen more pre-
dominantly at the linear system setup stages of the computation where an approximately
six-fold increase in computational time is seen. This is due to the following two reasons:

1. The linear basis function panel integration costs are approximately two times more
expensive than the constant basis counterpart.

2. The Galerkin testing involves a quadrature integration which at minimum comprises

the panel integration at 3 positions on the target panel. When compared with the
single evaluation point in a constant collocation approach, the Galerkin method is at
least three times more expensive.

118
Although the linear Galerkin approaches are more time consuming than constant colloca-
tion, the benefits of increased convergence and more refined solution representation out-
weigh the disadvantages.

6.7 Conclusions

in this chapter several example simulations highlighting the panel method formulations
implemented have been shown. The methods are all in good agreement with available
theoretical, computational and analytical results.

119
120
Chapter 7

Conclusions

This thesis presents a panel method which addresses the drawbacks with current panel
method implementations.

7.1 Contributions

Several contributions have been made in this thesis.

7.1.1 The general panel method framework

The panel method which has been developed in this thesis represents a synthesis of both
new and existing ideas, which when combined form a powerful tool for the rapid, hands off
analysis of steady and unsteady lifting body and wing-body flows. The tool is a contribution
in that it combines the most effective strategies currently available with contributions to the
field made in this thesis to provide an advanced potential flow analysis environment. The
tool which has been developed provides a framework in which potential flow simulations
can be easily performed requiring:
* Minimal user setup time due to automatic wake generation strategies.
* Minimal user expertise in wake modeling due to the use of automatic vortex particle
wake generation schemes.
* Minimal solution turnaround times due to the implementation of fast acceleration

121
approaches.
* Minimal user interference due to the hands off automatic wake strategies.
The approach relies on advanced acceleration algorithms (pFFT [109] and FMT [110])
to facilitate the rapid solution of a combined panel method-vortex-particle approach [63,
62].

7.1.2 Body Piercing Wake Formulation

The thesis presents a potential flow formulation which considers body piercing wakes
to handle the difficulties associated with wing-body geometries in both traditional panel
method approaches as well as combined panel method and vortex particle method ap-
proaches. The introduction of body piercing wakes permits wing-body simulations which
have not been possible in previous panel method-vortex particle approaches [62].

7.1.3 Higher Order Integration Approaches

This thesis presents the development of a conceptually simple, yet accurate, quadratic ba-
sis function, quadratic geometry self term panel integration approach for the double layer
integration. In addition, a desingularized projection based numerical integration approach
for the single layer self term has been implemented to form a quadratic basis, quadratic
geometry dense system BEM solver.

7.2 Future Work

Although the panel method presented in this thesis represents a complete framework, there
are several directions of future work which will be considered:

1. The continued development of the panel method to include:

(a) Viscous diffusion in the vortex particle wake. The addition of viscous diffu-
sion models will provide accurate vorticity wake models for flows which will
be considered (such as lower Reynold's number natural flapping flight applica-
tions).

122
 (b) Boundary layer corrections through the use of integral boundary layer methods.

(c) More panel integration schemes and orders into the panel method framework.
For example, implementing the curved panel, high order integration approaches
into the accelerated potential flow framework.

(d) Several wake models of differing complexity. Due to the generality of the
Helmholtz decomposition of the velocity field, many diverse wake models can
be considered. In certain situations, a filament or sheet based wake model might
be more appropriate than the vortex particle approach. As such, incorporating
those models into the general accelerated framework will prove useful.

2. The use of the panel method to analyze complex unsteady potential flow problems
particularly suited to hands-off, automatic wake handling, for example:

(a) Rotor aerodynamics: Due to the non-diffusive nature of the Lagrangian vortex
particle dynamics, accurate models of rotor aerodynamics are possible.

(b) Natural flapping propulsion and flight aerodynamics: Due to both the rapid so-
lution of unsteady flow problems as well as the automatic hands-off simulation
capabilities the method presented in this thesis represents an attractive approach
for natural flapping propulsion and flight analysis.

(c) Induced drag reduction in formation flight: Due to the automatic wake genera-
tion and advection, analyzing formation flight drag reduction is simplified when
compared to previous approaches.

123
124
Bibliography

[1] Noh, W.F. A time-dependent two-space dimensional coupled Euler-Lagrangian

code, Meth. in Comp. Physics, Vol. 3, 95(1):115-138, 1992.

[2] Franck, R.M., Lazarus, R.B., Mixed Eulerian-Lagrangian method. Meth. in Comp.
Physics, Vol. 3, 47-67, 1964.

[3] Belytschko, T., Kennedy, J.M., Schoeberle, D.F., Quasi-Eulerian finite element
formulation for fluid-structure interaction . Proceedings of Joint ASME/CSME
Presseure Vessels and Piping Conference. ASME: paper 78-PVP-60, 1978.

[4] Hamamoto, M., Ohta, Y., Hara, K, Hisada, T., Free-flight analysis offlapping flight
during turning by fluid structure interactionfinite element analysis based on ar-
bitrary lagrange-eulerian method, Intl. Conf. on robotics and Automation, IEEE,
2004.

[5] Peskin, C.S.,FlowPatternsAround Heart Valves:A Digital ComputerMethodfor

Solving the Equations of Motion, PhD. Thesis. Physiol., Albert Einstein Coll. Med.,
Univ. Microfilms, 378:72-30, 1972.

[6] Peskin, C.S., The immersed boundary method, Acta Numerica, 2002:479-517, Cam-
bridge University Press, 2002.

[7] Mittal, R., Iaccarino, G., Immersed boundary methods, Annu. Rev. of Fluid Mech.,
37:239-61, 2005.

[8] G.K.Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press,

Cambridge, 2000.

125
 [9] Anderson, J.D. Fundamentals of Aerodynamics, 3rd. Edition, McGraw-Hill Sci-
ence/Engineering/Math, 2001.

[10] Hess, J.L., Panel methods in computational fluid dynamics, Annu. Rev. of Fluid
Mech., 22:255-74, 1990.

[11] Hess, J.L. Higher-order numerical solution of the integral equation for the two di-
mensional Neumann problem. Comput. Methods Appl. Mesh. Eng., 2:1-15, 1973.

[12] Hess, J.L. The problem of three-dimensional lifting flow and its solution by means
of surface singularity distribution. Comput. Methods Appl. Mech. Eng. 4:283-319,
1974.

[13] Hess, J.L., Smith, A.M.O., Calculation of nonlifting potential flow about arbitrary
three-dimensional bodies. J. Ship Res., 8:22-44, 1964.

[14] Ahmadi, A., Widnall, S.E., Unsteady lifting line theory as a signular perturbation
problem, J. Fluid Mech., 153:59-81, 1985.

[15] Cenko, A., PAN AIR applictions to complex configurations, J. of Aircraft, 0021-
8669, vol 20, no. 10, 887-892, 1983.

[16] Park, M.A., Green, L.L., Montgomery, R.C., Raney, D.L., Determination of stability
and control derivativesusing computationalfluid dynamicsand automaticdifferen-
tiation, AIAA-99-3136, 1999.

[17] Nathman, J. Improvement of a panel method by including panel warp, AIAA-2004-

721, 42nd Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 2004.

[18] Tinoco, E.N., Ball, D.N., Rice, F.A., PANAIR analysis of a transport high-lift con-
figuration, J. Aircr., 24:181-86, 1987.

[19] T. R. Quackenbush, D. A. Washspress, A.H.Boschitsch, Rotor Aerodynamic Loads

Computation Using a Constant Vorticity Contour Free Wake Model, Juornal of Air-
craft, vol.32 no.5, pp 911-920, 1995.

126
[20] Keuthe, A.M., amd Chow, C.-Y,Foundations of Aerodynamics: Bases of Aerody-
namic Design, 4th ed., John Wiley and Sons, New York, 1986.

[21] Katz, J., Plotkin, A., Low-speed aerodynamics: from wing theory to panel methods,
McGraw-Hill, New York, 1991.

[22] Rosen, B.S., Laiosa, J.P., Davis, W.H., Stavetski, D., Splash Free-Surface Flow Code
Methodology for Hydrodynamic Design and Analysis of IACC Yachts, 11 Chesa-
peake Sailing Yacht Symposium, Annapolis, 1993.

[23] Sclavounos, P.D., Ship flow simulation in calm water and waves, SWAN 1 Version
3.1, User Manual, Boston Marine Consulting Inc., 1999.

[24] Lamb, H, Hydrodynamics., London, Cambridge University Press, 6th ed., 1932.

[25] Kellogg, O.D., Foundations of Potential Theory, Dover Publications, 1954.

[26] Brazier, J.G., and Smith, A.M.O., Development of Nose and Tail Shapes in Incom-
pressible, Irrotational Flow, Douglas Aircraft Co., Long Beach, CA, Rept. E.S.
0875, June 1947.

[27] Eppler, R., Airfoil Design and Data, Springer Verlag, Berlin, 1990.

[28] Selig, M.S., Maughmer, M.D., A multi-point inverse airfoil design method based on
conformal mapping, AIAA-1991-69, 29th Aerospace Sciences Meeting and Exhibit,
1991.

[29] Smith, A.M.O., Pierce, J., Exact solution of the Neumann problem. Calculation of
plane and axially symmetricflows aboutor within arbitraryboundaries,Proc. U.S.
Natl. Congr. Appl. Mech., 3rd, Providence, R.I., pp.807-15, 1958.

[30] Carmichael, R.L. and Woodward, F.A, An integrated approach to the analysis and
design of wings and wing-body combinations in supersonic flown NASA TN D-
3685, 1966.

[31] Morino, L., Luo, C.C., Subsonic potential aerodynamicsfor complex configurations.
A general theory. AIAA J., 12:191-97.

127
[32] Johnson, F.T., A general panel method for the analysis and design of arbitrary con-
figurations in incompressibleflows. NASA CR-3079, 1980.

[33] Oskam,B. Transonicpanel methodfor thefull potential equationappliedto multip-

component airfoils. AIAA J. 23:1327-34, 1985.

[34] Coopersmith, R.M., Youngren, H.H., Bouchard, E.E., Quadrilateral Element Panel
Method (QUADPAN), Theoretical Report (V. 3.0), Lockheed-California, LR, 1983.

[35] Maskew, B., USAERO,A TimesteppingAnalysis Methodfor the Flow About Multiple
Bodies in General Motions, User Manual, Technical Report, Analytical Methods,
Inc., Redmond, WA, 1990.

[36] High order panel emthod paper

[37] Fornasier, L., HISS-A high order subsonic/supersonic singularity method for calcu-
lating linearized potential flow, 17th AIAA Fluid dynamics, plasma dynamics and
lasers conference, 1984.

[38] Drela, M., XFOIL: An Analysis and Design System for Low Reynolds Number Air-
foils, Conference on Low Reynolds Number Aerodynamics, University of Notre
Dame, June 1989.

[39] Liebeck, R.H., Subsonic Airfoil Design, In, Applied Computational aerodynamics,
Vol 125 Progress in Astronautics and Aeronautics, AIAA, pp, 133-165, 1990.

[40] Mughal, B., Drela, M., A calculation method for the three-dimensional boundary
layer equations in integralform, AIAA Paper 93-0786, Reno, NV, 1993.

[41] Nishida, B., and Drela, M., Fully Simultaneous Couplingfor three dimensional vis-
cous/inviscidflows, AIAA-95-1806-CP, AIAA 13th Applied Aerodynamics Confer-
ence, June, 1995.

[42] Bulirsch, R., Stoer, J., The congugate-gradient method of Hestenes and Sti-
iefel,Section 8,7 in Introduction to Numerical Analysis., New York, Springer Verlag,
1991.

128
[43] Freund, R., Nachtigal, N. QMR: A quasi-minimal residual methodfor non-hermitian
linear systems, Numer. Math., 60:315-339, 1991.

[44] Nabors, K, White, J.K., FastCap: A Multipole Accelerated 3-D Extraction Program,
IEEE Trans. On Comp. Aided Design, Vol. 10, No. 11, 1991.

[45] Aluru, N.R., Nadkami, V., White, J.K., A parallel precorrected-FFT based capaci-
tance extraction program for signal integrity analysis, Proceedings of the 33rd De-
sign Automation Conference, Las Vegas, Jun, 1996.

[46] Board, J.A., Hakura, Z.S., Elliot, W.S., Gray, D.C., Blanke, W.J., Leathrum,
J.F.,Jr.Scalableimplementationsof multipole-acceleratedalgorithmsfor molecular
dynamics, Technical Report 94-002, Duke University, 1994.

[47] Willis, D.J., Lee, J., Coehlo, C.P.Bardhan, J, Hu, X., Ried, H., White, J.K., pFFT2+:
A Genericprecorrected-FFTalgorithmin C++: UsersManual,Work in progress.

[48] Frigo, M., Johnson, S.,FFTW Users Manual, Version 3.0.1, Massachusetts Institute
Of Technology, 2003.

[49] Vassberg, J, C. A fast surface-panel method capable of solving million-element prob-

lems, 35th AIAA Aerospace Sciences Meeting and Exhibit,AIAA-97-0168, 1997.

[50] Cooley, J.W., Tukey, O.W. An algorithm for the machine calculation of complex
fourier series., Math. Comput., 19:297-301, 1965.

[51] Jameson, A. Aerodynamics design via control theory, J. Sci. Comp. 3:233-60, 1988.

[52] Jameson, A. Solution of the Euler Equations by a Multigrid Method, Appplied math.
and Computation, Vol. 13, pp.327-356, 1983.

[53] L. Rosenhead, Theformation of vorticies from a surface of discontinuity, Proc. Roy.

Soc. London Ser. A 134, 170-192, 1931.

[54] Cottet, G-H, Koumoutsakos, P.D.,Vortex Methods:Theory and Applications, Cam-

bridge University Press, London, 2000.

129
[55] A. J. Chorin, Vortexmodels and boundary layer instability, SIAM J. Sci. Stat. Com-
put. 1:1-21, 1980.

[56] Cottet, G-H, Koumoutsakos, P., Ould-Salihi, M.L.,Vortex methods with spatially
varying cores, J. Comput. Phys. 162 164-185, 2000.

[57] J. T. Beale, A. Majda, Vortex methods I: convergence in three dimensions, Math.

Comput. 29 (159) 1-27, 1982.

[58] A. Leonard, Vortex methodsforflow simulation, J. Comput. Phys. 37,289-335, 1980.

[59] Leonard, A. Computing three-dimensional incompressible flows with vortex ele-

ments. Ann. Rev. Fluid Mech. 17:523-559, 1985.

[60] 0. M. Knio, A. F. Ghoniem, The three-dimensional structure of periodic vorticity

layers under non-symmetric conditions, J. Fluid Mech. 243, 353-392, 1992.

[61] Spalart, P.R., Two recent extensions of thevortex particle method, presented at the
22nd Aerospace Sciences Meeting and Exhibit, AIAA Paper 1984-343, Reno, Jan-
uary 1984.

[62] Voutsinas S.G., Belessis M.A., and Rados K.G., Investigation of the yawed opera-
tion of wind turbines by means of a vortex particle method. AGARD-CP-552 FDP
Symposium on Aerodynamics and Aeroacoustics of Rotorcraft, Berlin, Germany,
Paper 11, 1995.

[63] Rehbach, C., Calcul numerique d'ecoulement tridimensionels instationaires avec

nappes tourillonnaires, La Recherche Aerospatiale, pp.289-298, 1977.

[64] Huberson, S. Calcule d'ecoulements tridimensionnels instationnaires incomprress-

ibles par une methode particulaire. Journal de Mecaniques Theorique et Appliquee,
3(1): 805-819, 1984.

[65] Eldredge, J. D., Efficient tools for the simulation of flapping wing flows, presented
at 43rd AIAA Aerospace Sciences Meeting, AIAA Paper 2005-0085, Reno, January
2005.

130
[66] Y.M. Marzouk, A.F. Ghoniem, Mechanism of streamwise vorticity formation in a
transverse jet, 40th Aerospace Sciences Meeting and Exhibit, AIAA-2002-1063,
January 2002.

[67] A. Gharakhani, A. F. Ghoniem, Three-dimensional vortex simulation of time depen-

dent incompressible internal viscousflows, J. Comput. Phys. 134:75-95, 1997.

[68] Anderson, C., A method of local corrections for computing the velocity field due to
a distribution of vortex blobs. J. Comput. Phys., vol 62, 111-123, 1986.

[69] Ploumhans, P., Winckelmans, G.S., Vortex methods for high-resolution simulations
of viscous flow past bluff bodies of general geometry. J Comput Phys, 165:354-406,
2000.

[70] Appel, A.A., An efficient program for many body simulations', SIAM Journal of
scientific and statistical computing, vol. 16, n. 1, pp85-103, 1985.

[71] Maskew,B., Predictionof subsonic aerodynamiccharacteristics:a case for low-

order panel methods., J. Aircr., 19:157-63, 1982.

[72] Kerwin, J.E., Kinnas, S.A., Lee, J-T, Shih, W-ZA Surface Panel MEthod for the
Hydrodynamic Analysis of Ducted Propellers, SNAME Transactions, Vol. 95, 93-
122, 1987.

[73] Kinnas, S.A., Hsin, C.Y.,Boundary Element Methodfor theAnalysis of the Unsteady
Flow Around Extreme Propeller Geometries, AIAA Journal, Vol. 30, No. 3,688-696,
1992.

[74] Liu, P, Bose, N., Colbourne, B., A Broyden Numerical Kutta Condition for an Un-
steady Panel Method,

[75] Anderson, D.A., Tennehill, J.C., Pletcher, R.H.,Computational Fluid Mechanics and
Heat Transfer, Taylor and Francis, 1997.

[76] D.J.Willis, J.Peraire, and J.K.White To be presented, 44th AIAA Aerospace Sciences
Meeting and Exhibit, AIAA 2006-1253, 2006.

131
[77] Arfken, G. "Helmholtz's Theorem." Secl.15 in Mathematical Methods for Physi-
cists, 3rd ed. Orlando, FL: Academic Press, pp. 78-84, 1985.

[78] A. M. O. Smith, The Panel Method: Its Original Development, In, Applied Com-
putational aerodynamics, Vol 125 Progress in Astronautics and Aeronautics, AIAA,
pp, 3-20, 1990.

[79] D. L. Ashby, M. R. Dudley, and S. K. Iguchi, Development and Validation of an

Advanced Low Order Panel Method, NASA TM 101024, Oct 1988.

[80] B. Maskew,PROGRAM VSAERO: A Computer Program for Calculating the Non-

Linear Aerodynamic Characteristicsof Arbitrary Configurations", NASA CR-
166476, Nov 1982.

[81] B. Rosen, J. P. Laiosa, SPLASH Nonlinear and Unsteady Free-Surface Analysis

Code for Grand Prix Yacht Design, 13th Chesapeake Sailing Yacht Symposium,
Annapolis, Jan 1997.

[82] P.Koumoutsakos, Multiscale Flow Simulations Using Particles, Annu. Rev. Fluid
Mech., 37:457-487, 2005.

[83] WAMIT 6.2. WAMIT User Manual, Version 6.2, 6.2PC and 6. IS, 6.1S-PC, Dept. Of
Ocean Engineering, M.I.T.

[84] J. R. Philips and J. K. White, A Precorrected-FFT Method for Electrostatic Analysis

of Complicated 3-D Structures, IEEE Transactions On Computer-Aided Design of
Integrated Circuits and Systems, IEEE, Vol. 16, 1997.

[85] L. Greengard and V. Rohklin, A Fast Algorithm for Particle Simulations, J. Comp.
Phys., 73:pp325-384, 1987.

[86] J. Barnes and P. Hut, A Hierarchical O(N log N) Force Caluclation Algorithm, Na-
ture 324, pp 446-449, 1986.

[87] L. Morino,Steady, Oscillatory and Unsteady Subsonic and Supersonic Aerodynam-

ics - Production Version (SOUSSA) NASA CR-157130, 1980.

132
 [88] G. S. Winckelmansand A. Leonard, Contributions to VortexParticle Methods for the
Computation of Three-Dimensional Incompressible Unsteady Flows, J. Comp. Phys,
109, pp 247-273, 1993.

[89] Brebbia, C., Boundary Element Methods in Engineering: Proceedings, Springer,

New York, 1982.

[90] Chen, G., Zhou, J. Boundary Element Methods, 1st Ed., Academic Press, 1992.

[91] A. M. O. Smith and J. .L Hess, Calculation of potential fow about arbitrary bodies.,
Progress in Aeronautic Sciences, 8, 1960.

[92] R.Kress, Linear Integral Equations. Springer-Verlag, New York, 1989.

[93] K.E.Atkinson, and W.Han Theoretical Numerical Analysis: A Functional Analysis

Framework, Springer-Verlag, 2001.

[94] Willis, D.J. A pFFT accelerated, High Order Panel Method, Masters Thesis, MIT,
2003.

[95] Jackson, J.D., Classical Electrodynamics, Third Edition, John Wiley and Sons Inc.,
New York, 1998.

[96] Kythe, P.K. An introduction to boundary element methods., CRC Press, Boca Raton,
1995.

[97] E.Hairer, C.Lubich, G.Wanner, Geometric Numerical Integration, Springer Series

in Computational Mathematics, Springer, 2004.

[98] D.J.Willis, J. Peraire, J.K.White, Body Piercing Wakes For Automatic Wake Gener-
ation in the Potential Panel Method, Work in Progress.

[99] J. N. Newman,Distributionof Sourcesand Normal Dipoles Over a Quadrilateral

Panel, J. Eng. Math., 20, 1985.

[100] Y. Saad and M. Schultz, GMRES: A generalized Minimal Residual Algorithm for
Solving Non-Symmetrixc Linear Systems, SIAM, J. Sci. Stat. Comp., Vol 7, 1986.

133
[101] Z. Zhu, B. Song and J. K. White, Algorithms in FastIMP: A Fast Wideband
Impedance ExtractionProgramfor Complicated3D Geometries,Proceedings of
IEEE/ACM Design Automation Conference, Anaheim, California, June 2-6, 2003.

[102] R.W.Hockney, J.W.Eastwood, Computer Simulation Using Particles, Taylor and

Francis Inc., Bristol, 1988.

[103] J. Katz, A. Plotkin, Low-Speed Aerodynamics,Second Edition, Cambridge

Aerospace Series, Cambridge, 2001.

[104] T. Theodorsen, General Theory of Aerodynamic Instability and the Mechanism of

Flutter,NACA Report 496, 1935.

[105] J.Peraire, J.Piero, L.Formaggia, K.Morgan, and O.C.Zienkiewicz, Finite Element

Euler Computations in Three Dimensions, Int. J. Num. Meth. in Engn. 26, 2135-
2159, 1988.

[106] D. L. Ashby, M. R. Dudley, and S. K. Iguchi, Development and Validation of an

Advanced Low Order Panel Method, NASA TM 101024, Oct 1988.

[107] B. Maskew,PROGRAM VSAERO: A Computer Program for Calculating the Non-

Linear Aerodynamic Characteristicsof Arbitrary Configurations", NASA CR-
166476, Nov 1982.

[108] D.J. Willis, J. Peraire, J.K. White, A CombinedpFFT-Multipole Tree Code, Unsteady
Panel Method with Vortex Particle Wakes, 43rd AIAA Aerospace sciences Meeting
and Exhibit, AIAA-2005-0854, 10-13 January 2005.

[109] J. R. Philips and J. K. White, A Precorrected-FFT Methodfor Electrostatic Analysis

of Complicated 3-D Structures, IEEE Transactions On Computer-Aided Design of
Integrated Circuits and Systems, IEEE, Vol. 16, 1997.

[110] L. Greengard and V. Rohklin, A Fast Algorithm for Particle Simulations, J. Comp.
Phys., 73:pp325-384, 1987.

134
[111] X. Wang, J.N.Newman, and J.K.White, Robust Algorithms for Boundary Element
Integrals on Curved Surfaces , International Conference on Modeling and Simula-
tion of Micro systems, 2000.

[112] J. N. Newman,Distributionof Sources and Normal Dipoles Over a Quadrilateral

Panel, J. Eng. Math., 20, 1985.

[113] J. Katz, A. Plotkin, Low-Speed Aerodynamics, 2nd Edition, Cambridge university

Press, New York, 2001.

[114] O.D.Kellogg Foundations of Potential Theory,J. Springer, New York, 1929.

[115] C.Y. Hsin, J.E.Kerwin, and J.N.Newman, A high order panel method based on B-
Splines. Proceedings of Sixth International Conference on Numerical Ship Hydrody-
namics, eds. V.C. Patel and F. Stem, National Academy Press: Washington, pp133-
151, 1993.

[116] A.E. Magnus, and M.A.Epton, PAN AIR- A Computer Program for predicting Sub-
sonic or SupersonicLinear PotentialFlows about Arbitrary ConfigurationsUsing
a High Order Panel Method, Volume 1, Theory Document (Version 1.0), NASA
CR-3251, 1980.

[117] C.S. Lee and J.E. Kerwin, A B-Spline Higher-Order Panel Method Applied to Two
Dimensional Lifting Problems, Journal Of Ship Research, Volume 47, Number 4, pp
290-298, 2003.

[118] P. Ramachandran, S.C. Rajan, M. Ramakrishna, A Fast Multipole Methodfor Higher

Order VortexPanels in TwoDimensions, SIAM Journal on Scientific Computing, Vol
26, Issue 5, pp. 1620-1642, 2005.

[119] H. Maniar, A B-Spline Based Higher Order Method in 3D, Presented at the 10th
Workshop on Water Waves and Floating Bodies, Oxford, England, 1995.

[120] H. Maniar, A Three Dimensional High Order Panel Method Based on B-splines,
PhD. Thesis, Massachusetts Institute Of Technology, 1995.

135
[121] F.E. Ehlers, F.T.Johnson, and P.E. Rubbert, A Higher Order Panel Methodfor Lin-
earized Supersonic Flow. AIAA 76-381, July 1976.

[122] J. Hess, and A.M.O Smith, Calculation of Non-lifting potential flow about arbitrary
three dimensional bodies, J. Ship Research, 8, pp 22-44, 1994.

[123] U. Lemma, V.Marchese,and L. Morino, High-order BEM for potential transonic

flows, Computational Mechanics, Springer Verlag, Vol. 21, No. 3, pp. 243-252, 1998.

[124] J.H.M. Frijns, Improving the Accuracy of the Boundary Element Method by the use
of Second-Order Interpolation Functions., IEEE Transactions on Biomedical Engi-
neering, Vol. 47, No. 10, October 2000.

[125] A. Buchau, W. Rieger, and W.M.Rucker BEM Computations Using the Fast Multi-
pole Methodin Combinationwith HigherOrderElementsand the GalerkinMethod.,
IEEE Transactions on Magnetics, Vol. 37 No. 5, September 2001.

[126] K.J.Bathe Finite Element Procedures, Prentice Hall, 1995.

[127] K.E.Atkinson, and W.Han Theoretical Numerical Analysis: A Functional Analysis

Framework, Springer-Verlag, 2001.

[128] R. Klees, Numerical calculation of weakly singular surface integrals, Journal of

Geodesy, Vol 70, pp 781-797, 1996.

[129] J.-L. Guermond,NumericalQuadraturesfor layerpotentials over curved domains

in R3 , SIAM Jounral of Numercial Analysis, Vol. 29, No. 5, pp. 1347-1369, 1992.

[130] K. Hayami,C.A.BrebbiaQuadratureMethodfor Singularand Nearly SingularInte-

grals in 3-D Boundary Element Method, Boundary Element X, vol. 1, Mathematical
and Computational Aspects. 1987.

[131] Graglia, R.D. On the Numerical Integration of the Linear Shape Functions imes
the 3-D Green's Functionor its Gradienton a Plane Triangle,IEEE Transactions
on Antennas and Propagation, Vol. 41, No. 10, pp 1448-1456, 1993.

136
[132] L. Rossi, P.J.Cullen. On the fully numerical evaluation of the linear-shape func-
tion times the 3D Green's Function on a Plane Tirangle, IEEE Transactions on mi-
crowave theory and techniques, vol. 47, No. 4, pp. 398 - 403, 1999.

[133] P. Kim, and U.J.Choi, Two Trigonometric Quadrature Formulae for evaluating hy-
persingular integrals, International Journal For Numerical Methods in Engineering,
Vol. 56, pp 469-486, 2003.

[134] Cools, R., An encyclopedia of cubatureformulas, Journal of Complexity, Vol 19(3),

pages 445-453, 2003.

[135] Cools,R. Monomialcubaturerulse since "Stroud": a compilation-part 2. J. Com-

put. Appl. Math, 112(1-2): 21-27, 1999.

[136] Cools, R. and Rabinowitz Monomial cubature rules since "Stroud": a compilation,
J. Comput. Appl. Math, 48: 309-326, 1993.

[137] Kinnas, S.A., and Hsin, C.Y., Boundary element method for the analysis of the un-
steady flow around extreme propeller geometries, AIAA Journal, Vol. 30, No. 3,
pp.688-696, 1992.

[138] Kerwin, J.E., Kinnas, S.A., Lee, J-T and Shi, W-Z., A surface panel method for
the hydrodynamic analysis of ducted propellers, SNAME Transactions, Vol. 95, pp.
93-122, 1987.

[139] Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., Numerical Recipes in
C: The Art of Scientific Computing, 2nd Edition, Cambridge University Press, New
York, 1992.

[140] Broyden, C.G., Quasi-Newton Methods and Their Applications to Function Mini-
mization, Mathematics of Computation, Vol. 21, No. 99, pp. 368-381, 1967.

[141] Fletcher, R., Practical Methods of Optimization (2nd Ed.), John wiley, Chichester,
1987.

137
[142] Liu, P., Bose, N., and Colbourne, B.,A Broyden numerical Kutta condition for an
unsteady panel method Internaional Shipbuidling Progress, vol. 49, no. December,
pp. 263-274, 2002.

[143] Hall, K. C., Pigott, S.A., and Hall, S.R., Power Requirements for Large-Amplitude
Flapping Flight, Journal of Aircraft, Vol. 35, No. 3, pp.352-361, 1998.

138