PROPELLER PERFORMANCE ANALYSIS AND MULTIDISCIPLINARY
OPTIMIZATION USING A GENETIC ALGORITHM
Except where reference is made to the work of others, the work described in this
dissertation is my own or was done in collaboration with my advisory committee.
This dissertation does not include proprietary or classified information.
____________________________________
Christoph Burger
Certificate of Approval:
_______________________ _______________________
John E. Burkhalter Roy J. Hartfield, Jr., Chair
Professor Emeritus Professor
Aerospace Engineering Aerospace Engineering
_______________________ _______________________
Robert S. Gross Ronald M. Barrett
Associate Professor Associate Professor
Aerospace Engineering Aerospace Engineering
University of Kansa
_______________________
George T. Flowers
Interim Dean
Graduate School
ii
PROPELLER PERFORMANCE ANALYSIS AND MULTIDISCIPLINARY
OPTIMIZATION USING A GENETIC ALGORITHM
Christoph Burger
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of the
Requirements for the
Degree of
Doctor of Philosophy
Auburn Alabama
December 17, 2007
iii
PROPELLER PERFORMANCE ANALYSIS AND MULTIDISCIPLINARY
OPTIMIZATION USING A GENETIC ALGORITHM
Christoph Burger
Permission is granted to Auburn University to make copies of this dissertation at its
discretion, upon request of individuals or institutions at their expense. The author
reserves all publication rights.
_______________________
Signature of Author
_______________________
Date of Graduation
iv
VITA
Christoph Burger, son of Richard and Baerbel (Kundler) Burger, was born on
September 18, 1972 in Munich, Bavaria (Germany). Christoph graduated from High
School in St.Georgen, Freiburg (Germany) in June of 1993. From October of 1993 until
April of 1995 he served as an EMT in fulfillment of the mandatory military service. He
attended the University of Stuttgart in October of 1995 and graduated with a Diploma
Engineer (Master of Science) Degree in Aerospace Engineering, in 2001. In the fall of
2001, Christoph started work as a Research Scientist at the Aerospace Engineering
Department at Auburn University, Alabama. After about two years of research and
development work on unmanned aerial systems, he enrolled as a student at Auburn
University to complete his Doctor of Philosophy Degree in Aerospace Engineering.
v
DISSERTATION ABSTRACT
PROPELLER PERFORMANCE ANALYSIS AND MULTIDISCIPLINARY
OPTIMIZATION USING A GENETIC ALGORITHM
Christoph Burger
Doctor of Philosophy December 17, 2007
(Diploma Engineer (M.S.), University of Stuttgart 2001)
181 Typed Pages
Directed by Roy J. Hartfield
A propeller performance analysis program has been developed and integrated into
a Genetic Algorithm for design optimization. The design tool will produce optimal
propeller geometries for a given goal, which includes performance and/or acoustic
signature. A vortex lattice model is used for the propeller performance analysis and a
subsonic compact source model is used for the acoustic signature determination.
Compressibility effects are taken into account with the implementation of PrandtlGlauert
domain stretching. Viscous effects are considered with a simple Reynolds number based
model to account for the effects of viscosity in the spanwise direction. An empirical flow
vi
separation model developed from experimental lift and drag coefficient data of a NACA
0012 airfoil is included. The propeller geometry is generated using a recently introduced
Class/Shape function methodology to allow for efficient use of a wide design space.
Optimizing the angle of attack, the chord, the sweep and the local airfoil sections,
produced blades with favorable tradeoffs between single and multiple point optimizations
of propeller performance and acoustic noise signatures. Optimizations using a binary
encoded IMPROVE
?
Genetic Algorithm (GA) and a real encoded GA were obtained
after optimization runs with some premature convergence. The newly developed real
encoded GA was used to obtain the majority of the results which produced generally
better convergence characteristics when compared to the binary encoded GA. The
optimization tradeoffs show that single point optimized propellers have favorable
performance, but circulation distributions were less smooth when compared to dual point
or multiobjective optimizations. Some of the single point optimizations generated
propellers with proplets which show a loading shift to the blade tip region. When noise is
included into the objective functions some propellers indicate a circulation shift to the
inboard sections of the propeller as well as a reduction in propeller diameter. In addition
the propeller number was increased in some optimizations to reduce the acoustic blade
signature.
vii
ACKNOWLEDGMENTS
I would like to thank Dr. Roy J. Hartfield for his exceptional support and
guidance during the exertive time. Without his help and especially his patience it would
have been impossible to complete the Degree. I also would like to thank Dr. John E.
Burkhalter for his advice and help with the programming work.
viii
Style manual or journal used:
The American Institute of Aeronautics and Astronautics Journal
Computer software used:
Microsoft Word 2003, IMPROVE
?
3.1 Genetic Algorithm, Compaq Visual FORTRAN,
Tecplot 10, Microsoft Excel 2003
ix
TABLE OF CONTENTS
LIST OF FIGURES ........................................................................................................... xi
LIST OF TABLES............................................................................................................ xv
NOMENCLATURE ........................................................................................................ xvi
1 INTRODUCTION ...................................................................................................... 1
2 BASIC THEORY........................................................................................................ 5
2.1.1 LiftingLine Theory .................................................................................... 6
2.1.2 Lifting Surface Theory................................................................................ 9
2.2 Propeller Wake.................................................................................................. 13
2.3 Propeller Performance Prediction by VLM ...................................................... 16
2.4 Propeller Geometry........................................................................................... 30
2.4.1 2D Airfoil Geometry ............................................................................... 30
2.4.2 3D Propeller Geometry............................................................................ 35
2.5 Compressibility................................................................................................. 39
2.6 Viscous Model .................................................................................................. 46
2.7 Flow Separation Model..................................................................................... 48
2.8 Acoustic Signature Model................................................................................. 53
3 PROPELLER AERODYNAMIC PERFORMANCE PREDICTION PROGRAM. 63
4 CODE VERIFICATION........................................................................................... 66
4.1 Gauss Seidel Solver convergence ..................................................................... 66
4.2 Wake investigation............................................................................................ 70
4.3 Number of Chordwise and Spanwise Panels .................................................... 73
4.4 Validation of the numerical method ................................................................. 75
4.4.1 Comparison to lifting line and lifting surface method.............................. 75
4.4.2 NACA 109622 Propeller........................................................................... 78
5 GENETIC ALGORITHM ........................................................................................ 82
5.1 Optimization methods....................................................................................... 82
5.2 Genetic Algorithm ............................................................................................ 85
x
6 IMPLEMENTATION OF THE PROPELLER PERFORMANCE PROGRAM
INTO THE GENETIC ALGORITHM ............................................................................. 88
7 OPTIMIZATION RESULTS.................................................................................... 91
7.1 Comparison of propeller optimizations with and without viscous and
compressibility effects taken into account.................................................................... 92
7.1.1 Single point Optimization for Cruise Condition: Binary encoded GA..... 92
7.1.2 Single point Optimization for Cruise Condition: Real encoded GA ...... 100
7.2 Single point Optimization Launch Condition................................................. 107
7.3 Dual point Optimization Cruise and Launch Condition ................................. 111
7.4 Single Point Optimization in Cruise considering Propeller Noise.................. 116
7.5 Single/Dual Point Optimization in Cruise/Launch Condition considering
Propeller Noise and Performance ............................................................................... 119
8 CONCLUSION & RECOMMENDATIONS......................................................... 129
REFERENCES ............................................................................................................... 132
APPENDIX A: BIOTSAVART LAW.......................................................................... 142
APPENDIX B: PROPELLER GA VARIABLES .......................................................... 147
APPENDIX C: PROPELLER WAKE INVESTIGATION............................................ 149
APPENDIX D: PROPELLER PANEL DISCRETIZATION......................................... 150
APPENDIX E: NACA 109622 PROPELLER GEOMETRY ........................................ 151
APPENDIX F: NACA 0009 AIRFOIL DATA .............................................................. 154
APPENDIX G: PROPELLER GEOMETRY OF SKYHAWK 172............................... 156
APPENDIX H: WIND TUNNEL DATA CAM 13x7 PROPELLER ............................ 160
xi
LIST OF FIGURES
Figure 2.1 Liftingline model with single horseshoe vortex on finite wing ....................... 6
Figure 2.2 Distribution of horseshoe vortices, liftingline model....................................... 7
Figure 2.3 Induced drag due to induced downwash of the trailing vortices....................... 8
Figure 2.4 Lifting surface method by horseshoe elements ............................................... 10
Figure 2.5 Vortex ring element arrangement.................................................................... 11
Figure 2.6 Propeller blade grid discretization................................................................... 13
Figure 2.7 Propeller wake geometry................................................................................. 14
Figure 2.8 Propeller wake horseshoe trailer discretization............................................... 16
Figure 2.9 Induced velocity at point P(x,y,z) by a single vortex element ........................ 18
Figure 2.10 Scanning sequence on propeller blade panels ............................................... 21
Figure 2.11 Normal vector on propeller blade panel ........................................................ 22
Figure 2.12 Propeller boundary condition and orientation angles.................................... 23
Figure 2.13 Vortex circulation sequence on a two bladed propeller ................................ 24
Figure 2.14 Trailing vortex segments responsible for induced downwash ...................... 26
Figure 2.15 Total thrust and drag of the propeller............................................................ 28
Figure 2.16 Class function 2D design space..................................................................... 32
Figure 2.17 Bernstein polynomial representation of the unit shape function................... 33
Figure 2.18 Class/Shape function single and dual surface geometries............................. 35
Figure 2.19 Propeller sweep functions based on CST method ......................................... 38
Figure 2.20 Thin propeller blade geometry by Class/Shape function scheme ................. 39
Figure 2.21 Comparison of lift coefficient for incompressible and compressible................
flow with direct application of the PrandtlGlauert rule (Source .........................................
from Ref.[45] ).................................................................................................................. 40
Figure 2.22 PrandtlGlauert domain stretching in flow direction..................................... 41
Figure 2.23 PrandtlGlauert domain stretching of vortex ring elements .......................... 42
Figure 2.24 Spanwise lift coefficient normalized by blade radius.................................... 43
Figure 2.25 Generic two bladed propeller with wake....................................................... 44
Figure 2.26 Compressibility effects vs propeller tip Mach number ................................. 45
Figure 2.27 Spanwise lift coefficient normalized by blade radius
[23]
............................... 48
xii
Figure 2.28 Lift coefficient of NACA 0009 in post stall.................................................. 50
Figure 2.29 Drag coefficient of NACA 0009 in post stall................................................ 51
Figure 2.30 Post stall thrust behavior of a NACA109622 propeller with empirical ............
stall model......................................................................................................................... 52
Figure 2.31 Post stall power behavior of a NACA109622 propeller with empirical ...........
stall model......................................................................................................................... 52
Figure 2.32 Locations of noise sources............................................................................. 54
Figure 2.33 Acoustic signature analysis geometry ........................................................... 57
Figure 2.34 Linearized airfoil thickness distribution of a APC 14x7 propeller................ 59
Figure 2.35 1C160 quarter scale Cessna 172 propeller .................................................... 60
Figure 2.36 Pressure perturbation comparison experiment and theory from........................
Miller
[5]
............................................................................................................................. 61
Figure 2.37 Pressure perturbation of 1C160 quarter scale propeller ................................ 62
Figure 3.1 Flow chart propeller aerodynamic performance program............................... 65
Figure 4.1 Summation of LHS and RHS difference after matrix sweeps......................... 68
Figure 4.2 Accuracy of the solution after matrix sweeps by Gauss Seidel solver............ 69
Figure 4.24 Wake vortex sheet extensions for ?, 2 and 5 rotations downstream ........... 71
Figure 4.5 Propeller thrust variation due to wake vortex sheet extension........................ 72
Figure 4.6 Propeller Discretization schemes .................................................................... 73
Figure 4.7 Lift force convergence by discretization scheme and panel number .............. 74
Figure 4.8 Wing and wake discretization ......................................................................... 76
Figure 4.9 Lifting line and lifting surface method cl and cd of flat plate wing AR=6 ..... 77
Figure 4.10 NACA109622 propeller geometry with wake............................................... 79
Figure 4.11 Propeller efficiency comparison at different advance ratios ......................... 80
Figure 4.12 Propeller thrust coefficient at various advance ratios.................................... 81
Figure 4.13 Propeller power coefficients at various advance ratios ................................. 81
Figure 6.1.Flow chart propeller performance program implementation into GA............. 89
Figure 7.1 Optimized cruise propeller geometry from binary encoded GA; inviscid, .........
incompressible .................................................................................................................. 95
Figure 7.2 Optimized cruise propeller geometry from binary encoded GA; inviscid, .........
compressible ..................................................................................................................... 96
Figure 7.3 Optimized cruise propeller geometry from binary encoded GA; viscous,..........
incompressible .................................................................................................................. 97
Figure 7.4 Optimized propeller spanwise circulation distributions for cruise.................. 98
Figure 7.5 Objective function convergences binary encoded propeller ........................... 99
xiii
Figure 7.6 Objective function convergence inviscid, incompressible, for binary ................
encoded GA .................................................................................................................... 100
Figure 7.7 Optimized cruise propeller geometry from real encoded GA; inviscid, .............
incompressible ................................................................................................................ 102
Figure 7.8 Optimized cruise propeller geometry from real encoded GA; inviscid, .............
compressible ................................................................................................................... 103
Figure 7.9 Optimized cruise propeller geometry from real encoded GA; viscous, ..............
incompressible ................................................................................................................ 104
Figure 7.10 Optimized circulation distributions for cruise propeller, real encoded .............
GA ................................................................................................................................. 105
Figure 7.11 Objective function convergences real encoded GA, cruise condition......... 106
Figure 7.12 CAM 13x7 propeller efficiency at 19 [m/s] free stream velocity ............... 107
Figure 7.13 Optimized launch propeller geometry: viscous, incompressible case......... 109
Figure 7.14 Optimized propeller spanwise circulation distribution viscous,........................
incompressible case ........................................................................................................ 110
Figure 7.15 Objective function convergence from real encoded GA: viscous,....................
incompressible case ........................................................................................................ 111
Figure 7.16 Optimized cruise and launch propeller; viscous, incompressible case........ 114
Figure 7.17 Objective function convergences real encoded GA, cruise condition......... 115
Figure 7.18 Objective function convergences real encoded GA, cruise condition......... 116
Figure 7.19 Spanwise circulation distribution of the optimized propeller in cruise.............
condition ......................................................................................................................... 118
Figure 7.20 Optimized cruise propeller from real coded GA; viscous, incompressible.......
with noise considered...................................................................................................... 119
Figure 7.21 Optimized cruise propeller; viscous, incompressible case with noise and........
performance considered.................................................................................................. 122
Figure 7.22 Optimized launch propeller; viscous, incompressible case with noise and.......
performance considered.................................................................................................. 123
Figure 7.23 Optimized launch / cruise propeller; viscous, incompressible case with ..........
noise and performance considered.................................................................................. 124
Figure 7.24 Spanwise circulation distribution of the optimized propeller in cruise....... 125
Figure 7.25 Spanwise circulation distribution of the optimized propeller at launch...... 126
Figure 7.26 Spanwise circulation distribution of the optimized propeller in cruise.............
and launch condition....................................................................................................... 127
Figure 7.27 Objective function convergences, cruise, launch and combined case......... 128
xiv
Figure A.1 Velocity at point P due to a vortex distribution............................................ 144
Figure A.2 Induced velocity at point P by a vortex segment.......................................... 145
xv
LIST OF TABLES
Table 7.1 Fixed input parameters cruise condition........................................................... 92
Table 7.2 Optimized propellers performance parameter results cruise condition ............ 94
Table 7.3 Optimized propellers performance parameter results real encoded GA......... 101
Table 7.3 Fixed input parameters for launch condition .................................................. 108
Table 7.4 Optimized propellers performance parameter results for launch condition ... 108
Table 7.5 Fixed input parameters for dual point optimization........................................ 112
Table 7.6 Optimized propeller performance parameters: cruise and launch .................. 113
Table 7.7 Optimized propellers performance parameters cruise condition with noise.........
considered ....................................................................................................................... 117
Table 7.8 Propellers performance parameter results with noise optimization................ 121
xvi
NOMENCLATURE
Roman Letters
A
110
Bernstein polynomial coefficients, 2D airfoil section []
A
n
Shape function coefficients []
zyx AAA ,, Influence coefficient matrices [1/m]
AOA Angle of attack function [?]
aoa
n
AOA function coefficients []
Au, Al Coefficient functions for shape function []
au
n
, al
n
Shape function coefficients []
b Wing span [m]
b Unit propeller radius []
B Vector with right hand side of boundary condition []
zyx BBB ,, Influence coefficients horseshoe trailers [1/m]
c Section propeller chord [m]
c
d
Drag coefficient []
c
d_max
maximum drag coefficient []
xvii
c
f
Skin friction coefficient []
c Maximum chord [m]
1
2
N
N
C Class function with coefficients []
Chord Propeller chord function [m]
chord
n
Propeller chord function coefficients [m]
c
l
Section lift coefficient []
c
p
Propeller power coefficient []
c
t
Propeller thrust coefficient []
d Propeller diameter [m]
D Total induce drag [N]
D
i
Induced rag [N]
dl length of the vortex element [m]
D
visc
Viscous drag [N]
?D Panel drag force due to induced downwash [N]
?F Panel force perpendicular to the onset velocity [N]
?Q Panel torque [Nm]
?S Panel width, spanwise [m]
?T Panel thrust [N]
?V Velocity induced by single vortex element [m/s]
?y Panel width in the spanwise direction [m]
F
?
Propeller section drag force [N]
F
R
Propeller section radial force [N]
xviii
F
T
Propeller section thrust force [N]
h Rankine body diameter [m]
horeshoe_x,y,z x, y,z coordinates of the horseshoe elements [m]
i, j, k Program counters []
inf_coeff
x,y,z
Influence coefficients in x, y and z direction [1/m]
L Lift [N]
M Mach number []
m Wake gradient []
M
r
Mach number with respect to observer []
M
s
Mach number relative to propeller []
n Propeller revolutions [rad/sec]
n? Surface normal []
N
1
, N
2
Class function exponents []
nodal_x,y,z x, y,z coordinates of the propeller nodal elements []
nu
n
, nl
n
Class function coefficients []
nxbl Number of propeller blades []
P Propelr power [W]
pp
0
Presure perturbation [Pa]
P
visc
Power due to viscous torque [W]
Q Propeller torque [Nm]
qm Mas flux [kg/s]
Q
visc
Torque due to viscous drag [Nm]
r Distance collocation point to vortex element; observer [m]
xix
Re
c
Section chord Reynolds number []
r
i,j
Distance of center of vortex ring panel to origin [m]
S Surface area [m
2
]
S
r,n
Bernstein polynomial term []
Sweep Propeller sweep function [m]
sweep
n
Propeller sweep function coefficients []
T Propelr thrust [N]
t Airfoil thickness [m]
V
?
Free Stream velocity [m/s]
V
local
Local velocity [m/s]
V
onset
Product of free stream velocity and propeller RPM [m/s]
w
ind
Induced downwash velocity [m/s]
x, y, z Coordinate system x, y and z axis []
x Vector with circulations ? [m
2
/s]
Greek Letters
? Flow angle / Angle of source with y axis [?]
?
onset
Angle of V
onset
with xaxis [?]
?
i
Induced downwash angle [?]
? Wake angle [?]
? Circulation [m
2
/s]
xx
? Onset flow angle with x axis [?]
? Angle of panel normal with zaxis [?]
? Source/sink distance [m]
? Propeller efficiency []
? Angles in the x,yplane [?]
? Propeller panel dihedral angle [?]
? Viscosity coefficient []
? Density of air [kg/m
3
]
? Retarde time [s]
? Angle of observer with x axis [?]
? Fraction of local chord x/c []
? Angular velocity [1/sec]
Mathematical Expressions
? Increment []
t?
?
First order partial derivative in time []
Abbreviations
AR Aspect ratio
BP Bernstein polynomials
xxi
CST Class/Shape function transformation
GA Genetic Algorithm
LHS Left hand side
NACA National Advisory Committee for Aeronautics
PSO Particle Swarm Optimization
RHS Right hand side
RPV Remotely piloted vehicles
UAS Unmanned aerial system
VTL Vortex Lattice Method
1
1 INTRODUCTION
Even though aircraft propellers have been designed for over a century, investigations
of the performance and the design of propellers are as important as ever. The continued
interest in propeller propulsion is due to the fact that propellers are more efficient at low
and modest flight speeds when compared to turbofan propulsion. A relatively new field in
propeller design is the prediction and reduction of propeller noise, which is driven by the
public living close to airports or large wind turbines. Aircraft propeller noise is also of
special interest to the military and was first investigated in 1919 by Lynam and Webb.
[1]
This military interest in minimizing noise continues today with the recent fielding of low
flying Unmanned Aerial Systems (UAS) on reconnaissance missions.
This work describes an effort to optimize the propeller performance of small thin
airfoil propellers, which are found on electric powered UAS?s and Remotely Piloted
Vehicles (RPV). Previous work on 2D airfoil optimization using a Genetic Algorithm
(GA) has been investigated by Refs. [2] [4]. Fanjoy and Crossley
[3]
generated 2D airfoil
geometries using a Bspline method with 80 panels distributed over the airfoil surface
and used a panel method to do the aerodynamic analysis. Penalty functions were added to
provide a minimum thickness thus structural integrity is guaranteed. The two Genetic
Algorithm?s used are a binary tournament selection with uniform crossover and a real
encoded GA. The shortcoming is following: If there are cases for which the
2
analysis over predicts the desirable features of an airfoil, the population will move toward
those cases.
A 3D propeller analysis and optimization was done by Miller
[5]
. In his approach,
the Vortex Lattice Method (VLM) using a single panel in the chordwise direction was
employed to determine the aerodynamic propeller performance. The propeller geometries
were generated using an extended version of the classical method in which the blade
shank is the center of rotation. Airfoil camber is not considered due to the limitation of
using only a single chordwise panel. The work includes a subsonic noise analysis using
sectional forces located at the quarter chord to represent the blade loading and point mass
sources and sinks to represent blade thickness. The multivariable optimization method
used was developed by Powell
[6]
and is a conjugate directions technique. The propeller
optimization is limited because of constraints of the propeller geometry functions as well
as the constraints of the lifting line method to account for chordwise geometric shapes.
Olsen
[7]
described ship propeller optimizations using the VLM where the propeller blade
is replaced by a lattice of quadrilateral panels with constant circulation and the horseshoe
vortices follow regular helices. The propeller surface geometry is defined by a vector
which contains the rake, the skew of the midchord and the horizontal pitch angle. The
propeller optimization solves a variational problem to find the radial circulation
distribution which provides the lowest torque for a given thrust, thus generating
propellers with minimum losses. The propeller performance computation did not include
viscous effects and the propeller geometries were generated using simple cosine and sine
functions.
3
Hampsey
[1]
investigated the optimization of small wind turbines employing a three
dimensional panel method with distributions of source and doublet singularities over the
discretised blade surface. Airfoil geometries were generated by least squares Bspline
curve fitting through specified data points. The spanwise geometry was defined by
stringing a handful of scaled and rotated airfoils along the blade span and then skinning it
with a similar Bspline surface in the radial direction. The optimization scheme used was
a differential evolution algorithm which allows for multiobjective optimization. In
addition to the aerodynamic performance optimization, Hampsey included a wind turbine
blade structure model.
The main focus of the present work was:
? Develop a propeller program which accurately computes the propeller
aerodynamic performance parameters while being computationally
efficient;
? Implementation and development of a new way to generate propeller
geometries using the Class/Shape function Transformation (CST)
methodology;
? Increase computational accuracy by including viscosity, compressibility
and flow separation models;
? Implementation of an acoustic signature model.
The developed propeller performance program is attached to both a binary encoded GA
and a real coded GA.
A fundamentally new method to generate airfoil and propeller geometries based on a
Class/Shape function methodology is utilized. The VLM using lifting surfaces with
4
vortex ring elements has been used to determine the aerodynamic performance
parameters. In an effort to reduce computation times, a detailed wake extension
convergence model has been developed and a panel density investigation has been
completed. The focus of this work is the optimization of propellers with multiple
objectives, which has not been researched extensively. Finally the implementation of a
validated acoustic propeller signature model in combination with the VLM performance
analysis and the optimization scheme using a GA has been developed.
This investigation can be split up into four major areas: propeller geometry
determination, aerodynamic analysis, noise analysis and an optimization scheme. In each
Chapter, the background in the field begins with a discussion of the basic issues, followed
by a detailed description of the analysis method used, and concluded by the presentation
of the work done to verify the method. Optimization results include single and multipoint
performance and acoustic signature optimizations with performance tradeoff?s.
5
2 BASIC THEORY
The first and simplest method to predict propeller performance was developed by
Rankine
[9]
. The basic version is called axial momentum theory, in which the propeller is
replaced with a disc and the thrust is considered uniformly distributed. Even though the
theory was extended to account for swirl, friction losses
[10]
and compressibility
[11]
, it still
ignores many real world effects, which makes it useful only for predicting the upper
limits of propeller efficiency.
The blade element theory was formulated by Froude
[12]
and Drzewiecki.
[13]
Anderson
[14]
said ?Wilbur Wright was the first person to recognize ? that a propeller is
nothing more than a twisted wing.? Using information from their wind tunnel tests,
Wilbur Wright had invented ?blade element theory,? which is the idea that at each point
along the span, a propeller meets the air at a different angle and speed. Blade element
theory divides the propeller blade into elements, each of which is analyzed separately by
an analogy drawn between propeller and finite wing. The method takes axial and
associated rotational velocities at each spanwise section into account and uses 2D airfoil
coefficients to compute local and overall propeller performance. Even though this method
is more accurate than the axial momentum theory it does not consider the influence of the
propeller wake. Glauert
[15]
combined the momentum theory with the blade element theory
and assumed that the propeller is lightly loaded, hence slip stream contractions are small
6
and the radial components of the induced velocities are negligible. The Vortex Theory as
discussed above does not include tip losses. An improvement can be made using lifting
line theory which is discussed in the following section.
2.1.1 LiftingLine Theory
The liftingline theory was the first practical theory for incompressible, inviscid
flow to predict aerodynamic properties of unswept three dimensional wings.
[16]
The
theory was developed by Ludwig Prandtl and colleagues
[17]
between 1911 and 1918 and
is essentially an implicit potential flow solution. The wings are modeled as a single
bound vortex line, which is located at the quarter chord position and has a shed vortex
sheet extending to infinity to satisfy the Helmholtz theorem,
[16]
shown in Figure 2.1
below.
Figure 2.1 Liftingline model with single horseshoe vortex on finite wing
?
V
?
? c
b
?
?
?
.
Collocation point
7
To realistically simulate the downwash distribution of a wing, a large number of
horseshoe vortices are superimposed, each with a different bound vortex length located at
the quarter chord. Solution methods to the lifting line theory are described in numerous
texts including Katz & Plotkin
[16]
, Kuethe and Chow
[18]
and Anderson.
[14]
Figure 2.2 Distribution of horseshoe vortices, liftingline model
Figure 2.2 shows the geometric arrangement of the horseshoe elements over the wing
span with circulation ? changing along the liftingline.
The unknown circulation ? of the individual horseshoe vortices is obtained using
the BiotSavart Law
[16]
with the boundary condition, 0? =?nV
r
(zero flow normal to the
surface) at a predefined collocation point (also referred to as a control point). The
resulting form of the BiotSavart Law is:
3
4 r
rdl
dV
?
?
?
=
?
(2.1)
1
?
1
?
2
?
2
?
3
?
3
?
b/2
b/2
x
y
z
1
?
21
?+?
321
?+?+?
Lifting line
8
The local lift is obtained from the KuttaJoukowski Theorem
[18]
, which states that
lift, L per unit span, is directly proportional to the circulation ?, and is always
perpendicular to the local freestream velocity.
???=
?
VL ? (2.2)
The total lift is the summation of the local lift, L, which is perpendicular to the
free stream velocity. An illustration of the velocities and related forces is shown in Figure
2.3. To determine the induced drag, the resultant velocity, V, which has two components,
the free stream and the induced velocity, w
i
, must be determined. The induced angle, ?
i
,
is the angle between the free stream velocity and the resultant velocity. With the induced
angle of attack known, the induced drag is obtained as shown in Figure 2.3.
Figure 2.3 Induced drag due to induced downwash of the trailing vortices
The induced downwash velocity, w
i
, at any span location can be determined once
the circulation ? is known. This allows for the computation of the local flow angles on
V
?
x
y
V
?
w
i
?
i
???=
?
VL ?
( )
ii
VD ?? sin?????=
?
???= VL ?
V
9
the entire wing. A consequence of the downwash flow is that the direction of action of
each section?s lift vector is rotated relative to the freestream direction. The local lift
vectors are rotated backwards and hence give rise to a lift induced drag. By integrating
the component of section lift that acts parallel to the freestream across the span b, the
induced drag D
i
can be found as
()
?
?
????=
2/
2/
sin
b
b
ii
dyVD ?? (2.3)
2.1.2 Lifting Surface Theory
Due to the limitations of predicting accurately aerodynamic performance for only
straight medium to high aspect ratio wings, the LiftingLine Method was extended by
placing a series of lifting vortex lines along different wing chord stations.
[2], [14]
Each
lifting line has two associated trailing vortices which are parallel to the x axis. As the
location for lift calculations shifts downstream in the chordwise direction, starting at the
leading edge, additional superimposed trailing vortices are added. Thus vortex strength
dependency is gained in both the x and y direction. The vortex sheets cover both, the
spanwise and the chordwise direction and make up the lifting surface. Figure 2.4, a wing
with two chordwise vortex elements is shown.
10
Figure 2.4 Lifting surface method by horseshoe elements
As with the liftingline method, the unknown circulations ? must be determined
based on the flowtangency condition, ( )0? =?nV
r
, at all points on the wing, so that lift
and induced drag can be calculated. The central problem with lifting surface theory is to
solve for all the unknown circulations since circulation strengths are mutually dependent.
One numerical method is the Vortex Lattice Method (VLM) which was developed to
compute aerodynamic performance data of various wing geometries. Numerous
papers
[19] [22]
have been published which demonstrate the accuracy and applicability of
the VLM. In this method, the wing is covered by a net of bound and horseshoe elements
which is called a vortex lattice system. When the BiotSavart law and the flow tangency
condition is applied at all collocation points, a system of simultaneous equations is
obtained which then can be solved for the unknown circulation strengths, ?.
x
y
z
n
r
11
Lifting Surface Method by Vortex Ring or Quadrilateral Elements
In this method, vortex ring elements are distributed over the wing surface instead
of the horseshoe elements as previously shown. This simplifies the vortex model as well
as reduces computation times due to the lower number of vortex filaments. The boundary
condition on each ring vortex is the same as before: no flow through the wing at a control
point.
Since all vorticity based solutions are essentially implicit potential flow solutions,
flow separation, viscous and compressibility effects cannot be simulated and must be
taken into account separately, as discussed in Chapters 2.5 ? 2.7.
A single ring vortex consists of four straight bound vortex elements, which are
arranged in a geometric closed form. Positive circulation is defined according to the
righthand rotation rule, which is illustrated in Figure 2.5.
Figure 2.5 Vortex ring element arrangement
In previous efforts
[5]
the lifting surface solution by horseshoe elements was used,
which is computationally more expensive but yields similar results when compared to
vortex ring elements. The vortex geometry method used in this effort has been applied to
ship propellers
[7]
, aircraft propellers
[23]
and yacht sails
[24]
and is proven to give accurate
Element 2
Element 4
Element 1
?
Element 3
12
aerodynamic performance predictions. In the present work, a single propeller blade is
divided up into chordwise and spanwise panels. The corners of the panels are defined by
nodal points and the geometric center by a collocation point.
For thin propeller airfoils, a single layer of panels is placed on the mean camber
line. For thick airfoil propeller blades (>10%), the upper and lower blade surface is
covered with a single layer of vortex elements.
Vortex ring elements are placed on the propeller surface starting at the blade
leading edge and extending to the second last element in the chordwise direction. The last
chordwise panel elements are horseshoe elements with a spanwise bound vortex and two
trailers. The horseshoe trailers extend downstream to infinity following the pitch of the
last chordwise panel. Figure 2.6 illustrates the discretization of a propeller blade by
vortex elements.
13
Figure 2.6 Propeller blade grid discretization
The aerodynamic performance of the propeller is obtained, as described for the
lifting surface method with horseshoe elements by applying the boundary condition to
each panel that the flow normal to the panel surface is zero. For a detailed description of
the numerical computation of the aerodynamic performance of the propeller is given in
Section 2.3.
2.2 Propeller Wake
Generally there are two types of wake models, the prescribed wake
models
[23], [25], [26]
, which define the wake by a function and the free wake
methods
[24], [27], [28]
, which adjust the wake structure to account for wingtip roll up.
Blade trailing edge
Horseshoe element
Ring element
?
1
Collocation point
Nodal point
?
11
?
2
?
3
?
12
?
5
?
8
?
10
?
7
?
4
?
9
?
6
14
In this effort, the propeller wake geometry follows Goldstein?s
[25]
helical model, which
assumes lightly loaded propellers. Thus the interference flow of the vortex system is
small compared to the velocities of the blade. The pitch of the wake is constant and there
is no slipstream contraction. In Ref. [23], it is shown that aerodynamic propeller
performance can be accurately predicted using constant pitch wakes. To consider wing
tip rollup, a free wake model must be chosen. Most free wake models use an iterative
process to determine the exact wake geometry. The free wake models are not considered
within this optimization effort due to the increased computation time and minimal gains
in accuracy.
The wake of the propeller is simulated by horseshoe vortices which extend from
the last chordwise panel extending two revolutions downstream. More details on
propeller wake extensions are discussed in Section 4.2. The wake gradient follows the
trailing edge panel slope to closely approximate the Kutta condition.
Figure 2.7 Propeller wake geometry
z
x
Horseshoe trailer
sare divided up by
10
o
rotation
elements
?
First streamwise horseshoe
trailer
Collocation Points
Nodal Points
15
The trailers are divided up into 10 degree rotational increments which are
implemented in Equations 2.42. The y and z coordinates of the horseshoe trailers are
computed similarly but offset by 90 degrees. The x coordinates are determined by the
gradient of the last chordwise vortex elements multiplied by the distance traveled during
a 10 degree rotation. Once the gradient, m, of the last chordwise panels is obtained and
the angles ? of each last chordwise nodal point with the yaxis are determined, the x, y
and zcoordinates of the horseshoe vortex elements can be written as:
m
i
jncynodaljncxnodaljixhorseshoe ?
?
???+=
360
)(10
),(_2),(_),(_ ?
(2.4)
)
180
),()(10
cos()),(_),(_(),(_
5.022
jnci
jncznodaljncynodaljiyhorseshoe
?+?
?+=
(2.5)
)
180
),()(10
sin()),(_),(_(),(_
5.022
jnci
jncznodaljncynodaljizhorseshoe
?+?
?+= (2.6)
with m being the gradient of the last streamwise propeller panel
),1(_),(_
),1(_),(_
jncznodaljncznodal
jncxnodaljncxnodal
m
??
??
= (2.7)
and ? the angle in the xyplane associated with the last chordwise nodal point.
?
180
)),(_),(_(
),(_
sin),(
5.022
?
?
?
?
?
?
?
=?
jncznodaljncznodal
jncznodal
ajnc
(2.8)
16
The counter i runs to 72 to generate two full wake rotations and j from 1 to nr which are
the spanwise number of nodal points.
Figure 2.8 Propeller wake horseshoe trailer discretization
Initially, a single propeller blade including horseshoe trailers is generated. The geometry
is then rotated about the origin depending on the required number of propeller blades.
2.3 Propeller Performance Prediction by VLM
For the computation of the propeller loading, three different methods are available.
The methods differ in accuracy and computation time. One method is the Trefftz plane
analysis
[29]
in which it is assumed that the induced flow far downstream from the
propeller does not depend on the streamwise coordinate and thus is considered to be two
dimensional in all planes normal to the trailing vortex sheet. The induced velocities of the
X
Y
Z
10 degree rotation
increments
?
17
bound vortex elements are neglected since they are determined far downstream in a cross
plane normal to the wake. Thus, the analysis is essentially reduced to finding the flow
field of a two dimensional vortex sheet to satisfy certain boundary conditions. The second
method is the pressure summation method
[30]
which uses the velocity difference on the
upper and lower propeller surface to obtain pressure differences. The total force acting on
the propeller is obtained by surface integration. The third method is the Joukowski
method which determines the forces on each discrete vortex line in the lattice and through
summation obtains the total propeller loading. There are different ways of determining
the induced velocity which are described in Refs. [31] [33]. A comparison of the three
computational methods with respect to accuracy and computational effort, applied to
yacht sails is found in Ref. [24]. The analysis concludes that the Joukowski method has
good accuracy and provides the force distributions in both the spanwise and the
chordwise directions. The Joukowski method was also used by Cheung and generated
good accuracy in the performance prediction of propellers. Based on this analysis, the
Joukowski method was selected to determine propeller performance parameters.
LiftingSurface Solution by Vortex Ring Elements:
The Propeller blade is divided into chordwise and spanwise vortex ring elements.
The ring elements cover the entire blade except the last chordwise panels, which are
horseshoe elements. The horseshoe trailers leave the propeller blade with the same pitch
as the last chordwise panels as discussed in Section 2.1. The velocity induced by a vortex
element of the strength ?
with a length of dl on a collocation point P(x,y,z) is calculated
by:
18
( )
3
10
10
4
rr
rrdl
V
?
??
?
?
?
=?
?
(2.9)
with ( )
10
rr ? being the vector from the collocation point to the midpoint of the vortex
element. The circulation ? and the induced velocities are initially unknown and the last
part of Equation 2.9 contains the influence coefficients which are obtained with:
( )
3
10
10
4
1
inf_
rr
rrdl
coeff
?
??
?
?
=
?
(2.10)
The influence coefficients describe the geometric relation between the length of a vortex
element and the point at which the induced velocity is to be determined. Figure 2.9
illustrates the geometric relations of a vortex element and the point of velocity
evaluations.
Figure 2.9 Induced velocity at point P(x,y,z) by a single vortex element
dl
ds
?
P(x,y,z)
r
0
r
1
Origin
r
0
? r
1
?
?V
Vortex element
19
For a single quadrilateral vortex element the total influence coefficients in the x, y
and z direction are the summations of the influence coefficients in x, y and z direction of
the four straight bound vortex legs. Cheung
[23]
determines them by
()()()()[]
()()()()
()()()()Ixxyyyyxxcoeff
Izzxxxxzzcoeff
Iyyzzzzyycoeff
startendcolocationstartstartendcolocationstart
i
z
startendcolocationstartstartendcolocationstart
i
y
startendcolocationstartstartendcolocationstart
i
x
????????=
????????=
????????=
?
?
?
=
=
=
4
1
4
1
4
1
inf_
inf_
inf_
(2.11)
with
()
?
?
?
?
?
?
?
?
?
+?+
+
?
??
=
c
b
cba
ba
bca
I
2
1
2
(2.12)
and
()()( )
()( )()( )()( )
()()
222
222
ncollocatiostartncollocatiostartncollocatiostart
ncollocatiostartstartendncollocatiostartstartendncollocatiostartstartend
startendstartendstartend
zzyyxxc
zzzzyyyyxxxxb
zzyyxxa
?+?+?=
???+???+???=
?+?+?=
(2.13)
with I, a, b, c being variables which describing the geometric relations of a vortex
element with respect to a collocation point.
The same applies to the horseshoe elements, with the difference that the
summation of the influence coefficients consists of one bound vortex element and two
trailers which extend two revolutions downstream and which are divided into 10?
sections. It is of importance to follow the right hand rule when defining start and end
20
point of a single vortex element. The direction of the circulation within the panels cannot
change. Also once a direction of the circulation is selected, all following panel
circulations must be in the same rotational direction.
The sequence of calculating the influence coefficients with respect to one
collocation point is as follows: The first collocation point on the first panel of the first
propeller blade is selected and the influence coefficients in the x, y and z direction of the
first vortex ring element is computed. The obtained x, y and z values of the influence
coefficients are then stored in the first location of the matrices A
x 1,1
, A
y 1,1
, and A
z1,1
. The
next step is to keep the first collocation point and determine the influence coefficients of
the second vortex ring element, which is the second panel in the spanwise direction of the
propeller blade. Again the influence coefficients are stored in the matrices at the location
A
x 1,2
, A
y 1,2
, and A
z 1,2
. This process is repeated until all influence coefficients in x, y and
z direction are determined with respect to all vortex ring elements and horseshoe
elements of all propeller blades. Then the process is repeated with the collocation point of
the second vortex ring element selected. The individual influence coefficients are
determined and stored in the matrices A
x 2,n
, A
y 2,n
, and A
z 2,n
. The computations of all the
influence coefficients ends with the last collocation point on the last propeller blade and
the last horseshoe panel on the last propeller blade. Figure 2.10 the scanning sequence of
the propeller panels is shown.
21
Figure 2.10 Scanning sequence on propeller blade panels
All influence coefficients in x, y and z directions are stored in the matrices A
x
, A
y
,
and A
z
with the rows representing the dependencies with respect to a single collocation
point and the columns representing the dependencies of a single panel with respect to all
collocation points of all propeller blades. Thus the BiotSavart law becomes:
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
?
?
?
?
=?
MNMNMNMNMN
MN
MN
AAA
AAA
AAA
V
2
1
,2,1,
,22,21,2
,12,11,1
4
1
?
(2.14)
To calculate the unknown circulations, ?, the boundary condition that the velocity at each
collocation point normal to the individual panel is zero, or:
?
x
z
i
i+1
i+2
j
j+1
j+2
j+3
M
N
Counter k=M x N
Tip
Root
22
0? =?nV (2.15)
The surface normal of each quadrilateral and horseshoe element is determined as the
cross product of the diagonal vectors of each panel and is given as:
kk
kk
k
BA
BA
n
?
?
= (2.16)
Figure 2.11 Normal vector on propeller blade panel
The boundary condition applied to each collocation point is:
[]()()
kkkkonsetzijzijyijyijxijxij
M
i
k
N
j
VnAnAnA ??? cossin
11
???=?+?+???
??
==
(2.17)
with ?
k
being the angle of attack of the onset velocity V
onset_k
with the x axis, ?
k
the
angle of each panel normal with the zaxis and ? the orientation of the panels. Figure 2.12
illustrates the geometric relations in the propeller coordinate system.
i
i+1
i+2
j
j+1
j+2
A
k
B
k
n
k
23
Figure 2.12 Propeller boundary condition and orientation angles
The onset velocity V
onset
changes with each spanwise panel, since it is a combination of
the free stream velocity and the rotational speed of the propeller.
()( )
2
1
2
,
2
,
jionset
rVV
ji
?+=
?
? (2.18)
Equation 2.14 is a system of linearly dependent equations which are solved for the
unknown circulations ? with a Gauss Seidel solver. Section 4.1 describes in detail the
Gauss Seidel Solver.
[ ][ ][ ]BxA =? (2.19)
with A being the influence coefficients
[ ]
zijzijyijyijxijxij
M
i
N
j
nAnAnAA ?+?+?=
??
==11
(2.20)
?
onset
x
z
x
z
?
y
z
?
n
n
 +
V
onset
24
and x the unknown circulations:
[]
k
x ????= ,...,,
321
(2.21)
and B the boundary condition zero flow normal to the panel surface:
( ) ( )
()()
()()
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
=
kkkkonset
onset
onset
V
V
V
B
???
???
???
cossin
.....
cossin
cossin
2222
1111
(2.22)
The returned vector x contains the circulations ?
k
of each propeller blade panel as shown
in Figure 2.13 below.
Figure 2.13 Vortex circulation sequence on a two bladed propeller
With the returned circulations, the normal force of each quadrilateral and horseshoe
element is obtained by the use of the KuttaJoukowski theorem
[16]
. KuttaJoukowski
theorem states that the forces are always perpendicular to local onset flow, which is the
combination of the free steam velocity and the local rotational speed of the propeller. The
forces for the panels on the propeller blade, except for the forces on the leading edge
panels are:
?
?
1 ?
2
?
3 ?4 ?5 ?
6
?
7
?
8
25
( )
jjijijionsetji
yVF ???????=?
?? ,1,,_,
? (2.23)
For the leading edge panel (i=1) the forces are:
( ) yVF
jionsetji
?????=?
? ,,
? (2.24)
For the computation of the thrust and required power of the propeller, the induced
downwash at each collocation point must be determined. The influence coefficients from
only the trailing vortex segments with respect to each collocation point are computed in a
manner similar to the approach used for the quadrilateral panels. These coefficients are
stored in the matrices B
x
, B
y
, and B
z
. Figure 2.14 shows a propeller blade with related
horseshoe elements.
26
Figure 2.14 Trailing vortex segments responsible for induced downwash
Since the circulations ?
k
are known from the previous calculations, the induced
downwash on each collocation point is determined by:
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=
?
?
?
?
?
?
?
?
?
?
?
?
kkkkk
k
k
indk
ind
ind
BBB
BBB
BBB
w
w
w
2
1
,2,1,
,23,21,2
,12,11,1
2
1
...
............
...
...
...
(2.25)
The circulations used for the computations of the induced downwash are obtained from
the Gauss Seidel solution.
Trailing
vortex
segments
Propeller
geometry
10? horseshoe
wake segments
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Collocation
point
27
With the induced downwash, the induced drag on each panel is determined with:
() 1,
,,1,,_,_
>????????=?
?
iywD
jijijijiindjiind
? (2.26)
() 1,
,,,_,_
=??????=? iywD
jijijiindjiind
? (2.27)
The total induced drag is the summation of all individual panel contributions of one
propeller blade multiplied by the number of the propeller blades.
nxblDD
N
j
jiind
M
i
?
?
?
?
?
?
?
?=
??
== 1
,_
1
(2.28)
To compute the thrust and the required power of the propeller the forces and the induced
drag forces on each panel must be rotated in the propeller coordinate system as illustrated
in Figure 2.15.
28
Figure 2.15 Total thrust and drag of the propeller
The angle of rotation is defined by the onset flow angle ?
j
which changes in the spanwise
direction. Thus the propeller thrust contributions of the individual panels are:
( ) ( )
jjiindjjiji
DFT ?? cossin
,_,,
?????=? (2.29)
V
?
x
z
?
ind
??r
i
w
ind
?F
?
ind
V
onset
V
local
?D
ind
?T = ?F sin(?)  ?D
ind
?cos(?)
?D = ?F cos(?) + ?D
ind
sin(?)
?
?D
ind
= ?V
local
?sin(?
ind
)
?F = ?V
local
?cos(?
ind
) = ?V
onset
?
Rotor shaft
?
29
and the drag contributions are:
( ) ( )
jjiindjjiji
DFD ?? sincos
,_,,
??+??=? (2.30)
The total thrust of the propeller is obtained by the taking the summation of the individual
panel contributions multiplied by the number of the propeller blades:
nxblTT
N
j
ji
M
i
total
?
?
?
?
?
?
?
?=
??
== 1
,
1
(2.31)
To compute the required propeller power, the torque ?Q
i,j
generated by the individual
drag component of the panels, is multiplied by the distance from the collocation point to
the axis of rotation. The summation of the torque contributions of the individual panels
multiplied by the number of the propeller blades nxbl will determine the total propeller
torque Q
total
.
`
nxblQQ
rDQ
N
j
ji
M
i
total
jijiji
?
?
?
?
?
?
?
?=
??=?
??
== 1
,
1
,,,
(2.32)
The required propeller power is the total torque multiplied by the rotational speed ? of
the propeller.
??=
total
QP (2.33)
30
With the propeller thrust and power known, the thrust and power coefficients are:
42
dn
T
c
total
t
??
=
?
?
(2.34)
53
dn
P
c
p
??
=
?
?
(2.35)
with n being the revolutions in rad/sec and d the propeller diameter. The propeller section
lift coefficient, c
l
, and the propeller efficiency are other parameters which are commonly
used in the performance analysis.
()
cV
rc
il
?
??
=
?
2
(2.36)
P
TV
total
?
=
?
? (2.37)
2.4 Propeller Geometry
In general the geometric definition of propellers can be reduced to the airfoil type,
the chord length, the sweep and the angle of attack function. In this effort the CST
(Class/Shape Function Transformation) method for defining general aerodynamic shapes
has been applied to the generation of random propeller shapes which is discussed as
follows.
2.4.1 2D Airfoil Geometry
Several methods have been developed to represent general airfoil shapes for the
use in aerodynamic design optimizations. The goal is to define a simple analytic function
31
which efficiently describes the entire design space for airfoils. Some of the methods used
in previous efforts can be found in Ref. [34] [39].
The method chosen is based on the work of Ref. [40] and Ref. [41]. In this
approach a simple and well behaved analytic unit shape function, based on Bernstein
polynomials (BP) is introduced. This shape function directly controls key airfoil
parameters including leading edge radius, thickness and trailing edge angle. A unit class
function is added to the unit shape function to allow for a wider variety of general body
shapes.
The streamwise class function in the design space is defined as
( ) ( ) [ ]
211
2
1
NNN
N
C ??? ??= (2.38)
with ? being the fraction of the local chord. In the physical space the unit class function
is
21
1
NN
c
x
c
x
c
x
C
?
?
?
?
?
?
??
?
?
?
?
?
?
=
?
?
?
?
?
?
(2.39)
where c is the chord and x is the local chord fraction having a range from 01. The first
term of Equation 2.39 defines the shape of the airfoil leading edge and the second term
can be used to ensure a sharp trailing edge. If N
1
= 0.5 and N
2
= 1.0 a round airfoil
leading edge and a sharp trailing edge, is enforced. Figure 2.16 shows some examples for
different values of the class function coefficients N
1
and N
2
. Further geometry class
determinations due to N
1
and N
2
variation can be found in Ref. [40].
32
Figure 2.16 Class function 2D design space
The unit shape function is defined by a BP of the order n with the variable x ranging
from 01.0. The BP?s were chosen due to the mathematical property of ?Partition of
Unity? as described in Ref. [41]. The first term of the shape function defines the binomial
coefficients with increasing order n of the BP. For each order n of BP there exist n+1
terms which are defined by Equation 2.40.
()
()
?
=
?
?
?
?
?
?
?
??
?
=
n
r
rn
nr
c
x
rnr
n
xS
0
,
1
!!
!
(2.40)
The individual terms of the BP?s can be illustrated by means of a Pascal?s triangle as
shown in 2.17.
N
1
=1.0 , N
2
=1.0 N
1
=0.5 , N
2
=1.0 N
1
=0.5 , N
2
=0.5
c
x
33
Figure 2.17 Bernstein polynomial representation of the unit shape function
The first term of Equation 2.40, when written in polynomial form as shown in
Figure 2.17, defines the leading edge radius and the last term the trailing angle. The other
terms are shaping terms which do not influence airfoil leading or trailing edge shape.
The entire airfoil can be represented by one upper and one lower unit class function (if
thickness is included) multiplied by a unit shape function with Bernstein polynomials. In
Ref. [41] it was found that a BP of the order n of 69 matched the airfoil geometries and
aerodynamic forces. It was further suggested that for optimization purposes the order n
could be lowered to 4 which reduces the design variables and thus improves computation
times. Based on the findings in Ref. [41] the order of the BP is set to n=4 for all further
propeller optimizations.
Equations 2.41 and 2.42 define the upper and lower airfoil geometry. The first two terms
are the class function which sets the leading and trailing shape of the airfoil. The
1
x?1 x
2
)1( x?
3
)1( x?
4
)1( x?
5
)1( x?
xx)1(2 ?
2
x
xx
2
)1(3 ?
2
)1(3 xx?
3
x
xx
3
)1(4 ?
22
)1(6 xx?
2
)1(4 xx?
4
x
xx
4
)1(5 ?
23
)1(10 xx?
32
)1(10 xx?
4
)1(5 xx?
5
x
n=0
n=1
n=2
n=3
n=4
n=5
Leading Edge Radius Trailing Edge
Boattail Angle
?Partition of Unity?
() 1
0
,
=
?
=
n
r
nr
xS
?Shaping Terms? which do not affect
leading or trailing edge geometry
34
remaining terms are the BP for n=4 with the individual coefficients as shown in Figure
2.17. When thin airfoil propellers are considered only the upper airfoil geometry function
is used to define the propeller shape. A single layer of vortex elements is placed onto the
geometry function to represent the mean chord of the thin bladed propeller.
Upper surface definition:
()
[
]
4
5
23
4
22
3
3
2
4
1
114
1614
1
21
?
?
?
?
?
?
??+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
=
c
x
A
c
x
c
x
A
c
x
c
x
A
c
x
c
x
A
c
x
A
c
x
c
x
xy
NN
upper
(2.41)
Lower Surface definition:
()
[
]
4
10
23
9
22
8
3
7
4
6
114
1614
1
43
?
?
?
?
?
?
??+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?=
c
x
A
c
x
c
x
A
c
x
c
x
A
c
x
c
x
A
c
x
A
c
x
c
x
xy
NN
lower
(2.42)
with c being the maximum chord, x the local variable ranging from 0 ? 1 and A
110
the
coefficients (with a range of 01). This 2D class/shape function requires 14 variables for
a thick airfoil (upper and lower surface) and 7 variables for a thin propeller with a single
35
surface geometry function. Figure 2.18 illustrates two random airfoils from the design
space of the class/shape function, one thick and one thin.
Figure 2.18 Class/Shape function single and dual surface geometries
2.4.2 3D Propeller Geometry
The class/shape function methodology of representing 2D airfoils is extended here
to define general 3D propeller shapes. This approach is based on the work done by
Ref. [40] and [41] and then extended further to allow for a more open design space in the
creation of general propeller geometries. To accomplish spanwise geometry variation the
coefficients A
110
and N
14
for the top and the bottom airfoil side must be dependent on
the local spanwise position. Thus for each of the dependent variables A
110
a BP function
with n=3 describing the spanwise coefficient variation, is defined. For the coefficients of
the Class function, N
14
, a BP function with n=1 is introduced, to provide smooth
spanwise coefficient transitions. Thus the variables for the upper and lower airfoil
geometry are:
Class/Shape function single mean chord airfoil geometry
Class/Shape function upper lower airfoil geometry
36
[
]
[
]
?
?
?
?
?
?
?+
?
?
?
?
?
?
=
?
?
?
?
?
?
?
?
?
?
?
?
?+
?
?
?
?
?
?
=
?
?
?
?
?
?
?
?
?
?
?
?
??+
?
?
?
?
?
?
?
?
?
?
?
?
?
???
+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?=
?
?
?
?
?
?
?
?
?
?
?
?
??+
?
?
?
?
?
?
?
?
?
?
?
?
?
???
+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?=
?
?
?
?
?
?
+
+
++
+
++
+
b
y
nl
b
y
nl
b
y
Nu
b
y
nu
b
y
nu
b
y
Nu
b
y
al
b
y
b
y
al
b
y
b
y
al
b
y
al
b
y
Al
b
y
au
b
y
b
y
au
b
y
b
y
au
b
y
au
b
y
Au
jjj
jjj
ii
iii
ii
iii
1
1
113
13
113
13
2
2
3
15
2
10
2
5
3
3
15
2
10
2
5
3
(2.43)
b is the unit propeller radius, y is the local spanwise station, and i and j are the variable
counters from 15 and 12.
The total number of variables to describe the upper and lower airfoil geometry is
44. Twenty for each, the upper and lower airfoil shape function and 4 variables to define
the two class functions. The geometry of a single surface is described by the unit
class/unit shape function for the upper airfoil side as:
()
( ) ( )
[ () () ()
() () ]
4
5
23
4
22
3
3
2
4
1
114
1614
1,
21
?
?
?
?
?
?
??+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
=
c
x
yA
c
x
c
x
yA
c
x
c
x
yA
c
x
c
x
yA
c
x
yA
c
x
c
x
yxy
yNuyNu
upper
(2.44)
An angle of attack (AOA), chord length and sweep variation function is added to
the unit class/shape function to further extend the propeller design space. All three
functions are defined in the same manner as the unit shape function variables A
15
, with
37
the addition of a multiplication factor to define the max value of the function. Each
function requires 5 additional variables, four to define the BP with n=3 and one to set the
max value of the function.
The angle of attack, chord length and sweep function can be written as:
() [
]
3
4
2
3
2
2
3
15
113
13
?
?
?
?
?
?
??+
?
?
?
?
?
?
?
?
?
?
?
?
?
???
+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
??=
b
y
aoa
b
y
b
y
aoa
b
y
b
y
aoa
b
y
aoaaoayAOA
(2.45)
() [
]
3
4
2
3
2
2
3
15
113
13
?
?
?
?
?
?
??+
?
?
?
?
?
?
?
?
?
?
?
?
?
???
+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
??=
b
y
chord
b
y
b
y
chord
b
y
b
y
chord
b
y
chordchordyChord
(2.46)
with aoa
5
and chord
5
being the max angle of attack and max chord length.
() [
]
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??+
?
?
?
?
?
?
?
?
?
?
?
?
?
???
+
?
?
?
?
?
?
?
?
?
?
?
?
?
???+
?
?
?
?
?
?
??=
rootatchord
b
y
sweep
b
y
b
y
sweep
b
y
b
y
sweep
b
y
sweepsweepySweep
3
4
2
3
2
2
3
15
113
13
(2.47)
The variable, sweep
5
, defines the maximum sweep in multiples of the propeller root
chord. Figure 2.19 illustrates how the sweep function influences the propeller shape. This
brings the total number of variables to 66. Sixty three for the blade geometry, one for the
propeller radius, one for number of propeller blades and one for the position, which the
individual spanwise airfoil sections are rotated about.
38
Figure 2.19 Propeller sweep functions based on CST method
In this effort, propeller optimizations for small electric powered UAS are investigated
thus only thin airfoil propellers are considered. This eliminates one surface geometry
function and reduces the number of variables to 40. An example of a single surface
propeller blade with wake is shown in Figure 2.20.
Propeller root
Propeller tip
Sweep backwards
39
Figure 2.20 Thin propeller blade geometry by Class/Shape function scheme
2.5 Compressibility
Even though only small thin bladed propellers are considered in this effort, air
density variations, ?, due to Mach numbers higher than M > 0.3 can have noticeable
effects on the propeller performance. The general, PrandtlGlauert compressibility
correction uses the Mach number to correct pressure, lift and moment coefficients at
Mach numbers up to M=0.7 for thin airfoil and at small angle of attack. The derivation of
the compressibility correction is described in Ref. [14]. Gennaretti and Morino
[42]
showed
that at higher Mach numbers the PrandtlGlauert rule greatly over predicts the pressure,
moment and lift coefficients at helicopter rotor tip sections when applied directly to
incompressible results from lower Mach numbers as shown in Figure 2.21. Further
background information on compressibility is described in Ref. [31] [43] [44].
40
Figure 2.21 Comparison of lift coefficient for incompressible and compressible
flow with direct application of the PrandtlGlauert rule (Source
from Ref. [45] )
Szymendera
[45]
derived the PrandtlGlauert rule starting with the continuity equation
based on Ref. [31], [43], [44] and then applied it to the BiotSavart law. This resulted in a
domain stretching in the flow direction which is dependent on the Mach number and has
a direct effect on the induced velocity computation. Equation 2.48 and Figure 2.22 show
the relation between domain stretching and Mach number.
2
'
1 M
r
r
p
p
?
= (2.48)
As shown in Fig. 2.21 compressible results from Szymendera show better agreement with
the coefficients than the results computed directly by applying the PrandtlGlauert rule.
? Experiment
? ? Incompressible Method
 PrandtlGlauert Correction
____ Compressible Method
41
Figure 2.22 PrandtlGlauert domain stretching in flow direction
The implementation of compressibility into the propeller performance program is
done after the propeller geometry is determined based on the input parameters. Since the
tangential velocity increases along the propeller span, the stretching is more pronounced
at the tip where the velocities are the highest, as illustrated in Figure 2.23 below.
.
r
1
r
2
r
0
r
1
?
r
1
?
r
0
?
Stretched domain analogous to
incompressible flow ?
Actual domain analogous to
compressible flow
P
r
p
r
p
?
42
Figure 2.23 PrandtlGlauert domain stretching of vortex ring elements
The propeller performance code was run with and without compressibility taken into
account. The propeller diameter was 0.3 [m] and the propeller rpm was 10,000, which
results in a tip Mach number of 0.47. The section lift coefficient for both the
compressible and the incompressible solution is shown in Figure 2.24. This result agrees
with the findings of Gennaretti and Morino
[42]
. The section lift coefficient is higher for
the compressible case since the domain shrinking based on higher Mach numbers results
in a smaller chord which in return increases c
l
.
Stretched domain,
incompressible
Actual domain,
compressible
43
Radial position r/R []
S
e
ct
io
n
lif
t
c
o
e
f
f
i
c
i
e
n
t
C
l
[

]
0.2 0.4 0.6 0.8
0.06
0.08
0.1
0.12
0.14
Figure 2.24 Spanwise lift coefficient normalized by blade radius
To determine the influence of the compressibility on the overall performance of
small thin airfoil propellers at lower Mach numbers considered in this effort, a
second investigation has been done. A generic propeller was selected as shown in
Figure 2.25, with a diameter of 0.25 [m] operating at 36 [m/s] free stream velocity.
The propeller speed was raised in steps starting from 4000 to 10000 revolutions per
minute. The parameters of interest are the propeller efficiencies which are shown in
Figure 2.26 for the compressible and incompressible case. It is observable that at
Incompressible
Compressible
44
Mach numbers above 0.25 efficiencies are slightly higher in the compressible case.
The higher the tip Mach number the larger the discrepancy between incompressible
and compressible solution.
X
Z
Y
Figure 2.25 Generic two bladed propeller with wake
45
Mach number []
P
r
o
p
e
lle
r
E
f
f
ic
ie
n
c
y
?
[%
]
0.1 0.15 0.2 0.25 0.3 0.35 0.4
65
70
75
80
Compressible solution
Incompressible solution
Figure 2.26 Compressibility effects vs propeller tip Mach number
Most of today?s electric propulsion systems used on small UAS?s such as the US
Army?s Raven system (manufactured by AeroVironment) use propeller speeds in cruise
below 8000 revolutions per minute and propeller diameters of 0.20.4 [m]. This generates
a maximum propeller tip speed in the range of M=0.3 and small propeller efficiency
differences, of about 1% between incompressible and compressible solution, are
noticeable. Figure 2.26 illustrates the reduced efficiency when compressibility is
considered at different propeller speeds.
46
2.6 Viscous Model
Since the VLM is an inviscid solution technique the physically accurate treatment
of viscous flows is beyond its capability. Inviscid aerodynamic solvers predict lift and
drag well for high Reynolds numbers (Re) and moderate angle of attacks, but these
predictions become increasingly less accurate as the Re decreases and the sectional angle
of attack increases. Viscousinviscid interaction methods
[46] [47]
increase the range of the
prediction accuracy of the boundary element methods while maintaining computational
economy, but they are beyond the scope of this work. Most of the viscousinviscid
interaction methods use the solution obtained by a boundary method and apply a
numerical boundary layer model to account for viscous effects and flow separation. In
this effort, a simple Reynolds number based model is implemented in an attempt to
account for the effects of viscosity as it varies in the spanwise direction. Since the
propeller blades considered in this effort all have thin airfoils, a flat plate boundary layer
method
[14]
is used to determine the skin friction coefficients c
f
. It is assumed that laminar
flow exists over the first 15% of the chord and that the remaining 85% the boundary layer
flow is turbulent.
[48]
Thus the skin friction coefficients for laminar and turbulent flows
are:
c
arlaf
c
Re
328.1
min_
= (2.49)
5
1_
Re
074.0
c
turbulentf
c = (2.50
47
with Re
c
being the Reynolds number at the spanwise collocation points with respect to
the propeller section chord length c.
?
? cV
onset
c
??
=Re (2.51)
Once the Reynolds number is determined, the skin friction coefficient on each spanwise
collocation position is computed. The viscous drag is then obtained from:
[ ]
fonsetvisc
cSVD ????=
2
2
1
? (2.52)
with S being the summation of all sectional wetted areas at one spanwise station. With
the viscous drag known at each spanwise propeller blade position, the torque and the
power increase due to viscous effects are:
nxblQP
rDQ
viscvisc
iviscvisc
??=
?=
?
(2.53)
The viscous torque is multiplied by 2.0 to account for both the top and the bottom
side of the propeller surface. Finally the total power due to viscous effects is calculated
by multiplying the viscous torque by the number of propeller blades.
To gain confidence in the viscous flow modeling, a NACA 109622 straight blade
propeller was simulated and compared with results obtained from lifting line method with
and without viscous effects considered. To account for viscosity, Cheung
[23]
used drag
data obtained from Ref. [49] of a NACA16004 airfoil section. Figure 2.27 shows both the
results of from Cheung obtained with lifting line method and the flat plate viscous
48
computation method shown above. There is good agreement between the two inviscid
(lifting line and lifting surface) methods and small discrepancies in the results of the two
viscous methods at lower advance ratios.
0.1
0.15
0.2
0.25
0.3
0.35
1.7 1.8 1.9 2 2.1 2.2
Advance Ratio J []
Po
w
er C
o
efficien
t c
P
[]
Lifting Line Method no Thickness and cd=0
Lifting Line Method no Thickness with Drag Data
Lifting Surface Method no Thickness and cd=0
Lifting Surface Method viscous no thickness
Figure 2.27 Spanwise lift coefficient normalized by blade radius
[23]
2.7 Flow Separation Model
The flow field associated with propellers is highly three dimensional and unsteady
in nature and, as yet, not well understood.
[50] [53]
Of particular difficulty is the prediction
of propeller performance at high angles of attack when flow separation occurs. In recent
years, new efforts have been made to better predict post stall behavior since flow
separation on inboard regions of wind turbines are common during much of their
49
operational times.
[54] [59]
Most of the developed methods for predicting flow separation
use an empirical or semi empirical stall model due to the complexity of the flow. Tangler
and Kocurek
[60]
developed a global post stall method based on the Viterna?s
equations.
[61] [62]
The disadvantage is that the slope and the maximum value of the
aerodynamic coefficients must be known for the model to be used. Further if the
aerodynamic coefficients do not match the behavior of a flat plate, the method becomes
unusable up to 20? angle of attack. In Ref. [60], it was demonstrated that even when thick
airfoil sections are considered, post stall behavior of c
l
/c
d
of a rotor, beyond 20 degrees
angle of attack, follows flat plate theory. In this effort, the first part of the Viterna
equations is used, since it predicts the maximum section drag coefficient as a function of
the aspect ratio of the propeller. The maximum global drag coefficient prediction, shown
in Equation 2.54, was used by Ref. [60] for sectional post stall predictions and showed
good agreement when compared to experimental data obtained from a rotor with five
spanwise pressure sensors:
ARc
d
?+= 018.011.1
max_
(2.54)
For the computation of the stall and post stall behavior of the flow, data obtained
from a NACA 0009 symmetrical airfoil
[63]
were used to generate the empirical post stall
model. The NACA 0009 airfoil was selected because experimental data exist for it and
since it is still considered a thin airfoil which matches the thin propeller blade simulation
efforts of this work. This was done with a simple curve fitting of the lift and drag
coefficient slopes starting at ~11? and going up to 90? angle of attack. Figure 2.28 and
2.29 illustrates the curve fitting of the lift and drag coefficients. In Appendix F the figures
50
generated from the experimental data from Ref. [63] are shown. If a single vortex panel
exceeds 13? angle of attack, the flow separation model is applied to the corresponding
panel of the propeller blade and new lift and drag forces are computed. This allows for
performance computations of partially stalled propellers. During most of the
optimizations, the panel angles of attack stay below 13?, since at larger Angle of Attack
propeller efficiencies degrade. Thus the empirical flow separation model is mainly used
when fixed pitch propellers designed for high cruise speeds are analyzed at low free
stream velocities where local panel angles of attack exceed 13 degrees.
y = 0.0000000108x
5
+ 0.0000029014x
4
 0.0002898088x
3
+ 0.0127624897x
2
 0.2301008464x + 2.1352985456
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
10203040506070809010
Angle of Attack [deg]
Lift
C
o
ef
fici
en
t
C
l
[

]
Data experiment
Polynomial trend line
Figure 2.28 Lift coefficient of NACA 0009 in post stall
51
y = 0.00000290x
3
+ 0.00034614x
2
+ 0.00345460x  0.00587855
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
10203040506070809010
Angle of Attack [deg]
Dr
ag coef
f
i
ci
ent
Cd [

]
Data experiment
Polynomial trend line
Figure 2.29 Drag coefficient of NACA 0009 in post stall
The NACA109622 straight blade propeller was used to quantitatively verify the
effects of the flow separation model. Figure 2.30 shows that below an advance ratio of
1.3 certain propeller regions are above 13? angle of attack, which activates the flow
separation model, thus reducing the lift and increasing the required power. Figure 2.30
and 2.31 illustrate two runs with the propeller performance program where lift and power
with and without the flow separation model activated are investigated.
52
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
0 0.5 1 1.5 2 2.5
Advance Ratio J []
Thr
us
t T [N
]
Trust of NACA propeller no flow separation
Thrust of NACA propeller with flow separation
Figure 2.30 Post stall thrust behavior of a NACA109622 propeller with empirical
stall model
2.0E+05
2.5E+05
3.0E+05
3.5E+05
4.0E+05
4.5E+05
5.0E+05
00.511.522.5
Advance Ratio J []
Po
wer
P [W]
Power of NACA propeller no flow separation
Power of NACA propeller with flow separation
Figure 2.31 Post stall power behavior of a NACA109622 propeller with empirical
stall model
53
2.8 Acoustic Signature Model
The aerodynamic noise prediction of propellers dates back to 1919 when Lynam
and Webb first published their work, which was driven by the requirement to have
aircraft flying undetected over enemy territory. Since then, several different empirical,
semi empirical and theoretical methods of propeller noise prediction have been developed
and continuously improved. In recent years, propeller noise has again become of public
interest due to the rapid increase in air traffic and the popularity of wind turbines for
electric power generation. A summary of some of the different methods can be found in
Reference [64]. More recently, the introduction of UAS?s as reconnaissance platforms in
military operations, brought attention to the propulsion and propeller noise, since low
flying reconnaissance UAS?s have low detectability requirements. Engine related noise,
especially for small unmanned aircraft, is avoided by switching to an electric propulsion
system, which leaves the propeller as one of the main sources of noise.
[65]
In this effort, a subsonic steady analysis using point sources (point forces located at the
quarter chord for the elemental forces and point mass sources and sinks to represent
propeller thickness) is implemented. The location of the point sources is illustrated in
Figure 2.32
54
Mass sources
Point forces
Mass sinks
Figure 2.32 Locations of noise sources
The model is based on the work of Lowson
[66]
and was further developed by
Miller
[5]
to account for noise generated by radial forces. It is assumed that the forward
flight velocity is small M<0.1 which seems appropriate because only small low airspeed
aircraft are considered in this work, thus a static analysis technique is applicable. For
optimization efforts only the pressure perturbation, pp
0,
is considered and implemented
into the objective function to reduce propeller noise. If the perceived noise level dBA of a
human ear is of interest, the pressure perturbation must be converted to sound pressure
levels (SPL) and then adjusted based on the frequency to represent human ear sensitivity.
The pressure perturbation after Miller
[5]
is:
55
()
()
() ()()
() ()
() ()()
() ()
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+
?
?
??
?
?
+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+
?
?
??
?
???
+???
?
?
+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+
?
?
??
?
???
+???
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+
?
?
??
?
??
??
?
?
?
?
?
?
?
?
????
=?
r
MM
t
M
aM
r
aq
r
M
t
M
aM
rR
r
a
F
r
M
t
M
aM
r
r
a
F
r
M
t
M
aM
xx
F
rM
pp
srr
r
m
sr
r
obss
obsR
sr
r
obs
obs
sr
r
s
T
r
2
0
0
2
00
2
00
2
0
22
0
1
1
11
1
sinsin
cossin
11
1
sinsin
cossin
11
1
14
1
??
??
??
??
?
?
(2.55)
with
()()()()()[]
2
1
2
2
22
cossin2sincos
ssobsobssobs
RRrrxrr +?????+??= ???? (2.56)
( ) ( )
r
rM
M
obss
r
?? sinsin ???
= (2.57)
0
a
R
M
s
s
?
=
?
(2.58)
??? ?= (2.59)
()
() () () (){}????
??
sinsinsincos
sin
2
3
????+???
???
=
?
?
sobs
obssr
Rrr
r
rM
t
M
(2.60)
where M
s
is the Mach number of a point on the quarter chord in the spanwise direction, r
is the distance between the source and observer, r
obs
is the distance of the observer to the
origin, ? is the angular velocity, ? is the retarded time, R
s
is the distance between the
source and the origin, ? is the angle of the observer with the x axis, and ? is the angle of
the source with respect to the y axis. A geometric illustration is shown in Figure 2.33.
56
The volume displacement of the propeller is simulated by a Rankine
[5]
body with a single
source sink pair on each spanwise collocation position. The source strengths are
determined by solving for the Rankine body in the freestream with the major axis length
equal to the local propeller chord. Thus the maximum diameter, h, of the Rankine body is
determined from:
2
ht
c
t
S ?=??
?
?
?
?
?
?? ? (2.61)
with ?S the spanwise width of the panels and (t/c) the local thickness ratio.
57
Figure 2.33 Acoustic signature analysis geometry
Since the aerodynamic propeller analysis of this work does not consider airfoil thickness,
a generic thickness profile based on thin bladed propellers is applied, which is shown in
Figure 2.34.
.
?
?
x
y
z
x
s
x
y
observer, z=0
source
.
F
T
F
R
F
?
?
yxr ?=
()
() ()?? sin,cos,
0,,
??=
=
sss
RRxy
yxx
R
s
r
obs
58
The point mass flux, q
m
, is obtained by solving the source/sink position from
( ) 01
22
2
2
=+??? ???? (2.62)
where
h
c
=?
2
? (2.63)
and
chordmidfrompositionk
sourcec
____sin2
=?? (2.64)
Once ? is known, q
m
is obtained from
22
2
2
inf0
2
??
?
?
?? +?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
???=
c
Vq
m
(2.65)
To determine the overall pressure perturbation the time of each spanwise pressure
perturbation to reach the observer must be known as the propeller turns. This is done by
initially calculating the distance each blade pressure perturbation has to the observer and
dividing it by the speed of sound. The propeller blade is then rotated backwards in time in
1? increments and the observer time for each spanwise pressure perturbation is now
determined by computing the times the pressure perturbations take to reach the observer
minus the time the blade took to rotate 1?. Then, by matching the times of each pressure
perturbation, the total pressure perturbation is the sum of all the local pressure
perturbations with equal observer times.
59
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
012345678
Propeller Span [inch]
Th
i
c
k
n
e
s
s to Ch
ord
Ra
t
i
o t/
c
[]
Figure 2.34 Linearized airfoil thickness distribution of a APC 14x7 propeller
The test case for verifying the acoustic noise signature is the data obtained by
Succi
[67] [70]
. The simulated propeller is a ? scale model of a Cessna 172 Skyhawk
propeller with the type description 1C160. The geometry of the 1C160 quarter scale
propeller is taken from Succi who describes the chord and propeller pitch in the spanwise
direction. In Appendix G the propeller geometric description is given. Based on this
information, a flat plate propeller was generated with 18 panels spanwise and 3 panels
chordwise which is illustrated in Figure 2.35. In the spanwise direction are placed
sectional forces which are located at the quarter chord position. The point mass sources
and sinks are located with equal distance from the middle chord of each spanwise panel
center.
60
Figure 2.35 1C160 quarter scale Cessna 172 propeller
For the test case verification, the observer was located 0.5 meters from the propeller. The
wave form is for an observer in the plane of rotation behind the propeller ? = 130?. The
propeller speed was 10000 rpm and the free stream velocity 29 [m/s]. In Figure 2.36 the
experimental and theoretical results of Succi as compared to the results of Miller. The
implemented acoustic signature model of his work is shown in Figure 2.37 and agrees
well when compared to the results of Succi and Miller.
61
Figure 2.36 Pressure perturbation comparison experiment and theory from
Miller
[5]
62
Time [ms]
P
r
es
s
u
r
e
p
e
r
t
u
b
a
t
i
o
n
[
P
a]
1.5 1 0.5 0 0.5 1
40
30
20
10
0
10
20
Figure 2.37 Pressure perturbation of 1C160 quarter scale propeller
63
3 PROPELLER AERODYNAMIC PERFORMANCE PREDICTION
PROGRAM
The propeller performance program is the principal portion of the objective
function for the GA. The program consists of nine major routines.
? Airfoil and propeller geometry generation
? Influence coefficient determination
? Gauss Seidel solver subroutine to solve for the unknown circulations
? Main program to connect the subroutines and to compute the lift, drag and
power of the propeller
? Compressibility correction subroutine
? Viscous model analysis
? Pressure perturbation determination
? Flow separation model
The geometry of the propeller is defined by 40 variables:
? 22 variables describe the airfoil geometries
? 5 variables describe the spanwise angle of attack function
? 5 variables describe the spanwise chord length variation function
? 5 describe the propeller sweep function
? 1 for the number of propeller blades
? 1 for the propeller radius
? 1 is related to the position about which the pitch of the propeller is changed
64
The 40 propeller geometry design variables are allowed to change during the
optimization process and other fixed parameters such as the free stream velocity,
rotational speed, number of panels and propeller hub diameter are constants. Figure 3.1
illustrates the program structure. A description of the variables can be found in Section 2
and Appendix B. The open program structure allows for optional activation of certain
individual subroutines to account for compressibility effects, viscous and acoustic
signature determination. In addition to single (design and operating) point analysis the
program is capable of analyzing a propeller over a range of free stream velocities and
propeller rotational speeds. Additionally variable pitch propellers can be designed by
adding variables which offset the angle of attack function of the propeller by a constant
value.
The program starts by loading the propeller geometry input parameters from the
file ?Snglerun.dat.? The subroutine ?coordgen? generates the propeller shape based on
the input parameters. The influence coefficients are determined and the unknown
circulations are returned by a GaussSeidel solver. The accuracy of this basic inviscid
solution can be increased by activating a viscous, compressibility model. Additionally an
acoustic signature model can be used to determine propeller pressure perturbations.
Based on the optimization requirements either thrust, power input, efficiency or acoustic
signatures can be returned and used as an optimization objective.
65
Figure 3.1 Flow chart propeller aerodynamic performance program
Input Parameters:
pr = (x(1)) nxbl = (x(2)) ua_1 = (x(3))
ua_2 = (x(4)) ua_3 = (x(5)) ua_4 = (x(6))
ua_5 = (x(7)) ua_6 = (x(8)) ua_7 = (x(9))
ua_8 = (x(10)) ua_9 = (x(11)) ua_10 = (x(12))
ua_11 = (x(13)) ua_12 = (x(14)) ua_13 = (x(15))
ua_14 = (x(16)) ua_15 = (x(17)) ua_16 = (x(18))
ua_17 = (x(19)) ua_18 = (x(20)) ua_19 = (x(21))
ua_20 = (x(22)) fau_1 = (x(23)) fau_2 = (x(24))
? Compressibility Model
? Viscous Model
? Flow Separation Model
? Noise Model
?Propeller airfoil and Geometry
?Grid Discretization
?Compressibility Domain stretching
Influence Coefficient Determination
Gauss Seidel Solver
Write Output
66
4 CODE VERIFICATION
The aerodynamic propeller performance code was developed without significant use
of existing programs. Thus the numerical scheme was verified using experimental
published and existing numerical data. In an effort to reduce optimization computation
times the Gauss Seidel subroutine which solves for the unknown circulations was
investigated to determine the minimum requirement for the convergence of solution. The
minimum propeller wake extension and the minimum number of chordwise and spanwise
panels which still produce accurate results were also investigated to further reduce
computation times.
4.1 Gauss Seidel Solver convergence
After the determination of the influence coefficients in the x, y and zdirection and
the application of the BC that the flow is parallel to the propeller surface 0? =?nV
r
, a
system of linear dependent equations is obtained.
( )()
()()
()()
?
?
?
?
?
?
?
?
?
?
?
?
?+??
?+??
?+??
=
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
mmmmmm
m
Vr
Vr
Vr
b
bb
bbb
???
???
???
?
cossin
...
cossin
cossin
......
...
4
1
222
111
2
1
2221
11211
(4.1)
67
with
),(),(inf_),(),(inf_),(),(inf_ jinjicoeffjinjicoeffjinjicoeffb
zzyyxxij
?+?+?= (4.2)
To solve for the unknown circulations ?, an implicit or explicit solver can be
used. Due to the large number of panels on the propeller (up to 5000) an explicit Gauss
Seidel solver
[71]
is used to determine the individual circulations of each quadrilateral and
horseshoe element. The advantage of the explicit solver is that it is much faster when
solving a system with a large number of equations.
[71]
The Gauss Seidel solver begins with the first linear equation, solves it for the first
unknown ? as shown in Equation 4.3. The second linear equation is then solved after the
circulation for the second panel and is updated with the vortex strength of the first panel.
This procedure is repeated until all equations are solved once. The process is then
reversed and all linear equations are solved again starting at the bottom of the matrix.
[ ]
[]
[]
mm
m
n
mm
n
m
n
mn
m
n
mm
nn
n
n
mm
nn
n
b
RHSbbb
b
RHSbbb
b
RHSbbb
+??+??+??
=?
+??+??+??
=?
+??+??+??
=?
+
??
++
+
+
+
+
1
11
1
22
1
111
22
22323
1
1211
2
11
113132121
1
...
..........
..........
...
...
(4.3)
The convergence of the solution is measured by the summation of the errors
between the left hand sides (LHS) and the right hand sides (RHS) of the linear equation
after each matrix sweep. Equation 4.4 defines the determination of the related errors.
68
()
?
=
+++
+??+??+??=
m
i
i
n
mim
n
i
n
i
RHSbbberror
1
11
22
1
11
... (4.4)
A two bladed propeller with panel numbers ranging from 32 to 1568 has been
investigated to determine the minimum number of required matrix sweeps which still
produce accurate results. Figure 4.1 shows the convergence of the solution over the
number of matrix sweeps. Noticeable is that after about 250 matrix sweeps the
summation of the errors for all panel densities reach the numerical accuracy of the
computer.
Number of Matrix Sweeps []
S
u
mma
t
i
o
n
o
f
E
r
r
o
r
[
m/
s
]
200 400 600
10
16
10
14
10
12
10
10
10
8
10
6
10
4
4x4
8x8
12 x 12
16 x 16
20 x 20
24 x 24
28 x 28
Figure 4.1 Summation of LHS and RHS difference after matrix sweeps
69
To evaluate the convergence of the solution the thrust of a two bladed propeller
using 512 and 1565 panels respectively is plotted over the number of matrix sweeps.
These results of this study are shown in Figure 4.2 below.
1
10
100
1000
10000
11010
Number of Sweeps []
Thr
ust
[
N
]
0.00001
0.0001
0.001
0.01
0.1
1
10
100
% o
f
th
r
u
s
t
di
ffer
e
nce
w
i
th
r
esp
. to
1000 sw
eeps
Thrust [N] 16 x16
Thrust [N] 28 x 28
% difference to 1000 16 x 16
% difference to 1000 28 x 28
Figure 4.2 Accuracy of the solution after matrix sweeps by Gauss Seidel solver
A single matrix sweep is defined by solving all equations in the matrix once
downstream and once upstream while updating the circulations continuously. The results
shown in Figure 4.2 are referenced against the solution obtained after 1000 matrix
sweeps. It is observed that with increasing panel numbers the number of matrix sweeps
must be increased to maintain a prescribed numerical accuracy. Figure 4.2 shows that at
about 70 sweeps, the solutions are within 1%. Thus for all further propeller performance
70
computations the number of matrix sweeps is set to 100, which will provide an accuracy
of 0.01 in this case while maintaining computational efficiency.
4.2 Wake investigation
To minimize the computation times during the optimization process, a wake
analysis was done to determine the convergence behavior of the aerodynamic
performance parameters. The goal is to reduce the number of downstream rotations of the
propellers rigid helical vortex sheet, while maintaining accuracy of the solution.
The wake analysis was done on optimized two and four bladed propellers with a
diameter of 0.32 meters, a freestream velocity of 21 [m/s] and a propeller speed of 5800
RPM. The wake vortex sheet was extended in steps from 1/8 to 5 rotations downstream
and propeller performance parameters were noted. An illustration of the wakes is shown
in Figures 4.24.4:
71
X
Z
Y
X
Z
Y
Figure 4.24 Wake vortex sheet extensions for ?, 2 and 5 rotations downstream
72
The investigated performance parameters are thrust and power since in most
propeller optimizations these are of special interest. Figure 4.5 shows the thrust generated
form two and four bladed propellers, depending on the wake extension in numbers of
360
o
rotations. Also shown is the percent difference of the propeller thrust with respect to
the max vortex sheet (5 rotations downstream). More details of the wake convergence
investigation are given in Appendix C.
2.5
3.5
4.5
5.5
6.5
7.5
0123456
Number of wake rotations []
Thr
us
t
T [
N
]
0
1
2
3
4
5
6
7
8
9
10
Di
f
f
e
r
e
n
c
e
i
n
Thr
us
t
[
%
]
Thrust 2 bladed propeller
Thrust 4 bladed propeller
% difference in thrust 2 bladed propeller
% difference in thrust 4 bladed propeller
Figure 4.5 Propeller thrust variation due to wake vortex sheet extension
In both, the two and four bladed propeller cases, the thrust converges within 0.3%
after only 2 full downstream rotations of the wake vortex sheet. This result is confirmed
by Cheung
[23]
who showed that, even with a limited wake, good agreements with
73
experimental data is obtained. Thus in all further optimizations the propeller wake is
extended only to two 360? rotations downstream.
4.3 Number of Chordwise and Spanwise Panels
To minimize the optimization computation times, a grid convergence study has been
done. The goal is to reduce the numbers of chordwise and spanwise panels while
maintaining solution accuracy. Two different propeller grid discretizations have been
investigated as shown in Figure 4.6 below. The chosen schemes are based on the work
done by S. Werner
[24]
, Fiddes & Gaydon
[28]
and Labrujere & Zandbergen
[72]
. The
investigated schemes are
1. Full cosine distribution in the chordwise direction and even spacing in the
spanwise direction
2. Full cosine distribution in the chordwise direction and half cosine in spanwise
direction
Figure 4.6 Propeller Discretization schemes
Propeller root
Propeller tip
1 2
74
Figure 4.7 shows that the discretization scheme with a full cosine distribution in the
chordwise direction and the half cosine in the spanwise direction converge the fastest.
The accuracy when compared to the converged thrust value is within 0.5% with a 7x14
(7 panels in chordwise and 14 panels in spanwise direction) panel arrangement. This
agrees with the findings form Cheung
[23]
, who found good agreement with experimental
data with a blade discretization of 7 chordwise and 10 spanwise panels.
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0 25 50 75 100 125 150 175 200
Number of Panels []
Thr
u
st
T
[
N
]
Full cosine chordwise 3x20
Full cosine chordwise 4x20
Full cosine chordwise 5x20
Full cosine chordwise 6x20
Full cosine chordwise 7x20
Full cosine chordwise 7x5
Full cosine chordwise 7x10
Full cosine chordwise 7x25
Full cosine chordwise 7x30
Full cosine chordwise 7x n
half cosine spanwise
Figure 4.7 Lift force convergence by discretization scheme and panel number
75
4.4 Validation of the numerical method
4.4.1 Comparison to lifting line and lifting surface method
The VLM by lifting surface was applied to a simple test case to validate the
numerical method. The validation includes the accuracy of lift coefficient and drag
coefficient for a flat plate wing with zero thickness and an aspect ratio of AR=6 at
different angles off attack. The wing is simulated with a single layer of vortex ring
elements and is discretized in four even spaced panels spanwise and two even spaced
panels chordwise. Figure 4.8 illustrates the wing geometry and the horseshoe trailers.
Even though the lifting line method agrees well with results obtained by the lifting
surface solution when compared to flat plates, the lifting line method is not capable of
predicting more complex geometries accurately.
76
X
Y
Z
Figure 4.8 Wing and wake discretization
The results are compared to the solutions obtained from lifting line method by
Kuethe and Chow
[18]
, which is illustrated in Figure 4.9.
77
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0 2 4 6 8 10 12 14
Angle of Attack [deg]
L
i
ft
Coeffi
c
i
en
t cl
[]
cl from Kuethe & Chow (lifting line method)
cl from Propeller Program (lifting surface method)
0
0.01
0.02
0.03
0.04
0.05
0.06
2 4 6 8 101214
Angle of Attack [deg]
Dr
a
g
Coef
fic
i
e
n
t c
d
[]
cd from Propeller program (lifting surface method)
cd from Kuethe & Chow (lifting line method)
Figure 4.9 Lifting line and lifting surface method cl and cd of flat plate wing AR=6
78
Good agreement between the lifting surface solution and the lifting line method is
obtained for both, the lift and the drag coefficient for the investigated flat plate wing with
AR=6.
4.4.2 NACA 109622 Propeller
In a second validation case the lifting surface method was used to predict the
performance of a NACA109622 straight blade propeller with three blades. The blades are
modeled by ten spanwise and two chordwise panels with a blade angle ?=45.4? at the
75% propeller radius. The geometry of the propeller is given in Appendix E. Figure 4.10
shows the propeller geometry with the horseshoe trailers extending two full propeller
rotations downstream.
79
X
Z
Y
Figure 4.10 NACA109622 propeller geometry with wake
The results are plotted in Figure 4.1113 and incorporate experimental data as well
as data obtained from lifting line method by Cheung
[23]
with and without 2D drag data.
The lifting line method with no drag data and the lifting surface method agree very well.
This is primarily due to the fact that the propeller blade is not cambered, thus the lifting
surface method does not improve the accuracy of the solution.
80
0
0.2
0.4
0.6
0.8
1
1.2
1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5
Advance Ratio J []
E
f
ficie
n
c
y
?
[
]
Experiment
Lifing Line Method no Thickness and cd=0
Lifting Line Method no Thickness with Drag Data
Lifting Surface Method no Thickness and cd=0
Figure 4.11 Propeller efficiency comparison at different advance ratios
Figures 4.12 and 4.13 show the thrust and power coefficients obtained from experimental
data from a lifting line solution with and without drag data included and from the lifting
surface solution developed in this effort. The lifting line and the lifting surface solution
agree very well for both the thrust and the power coefficient. Discrepancies between the
experimental data are mainly due to the simplified model which does not include
thickness or Reynolds number effects.
81
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5
Advance Ratio J []
Th
rust Coeffici
ent c
T
[
]
Experiment
Lifting Line Method no Thickness and cd=0
Lifting Line Method no Thickness with Drag Data
Lifting Surface Method no Thickness and cd=0
Figure 4.12 Propeller thrust coefficient at various advance ratios
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5
Advance Ratio J []
P
o
wer Coefficien
t
c
P
[]
Experiment
Lifting Line Method no Thickness and cd=0
Lifting Line Method no Thickness with Drag Data
Lifting Surface Method no Thickness and cd=0
Figure 4.13 Propeller power coefficients at various advance ratios
82
5 GENETIC ALGORITHM
The optimization of designs in early stages of product development has become
more popular in recent years. The goal is to find an optimal solution for any given design
problem. This effort is driven by cost savings in the development stage and product
competitiveness. Real world applications are found today in every sector of life, which is
mainly due the introduction of very powerful and cheap computers. This lead to the
development of various mathematical optimization theories which range from simple
methods using first derivatives to advanced evolutionary algorithms.
5.1 Optimization methods
Even though optimizations existed before the 17
th
century, with the development of
Newton?s calculus, initially differentiation was used with simple functions to find local
minimums and maximums. Problems emerged when multiple functions were combined
for a system description and nonlinear algebraic equations were produced which were not
easy to solve. In the 19
th
century gradient based methods emerged, which used the
steepest ascent or descent (steps are proportional to the gradient) to find the functions
maximum or minimum. The disadvantages of this method are that the function must be
differentiable in the design space with respect to all independent variables and only local
minimums or maximums can be found, which requires a good initial guess of the solution
83
before starting the optimization. Linear programming as stated by Gabasov and
Kirillova
[73]
was first formulated by a Soviet mathematician L.V. Kantorovich in the
1930?s in which an objective function with linear constraints was used. The disadvantage
of this method is that it cannot handle complex problems with multiple variables since the
objective function and the constraints must be linear. Another method developed around
1960 by Hooke and Jeeves is pattern search optimizations which looks at the behavior of
a function to be analyzed and based on this information generates new points which are
closer to the optimum. The disadvantages of pattern methods are that they may find a
local minimum or maximum depending on the initially guessed starting solution. The
concept of ?Design of Experiments? was initially developed in the 1920s with the
requirement to analyze vast amounts of experimental data. Weber and Skillings
[74]
defined steps for the experimenter to define the structure of the experiment and to find
the proper testing setup. In this method, goals are defined, experiment variables and
ranges are set, the experiment is run and the results are analyzed. A linear statistical
model is then used to describe the structure of the data resulting from the experiment.
Monte Carlo methods are also popular methods since they are able to solve
complex physical and mathematical problems with many variables. They are
computationally expensive since they randomly ?walk? through the design space by
changing variable numbers to find a minimum or maximum. An incorporated gradient
based method can facilitate the determination of the optimum.
A relatively new method is the Particle Swarm Optimization method which was
developed by Kennedy and Eberhart
[75]
in the mid 1990s. It is a population based
stochastic optimization technique inspired by social behavior of bird flocking or fish
84
schooling. PSO shares many similarities with evolutionary computational techniques
such as Genetic Algorithms (GA?s). The system is initialized with a population of
random solutions and searches for optima by updating generations. However, unlike
GA?s, PSO?s have no evolutionary operators such as crossover and mutation. In PSO, the
potential solutions, called particles, fly through the problem space by following the
current optimum particles. Each particle keeps track of its coordinates in the problem
space which are associated with the best solution (fitness) it has achieved so far and
stores the fitness value. Another "best" value that is tracked by the particle swarm
optimizer is the best value, obtained ?to date? by any particle in the neighborhood of the
particle. If a particle takes all the population as its topological neighbors, the best value is
a global best. The particle swarm optimization concept consists of an evaluation of the
performance of each particle at each time step and based on the particles and the global
best solution a change in the velocity of (accelerating) each particle is determined.
Acceleration is weighted by a random term, with separate random numbers being
generated for acceleration toward best (individual particle and global) locations. In the
past several years, PSO has been successfully applied in many research and application
areas Ref [76] [81]. It has been demonstrated that PSO gets results in a faster, cheaper
way as compared with Genetic Algorithms
[82] [84]
. Another reason that PSO is attractive is
that there are few parameters to adjust. One version, with slight variations, works well in
a wide variety of applications. The drawbacks of PSO?s are that they may prematurely
converge and are not able to improve on the solution. The second drawback is that
stochastic approaches have problem dependent performance. This dependency normally
results from the parameter settings of each algorithm.
85
5.2 Genetic Algorithm
Mitchell
[85]
states that ?Genetic Algorithms were invented by John Holland in the
1960?s and were developed by Holland and his colleagues at the University of Michigan
in the 1960?s and 1970?s?. Holland wanted to take the evolutionary processes that are
hypothesized to occur in nature (adaptation, survival of the fittest, etc) and incorporate
them into a computer system. Holland?s classic book ?Adaptation in Nature and Artificial
Systems? formulated evolutionary/population based algorithms that can be used to
optimize a variety of real world systems.
According to Coley
[86]
, Genetic Algorithms are numerical optimization algorithms
built around natural selection and natural genetics. Sometimes GA?s are referred to as
evolutionary optimizers because they mimic some of the process in Darwinian evolution
theory. Darwinian evolution theorizes that totally new species can be produced via
mutation and random chance. However, unlike evolution, GAs operate more on the
principal of improvement of an initial population of solutions to a design problem rather
than pure optimization. Coley
[86]
describes the makeup of a typical GA in a number of
ways. First, a GA can be described as a number, or population, of guesses of the solution
to the problem. Second, The GA employs a method of calculating how good or bad the
individual solutions within the population are. This method is called determining the
fitness of the solution. Additionally, a method for mixing fragments of the better
solutions to form new and on average, even better solutions can be used. This method is
called crossover. Finally, the GA is able to use a mutation operator to avoid permanent
loss of diversity within the solutions.
86
One of the main benefits of using a GA is the fact that the algorithm can start
without a single point, or initial guess, to get the optimization running. The previously
described direct methods all, except the design of experiments and PSO, need an initial
guessed solution to ?march? toward the desired optimized result, which can be a
maximum or a minimum value. The GA uses a population of guesses that are random and
spread throughout the search space. Powerful operators such as selection, crossover and
mutation help direct members of each population toward the desired goal(s) of the
problem. A binary encoding system allows for a host of variables to be manipulated by
the GA and then used in a suite of performance codes. These codes analyze the
performance of each member of the population and the GA ranks each one according to
how well the member of the population meets the desired goal(s). In GA terminology, the
objective function is the function that determines the performance of a particular
chromosome (i.e. member of the population).
[87]
Thus, in this study, the suite of
performance codes, grouped together as a whole, represent the objective function.
At the same time the GA has some disadvantages. In using the GA there is a greater
likelihood that a global optimum solution will be found. However, finding this global
optimum is not guaranteed. Even if the GA is in the neighborhood of the global optimum,
there is a possibility, that through crossover and mutation the global optimum may not be
selected. Also the GA does not address the robustness of the individual design solutions it
creates. The GA simply attempts to meet the desired goals and will adjust the design
parameters accordingly. Thus, it is up to the user to ensure the proper operation of the GA
and to verify that the results are genuine. Finally, the satisfactory operation of the GA
87
relies on the accuracy of the system models that make up the objective function. The
accuracy of the system and overall objective function is described in Chapter 2 and 3.
In this effort a binary encode tournament based GA which was developed by
Anderson
[87]
was used to drive the design optimization. The population size of each
generation was 400 and the optimizations were run for about 700 generations. Results
from initial GA propeller optimizations showed premature convergence behavior. A
second real coded GA, which was recently developed by John
[88]
, was implemented and
showed better convergence to an optimum solution when compared to the binary encoded
GA.
88
6 IMPLEMENTATION OF THE PROPELLER PERFORMANCE PROGRAM
INTO THE GENETIC ALGORITHM
The propeller aerodynamic performance program is integrated in the GA as a
subroutine which, when required, is called by the GA. The GA design variables are listed
in a separate file with min and max value limitations as well as step resolution. This input
file also contains the GA parameters, population size, number of generations and
advanced features such as elitism, niching and pareto.
The GA starts by creating a population of different propellers based on a random
selection of the blade geometry variables. All members are then sent to the propeller
program where the performance parameters lift, drag, power and efficiency are
determined. In the final step the value of the objective parameter or function is send back
to the GA for evaluation. The best performing members of the first generation are
selected and, based on the values of these variables, a new generation of propellers is
created. Figure 6.1 shows the optimization structure.
89
Figure 6.1.Flow chart propeller performance program implementation into GA
The objective function returns a parameter which can be minimized or maximized
by the GA. Based on the requirements of the user, different propeller performance
parameters can be integrated into the objection function. When a combination of
parameters is included in one single objective function, close attention must be given to
how each parameter is weighted in the function. To avoid multiple design variables
Genetic Algorithm (binary)
Generation with 400propeller members
Best performers are chosen and new population is
generated based on
Gannl.dat
Optimization parameters
Population size
Number of generations
Optimization characteristic (niche,
pareto, elitist?)
Propeller Parameters
(max, min, resolution)
Number of propeller blades
Propeller span
Propeller 3D geometry function (22)
AOA function (5)
Sweep function (5)
Chord function (5)
Propeller hinge point
Propeller Performance
3D propeller geometry (with Chord,
AOA and Sweep function)
Grid generation
Results:
Vortex strength
Normal force
Viscous effects
Compressibility
Flow separation model
Propeller acoustic signature
Aerodynamic coefficients, power,
efficiency
Objective function
90
within a function, the GA allows for simultaneous optimization of multiple objective
functions. This simplifies the goal function and leads to more diverse goal options.
91
7 OPTIMIZATION RESULTS
The motivation behind this effort was to investigate the feasibility of designing
propellers using a binary coded GA and a real encoded GA. The goal has been to speed
up the design process while optimizing the propeller for a specific objective composed of
a single or multiple goals. For the purpose of finding an optimal propeller the design
space of available geometries has been developed with a wide range of design
parameters. Throughout all optimizations, the number of chordwise and spanwise panel
elements (7 and 14) has been fixed based on the findings of Section 4.3. In addition, the
number of propeller blades was limited in the program to a maximum of four blades.
The cases are broken down into single and multiple point optimizations as well as
inclusions of compressibility, viscous effects and propeller acoustic signature. The
objective was to design a propeller with better performance than existing propellers. The
propeller operating conditions are defined by the cruise power, the flight speed and the
maximum propeller diameter. These propeller operating conditions were set to match a
small generic UAS like the Raven which is currently used by the US Army and which
was developed by AeroVironment. Therefore the power was set to 200300 Watts in
cruise, the maximum propeller diameter was 0.25 meters and the flight speed 36 [m/s].
The test cases considered are described in the following.
92
7.1 Comparison of propeller optimizations with and without viscous and
compressibility effects taken into account
For comparison of viscous and compressibility effects, three different propeller
optimization cases at different operational points have been investigated. In the first case
a fixed pitch propeller was optimized for a cruise airspeed of 36 [m/s]; in the second case
a fixed pitch propeller is optimized for a launch condition of 9 [m/s]; and in the third case
a fixed pitch propeller is optimized using both cruise and launch conditions as the
objective. The parameters of interest are the propeller thrust (the first goal to be matched)
and power (second goal to be minimized to maximize propeller efficiency).
7.1.1 Single point Optimization for Cruise Condition: Binary encoded GA
Three different single point optimization runs have been conducted to determine
the influence of viscous and compressibility effects. The first case did not include
compressibility or viscous effects, the second one considered only compressibility, and
the last took only viscous effects into account. All of the cases had the same fixed input
parameters as shown in Table 7.1.
Table 7.1 Fixed input parameters cruise condition
Free stream velocity 35 [m/s]
Propeller speed 7500 [rpm]
Propeller max diameter 0.25 [m]
Propeller hub diameter 0.06 [m]
93
The binary encoded GA was set up to drive the design using two objective functions
which must be minimized. Within the GA the two objective functions are combined to
one single function which is stored in the ?GAout1? file and used to determine
convergence behavior. The first objective function contains the required thrust in cruise,
which was set in this design effort to 7.0 [N]. The cruise thrust selected was a guess,
based on available data for small, UAS and does not represent a particular system.
Subtracted from the required thrust is the generated thrust determined by the aerodynamic
prediction program. Equation 7.1 shows the function which is minimized when the
propeller thrust reaches the required thrust.
Thrustobj ?= 0.71_ (7.1)
In the second objective function the power was set as a goal to be minimized which
drives the propeller efficiency. The power was divided by 300 to equalize the weighting
of the objective functions.
300
2_
Power
obj = (7.2)
Table 7.2 shows the results of the three optimization cases which have been run for about
200 generations with each generation consisting of 400 members. In all of the cases the
thrust required was matched. The propeller efficiencies range between 78 and 87%
depending if viscosity or compressibility are taken into account. A variable pitch
propeller with similar diameter was tested in the wind tunnel at 18 [m/s] free stream
94
velocity and was rated with a max efficiency of 62 %. Wind tunnel data are given in
Appendix H.
Table 7.2 Optimized propellers performance parameter results cruise condition
Case 1
Inviscid
Incompressible
Case 2
Inviscid
Compressible
Case 3
Viscous
Incompressible
Thrust [N] 7.0 7.0 7.0
Power [W] 306.9 290.2 321.3
Efficiency [%] 82.0 86.8 78.4
In Figures 7.1, 7.2 and 7.3 the optimized propellers for inviscid, compressible and
viscous solutions are illustrated. All blades are two bladed with the inviscid and the
viscous solutions showing some form of a proplet. Noticeable in Figures 7.17.3 are also
the leading edge panels of the propeller blades are bent downwards in the tip region
which would explain why the efficiencies in some cases are lower.
95
Figure 7.1 Optimized cruise propeller geometry from binary encoded GA; inviscid,
incompressible
96
Figure 7.2 Optimized cruise propeller geometry from binary encoded GA; inviscid,
compressible
97
Figure 7.3 Optimized cruise propeller geometry from binary encoded GA; viscous,
incompressible
98
Nondimensional Propeller Span []
Ci
r
c
u
l
a
t
i
o
n
[
m
2
/s
]
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
inviscid, incompressible
inviscid, compressible
viscous, incompressible
Figure 7.4 Optimized propeller spanwise circulation distributions for cruise
The spanwise circulations for the propellers considered in this section are presented in
Figure 7.4. All the cases have discontinuities in the spanwise circulation, which reduces
propeller efficiencies and is a sign that only a local optimum has likely been found. The
viscous incompressible case in Figure 7.3 developed a noticeable large proplet at the tip,
which allows for a higher tip loading which is shown in Figure 7.4 by an increased
circulation.
99
Number of Generations []
Number of Evaluations []
O
b
je
c
t
i
v
e
F
u
n
c
t
io
n
[

]
0 50 100 150 200 250
0 50000 100000
1
1.5
2
inviscid, incompressible
inviscid, compressible
viscous, incompressible
Figure 7.5 Objective function convergences binary encoded propeller
Since the three optimizations have been run for only 200 generations a separate run has
been done to investigate the long term convergence behavior. Therefore the inviscid
incompressible case has bee run for 1300 generations which relates to 520,000 propeller
evaluations and took about 6 days of computation. The convergence result is illustrated in
Figure 7.6 and demonstrates that after about 100 generations no significant improvements
are achieved even though a gain of about 34 % in propeller efficiency is possible.
100
Number of Generations
O
b
j
e
c
t
i
v
e
F
unc
t
i
on
0 500 1000
1
1.2
1.4
1.6
1.8
2
2.2
inviscid, incompressible
Figure 7.6 Objective function convergence inviscid, incompressible, for binary
encoded GA
7.1.2 Single point Optimization for Cruise Condition: Real encoded GA
The three different cases presented in section 7.1, the inviscid, incompressible, the
inviscid, compressible and the viscous incompressible were run again using a real
encoded GA to investigate convergence behavior to determine whether on not the
objective functions are minimized to an lower value. All input parameters and objective
functions are identical to the cases shown in Section 7.1.1.
In table 7.3 are illustrated the results obtained from the optimizations of the real coded
GA. The thrust objective in all cases is satisfied exactly with the inviscid incompressible
101
case having the highest efficiencies. The compressible case lags the efficiency by about 7
% when compared to the results obtained be the binary encoded GA, which indicates a
premature convergence thus a global optimum is not reached.
Table 7.3 Optimized propellers performance parameter results real encoded GA
Case 1
Inviscid
Incompressible
Case 2
Inviscid
Compressible
Case 3
Viscous
Incompressible
Thrust [N] 7.0 7.0 7.03
Power [W] 288.9 320.5 320.4
Efficiency [%] 87.2 79.2 79.0
The propeller geometries generated by the real coded GA are shown in the Figures 7.7 ?
7.9. It is again characteristic that two propeller geometries have proplets which would
indicate a loading shift to the tip. Though the circulation distributions in Figure 7.10
shows that the loading is actually more spread out when compared with the results
obtained by the binary encoded GA. The lag in efficiency of the viscous incompressible
case is also seen in the discontinuous circulation distribution illustrated in Figure 7.10.
102
Figure 7.7 Optimized cruise propeller geometry from real encoded GA; inviscid,
incompressible
103
Figure 7.8 Optimized cruise propeller geometry from real encoded GA; inviscid,
compressible
104
Figure 7.9 Optimized cruise propeller geometry from real encoded GA; viscous,
incompressible
105
Nondimensional Propeller Span []
Ci
r
c
u
l
a
t
i
o
n
[
m
2
/s
]
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0.2
0.4
0.6
0.8
1
1.2
inviscid, incompressible
inviscid, compressible
viscous, incompressible
Figure 7.10 Optimized circulation distributions for cruise propeller, real encoded
GA
When comparing Figure 7.5 and 7.10 it is noticeable that even though both GA versions
seem to level off in the convergence behavior around 100 generations, the real encoded
GA shows a higher gradient in the objective functions after the first 100 generations.
Thus for all the following optimizations the real coded GA is used. The case where
viscous effects are taken into account is the only one used in the remaining optimizations
since the subroutine for compressibility effects should only be used in determining the
performance of a propeller when compared to an inviscid or viscous solution.
106
Number of Generations []
Number of Evaluations []
O
b
j
e
c
t
i
v
e
F
unc
t
i
on
[

]
0 50 100 150 200 250
0 50000 100000
0.8
1
1.2
1.4
1.6
1.8
2
inviscid, incompressible
inviscid, compressible
viscous, incompressible
Figure 7.11 Objective function convergences real encoded GA, cruise condition
To verify that the optimized propeller has better performance a CAM 13x7 propeller, set
up with a variable pitch mechanism, was tested in a wind tunnel at about 19 [m/s] free
stream velocity. The power input was kept constant while the propeller pitch was
changed gradually. Figure 7.12 shows that the maximum propeller efficiency of the CAM
13x7, at 19 [m/s], is 63% while the optimized propeller operates at about 80% efficiency.
107
0
10
20
30
40
50
60
70
0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
Propeller Torque [Nm]
Pro
p
e
l
l
e
r
Ef
f
i
ci
en
cy
[
%
]
CAM 13x7 Propeller
Figure 7.12 CAM 13x7 propeller efficiency at 19 [m/s] free stream velocity
7.2 Single point Optimization Launch Condition
A single point optimization was conducted to design a propeller for a launch condition of
9 [m/s]. The investigation included the case where viscous effects are taken into account
as described in section 7.1.1. The only difference in the fixed input parameters is the
launch airspeed which simulates a hand launch of a UAS by a human and the thrust
requirement of 14 [N]. The input parameters for this case are shown in Table 7.3. The
thrust required for the launch was set to a high value since the operational mode of this
108
system does not include a mechanical launch by catapult or bungee cord. Thus the
optimization objective functions are
Thrustobj ?= 0.141_ (7.3)
300
2_
Power
obj = (7.4)
Table 7.3 Fixed input parameters for launch condition
Free stream velocity 9 [m/s]
Propeller speed 7500 [rpm]
Propeller max diameter 0.25 [m]
Propeller hub diameter 0.06 [m]
The optimized propeller satisfies the thrust requirement by generating exactly 14 [N] of
thrust. The efficiency of 65 percent seems initially low, but this is due to the
characteristic of the efficiency equation which only produces high efficiencies at higher
free stream velocities when propeller size is limited. Results from wind tunnel tests of a
CAM 13x7 propeller reached an efficiency of 58 %. The data of the wind tunnel test are
shown in Appendix H.
Table 7.4 Optimized propellers performance parameter results for launch condition
Case:
Viscous, Compressible
Thrust [N] 14.0
Power [W] 193.6
Efficiency [%] 65.1
109
Figure 7.13 Optimized launch propeller geometry: viscous, incompressible case
The optimized propeller has an efficiency of 65 % which is 7% better than the propeller
tested during wind tunnel experiments. The propeller which has four blades due to the
thrust requirement of 14 [N] is shown in Figure 7.13.
110
Nondimensioal Propeller Span []
Ci
r
c
u
l
a
t
i
o
n
[
m
2
/s
]
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0.2
0.4
0.6
0.8
1
1.2
viscous, incompressible
Figure 7.14 Optimized propeller spanwise circulation distribution viscous,
incompressible case
The spanwise circulation of the launch propeller shown in Figure 7.14 is smooth in the tip
region, but shows breaks toward the propeller root. The reason is due to either premature
convergence or convergence that has not yet been reached.
111
Number of Evaluations []
O
b
je
c
t
iv
e
F
u
n
c
t
io
n
[

]
0 10000 20000 30000 40000
0.5
1
1.5
2
2.5
viscous, incompressible
Figure 7.15 Objective function convergence from real encoded GA: viscous,
incompressible case
The convergence of the objective function in Figure 7.15 is obtained after about 32,000
evaluations.
7.3 Dual point Optimization Cruise and Launch Condition
In this case a propeller is optimized for both, the launch and the cruise conditions. The
only parameter which changes is the free stream velocity which, for the launch condition
is equal to the single point optimization of section 7.1.2 and, for the cruise condition is 36
[m/s] as described section 7.1.1. The first objective function (Equation 7.5) includes the
112
absolute value of required thrust for launch minus the generated propeller thrust and the
second objective function (Equation 7.6) contains the cruise thrust requirement as well as
the actual cruise power. The weighting of the propeller thrust and power in the second
objective function (Equation 7.6) was normalized by dividing the power by 300. The first
objective function did not contain any power term since it is assumed that the launch
power setting is only a very short time period thus it is negligible in the overall flight
endurance.
launch
Thrustobj ?= 0.141_ (7.5)
300
0.72_
cruise
cruise
Power
Thrustobj +?= (7.6)
The fixed input parameters shown in Table 7.5 are the same as from section 7.1.1 and
7.1.2.
Table 7.5 Fixed input parameters for dual point optimization
Free stream velocity 9 / 36 [m/s]
Propeller speed 7500 [rpm]
Propeller max diameter 0.25 [m]
Propeller hub diameter 0.06 [m]
The thrust required in both the launch and cruise condition is well matched. When
compared to the single point optimizations the efficiencies in the dual point optimizations
are lower which is expected since the propeller has to satisfy both operating conditions.
Noticeable is that the geometry of the optimized two bladed propeller appears to be very
113
smooth even in the tip region where the single point optimized propellers exhibits small
ripples in the geometry.
Table 7.6 Optimized propeller performance parameters: cruise and launch
Case:
Viscous, Incompressible
Launch Cruise
Thrust [N] 14.01 6.97
Power [W] 373.9 396.6
Efficiency [%] 33.7 63.7
114
Figure 7.16 Optimized cruise and launch propeller; viscous, incompressible case
The circulation distributions for both the launch and cruise conditions shown in Figure
7.17 have smooth distributions, with the launch case showing a higher circulation due to
the larger thrust requirement.
115
Nondimensional Propeller Span []
Ci
r
c
u
l
a
t
i
o
n
[
m
2
/s
]
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
launch condition v=9 [m/s]
cruise condition v=36 [m/s]
Figure 7.17 Objective function convergences real encoded GA, cruise condition
The convergence behavior of the goal for the two advance ratio optimization is shown in
Figure 7.18. The objective function convergence appears to occur at about 60,000 case
evaluations which are about 20,000 evaluations higher than the single point optimizations
done before.
116
Number of Evaluations []
O
b
j
e
c
t
i
v
e
F
unc
ti
on
[

]
0 50000 100000
0.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
3
Figure 7.18 Objective function convergences real encoded GA, cruise condition
7.4 Single Point Optimization in Cruise considering Propeller Noise
The first case was set to optimize the cruise propeller as described in 7.1.1 with the
same fixed input parameters. The objective functions contained the total pressure
perturbation (divided by 5 to obtain about equal weighting) and the thrust required. Thus
the two objective functions are
Thrustobj ?= 0.71_ (7.7)
117
( )
0.5
2_
0 total
pp
obj
?
= (7.8)
The propeller power was not part of the optimization goal since the propeller efficiency
was not of interest in this optimization.
The optimized propeller satisfies the thrust requirement by producing exactly the
required thrust as shown in Table 7.7. Figure 7.19 illustrates the loading shift of the
circulation to the inboard section which relates to a noise reduction. Also the span of the
propeller is reduced, which results in lower noise due to lower tip Mach numbers. Finally
the optimization lead to an increase from two to four blades, shown in Figure 7.20, which
is known to lower the overall propeller noise. All of these findings follow the conclusions
made by Miller
[5]
who investigated propeller noise optimization. The optimized propeller
efficiency, shown in Table 7.7, is low since the power minimization was not an
optimization objective.
Table 7.7 Optimized propellers performance parameters cruise condition with noise
considered
Case:
Viscous, Incompressible
Thrust [N] 7.0
Power [W] 509.4
Efficiency [%] 49.5
118
Nondimensioal Propeller Span []
Ci
r
c
u
l
at
i
o
n
[
m
2
/s
]
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
Noise included in objective function
Noise not included in the objective function
Figure 7.19 Spanwise circulation distribution of the optimized propeller in cruise
condition
119
Figure 7.20 Optimized cruise propeller from real coded GA; viscous, incompressible
with noise considered
7.5 Single/Dual Point Optimization in Cruise/Launch Condition considering
Propeller Noise and Performance
The optimizations of three different propeller operational points are considered in this
section. A cruise propeller at 36 [m/s] free stream velocity as described in section 7.1.1, a
launch propeller at 9 [m/s] free stream velocity from section 7.1.2 and the combined case
with cruise and launch condition described in section 7.1.3. For this analysis the propeller
acoustic signature is included in the objective functions to reduce the overall propeller
noise. The acoustic computation is done with the model described in section 2.8 without
120
the thickness noise taken into account. This is done because the thickness distribution
cannot be changed during the optimization and acts only as a fixed source of noise.
Additionally the thickness noise model only applies to propellers with a high aspect ratio
and does not predict the thickness noise well for propellers with large chords. Thus the
objective functions for the three different propeller operational points are as follows:
For cruise condition:
Thrustobj ?= 0.71_ (7.9)
( )
0.3000.10
2_
_0
Power
pp
obj
cruisetotal
+
?
= (7.10)
For launch condition:
Thrustobj ?= 0.141_ (7.11)
( )
0.3000.10
2_
_0
Power
pp
obj
launchtotal
+
?
= (7.12)
For launch / cruise condition:
launchcruise
ThrustThrustobj ?+?= 0.140.71_ (7.13)
( ) ( )
0.100.3000.10
2_
_0_0 cruisetotal
cruise
launchtotal
pp
Power
pp
obj
?
++
?
= (7.14)
121
The results from the three different optimizations are included in Table 7.8. The thrust
requirements are all matched exactly and the propeller efficiencies seemed close to the
ones where noise was not considered.
Table 7.8 Propellers performance parameter results with noise optimization
Case 1
Cruise
36 [m/s]
Case 2
Launch
9 [m/s]
Case 3
Launch / Cruise
9 / 36 [m/s]
Thrust [N] 7.0 14.0 13.7 /6.62
Power [W] 322.9 267.8 366.8 / 350.9
Efficiency [%] 78.0 47.0 30.0 / 70.2
The geometries of the optimized propeller blades, illustrated in Figures 7.217.23 share
similar shapes. The root chord is small and increases to the propeller mid section with a
decrease in chord length towards the tip. The propeller blade geometries have a smooth
shape with no geometric perturbation visible.
122
Figure 7.21 Optimized cruise propeller; viscous, incompressible case with noise and
performance considered
123
Figure 7.22 Optimized launch propeller; viscous, incompressible case with noise and
performance considered
124
Figure 7.23 Optimized launch / cruise propeller; viscous, incompressible case with
noise and performance considered
125
Nondimensional Propeller Span []
Ci
r
c
u
l
a
t
i
o
n
[
m
2
/s
]
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0.2
0.4
0.6
0.8
1
1.2
1.4
36 [m/s] Free Stream Velocity, Noise considered
36 [m/s] Free Stream Velocity, No noise considered
Figure 7.24 Spanwise circulation distribution of the optimized propeller in cruise
The optimization of the cruise propeller produced a significant reduction in the
circulation strength and a shift towards the propeller mid section, which is illustrated in
Figure 7.24. The noise optimized propeller shown in Figure 7.24 is normalized with
respect to the propeller span of the optimized propeller where no noise was considered.
Besides the reduction in circulation the noise optimized propeller shows also a smaller
span which lowers the tip Mach number and as a result reduces the acoustic signature.
[5]
Even though the circulation distribution in the case of the noise optimized propeller is
more uniform, the propeller efficiency is a little lower. The reasons for that are the
126
smaller propeller diameter which lowers the efficiency and the wider blade chord which
generates higher viscous drag.
Nondimensional Propeller Span []
Ci
r
c
u
l
a
t
i
o
n
[
m
2
/s
]
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0.2
0.4
0.6
0.8
1
9 [m/s] Free stream Velocity, Noise considered
9 [m/s] Free Stream Velocity, No noise consiedered
Figure 7.25 Spanwise circulation distribution of the optimized propeller at launch
The circulation distribution of the launch condition propeller does not show any
significant shift or reduction in the circulation. This may be explained by the significant
objective function improvement toward the end of the optimization, where the design was
driven by improving the efficiency and not by the acoustic propeller signature.
127
Nondimensional Propeller Span []
Ci
r
c
u
l
a
t
i
o
n
[
m
2
/s
]
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0.5
1
1.5
2
36 [m/s] Free Stream Velocity, Noise considered
9 [m/s] Free Stream Velocity, Noise considered
36 [m/s] Free Stream Velocity, No noise considered
9 [m/s] Free Stream Velocity, No noise considered
Figure 7.26 Spanwise circulation distribution of the optimized propeller in cruise
and launch condition
The last case considered both acoustic signature at launch and the acoustic signature at
cruise condition. The noise optimized propeller shows a reduction in the circulation
strength at both operating points as shown in Figure 7.26. The loading seemed to shift
towards the tip region which is an indication that power and thrust required in the
objective function are dominant.
128
Number of Generations []
Number of Evaluations []
O
b
j
e
c
t
i
v
e
F
unc
ti
on
[

]
0 50 100 150 200 250
0 50000 100000
1
1.5
2
2.5
Free Stream Velocity 36 [m/s]
Free Stream Velocity 9 [m/s]
Free Stream Velocity 9 / 36 [m/s]
Figure 7.27 Objective function convergences, cruise, launch and combined case
The convergence of the objective functions of the three investigated propeller operating
points is illustrated in Figure 7.27. The single point optimizations for launch and the
cruise conditions show a convergence after about 30,000 evaluations, but the case
considering the launch and the cruise condition together shows after an initial
convergence tendency a steep gradient in the objective function around 90,000
evaluations which indicates that an optimum has not yet been achieved.
129
8 CONCLUSION & RECOMMENDATIONS
The goal of this research was to develop a tool for performance and acoustic noise
signature optimization applicable to small thinairfoil propellers with complex
geometries. The optimization scheme consists of a propeller performance program, an
acoustic signature subroutine and a GA for the optimization. A major thrust of this work
was the development of a general geometric description of propeller blades which allows
for a wide design space while generating smooth, viable blade geometries. The VLM
using a lifting surface scheme with vortex ring elements was used to determine the
aerodynamic propeller performance parameters. This allows for the computation of swept
propellers including dihedral and camber variation. All propeller blades were considered
to be thin and modeled using a single layer of vortex elements to describe the geometry.
The research includes an investigation to reduce computational times by finding the
minimum number of required panel elements and the minimal wake extension which
produces accurate propeller performance results. The performance analysis scheme was
verified by comparing it with flat plate wing results as well as a straightblade
NACA109622 propeller. The implemented acoustic signature model to determine
propeller noise is based on the work of Miller,
[5]
and during the optimization, only
loading noise was considered.
130
The optimizations were performed to generate a propeller for a small electric powered
UAS. Fixed input parameters such as the maximum chord, the maximum propeller span
and the propeller speed were kept constant to allow for comparisons of the optimizations.
Three different propeller operational conditions were considered. One launch condition
with a free stream velocity of 9 [m/s], one cruise condition with 36 [m/s] free stream
velocity and one case were launch and cruise condition were considered as the same time
to demonstrate multiobjective optimization. The analysis included inviscid
incompressible, compressible and viscous flows. Additionally, a binary and a real
encoded GA have been used to investigate the method that converges quicker and which
method produces better performing propellers. Results showed that even though both
GA?s converged within the same number of case evaluations, the real encoded GA was
generally more successful over longer optimization runs. Some optimizations show
premature convergence which can be improved upon by changing ranges and limitations
on GA variables. The method produced tradeoffs between best performing single run
and dualpoint optimizations. Some of the propellers showed the development of a
proplet with a loading shift to the blade tip section. Since the optimization did not include
a propeller structure model which would eliminate propellers with unrealistic high
bending moments, propellers with proplets are competitive candidates. Noticeable is that
the circulation distributions for the single point optimizations showed unevenness
whereas the case with multiple objective goals generated a smoother circulation
distribution. The inclusion of noise in the objective function resulted in a reduction of the
propeller radius for some cases and a loading shift of the circulation towards the propeller
mid section. When multiple parameters were used in a single objective function the
131
weighting of the individual goals within the objective functions must be carried out
judiciously to achieve a propeller design which satisfies all goals. Further research on the
effects of weighting of individual propeller parameters within an objective function
should be investigated.
Recommendations
? The recommendations for improvement of this design tool are concerned mainly
with premature convergence during the optimizations. Further investigations must
be completed to fully understand why in some cases, the GA converges to a local
and not a global minimum.
? The implemented empirical flow separation and viscous model could be replaced
by a viscous/inviscid interaction method to describe flow separation and viscous
effects better.
? To improve the noise prediction for blades with large chords, multiple chordwise
noise sources are required.
? To allow for the simulation of thick airfoil propellers, a second vortex layer must
be added to describe the bottom airfoil geometry.
? The implementation of a structure model based on the two vortex layers to
describe the propeller geometry would generate a realistic thickness distribution.
? The extension of the aerodynamic propeller performance program to include a
duct around the propeller would allow for the optimization of ducted propeller/fan
combinations.
? Finally, for highly loaded and low advance ratio conditions, a free wake model
can improve the aerodynamic prediction.
132
REFERENCES
[1] Lynam, F. J. H. and Webb, H. A. ?the Emission of Sound by Airscrews.?
Technical Report of the Advisory Committee if Aeronautics for the Year
19181919 Vol. 2 (London) Technical Reports and Memoranda No. 624,
March 1919, pp. 792801.
[2] J. Grasmeyer, ?Application of a Genetic Algorithm with Adaptive Penalty
Functions to Airfoil Design?, AIAA 19977, Aerospace Science Meeting
Exhibit, 35
th
, Reno, NV, Jan. 69, 1997
[3] Fanjoy, D. W., Crossley W. A., ?Aerodynamic Shape Design for Rotor Airfoils
via Genetic Algorithm?. American Helicopter Society 53
rd
Annual Forum,
Virginia Beach, VA, April 19may 1, 1998.
[4] Gardner A. B., Selig M. S., ?Airfoil Design using a Genetic Algorithm and an
Inverse Method?. AIAA 200343, 41
st
Aerospace Science Meeting and
Exhibit, 69 January, 2003.
[5] Miller C. J., ?Optimally Designed Propellers Constrained by Noise?, Purdue
University, Dissertation, December 1984.
[6] Powell, M. J. D., ?An Efficient Method for Finding the Minimum of a Function of
Several Variables without Calculating Derivatives,? Computer Journal,
Volume 7, 1964, pp. 155162.
[7] Anders Smaerup Olsen, ?Optimization of Propellers using the Vortex Lattice
Method?, Technical University of Denmark, Dec 2001.
133
[8] Hampsey M., ?Multiobjective Evolutionary Optimization of Small Wind Turbine
Blades.? Dissertation, Department of mechanical Engineering, University of
Newcastle, 2002.
[9] Rankine, W. J. M., ?On the Mechanical Principles of the Action of propellers?,
Transactions of the Institute of Naval Architects, Volume 6, 1865.
[10] Froude R. E., ?On the Part Played in Propulsion by Difference in Fluid Pressure?,
Transactions of the Institute of Naval Architects, Volume 30, 1889.
[11] Vogeley A. W., ?Axial Momentum Theory for Propellers in Compressible Flow?,
NACA TN No. 2164, August, 1950.
[12] Froude W., ?On the elementary Relation between Pitch, Slip and Propulsive
Efficiency?, Transaction of the Institute of Naval Architects, Volume 19,
1878.
[13] Drzewiecki S., ?Theorie Generale de l?H?eleice?, Paris, 1920.
[14] Anderson J.D., Jr. ?Fundamentals of Aerodynamics?, 3rd Edition, McGrawHill,
New York 2001.
[15] Glauert H., ?An Aerodynamic Theory of the Airscrew?, Grear Britain
Aeronautical Research Committee, Reports & Memoranda No. 786, 1922.
[16] Katz J., Plotkin A., ?Low Speed Aerodynamics? 2nd Edition, Cambridge
University Press, New York, 2005.
[17] Prandtl L., ?Tragfluegel Theorie? 1 und 2 Mitteilung. Nachrichten von der
Koeniglichen Gesells. D. Wissens. Math.Phys. Klasse 451 (1918) und 107
(1919). Reprinted in Abhandlungen zur Hydrodynamik und Aerodynamik.
Goettingen 9(1927).
[18] Kuethe A. M., Chow C. Y., ? Foundations of Aerodynamics, Bases of
Aerodynamic Design? 5th Edition, John Wiley & Sons, Inc., New York, 1998
134
[19] Purvis J. W. and Burkhalter J. E., ?Simplified Solution of the Compressible
Lifting Surface Problem?, AIAA Journal, Volume 20, No. 5, Page 589, May
1982.
[20] Burkhalter J. E., Harris M., ?Nonlinear aerodynamic analysis of grid fin
configurations?, AIAA Applied Aerodynamics Conference, AIAA1995
1894, 13th, San Diego, CA, June 1922, 1995.
[21] Brooks R. A., Burkhalter J. E., ?Experimental and analytical analysis of grid fin
configurations?, Journal of Aircraft, 00218669, vol.26, No.9 (885887),
1989.
[22] Burkhalter J. E. (U.S. Air Force Academy, Colorado Springs, CO), ?Downwash
measurements on a pitching canardwing configuration?, Journal of Aircraft,
00218669, vol.30, No.6 (10051008), 1993.
[23] Shiu Wing Cheung, R. ?Numerical Prediction of Propeller Performance by Vortex
Lattice Method?, UTIAS Technical Note No.265 Nov 1987.
[24] Sofia Werner, ?Application of the Vortex lattice Method to Yacht Sails?,
University of Aukland, Master Thesis, July 2001.
[25] Goldstein, S., ?On the Vortex Theory of Screw Propellers?, Proceedings of the
Royal Society of London, Series A, Vol. 123, 1929.
[26] Christopher J. Fisichella, ?A case for improving prescribed wake analysis?, AIAA
ASME Wind Energy Symposium, 19th, Aerospace Science Meeting and
Exhibit, 38th NV, Jan 1013 2000.
[27] A. J. Chorin and P.S. Bernard, ?Discretization of a vortex sheet, with an example
of rollup, J. Comp. phys, (13):423, 1973.
[28] S.P. Fiddes, and J. Gaydon, ?A vortex lattice method for calculating the flow past
yacht sails?, Journal of Wind Engineering and Industrial Aerodynamics,
63:3559, 1996.
135
[29] Karamcheti K. ?Principles of IdealFluid Aerodynamics?, Rep. Edition, John
Wiley & Sons, Inc., New York, 1980.
[30] D. F. Thrasher, D. T. Mook, and A. H. Nayfeh, ?A computer based method for
analyzing the flow over sails?, Chesapeake Sailing Yacht Symposium,
pages119127, SNAME, 1983.
[31] J. J. Bertin and M. L. Smith, ?Aerodynamics for Engineers?, PrenticeHall, Inc.,
Englewood Cliffs, 2nd Edition, 1989.
[32] A. H. Day, ?Steps toward an optimal yacht sail?, Technical Report, The Royal
Institution of Naval Architecture, 1993.
[33] J. L. Guermond, ?Collocation methods and lifting surfaces?, European Journal of
mechanics, B/Fluids, 8(4):283305, 1989.
[34] Sobieczky, H., ?Parametric Airfoils and Wings?, Notes on Numerical Fluid
Mechanics, Vol. 68, pp.7188, Vieweg Verlag, 1988.
[35] Samareh, J.A., ?Survey of Shape Parameterization Techniques for HighFidelity
Multidisciplinary Shape Optimization?, AIAA Journal Vol. 39, No.5, May
2001.
[36] Robinson, G. M., Keane, A. J., ?Concise Orthogonal Representation of
Supercritical Airfoils?, Journal of Aircraft, Vol. 38, No.3.
[37] Song, W., Keane, A. J., ?A Study of Shape Parameterization Airfoil
Optimization?, AIAA2004448 10th AIAA, ISSMO Multidisciplinary
Analysis and Optimization Conference, Albany, New York, Aug. 3001,
2004.
[38] Padula, S., Li, W., ?Options for Robust Airfoil Optimization Under Uncertainty?,
9th AIAA Multidisciplinary Analysis and Optimization Symposium, 46
September2002.
136
[39] Hicks, R. M., Henne, P.A.,? Wing design by numerical optimization?, Journal of
Aircraft, Vol.15, pp. 407412, 1978.
[40] Kulfan B., M., Bessoletti J., E., ?Fundamental Parametric Geometry
Representation for Aircraft Component Shapes?, 11th AIAA/ISSMO
Multidisciplinary Analysis and Optimization Conference, 66 September
2006, Portsmouth, Virginia.
[41] Kulfan B., M., ?A Universal Parametric Geometry Representation Method ?
CST?, 45th AIAA Aerospace Sciences Meeting and Exhibit, 811 January,
Reno, Nevada, 2007.
[42] Gennaretti M. and Morino L., ?A Boundary Element Method for the Potential,
Compressible Aerodynamics of Bodies in Arbitrary Motion?, Aeronautical
Journal, Vol. 96, Jan. 1992, pp1519.
[43] Curri I. G., ?Fundamental Mechanics of Fluids?, McGrawHill, 1993.
[44] Fung Y. C., ?An Introduction to the Theory of Aeroelasticity?, John Wiley &
Sons, Inc., New York, NY, 1955.
[45] C. J. Szymendera, ?Computational Free Wake Analysis of a Helicopter Rotor?,
Master Thesis, Pennsylvania State University, 2002.
[46] Tromp, Jeffrey C., ?Numerical Study of Three Viscous/Inviscid Interaction
Methods?, ADA202575, Defense Technical Infromation Center, Sep. 1988.
[47] Cebeci, T. and Besnard, E. ?An efficient and accurate approach for analysis and
design of high lift configurations.? Canadian Aeronautics and Space Journal,
44(4):117 1998.
[48] John Dreese, http://www.dreesecode.com/primer/airfoil5.html, 2004
137
[49] Borst, H. V. and Associates, ?Summary of Propeller Design.? Summary of
Propeller Design Procedures and Data, Volume 1 Aerodynamic Design and
Installation, AD774 831, NTIS, 1973.
[50] Leishman, J. G., ?Challenges in Modeling the Unsteady Aerodynamics of Wind
Turbines,? Wind Energy, Vol. 5, 2002, pp. 86?132.
[51] Snel, H., ?Review of the Wind Turbine Aerodynamics,? Wind Energy, Vol. 6, No.
3, July/September, 2003, pp. 203?211.
[52] Hansen A. C. and Butterfield C. P., ?Aerodynamics of HorizontalAxis Wind
Turbines,? Annual Review of Fluid Mechanics Vol. 25, 1993, pp. 115?149.
[53] Fingersh, L. J., Simms, D., Hand, M., Jager, D., Cotrell, J., Robinson, M.,
Schreck, S. and Larwood, S. ?Wind Tunnel Testing of NREL?s Unsteady
Aerodynamics Experiment,? In proceedings of 39th Aerospace Sciences
Meeting, AIAA Paper 2001?0035, Reno, NV, Jan. 8?11, 2001.
[54] Corrigan, J. J., and Schilings, J. J., ?Empirical Model for Stall Delay Due to
Rotation,? American Helicopter Society Aeromechanics Specialists
Conference, San Francisco, CA, 1994.
[55] Du, Z. H. and Selig, M. S., ?A 3D Stall Delay Model for Horizontal Axis Wind
Turbine Performance Prediction, ?ASME/AIAA Joint Wind Energy
Symposium, AIAA Paper 1998?0021, Reno, NV, Jan. 1998, pp. 9?19.
[56] Pierce, K. and Hansen, A. C., ?Prediction of Wind Turbine Rotor Loads Using the
Beddoes?Leishman Model for Dynamic Stall,? Journal of Solar Energy
Engineering, Vol. 117, 1995, pp. 200?204.
[57] Beddoes, T. S., ?Representation of Aerofoil Behaviour,? Vertica, Vol. 7, 1983,
pp. 183?197.
[58] Leishman, J. G., and Beddoes, T. S., ?A SemiEmpirical Model for Dynamic
Stall,? Journal of American Helicopter Society, Vol. 34, No. 3, 1989, pp. 3?
17.
138
[59] Schreck, S. J., and Robinson, M. C., ?Dynamic Stall and Rotational
Augmentation in Recent Wind Turbine Aerodynamic Experiments,? In
Proceedings of 32nd AIAA Fluid Dynamics Conference and Exhibit, St. Louis,
MS, June 24?26, 2002.
[60] James L. Tangler, J. David Kocurek, ?Wind Turbine PostStall Airfoil
Performance Characteristics Guidelines for BladeElement Momentum
Methods?, 43rd AIAA Aerospace Sciences Meeting and Exhibit 1013 January
2005, Reno, Nevada.
[61] Viterna, L .A., and Corrigan, R. D., ?Fixed Pitch Rotor Performance of Large
Horizontal Axis Wind Turbines, ?DOE/NASA Workshop on Large Horizontal
Axis Wind Turbines, Cleveland, Ohio, July 1981.
[62] Viterna, L .A., and Janetzke, D. C., ?Theoretical and Experimental Power from
Large HorizontalAxis Wind Turbines,? NASA TM82944, Sept 1982.
[63] Robert E. Sheidahl R., E., Klimas P., C., ?Aerodynamic Characteristics of Seven
Symmetrical Airfoil Sections through 180Degree Angle of Attack for use in
Aerodynamic Analysis of Vertical Axis Wind Turbines?. Sandia National
Laboratories, SAND802114, 1981.
[64] F. B. Metzger, ?A Review of Propeller Noise Prediction Methodology 1919
1994?, NASA Contractor Report 198156, National Aeronautics and Space
Administration, Langley Research Center, June 1995.
[65] Kenneth R. Fidler, ?Subsystem acoustic testing of a Vertical Takeoff and Landing
Ducted Propeller Unmanned Aerial Vehicle? Technical Report, AMRSS04
05.
[66] Lowson, M. V., ?The sound field for singularities in motion.? Proceedings of the
Royal Society (London), Volume A286, 1965, pp. 559572.
[67] Succi G. P., ?Noise and Performance of General Aviation Aircraft: A Review of
the MIT Study,? SAE Technical Paper No. 810586, April 1981.
139
[68] Succi G. P., Munro, D. H., Zimmer, J. A., Dunbeck, P. D., Lasabee E. E., K. U.
and Kerrenbrock J. L., ?Noise and Performance of Propellers for Light
Aircraft,? NASA CR165732, 1981.
[69] Succi G. P., Munro, D. H. and Zimmer J. A., ?Experimental verification of
Propeller Noise Prediction,? AIAA Journal, Volume 2, No. 11, pp.1483
1491., November 1982.
[70] Sullivan J. P., ?The Effect of Blade Sweep on Propeller Performance,? AIAA
paper No. 77716, 1977.
[71] Anderson, Jr. J. D., ?Computational Fluid Dynamics, The basics with
Applications?. McGrawHill, Inc., New York, 1995.
[72] Th. E. Labrujere and P.J. Zandenbergen. ?On the application of a new version of
lifting surface theory to nonslender and kinked wings?. Journal of
Engineering mathematics, 7:8596, 1973.
[73] Gabasov, R.F. and Kirillova, F.M., Methods of Optimization, Optimization
Software Inc., New York, NY, 1988.
[74] Weber, D.C. and Skillings, J.H., A first Course in the Design of Experiments: A
linear Models Approach, CRC Press, Boca Raton, FL, 2000.
[75] Kennedy, J. Eberhart, R. ?Particle Swarm Optimization?, Proceedings of the
IEEE International Conference on Neural Networks, Perth, Australia 1995, pp
19421945.
[76] Fukuyama, Y., Takayama, S., Nakanishi, Y., and Yoshida, H. ?A particle swarm
optimization for reactive power and voltage control in electric power systems.
Proceedings of the Genetic and Evolutionary Computation?, Conference 1999
(GECCO 1999), Orlando, Florida, USA. pp. 15231528, 1999.
140
[77] Carlisle, A. and Dozier, G. ?Adapting particle swarm optimization to dynamic
environments?. Proceedings of International Conference on Artificial
Intelligence (ICAI 2000), Las Vegas, Nevada, USA. pp. 429434, 2000.
[78] Venter, G. and SobieszczanskiSobieski, J., "Multidisciplinary optimization of a
transport aircraft wing using particle swarm optimization," Structural and
Multidisciplinary Optimization, pp. preprint, 2004.
[79] Onwubolu, G. C. and Clerc, M., "Optimal path for automated drilling operations
by a new heuristic approach using particle swarm optimization," International
Journal of Production Research, vol. 4 pp. 473491, 2004.
[80] Zhou, Y., Zeng, G., and Yu, F., "Particle swarm optimizationbased approach for
optical finite impulse response filter design," Applied Optics, vol. 42, no. 8,
pp. 15031507, Mar.2003.
[81] Zheng, Y., Ma, L., Zhang, L., and Qian, J. Robust PID controller design using
particle swarm optimizer. Proceedings of IEEE International Symposium on
Intelligence Control 2003, pp. 974979, 2003.
[82] Jenigiri, S. ?Comparative study of efficiency of genetic algorithms and particle
swarm optimization technique to solve permutation problems?. 1996.
[83] Angeline, P. J. Evolutionary optimization versus particle swarm optimization:
philosophy and performance differences. Evolutionary Programming VII:
Proceedings of the Seventh Annual Conference on Evolutionary
Programming, 1998.
[84] Eberhart, R. C. and Shi, Y. Comparison between genetic algorithms and particle
swarm optimization. Evolutionary Programming VII: Proceedings of the
Seventh Annual Conference on Evolutionary Programming, San Diego, CA.
1998.
[85] Mitchell, M., An Introduction to Genetic Algorithms, MIT Press, Cambridge,
MA, 1998.
141
[86] Coley, D.A., An Introduction to Genetic Algorithms for Scientists and Engineers,
World Scientific Publishing Co., Singapore, 1999.
[87] Anderson, M.B., ?Design of a Missile Interceptor Using a Genetic Algorithm?,
Doctoral Dissertation, Auburn University, 1998.
[88] Dyer John, personal communication, Auburn University, 2007
142
APPENDIX A: BIOTSAVART LAW
The following describes the derivation of the BiotSavart law.
BiotSavart Law:
The continuity equation in nonconservation form is
0
q
=???+?
?
Dt
D
(A.1)
with ? being the density and q the velocity.
For incompressible fluids (?=const.) the continuity equation reduces to
0q
0q
=
?
?
+
?
?
+
?
?
=??
=??
z
w
y
v
x
u (A.2)
Angular velocity and vorticity
The motion of a fluid element consists of translation, rotation, and deformation. The
angular velocity of a fluid element is due to the velocity variation within the element. As
a result the fluid element may deform and rotate. The derivation of the angular velocity is
described in Ref. [16]. Thus the vector form of the angular velocity is
143
V???=
2
1
? (A.3)
Vorticity ? is defined as twice the angular velocity.
curlVV =??=?? ?2? (A.4)
To obtain the velocity field as a function of the vorticity distribution, Equation (A.4) must
be inverted. As described in Ref. [16] the velocity field may be expressed as the curl of a
vector field B, such that
Bq ??= (A.5)
Since the curl of a gradient vector field is zero, B is indeterminant to within the gradient
of a scalar function of position and time, and B can be selected such that
0=?? B (A.6)
The vorticity then becomes
( ) ( ) BBB
2
q? ?????=????=??= (A.7)
The application of Equation A.6 reduces this to Poisson?s equation for the vector
potential B:
B
2
? ??= (A.8)
The solution of this equation, using Green?s theorem (see Ref. [29]) is
144
dVB
V
?
?
=
10
rr
?
4
1
?
(A.9)
Here B is evaluated at point P (which is a distance r
0
from the origin, shown in Figure
A.1) and is a result of integrating the vorticity ? (at point r
1
) within the volume V. The
velocity field is then the curl of B:
dV
V
?
?
??=
10
rr
?
4
1
q
?
(A.10)
Figure A.1 Velocity at point P due to a vortex distribution
Before proceeding with the integration, an infinitesimal piece of the vorticity filament ?,
as shown in Figure A.2 is considered. The crosssectional area dS is selected such that it
P
.
r
0
r
1
?
dV
V
r
1
r
0
145
is normal to ?, and the direction dl on the filament is
dld
?
?
l = (A.11)
Also, the circulation ? is
dS?=? (A.12)
and
dldSdV = (A.13)
Figure A.2 Induced velocity at point P by a vortex segment
so that
dl
ds
?
P
r
0
r
1
Origin
r
0
? r
1
?
146
1010
rr
l
rr
?
?
???=
?
??
d
dV (A.14)
Carrying out the curl operation while keeping r
1
and dl fixed we get
( )
3
10
10
10
rr
rrl
rr
l
?
??
?=
?
???
d
d
(A.15)
Substitution of this result back into Equation A.6 results in the BiotSavart law, which
states
( )
?
?
??
?
=
3
10
10
rr
rrl
4
q
d
?
(A.16)
or in differential form
( )
3
10
10
rr
rrl
4
q
?
??
?
=?
d
?
(A.17)
147
APPENDIX B: PROPELLER GA VARIABLES
Shown below are the propeller variables which define the propeller geometry, the number
of propeller blades as well as the max propeller span. If a variable pitch propeller is
considered, the optional variables, which offset the angle of attack by a constant value are
activated.
Propeller GA Variables
pr xray(1) Propeller span
nxbl xray(2) Number of propeller blades
ua_1 xray(3) BP Shape Function Parameter airfoil upper side
ua_2 xray(4) BP Shape Function Parameter airfoil upper side
ua_3 xray(5) BP Shape Function Parameter airfoil upper side
ua_4 xray(6) BP Shape Function Parameter airfoil upper side
ua_5 xray(7) BP Shape Function Parameter airfoil upper side
ua_6 xray(8) BP Shape Function Parameter airfoil upper side
ua_7 xray(9) BP Shape Function Parameter airfoil upper side
ua_8 xray(10) BP Shape Function Parameter airfoil upper side
ua_9 xray(11) BP Shape Function Parameter airfoil upper side
ua_10 xray(12) BP Shape Function Parameter airfoil upper side
ua_11 xray(13) BP Shape Function Parameter airfoil upper sid
ua_12 xray(14) BP Shape Function Parameter airfoil upper side
ua_13 xray(15) BP Shape Function Parameter airfoil upper side
ua_14 xray(16) BP Shape Function Parameter airfoil upper side
ua_15 xray(17) BP Shape Function Parameter airfoil upper side
ua_16 xray(18) BP Shape Function Parameter airfoil upper side
ua_17 xray(19) BP Shape Function Parameter airfoil upper side
ua_18 xray(20) BP Shape Function Parameter airfoil upper side
ua_19 xray(21) BP Shape Function Parameter airfoil upper side
ua_20 xray(22) BP Shape Function Parameter airfoil upper side
fau_1 xray(23) Class Function upper Airfoil Parameter
fau_2 xray(24) Class Function upper Airfoil Parameter
aoa_1 xray(25) BP Angle of Attack Function Parameter
aoa_2 xray(26) BP Angle of Attack Function Parameter
148
aoa_3 xray(27) BP Angle of Attack Function Parameter
aoa_4 xray(28) BP Angle of Attack Function Parameter
aoa_5 xray(29) Angle of Attack max min value
chord_1 xray(30) BP Chord Length Function Parameter
chord_2 xray(31) BP Chord Length Function Parameter
chord_3 xray(32) BP Chord Length Function Parameter
chord_4 xray(33) BP Chord Length Function Parameter
chord_5 xray(34) Chord Length max min
sweep_1 xray(35) BP Propeller Sweep Function Parameter
sweep_2 xray(36) BP Propeller Sweep Function Parameter
sweep_3 xray(37) BP Propeller Sweep Function Parameter
sweep_4 xray(38) BP Propeller Sweep Function Parameter
sweep_5 xray(39) Propeller Sweep in multiples of max chord length
cprp xray(40) Propeller rotational point
149
APPENDIX C: PROPELLER WAKE INVESTIGATION
Propeller Wake Investigation of a two and four Bladed Propeller
Below are shown the results of the thrust and power obtained from the propeller
performance program based on the wake extension. The propeller speed was 5800 rpm,
the free stream velocity V
inf
= 21.0 [m/s] and propeller diameter was d=0.32 [m].
Wake extenstion 2 bladed propeller 4 bladed propeller
in revolutions Thrust [N] Power [W] Thrust [N] Power [W]
0.125 3.97 104.07 7.70 203.43
0.250 3.93 103.22 7.45 198.42
0.500 3.87 102.01 7.22 193.83
1.000 3.83 101.28 7.10 191.36
2.000 3.82 100.99 7.05 190.48
3.000 3.81 100.94 7.04 190.23
4.000 3.81 100.92 7.03 190.16
5.000 3.81 100.91 7.03 190.12
150
APPENDIX D: PROPELLER PANEL DISCRETIZATION
Propeller Panel Discretization and Panel Density
Below are the results of the panel discretization and panel density investigation. This was
done to determine the least number of panels which still produce accurate aerodynamic
performance prediction in an effort to reduce computation times during the optimization
process. The investigation contains two panel arrangement types. The first one uses a full
cosine distribution of the panels in the chordwise direction. This assures that the panels
are denser at the leading and training edge. The second panel arrangement has besides the
full cosine distribution of the panels in the chordwise direction a half cosine distribution
in the spanwise direction. In addition to that the Table below shows also the effects on
the thrust when the number of spanwise and chordwise panels are changed.
Full cosine panel distribution
in the chordwise direction
Full cosine panel distribution in the chordwise and
half cosine panel distribution in the spanwise direction
Panel Numbers
Chordwise x Spanwise Thrust [N]
Panel Numbers
Chordwise x Spanwise Thrust [N]
3x20 2.8144 7x7 3.9418
4x20 2.9766 7x10 3.8095
5x20 3.3385 7x12 3.7544
6x20 3.6052 7x14 3.7281
7x5 3.936 7x16 3.7195
7x10 3.8262 7x18 3.7183
7x20 3.8162 7x20 3.7203
7x25 3.8152 7x24 3.7269
7x30 3.8145 7x26 3.7316
151
APPENDIX E: NACA 109622 PROPELLER GEOMETRY
NACA 109622 Propeller Geometry Data
The propeller performance program was validated using a NACA 109622 straight blade
propeller. The data below are from Ref. [23], where a three bladed propeller was
simulated with a blade angle of ?
.75
=45.4?.at the 75% spanwise position. The x, y and z
coordinates are given of the spanwise geometric description with C/r being the section
chord to span ratio.
x y z c/r ?
0.0 0.275 0.0 0.2415 1.15581
0.0 0.3475 0.0 0.2415 1.07155
0.0 0.42 0.0 0.2415 1.0034
0.0 0.4925 0.0 0.2415 0.94246
0.0 0.565 0.0 0.2415 0.88481
0.0 0.6375 0.0 0.2415 0.83769
0.0 0.71 0.0 0.2415 0.77661
0.0 0.7825 0.0 0.2415 0.74869
0.0 0.855 0.0 0.2415 0.70855
0.0 0.9275 0.0 0.2415 0.67713
0.0 1 0.0 0 0
152
The propeller geometry is generated based on the data from Table above. The propeller
has 10 spanwise and two chordwise panels. Thus the individual nodal and collocation
points of the first propeller blade are
Nodal Points
Leading Edge Center Trailing Edge
x y z x y z x y z
0.0000 0.2750 0.0000 0.1113 0.2750 0.0449 0.2225 0.2750 0.0899
0.0000 0.3475 0.0000 0.1071 0.3475 0.0542 0.2142 0.3475 0.1083
0.0000 0.4200 0.0000 0.1032 0.4200 0.0613 0.2063 0.4200 0.1226
0.0000 0.4925 0.0000 0.0992 0.4925 0.0675 0.1984 0.4925 0.1350
0.0000 0.5650 0.0000 0.0952 0.5650 0.0731 0.1903 0.5650 0.1462
0.0000 0.6375 0.0000 0.0916 0.6375 0.0775 0.1832 0.6375 0.1550
0.0000 0.7100 0.0000 0.0867 0.7100 0.0829 0.1734 0.7100 0.1659
0.0000 0.7825 0.0000 0.0844 0.7825 0.0853 0.1687 0.7825 0.1707
0.0000 0.8550 0.0000 0.0809 0.8550 0.0886 0.1618 0.8550 0.1773
0.0000 0.9275 0.0000 0.0781 0.9275 0.0911 0.1561 0.9275 0.1823
0.0710 1.0000 0.0700 0.0779 1.0000 0.0913 0.0860 1.0000 0.1000
153
Collocation Points
Leading Edge Panels Trailing Edge Panels
xc yc zc xc yc zc
0.0556 0.3113 0.0225 0.1669 0.3113 0.0674
0.0535 0.3838 0.0271 0.1606 0.3838 0.0812
0.0516 0.4563 0.0307 0.1547 0.4563 0.0920
0.0496 0.5288 0.0337 0.1488 0.5288 0.1012
0.0476 0.6013 0.0365 0.1428 0.6013 0.1096
0.0458 0.6738 0.0387 0.1374 0.6738 0.1162
0.0434 0.7463 0.0415 0.1301 0.7463 0.1244
0.0422 0.8188 0.0427 0.1266 0.8188 0.1280
0.0404 0.8913 0.0443 0.1213 0.8913 0.1330
0.0390 0.9450 0.0456 0.1171 0.9450 0.1367
154
APPENDIX F: NACA 0009 AIRFOIL DATA
The lift and drag coefficient data of a NACA 0009 airfoil section.
[63]
are taken to generate
the empirical stall model for the propeller.
155
156
APPENDIX G: PROPELLER GEOMETRY OF SKYHAWK 172
The geometric description of the ? scale Cessna 172 Skyhawk propeller as described by
Succi
[69]
and the generated geometry for the numerical verification of the test case is
given in the following table.
Geometric propeller description by Succi:
Radius r/R Chord C/R Thickness t/C Pitch ?
0.133 0.149 0.43 31.5
0.244 0.153 0.204 30.9
0.32 0.153 0.156 29.5
0.4 0.154 0.127 27.1
0.48 0.153 0.109 24.6
0.64 0.141 0.092 20.3
0.8 0.119 0.092 17.3
0.88 0.103 0.08 16.1
0.96 0.079 0.08 15.1
1 0 0 14.3
157
The collocation pints for the for the Cessna 172 1C160 blade model are
xControl Point yControl Point zControl Point
0.00471 0.03879 0.00773
0.00473 0.05218 0.00786
0.00469 0.06346 0.00795
0.00459 0.07263 0.00801
0.00447 0.08204 0.00809
0.00431 0.09169 0.00821
0.00413 0.10135 0.00830
0.00394 0.11100 0.00836
0.00361 0.12548 0.00829
0.00317 0.14478 0.00809
0.00274 0.16408 0.00770
0.00233 0.18339 0.00714
0.00203 0.19787 0.00663
0.00182 0.20752 0.00619
0.00160 0.21717 0.00563
0.00136 0.22682 0.00494
0.00104 0.23406 0.00389
0.00064 0.23889 0.00247
0.01412 0.03879 0.02318
0.01419 0.05218 0.02357
0.01407 0.06346 0.02385
0.01378 0.07263 0.02402
0.01340 0.08204 0.02428
0.01293 0.09169 0.02463
0.01240 0.10135 0.02491
0.01182 0.11100 0.02509
0.01084 0.12548 0.02488
0.00951 0.14478 0.02426
0.00822 0.16408 0.02310
0.00699 0.18339 0.02142
0.00609 0.19787 0.01990
0.00547 0.20752 0.01858
0.00480 0.21717 0.01689
0.00408 0.22682 0.01483
0.00311 0.23406 0.01167
0.00192 0.23889 0.00740
158
The nodal pints for the for the Cessna 172 1C160 blade model are
xNodal Point yNodal Point zNodal Point
0.00000 0.03209 0.00000
0.00000 0.04549 0.00000
0.00000 0.05888 0.00000
0.00000 0.06805 0.00000
0.00000 0.07722 0.00000
0.00000 0.08687 0.00000
0.00000 0.09652 0.00000
0.00000 0.10617 0.00000
0.00000 0.11582 0.00000
0.00000 0.13513 0.00000
0.00000 0.15443 0.00000
0.00000 0.17374 0.00000
0.00000 0.19304 0.00000
0.00000 0.20269 0.00000
0.00000 0.21234 0.00000
0.00000 0.22200 0.00000
0.00000 0.23165 0.00000
0.00000 0.23647 0.00000
0.00000 0.24130 0.00000
0.00939 0.03209 0.01533
0.00944 0.04549 0.01558
0.00948 0.05888 0.01584
0.00929 0.06805 0.01595
0.00909 0.07722 0.01607
0.00878 0.08687 0.01631
0.00846 0.09652 0.01654
0.00807 0.10617 0.01667
0.00768 0.11582 0.01678
0.00677 0.13513 0.01639
0.00590 0.15443 0.01596
0.00505 0.17374 0.01485
0.00427 0.19304 0.01371
0.00385 0.20269 0.01283
0.00345 0.21234 0.01194
0.00295 0.22200 0.01057
0.00248 0.23165 0.00920
0.00167 0.23647 0.00636
0.00089 0.24130 0.00351
0.01879 0.03209 0.03066
0.01887 0.04549 0.03117
0.01896 0.05888 0.03168
0.01857 0.06805 0.03191
0.01818 0.07722 0.03213
0.01756 0.08687 0.03261
159
0.01693 0.09652 0.03308
0.01615 0.10617 0.03333
0.01537 0.11582 0.03357
0.01355 0.13513 0.03278
0.01180 0.15443 0.03191
0.01011 0.17374 0.02970
0.00854 0.19304 0.02742
0.00770 0.20269 0.02565
0.00689 0.21234 0.02388
0.00591 0.22200 0.02115
0.00497 0.23165 0.01840
0.00334 0.23647 0.01272
0.00179 0.24130 0.00701
160
APPENDIX H: WIND TUNNEL DATA CAM 13x7 PROPELLER
Propeller wind tunnel data for a CAM 13x7 propeller. The propeller hub allowed for
blade pitch changes to maximize the thrust. Motor efficiencies are about 72% for a power
input of 67 [W]. The first Table below contains thrust and torque data of the CAM 13x7
propeller at 60 feet free stream velocity and a constant power input of 67 Watts, while
changing blade pitch.
Pitch change 60 feet/sec free stream velocity
Thrust [lbf] Torque [lb inch]
0.032 0.685
0.282 0.882
0.335 0.987
0.345 1.035
0.348 1.072
0.37 1.24
0.362 1.308
0.345 1.369
0.339 1.534
0.27 1.565
The Table below contains the thrust data of the 13x7 propeller in a fixed pitch mode as it
changes with increased free stream velocity and with a constant power input of 67 Watts.
Free Stream Velocity [ft/sec] Fixed pitch propeller thrust [lbf]
0 1.07
20 0.92
30 0.7
40 0.524
50 0.26
60 0.22