Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Aerodynamic optimization using high-order finite-volume CFD simulations Azab, Mohammad Baher 2011

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2011_fall_azab_mohammad.pdf [ 2.42MB ]
Metadata
JSON: 24-1.0072307.json
JSON-LD: 24-1.0072307-ld.json
RDF/XML (Pretty): 24-1.0072307-rdf.xml
RDF/JSON: 24-1.0072307-rdf.json
Turtle: 24-1.0072307-turtle.txt
N-Triples: 24-1.0072307-rdf-ntriples.txt
Original Record: 24-1.0072307-source.json
Full Text
24-1.0072307-fulltext.txt
Citation
24-1.0072307.ris

Full Text

Aerodynami
 Optimization Using High-Order Finite-Volume CFD Simulations by Mohammad Baher Azab B.S
., Aerospa
e Engineering Cairo University (Egypt), 1999 M.S
., Aerospa
e Engineering Cairo University (Egypt), 2002 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Fa
ulty of Graduate Studies (Me
hani
al Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Van
ouver) O
tober 2011 © Mohammad Baher Azab 2011 Abstra
t The growth of 
omputer power and storage 
apa
ity allowed engineers to ta
kle engineering design as an optimization problem. For transport air
raft, drag mini- mization is 
riti
al to in
rease range and redu
e operating 
osts. Lift and geometri 
onstraints are added to the optimization problem to meet payload and rigidity 
on- straints of the air
raft. Higher order methods in CFD simulations have proved to be a valuable tool and are expe
ted to repla
e 
urrent se
ond order CFD methods in the near future; therefore, exploring the use of higher order CFD methods in aerodynami
 optimization is of great resear
h interest and is one goal of this thesis. Gradient-based optimization te
hniques are well known for fast 
onvergen
e, but they are only lo
al minimizers; therefore their results depend on the starting point in the design spa
e. The gradient-independent optimization te
hniques 
an nd the global minimum of an obje
tive fun
tion but require vast 
omputational eort; therefore, for global optimization with reasonable 
omputational 
ost, a hybrid op- timization strategy is needed. A new least-squares based geometry parametrization is used to des
ribe airfoil shapes and a semi-torsional spring analogy mesh morphing tool updates the grid everywhere when the airfoil geometry 
hanges during shape optimization. For the gradient based optimization s
heme, both se
ond and fourth order sim- ulations have been used to 
ompute the obje
tive fun
tion; the adjoint approa
h, well known for its low 
omputational 
ost, has been used for gradient 
omputa- tion and mat
hes well with nite dieren
e gradient. The gradient based optimizer have been tested for subsoni
 and transoni
 inverse design problems and for drag minimization without and with lift 
onstraint to validate the developed optimizer. The optimization s
heme used is Sequential Quadrati
 Programming (SQP) with the BFGS approximation of the Hessian matrix. A mesh renement study is presented for an aerodynami
ally 
onstrained drag minimization problem to show how se
ond ii and fourth order optimal results behave with mesh renement. A hybrid parti
le swarm / BFGS s
heme has been developed for use as a global optimizer. It has been tested on a drag minimization problem with a lift 
onstraint; the hybrid s
heme obtained a sho
k free proles, while gradient-based optimization 
ould not in general. iii Prefa
e The resear
h ideas and methods explored in the three 
o-authored manus
ripts of this thesis are the fruits of a 
lose working relationship between Dr Carl Ollivier-Goo
h and Mohammad Azab. The implementation of the methods, the data analysis, and the manus
ript preparation were done by Mohammad Azab with invaluable guidan
e from Carl Ollivier-Goo
h throughout the pro
ess. iv Table of Contents Abstra
t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv Nomen
lature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Finite Volume Flow Solver . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Numeri
al Aerodynami
 Optimization . . . . . . . . . . . . . . . . 5 1.2.1 Gradient-based aerodynami
 optimization . . . . . . . . . . 6 1.2.2 Gradient free optimization . . . . . . . . . . . . . . . . . . . 8 1.3 Contributions of the Thesis . . . . . . . . . . . . . . . . . . . . . . . 10 2 GEOMETRY PARAMETRIZATION . . . . . . . . . . . . . . . . . 12 2.1 Survey of Geometri
 Parametrization Te
hniques . . . . . . . . . . . 12 2.1.1 Analyti
al parametrization . . . . . . . . . . . . . . . . . . . 13 2.1.2 Pie
e-wise spline parametrization . . . . . . . . . . . . . . . 14 2.1.3 CAD parametrization . . . . . . . . . . . . . . . . . . . . . . 14 2.1.4 Free form deformation (FFD) . . . . . . . . . . . . . . . . . 15 2.1.5 Multidis
iplinary aerodynami
/stru
tural shape optimization using deformations (MASSOUD) . . . . . . . . . . . . . . . . 15 2.2 New Least Squares Parametrization Te
hnique . . . . . . . . . . . . 16 v 2.2.1 Airfoil geometry parametrization . . . . . . . . . . . . . . . . 16 2.2.2 Thi
kness 
onstraint . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.3 Validation 
ases . . . . . . . . . . . . . . . . . . . . . . . . . 23 3 MESH MORPHING AND MESH SENSITIVITY . . . . . . . . . 30 3.1 Mesh Morphing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2 Testing Mesh Morphing . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3 Mesh Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4 GRADIENT CALCULATION USING ADJOINT APPROACH 35 4.1 Forward and Adjoint Formulations . . . . . . . . . . . . . . . . . . 35 4.2 Element and Fa
e Geometri
 Properties Dependen
y on Mesh Coor- dinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.2.1 Element area mesh dependen
y . . . . . . . . . . . . . . . . 45 4.2.2 Fa
e length mesh dependen
y . . . . . . . . . . . . . . . . . 46 4.2.3 Fa
e normal mesh dependen
y . . . . . . . . . . . . . . . . . 49 4.3 Fa
e Flow Properties Dependen
y on The Mesh . . . . . . . . . . . 50 4.3.1 Fa
e ow properties re
onstru
tion . . . . . . . . . . . . . . 50 4.3.2 Mesh dependen
e of the fa
e ow property re
onstru
tion . . 57 5 GRADIENT VALIDATION . . . . . . . . . . . . . . . . . . . . . . . 60 5.1 Subsoni
 Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.2 Transoni
 Test Case with No Limiter . . . . . . . . . . . . . . . . . 61 5.3 Transoni
 Test Case with Limiter . . . . . . . . . . . . . . . . . . . 65 5.4 Sensitivity of nite dieren
e gradient to perturbation magnitude . 67 6 GRADIENT BASED OPTIMIZATION TEST CASES . . . . . . 70 6.1 Subsoni
 Inverse Design . . . . . . . . . . . . . . . . . . . . . . . . . 70 6.2 Transoni
 Inverse Design . . . . . . . . . . . . . . . . . . . . . . . . 71 6.3 Drag Minimization without Lift Constraint . . . . . . . . . . . . . . 74 6.4 Drag Minimization with Lift Constraint . . . . . . . . . . . . . . . . 80 6.5 Mesh Renement Study of a Drag Minimization with Lift Constraint 86 vi 7 PARTICLE SWARM OPTIMIZATION AND A NEW HYBRID OPTIMIZATION METHOD . . . . . . . . . . . . . . . . . . . . . . 88 7.1 Introdu
tion to Gradient Free Optimization . . . . . . . . . . . . . . 88 7.2 Swarm Intelligen
e . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 7.2.1 Premature 
onvergen
e dete
tion . . . . . . . . . . . . . . . . 93 7.2.2 Redening de
ision variables range . . . . . . . . . . . . . . . 94 7.2.3 Regrouping and position 
lipping . . . . . . . . . . . . . . . 94 7.3 Hybrid SQP-RegPSO Te
hnique . . . . . . . . . . . . . . . . . . . . 95 7.4 Optimization Test Cases . . . . . . . . . . . . . . . . . . . . . . . . 95 7.4.1 Constrained drag optimization of the NACA 0012 . . . . . . 97 7.4.2 Constrained drag minimization of NACA 00083 . . . . . . . 102 7.4.3 Aerodynami
 and thi
kness 
onstraint drag minimization of NACA 0012 . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.4.4 Aerodynami
 and thi
kness 
onstraint drag minimization of NACA 00083 . . . . . . . . . . . . . . . . . . . . . . . . . . 108 7.5 Comparing SQP-RegPSO with RegPSO-SQP optimization strategies 112 8 CONCLUSIONS AND RECOMMENDATIONS . . . . . . . . . . 119 8.1 Contributions and Con
lusions . . . . . . . . . . . . . . . . . . . . . 119 8.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 vii List of Tables 2.1 RMS error in dierent parametrized airfoil geometries . . . . . . . . 29 4.1 Two point Gauss quadrature rule . . . . . . . . . . . . . . . . . . . . 48 5.1 The magnitude of se
ond and fourth order CL gradients and angles between the evaluated gradients for NACA 0012 in subsoni
 ow. . . 62 5.2 The magnitude of se
ond and fourth order CD gradients and angles between the evaluated gradients for NACA 0012 in an unlimited tran- soni
 ow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.3 Magnitudes of se
ond and fourth order CD gradients and angles be- tween the evaluated gradients for NACA 0012 in Venkatakrishnan limited transoni
 ow . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4 The magnitude of se
ond and fourth order CD gradients and angles between the evaluated gradients for NACA 0012 using higher order limiter in transoni
 ow. . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.5 The magnitudes and angles between the evaluated modied fourth order CD gradients using adjoint, and nite dieren
e for NACA 0012 using Venkatakrishnan and higher order limiters in transoni
 ow. . 67 5.6 Finite dieren
e drag gradient sensitivity with respe
t to perturbation amplitude ǫ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.1 Aerodynami
 
oe
ients of original and optimized RAE 2822 airfoil at transoni
 
onditions. . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.2 Aerodynami
 
oe
ients of original and optimized RAE 2822 airfoil at transoni
 
onditions. . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.3 Lift penalty weight ee
t on Drag minimization of RAE 2822 with lift 
onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 viii 7.1 Angles between adjoint, FD., and SQP-hybrid ve
tors . . . . . . . . 101 7.2 Angles between adjoint, FD, and BFGS-hybrid best solutions ve
tors 103 7.3 Thi
kness penalty weight impa
t on the optimization results of NACA 00083 using hybrid s
heme . . . . . . . . . . . . . . . . . . . . . . . . 110 7.4 Optimization results of the hybrid s
heme in Phase I, II . . . . . . . 113 7.5 Optimization results of the hybrid T s
heme in Phase I, II . . . . . . 113 ix List of Figures 1.1 Airfoil aerodynami
 optimization 
y
le using gradient-based optimiza- tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Least square surfa
e presentation of RAE 2822 airfoil using two poly- nomials P1 (x) & P2 (x) tted using nine 
ontrol points. . . . . . . . 18 2.2 Parametrized airfoil using four polynomials . . . . . . . . . . . . . . 21 2.3 Parametrized NACA 0011SC, 'O' are the original airfoil ordinates . . 24 2.4 Parametrized NACA 0012, 'O' are the original airfoil ordinates . . . 24 2.5 Parametrized NACA 6509, 'O' are the original airfoil ordinates . . . 25 2.6 Parametrized NACA 16006, 'O' are the original airfoil ordinates . . . 25 2.7 Parametrized NACA 63412, 'O' are the original airfoil ordinates . . . 26 2.8 Parametrized NACA 644421, 'O' are the original airfoil ordinates . . 26 2.9 Parametrized RAE2822, 'O' are the original airfoil ordinates . . . . . 27 2.10 Parametrized LV2 laminar airfoil, parametrized using 5 polynomials per surfa
e, 'O' are the original airfoil ordinates . . . . . . . . . . . . 28 2.11 LV2 airfoil parametrized using dierent methods with 20 design vari- ables, presented with permission of Brezillon. . . . . . . . . . . . . . 28 3.1 S
hemati
 drawing of an edge pq and its fa
ing angles a and b. . . . 31 3.2 Mesh movement s
heme results of doubling the thi
kness of NACA 0012 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Mesh movement results, large outer boundary deformation of a re
t- angular domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.1 S
hemati
 drawing of an element fa
e κ, and illustration of its left and right sides . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2 General triangular element with unit normal n̂ on one of its fa
es κ. 44 x 4.3 Boundary element with 
urved fa
e for high order integration s
heme 47 4.4 S
hemati
 drawing of a 
urved fa
e with two Gauss points used . . . 48 4.5 S
hemati
 drawing of rst, se
ond and third neighbor layers of trian- gular element k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.1 The pressure sensitivity with respe
t to one of the design 
ontrol points 
omputed for subsoni
 ow over NACA 0012, 
omparing the sensitivity 
al
ulation of Eq. 4.6 with nite dieren
e results. . . . . 61 5.2 CL gradient error in se
ond and fourth order s
hemes with respe
t to the design points normalized by gradient magnitude. . . . . . . . . . 62 5.3 The pressure sensitivity with respe
t to one of the design 
ontrol points 
omputed for an unlimited transoni
 ow around NACA 0012, 
omparing the sensitivity 
al
ulation of Eq. 4.6 with nite dieren
e results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.4 CD gradient error in se
ond and fourth order s
hemes with respe
t to the design points, normalized by gradient magnitude. . . . . . . . . . 64 5.5 The pressure sensitivity with respe
t to one of the design 
ontrol points 
omputed for an unlimited transoni
 ow around NACA 0012. 65 5.6 The normalized CD gradient error in se
ond, fourth, and modied fourth order s
hemes with respe
t to the design points in a limited transoni
 ow (Venkatakrishnan limiter) . . . . . . . . . . . . . . . . 68 5.7 The normalized CD gradient error in se
ond, fourth, and modied fourth order s
hemes with respe
t to the design points in a limited transoni
 ow (higher order limiter). . . . . . . . . . . . . . . . . . . 68 6.1 Subsoni
 NACA 2412 inverse design pressure distributions for the ini- tial, target, and optimized airfoil proles . . . . . . . . . . . . . . . . 72 6.2 Se
ond and fourth order optimization 
onvergen
e history. . . . . . . 72 6.3 Gradient norm history . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.4 Subsoni
 inverse design optimal airfoil shapes . . . . . . . . . . . . . 73 6.5 The dieren
e between the target prole and the optimized proles, se
ond and fourth order . . . . . . . . . . . . . . . . . . . . . . . . . 74 xi 6.6 Subsoni
 NACA 2412 inverse design pressure distributions for the ini- tial, target, and optimized airfoil proles . . . . . . . . . . . . . . . . 75 6.7 Se
ond and fourth order optimization 
onvergen
e history. . . . . . . 75 6.8 Gradient norm history . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.9 The dieren
e between the target prole and the optimized proles, se
ond and fourth order . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.10 Transoni
 inverse design optimal airfoil shapes. . . . . . . . . . . . . 77 6.11 Pressure 
ontours of RAE 2822 at Ma
h 0.73 and angle of atta
k 2 . 78 6.12 Optimized RAE 2822. . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6.13 Optimization surfa
e displa
ements of the original RAE2822 surfa
es. 79 6.14 Pressure distribution 
omparison for the original and the optimized geometries. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 6.15 Obje
tive fun
tion response surfa
e along positive and negative gra- dient dire
tions 
entered at the optimal prole of RAE 2822 transoni drag minimization without lift 
onstraint. . . . . . . . . . . . . . . . 80 6.16 Se
ond and fourth order optimization 
onvergen
e history. . . . . . . 81 6.17 Surfa
e pressure distribution, of the initial and optimized geometries 82 6.18 Se
ond and fourth order optimization 
onvergen
e history. . . . . . . 82 6.19 Optimal shapes 
omparison: se
ond order, fourth order, and opti- mized prole by Brezillon and Gauger 
ompared with the original RAE2822. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.20 Se
ond and fourth order optimized prole surfa
e displa
ements from the initial shape. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.21 Obje
tive fun
tion response surfa
e along positive and negative gra- dient dire
tions 
entered at the optimal prole of RAE 2822 transoni drag minimization without lift 
onstraint. . . . . . . . . . . . . . . . 84 6.22 initial and optimal pressure distribution obtained by Brezillon and Gauger [11℄ (presented with permission) . . . . . . . . . . . . . . . . 85 6.23 Optimal CD value with mesh renement . . . . . . . . . . . . . . . . 87 7.1 Geneti
 Algorithm optimization . . . . . . . . . . . . . . . . . . . . . 90 7.2 Velo
ity 
omponents and position update of a parti
le . . . . . . . . 93 7.3 Pressure Field of the original NACA 0012 and the optimized proles 98 xii 7.4 Surfa
e pressure distribution of original NACA 0012 and the opti- mized proles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 7.5 Airfoil shapes of original NACA 0012 and the optimized proles . . . 99 7.6 Fun
tion minimization iterations in the se
ond phase of the hybrid s
heme (RegPSO) with a NACA 0012 airfoil as the starting point. . 100 7.7 Obje
tive fun
tion values between SQP and hybrid optimal points (starting geometry: NACA 0012). . . . . . . . . . . . . . . . . . . . . 101 7.8 Pressure surfa
e of NACA 00083 and the optimized proles . . . . . 103 7.9 Surfa
e pressure of NACA 00083 and the optimized proles . . . . . 104 7.10 Prole 
omparison of NACA 00083, BFGS optimal, and hybrid opti- mal airfoils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.11 Fun
tion minimization iterations in the se
ond phase of the hybrid s
heme (RegPSO) of NACA 00083 start . . . . . . . . . . . . . . . . 105 7.12 Obje
tive fun
tion values between BFGS and hybrid optimal points, starting geometry NACA 00083 . . . . . . . . . . . . . . . . . . . . . 105 7.13 Pressure surfa
e of NACA 0012 and the optimized proles with thi
k- ness 
onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.14 Surfa
e pressure of NACA 0012 and the optimized prole with thi
k- ness 
onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.15 Pressure surfa
e of NACA 00083 and the optimized proles with thi
k- ness 
onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 7.16 Surfa
e pressure of NACA 00083 and the optimized prole with thi
k- ness 
onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7.17 Prole 
omparison of the hybrid s
heme optimal proles, starting ge- ometries NACA 0012 and NACA 00083 . . . . . . . . . . . . . . . . 110 7.18 Thi
kness penalty weight impa
t on the Optimal pressure distribution of NACA 00083 using hybrid optimization s
heme. . . . . . . . . . . 111 7.19 Thi
kness penalty weight impa
t on the Optimal prole with NACA 00083 starting geometry using hybrid optimization s
heme. . . . . . 111 7.20 Hybrid and hybrid T optimal pressure distribution for NACA 0012 (Case I). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 7.21 Hybrid and hybrid T optimal pressure distribution for NACA 00083 (Case II). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 xiii 7.22 Hybrid and hybrid T optimal pressure distribution for NACA 0012 with 10% thi
kness 
onstraint and unit thi
kness penalty weight (Case III). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 7.23 Hybrid and hybrid T optimal pressure distribution for NACA 00083 with 10% thi
kness 
onstraint and unit thi
kness penalty weight (Case IV). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.24 Hybrid and hybrid T optimal pressure distribution for NACA 00083 with 10% thi
kness 
onstraint and ten thi
kness penalty weight (Case V). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 xiv List of Algorithms 7.1 Parti
le swarm pseudo 
ode . . . . . . . . . . . . . . . . . . . . . . . 91 7.2 Hybrid SQP-RegPSO pseudo 
ode . . . . . . . . . . . . . . . . . . . 96 7.3 Pseudo 
ode of RegPSO-SQP hybrid T  optimization s
heme . . . . 112 xv Nomen
lature α Angle of atta
k △Q Change in 
onservative ow properties △t Time step δ norm Normalized swarm radius ∂R ∂Q ja
obian matrix ∂U ∂M Dependen
y of premitive ow variables solution on mesh. dCL,D,M dD Gradient of lift, drag. and moment 
oe
ient dF dD Obje
tive fun
tion gradient dR dD Residual dependen
y of design variable γ Spe
i
 heats ratio n̂ Fa
e unit normal Ωi Control volume i Uk Control volume avaraged primitive ow properties −→ F Flux ve
tor a
ross 
ontrol volume boundary −→g Swarm best known position in the design spa
e −→pi Parti
le i best known position in the design spa
e −→rb Position ve
tor of point b xvi −→ Vi Parti
le i velo
ity −→xi(k) Parti
le position in iteration k −→xi(k + 1) Parti
le position in iteration k+1 ∂M/∂D Mesh sensitivity ∂R/∂D Residual senstivity with respe
t to design variables ∂−→rb/∂Di Dependen
y of boundary point position −→rbve
tor on design variable D ΨL,D,M Adjoint ve
tor of of lift, drag. and moment 
oe
ient ρ Density ρ̃ Roe avaraged density à Roe averaged ux ja
obian matrix ã Roe avaraged speed of sound h̃ Roe avaraged enthalpy ũ Roe avaraged x-velo
ity ṽ Roe avaraged y-velo
ity ~F (QR) , ~F (QL) Right and left fa
e uxes CD Drag 
oe
ient CL Lift 
oe
ient CM Moment 
oe
ient Cp, Cv Spe
i
 heats at 
onstant pressure and volume Di Desigen variable i et Spe
i
 total energy Kib Stiness matrix intries relates interior nodes with boundary points xvii Kii Stiness matrix intries relates interior points PRec Flow property re
onstru
tion polynomial Ps Surfa
e pressure Q Conservative ow variables ve
tor QR, QL 
onserved ow properties at right and left of the 
ell fa
e Ubg Boundary Gauss point ow premitive properties Ufi Re
onstru
ted primitive ow properties at fa
e Gaussian integration point i Ui, Ub Interior and boundary mesh point displa
ements respe
tively Un Velo
ity in the dire
tion normal to the 
ontrol volume boundary ws Gaussian integration weight yc Geometry 
ontrol point variable y ordinate (design veriable) a Speed of sound E Energy per unit volume e Spe
i
 internal energy F Aerodynami
 obje
tive fun
tion GMRES Generalized minimal residual h Spe
i
 enthalpy M Ma
h number P Pressure R Residual RMS Root mean square T Temperature xviii A
knowledgments I would like to express my deepest gratitude to my supervisor Dr Carl Ollivier- Goo
h. You are su
h a great resear
h supervisor and gentleman. I would like also to a
knowledge my supervision 
ommittee, Dr Mi
hael Friedlander, Dr Philip Loewen, and Dr James Olson for their valuable remarks that helped me to justify my resear
h strategy. I would like also to thank my fellows in the resear
h and development department of the Egyptian Air For
e espe
ially Col. Dr Eng. Mohamed Ibrahim Mustafa and Major General Dr Abdalla El-Ramsisy for their great support during my resear
h in Egypt whi
h has had a great impa
t on my PhD resear
h progress. Great thanks goes also to my dear professor Dr Chen Grief for being there when I needed his advi
e. Finally, I would like to thank my wife and the light of my life, Deena, for her support and bringing happiness to my life. Mohammad Baher Azab xix Dedi
ation To my home 
ountry, the land of Pharaohs Egypt To the land of ho
key Canada To the souls of my parents Baher and Laila xx Chapter 1 INTRODUCTION Drag redu
tion of transport air
raft is of great importan
e be
ause it redu
es air
raft fuel 
onsumption and thereby redu
es operating 
osts and environmental impa
t in the form of pollution and global warming. Drag, lift and other aerodynami
 for
es 
an be predi
ted using CFD simulations, whi
h have be
ome an essential tool for aerodynami
 analysis and design. CFD simulations 
arried out using unstru
tured grids give a

urate aerodynami
 for
e predi
tions, and unstru
tured grids have the advantage of easily representing any 
omplex shape. As most transport air
raft travel at transoni
 speed, it is of great importan
e to redu
e the wave drag by weakening or entirely eliminating sho
k waves on the wing. Aerospa
e engineers use CFD simulations with numeri
al optimization te
hniques for aerodynami
 design; the optimization problem is to minimize an aerodynami
 obje
tive (often drag) by 
hanging geometri
 design variables, given an initial aerodynami
 shape and subje
t to some geometri
 and aerodynami
 
onstraints. This pro
ess requires a

urate assessment of the aerodynami
 
hara
teristi
s of a given geometry. The ow model used in this resear
h is Euler's ow model; although this model negle
ts vis
ous ow ee
ts, redu
ing the drag using this invis
id ow model will ultimately redu
e the drag on real 
ongurations as the invis
id drag is about 40% of the total drag [47℄. A higher order nite volume CFD solver developed by C. Ollivier-Goo
h and 
o- workers [68, 61, 56℄ is used in this resear
h. Optimization s
hemes 
an be divided into two main 
ategories: gradient-based and non-gradient-based s
hemes. Gradient-based s
hemes require 
omputing the ob- je
tive fun
tion value and its gradient with respe
t to the design variables, while the non-gradient-based te
hniques require only 
omputing the obje
tive fun
tion value. Gradient-based s
hemes are fast to 
onverge to a lo
al optimal point in the design spa
e 
ompared to the non-gradient-based methods, but the optimal solution found by gradient-based optimization depends on the starting point [33℄. Non-gradient- 1 based s
hemes 
an nd global minimum solution regardless of the starting point, but with larger 
omputational 
ost 
ompared to gradient-based s
heme. 1.1 Finite Volume Flow Solver The two dimensional integral form of Euler's equation 
an be written for a 
ontrol volume Ωi as ∂ ∂t ¨ Ωi QdV + ˛ dΩi −→ F · n̂dl = 0 (1.1) where n̂ = nxi + nyj is the outward pointing normal to the 
ontrol volume fa
es; Q and −→ F are respe
tively the 
onserved variable ve
tor and the ux a
ross the boundary of 
ontrol volume Ωi boundaries. These 
an be expressed as Q =  ρ ρu ρv E  , −→F · n̂ =  ρUn ρUnu+ nxp ρUnv + nyp (E + p)Un  p = (γ − 1) [ E − ρ ( u2 + v2 ) 2 ] Un = nxu+ nyv where is the uid density, u and v are the Cartesian velo
ity 
omponents, p is the pressure, and E is the energy per unit volume. Un is the velo
ity in the dire
tion normal to the 
ontrol volume boundary. Other thermodynami
 relations like the spe
i
 heats at 
onstant pressure (Cp) and 
onstant volume (Cv), the ratio of spe
i heats γ 
an be expressed in terms of the gas 
onstant for air, R: γ = Cp Cν , Cp = γR γ − 1 , Cν = R γ − 1 2 For a thermally and 
alori
ally perfe
t gas, thermodynami
 properties 
an be related by P = ρRT e = CvT et = e+ 1 2 ( u2 + v2 ) E = ρet h = e+ P ρ a = √ γRT where T is the temperature, e is the spe
i
 internal energy, et is the spe
i
 total energy, h is the spe
i
 enthalpy, and a is the speed of sound. For dis
rete solution of the Euler equations, ow properties are normalized by some referen
e values in order to redu
e the round-o errors in the dis
retized linear system; for external aerodynami
s, the normalized ow properties are as follows ρ̄ = ρ ρ∞ , ā = a a∞ , P̄ = P ρ∞a2∞ ū = u a∞ , v̄ = v a∞ , Ē = P̄ γ − 1 + 1 2 ρ̄ ( ū2 + v̄2 ) x̄ = x L , ȳ = y L , t̄ = ta∞ L With this normalization, the normalized Euler equations are identi
al to their di- mensional form with the addition of [̄·] to every variable. From this point on, the normalized ow properties are used and therefore [̄·] will be dropped. The ux is evaluated using Roe ux dieren
e splitting te
hnique [74℄ and eval- uated at ea
h 
ell fa
e κ using the following formula ~F = 1 2 [ ~F (QR) + ~F (QL)− ∣∣∣Ã∣∣∣ (QR −QL)] κ (1.2) where QR, QL are the 
onserved ow properties at right and left of the 
ell fa
e κ and the Roe averaged matrix à is the ux Ja
obian ∂ ~F/∂Q at Roe-averaged quantities 3 as follows ρ̃ = √ ρLρR ũ = uL + uR √ ρR ρL 1 + √ ρR ρL , ṽ = vL + vR √ ρR ρL 1 + √ ρR ρL h̃ = hL + hR √ ρR ρL 1 + √ ρR ρL , ã2 = (γ − 1) ( h̃− ( ũ2 + ṽ2 ) 2 ) (1.3) The Roe-averaged Ja
obian matrix à has four eigenvalues. Three of these are dis- tin
t, but the eigensystem is 
omplete. The Roe dissipation matrix 
an be written in terms the eigenvalues and eigenve
tors of the Ja
obian as ∣∣∣Ã∣∣∣ = X̃  ∣∣∣Ũn∣∣∣ 0 0 0 0 ∣∣∣Ũn∣∣∣ 0 0 0 0 ∣∣∣Ũn + ã∣∣∣ 0 0 0 0 ∣∣∣Ũn − ã∣∣∣  X̃ −1 (1.4) where X̃ is the eigenve
tor matrix evaluated using Roe-averaged ow quantities; the term ∣∣∣Ã∣∣∣ (QR −QL) 
an be written a

ording to Frink as follows [27, 26℄∣∣∣Ã∣∣∣ (QR −QL) = ∣∣∣Ã∣∣∣△Q = ∣∣∣△F̃1∣∣∣+ ∣∣∣△F̃4∣∣∣+ ∣∣∣△F̃5∣∣∣ (1.5) 4 where ∣∣∣△F̃1∣∣∣ = ∣∣∣Ũn∣∣∣  ( △ρ− △P ã2 ) 1 ũ ṽ ũ2+ṽ2 2 + ρ̃  0 △u− nx△U △v − ny△U ũ△u+ ṽ△v − Ũn△U   ∣∣∣△F̃4,5∣∣∣ = ∣∣∣Ũn ± ã∣∣∣ (△P ± ρ̃ã△U 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã  △ρ = ρR − ρL, △u = uR − uL, △v = vR − vL, △P = PR − PL, △U = nx△u+ ny△v Higher-order a

ura
y is obtained by least-squares re
onstru
tion of the non-
onserved variables U = [ ρ u v p ]T and Gauss quadrature in ux integration [68, 61℄. After integrating uxes around ea
h 
ontrol volume in the mesh, an impli
it time dis
retization leads to a sparse system of linear equations, whi
h for the simplest 
ase of a global time step 
an be written as,[ I △t + ∂R ∂Q ]n △Qn = Rn (1.6) where △Q = [ △ρ △ρu △ρv △E ]T , ∂R ∂Q is the global Ja
obian matrix, and R is the residual. The steady state solution is obtained iteratively when △Q → 0. In pra
ti
e, we use a quasi-Newton generalization of Eq. 1.6 that in
ludes residual-based lo
al time stepping [57℄ and solve the system using GMRES [75℄. 1.2 Numeri
al Aerodynami
 Optimization Aerodynami
 design used to rely on CFD simulations in 
onjun
tion with experimen- tal testing and engineering intuition of the designer. With the growth of high speed 
omputers, integrating numeri
al optimization s
hemes with CFD simulations has 5 be
ome possible and is now used for aerodynami
 design and optimization. Gradient- based optimization te
hniques are widely used be
ause they rea
h an optimized shape after a reasonable exe
ution time; however, the nal optimal shape is the lo
al min- imum lo
ated downhill from the optimization starting point. Non-gradient-based methods like geneti
 algorithm (GA) or parti
le swarm (PS) are slower to nd an optimum but 
an nd the global minimum regardless of the starting point; their drawba
k is the large number of iterations required to rea
h this minimum 
om- pared to gradient-based s
hemes. 1.2.1 Gradient-based aerodynami
 optimization Gradient-based optimization depends on evaluating the gradient of the obje
tive fun
tion with respe
t to the design variables and using the gradient in a linear model (steepest des
ent) or a quadrati
 model (Newton or quasi-Newton model) to nd a sear
h dire
tion; this sear
h dire
tion is the dire
tion in whi
h the de- sign variables should 
hange their values to minimize the obje
tive fun
tion [66℄. Gradient-based te
hniques have been widely used in aerodynami
 optimization due their fast 
onvergen
e to an optimal solution. The obtained optimal shape is bi- ased by the optimization starting point (initial aerodynami
 shape), and there is no guarantee that gradient-based methods 
an nd the global best optimal shape in the design spa
e [2℄. Hi
ks and Henne where among the rst to apply gradient- based optimization te
hniques in aerodynami
 design in the late 1970's [35℄. Sin
e then many resear
hers have investigated the use of gradient-based optimization te
h- niques like steepest des
ent and quadrati
 programming in aerodynami
 optimiza- tion [17, 31, 28℄. The most expensive part of gradient-based optimization is the gradient 
al
ulation. Hi
ks and Henne used nite dieren
e rules to 
al
ulate the obje
tive fun
tion gradient. This means that two CFD simulations were required for ea
h design variable to 
ompute the gradient (using a 
entral nite dieren
e rule), whi
h is 
omputationally expensive. The same strategy was applied by Consentino and Holst to optimize transoni
 wings [17℄. The use of the adjoint method, whi
h was originally applied to aerodynami
s problems by Jameson, to 
ompute the gra- dient redu
ed the 
omputational 
ost of gradient 
al
ulation to the 
ost of one ow simulation regardless of the number of design variables [41, 43, 60, 38, 39, 42, 11, 62℄. 6 Figure 1.1: Airfoil aerodynami
 optimization 
y
le using gradient-based optimization Numerous resear
hers have applied gradient-based optimization using the ad- joint approa
h in aerodynami
 optimization sin
e the early 1990's. Reuther and his 
o-workers applied the adjoint approa
h for aerodynami
 optimization of air
raft 
onguration using Euler's ow model [73℄. Jameson developed an adjoint formula- tion for the Navier-Stokes equations and applied it to transoni
 wing optimization [43℄; in this work, Jameson suggested that an optimal pressure distribution rst be obtained using Euler's ow model then used as a target pressure for an inverse design optimization problem using a Navier-Stokes adjoint optimizer to redu
e the overall 
omputational 
ost. However, Jameson treated the eddy vis
osity as 
onstant whi
h was later shown to be a bad assumption. Anderson and Bonhaus examined the ee
t of the strength of 
oupling of the turbulen
e model to the ow equations. They 
om- pared the adjoint gradient of the ow equations, with the eddy vis
osity frozen, with the nite dieren
e gradient of the 
ombined ow and turbulen
e equations, and found that freezing the eddy vis
osity 
an lead to signi
ant error in the 
omputed gradient. Therefore, they developed a ow solver that 
ouples the Spalart-Allmaras one-equation turbulen
e model with the ow equations; their 
oupled solver used 5× 5 blo
ks for two dimensional ows and 6× 6 blo
ks for three dimensional 
ases. They found that this 
oupling of the turbulen
e equations with the ow equations led 7 to an a

urate adjoint gradient [3℄. The same observation was veried by Nielsen and Kleb who extended their adjoint solver to deal with 
hemi
ally rea
ting ows [65℄. Zymaris et al developed a 
ontinuous adjoint optimizer for turbulent ow using the k − ε turbulen
e model and applied it to du
t optimization; they showed also that the assumption of 
onstant eddy vis
osity leads to great ina

ura
y in the 
omputed gradient and this leads to a poor sear
h dire
tion [88℄. Regarding the optimization te
hnique used, early resear
hers used the steepest des
ent s
heme, in whi
h the design variables are updated on a sear
h dire
tion exa
tly opposite to the gradient of the obje
tive fun
tion [66, 4℄. This s
heme has been implemented by Jameson and other resear
hers [37, 41, 73, 3℄, but as the steepest des
ent s
heme requires a large number of iterations to 
onverge to minimal point, the sequential quadrati
 programming (SQP) s
heme seems to be an attra
tive 
andidate as an optimization te
hnique. The use of SQP requires 
omputing the Hessian of the obje
tive fun
tion with respe
t to the design variables. The exa
t Hessian is expensive to 
ompute and may not be positive denite; therefore, the BroydenFlet
herGoldfarbShanno (BFGS) approximate formula is often adopted to approximate the Hessian using obje
tive fun
tion gradient history [12, 13, 87℄. The BFGS approximation always gives a positive denite approximate Hessian, therefore a real optimization sear
h dire
tion is guaranteed. Dadone et. al used BFGS in aerodynami
 optimization for transoni
 and supersoni
 wings and 
ompared BFGS and steepest des
ent; there results showed that the BFGS method is more e
ient in nding the optimal solution and is less sensitive to any ina

ura
y 
aused by approximation in gradient 
omputations [19℄. BFGS optimization has also been used by Neme
 and Zingg in subsoni
 and transoni
 turbulent aerodynami
 optimization of two dimensional airfoil [63℄. 1.2.2 Gradient free optimization Gradient free or gradient independent optimization methods, also known as heuris- ti
 optimization methods, are optimization te
hniques that do not require obje
tive fun
tion gradient 
omputation and therefore 
an be applied to non-dierentiable problems. They 
an be 
ategorized as evolutionary s
hemes (in
luding geneti
 algo- rithms) and random sear
h s
hemes (like the parti
le swarm te
hnique). 8 Geneti
 algorithm optimization is an evolutionary optimization algorithm that is inspired by Darwin's theory of evolution and natural sele
tion [30℄. The design variables are treated as 
hromosomes and optimization is 
arried out by 
rossing and mutating these 
hromosome to nd a better solution that minimizes the obje
tive fun
tion. The initial population is randomly generated, the obje
tive fun
tion value is 
omputed using CFD for ea
h population member and a tness value is 
omputed based on that; some members are sele
ted based on their tness to be the parents of the next generation and are used to generate the 
hromosomes of the ospring of the next optimization iteration. Transoni
 wing optimization using a geneti
 algorithm was explored by Gregg and Misegades [32℄ and by Gage and Kroo [29℄ in the late 1980's to minimize drag with lift 
onstraint. Somewhat later, Anderson applied a geneti
 algorithm in subsoni
 wing optimization with stru
tural 
onstraints [1℄; he added the geometri
 and aerodynami
 
onstraints to the obje
tive fun
tion as penalty terms. Jang and Lee applied a geneti
 algorithm in subsoni
 and transoni invis
id airfoil optimization; their obje
tive was to maximize the lift-to-drag ratio of an airfoil, beginning from the NACA 0012 1 [44℄. Oyama et al applied a geneti algorithm with a Navier-Stokes solver for transoni
 wing optimization [70℄. They also explored the use of fra
tal analysis in GA aerodynami
 optimization [71℄. The parti
le swarm method (PSOpt) is a sto
hasti
 optimization sear
h method developed by Eberhart and Kennedy in 1995, inspired by the so
ial behavior of bird o
ks [46, 21℄. The general idea of the parti
le swarm optimization is to randomly generate a swarm of parti
les in the design spa
e. For ea
h parti
le a tness value is 
al
ulated based on CFD simulation. Then parti
les y in the design spa
e a

ording to a simple formula that takes into a

ount that parti
le's own best t- ness position and the swarm's overall best tness position [22℄. PSOpt is known to suer from premature 
onvergen
e prior to dis
overing the true global minimizer; Evers suggests an automati
 regrouping PSO (RegPSO) that automati
ally triggers swarm regrouping when premature 
onvergen
e is dete
ted. The suggested regroup- ing strategy aims to liberate parti
les from sub-optimal solutions and enables nding the global minimum [24℄. Although the PSO algorithm has been applied to a wide 1 The National Advisory Committee for Aeronauti
s (NACA) airfoil family geometri
 
oordinates 
an be found at University of Illinois Urbana-Champaign website http://www.ae.illinois.edu/m- selig/ads/
oord_database.html. An explanation of the meaning of the digits in the NACA airfoil naming s
heme 
an be found in [36℄. 9 range of engineering problems in the literature, very few aerodynami
 optimization appli
ations are known. Venter and Sobiesz
zanski applied the parti
le swarm op- timization te
hnique in the multidis
iplinary optimization of a wing; the obje
tive was to maximize the air
raft range by maximizing lift-to-drag ratio and redu
ing the wing weight subje
t to geometri
 
onstraints [85℄. Chandrashekarappa and Du- vigneau [15℄ used parti
le swarm optimization s
heme to aerodynami
ally optimize wings in supersoni
 
onditions; Duvigneau also applied parti
le swarm optimization to aerodynami
 optimization of wings at transoni
 speed with free stream Ma
h number un
ertainty [20℄. 1.3 Contributions of the Thesis High order CFD methods 
an 
ompute an a

urate value of an aerodynami
 obje
tive fun
tion at lower 
omputational 
ost than required when a se
ond order method is used. The rst major 
ontribution of this thesis is a study of whether the e
ien
y of high order methods for CFD analysis translates into improved e
ien
y for gradient- based aerodynami
 optimization. Gradient based optimization te
hniques are known to be lo
al minimizers, with results depending on the starting airfoil geometry. The se
ond major 
ontribution of this thesis is the development and study of an optimization s
heme that 
an rea
h a true global optimum (for invis
id transoni
 aerodynami
s, a sho
k free airfoil subje
t to aerodynami
 and geometri
 
onstraints) after a reasonable number of CFD simulations. This s
heme is a hybrid (BFGS + regrouped parti
le swarm) s
heme that takes the advantages of gradient-based and gradient-free optimization s
hemes. To a
hieve these goals, the following 
omponents were needed, in addition to the pre-existing high-order a

urate ow solver: An e
ient geometry parametrization method. For optimization purposes, the airfoil shape must be represented by a nite number of design variables. Chap- ter 2 des
ribes the requirements on su
h a parametrization and presents a new least-squares spline parametrization method developed for this work. Robust mesh movement. As the airfoil shape 
hanges during optimization, the 
omputational mesh must be updated. Be
ause mesh regeneration will intro- 10 du
e una

eptably large 
hanges in the dis
retization error, mesh movement is strongly preferred. Chapter 3 des
ribes the semi-torsional spring mesh move- ment s
heme used in this thesis. Mesh sensitivity 
al
ulation. Cal
ulating the gradient of the obje
tive fun
tion a

urately and e
iently requires information about the movement of the mesh with 
hanges in design variables. Se
tion 3.3 des
ribes this pro
ess, whi
h 
om- bines aspe
ts of the geometri
 parametrization and mesh movement s
hemes. Obje
tive fun
tion gradient 
al
ulation. E
ient gradient 
al
ulation requires the solution of the adjoint to the governing ow equations. Chapter 4 des
ribes how this is done. A key innovation of the thesis is e
ient solution of the dis
rete adjoint problem for a high-order a

urate ow solution s
heme. In addition to showing the formulation of the gradient 
al
ulation, this 
hapter 
ompares nite dieren
e and adjoint gradients for subsoni
 and transoni
 ows for se
ond and fourth order nite volume s
hemes. For transoni
 ow, results are presented both with and without limiting of the 
omputational solution to prevent overshoots. Optimization drivers. The gradient-based optimizer using in this thesis is the BFGS-based quasi-Newton solver in Matlab's optimization toolbox. Gradient based optimization test 
ases for inverse design problems and for drag mini- mization with and without a lift 
onstraint is shown in Chapter 6. Also, the impa
t of mesh renement on the se
ond and fourth order optimal results is studied. The regrouped parti
le swarm optimization 
ode was written by the author. This s
heme and its hybridization with the BFGS s
heme is presented in Chap- ter 7. Examples are given to show the ee
tiveness of the hybrid s
heme in nding global optima when the gradient-based s
heme is unable to. 11 Chapter 2 GEOMETRY PARAMETRIZATION The geometry of engineered obje
ts is dened mathemati
ally in 
omputer-aided design (CAD) software, then exported as a group of points or polygons whi
h ap- proximates the original geometry. This dis
rete data provides input to mesh gener- ation software that 
reates a dis
rete representation of the 
omputational domain (a mesh), whi
h in turn is used as input for CFD analysis. Aerodynami
 design and optimization modies the aerodynami
 shape by 
hanging a set of geometri
 design variables. An obvious 
hoi
e is to use all the surfa
e grid points of a wing, but this approa
h 
auses two problems. First, it makes the design spa
e very large and this may lead to a highly expensive optimization. Se
ond, it may lead to a non-smooth geometry due to the displa
ement independen
e of a surfa
e mesh point; this prob- lem 
an be solved by the use of a smoothing fun
tion, as des
ribed by Jameson [41℄. To avoid these problems, most optimization approa
hes rely on some form of geo- metri
 parametrization. Geometri
 parametrization te
hniques 
an be 
lassied as analyti
al; pie
e-wise spline tting; CAD based; and free form deformation (FFD) approa
hes. In Se
tion 2.1, a review of various geometry parametrization te
hniques is presented. Se
tion 2.2 des
ribes in detail the geometry parametrization te
hnique used in this resear
h, a novel pie
e-wise least-squares tting te
hnique [5, 6℄. 2.1 Survey of Geometri
 Parametrization Te
hniques In this review, various te
hniques applied to aerodynami
 optimization are presented; the basi
s of these te
hniques and their appli
ation limitations are dis
ussed, in
lud- ing referen
e to the original papers des
ribing them in more detail . 12 2.1.1 Analyti
al parametrization The analyti
al approa
h was rst applied by Hi
ks and Henne for airfoil optimiza- tion [35℄. They suggest that weighted sinusoidal displa
ement shape fun
tions be added to the base geometry to modify the airfoil shape; the weights are the opti- mization design variables. The sinusoidal displa
ement shape fun
tion is expressed as h(x) = ( sin ( πx ln(0.5) ln(b) ))a 0 < x < 1 (2.1) where a and b are 
onstants to 
ontrol the peak lo
ation and the width of the sinusoidal displa
ement shape fun
tion. a = 4 is re
ommended for most 
ases, while b must be between 0 and 1 [45℄. The Class-Shape fun
tion-Transformation method (CST) is another analyti
al geometry parametrization te
hnique presented by Kulfan [49, 48℄. CST parametrizes the airfoil geometry using the following formula y = x0.5(1− x)S (x) + x△ZTE (2.2) where S (x) is the shape fun
tion and △ZTE is the nite thi
kness of the trailing edge. Kulfan and Bussoletti re
ommended using a weighted Bernstein binomial of order n as a shape fun
tion S (x) = n∑ i=0 bi n! i! (n− i)!x i (1− x)n−i (2.3) where the weights bi are used as the design variables. Although CST gives non-wavy proles, it is not 
apable of representing 
omplex geometries. Mousavi and Nadiragah 
ompared the impa
t of using dierent geometry parametrization te
hniques on the optimal wing shape; using CST parametrization gave drag 
oe
ients higher by about 15% than the optimal geometry from a B-spline parametrization [58℄ for a three-dimensional lift-
onstrained transoni
 drag minimization problem. 13 2.1.2 Pie
e-wise spline parametrization Bezier 
urves 
an be used for airfoil shape parametrization. Obayashi used Bezier 
urve for aerodynami
 optimization using geneti
 algorithm [67℄; he noti
ed that the Bezier 
urve representation fails to represent geometries that gives rooftop pres- sure distributions 2 be
ause Bezier 
urves are always 
onvex. B-splines, a general- ization of Bezier 
urves, were found to be more suitable. A 
ubi
 B-spline rep- resentation is a very good geometry parametrization te
hnique. Better 
ontrol of the 
ubi
 spline representation 
an be obtained by in
reasing the number of splines that represent the airfoil; Li et al. optimized the NACA 0012 at single point and multi-point operating 
onditions with a lift 
onstraint using spline representation as geometry parametrization te
hnique [50℄. CAD systems typi
ally use non-uniform rational B-spline (NURBS) representation for geometry modeling, allowing them to represent any 
omplex geometry; a detailed dis
ussion of NURBS 
an be found in The NURBS Book [72℄. Mengistu and Ghaly applied su

essfully a NURBS parametrization s
heme to turboma
hine blade aerodynami
 optimization using a gradient-free method [53℄. While pie
e-wise spline parametrization is well-suited for two-dimensional shapes and simple three-dimensional geometries, 
omplex shapes require a large number of 
ontrol points, redu
ing the ee
tiveness of gradient-based optimization te
hniques [76℄. 2.1.3 CAD parametrization Computer aided design pa
kages have evolved to implement NURBS for geometry representation due the ex
ellent properties of NURBS. Linking the CAD and grid generation software 
an be done using an API that allows a

ess to the CAD system's internal interfa
e [81, 8℄. However, imposing geometri
 
onstraints is still an obsta- 
le. Mesh sensitivity 
al
ulation is another obsta
le for gradient-based optimization te
hniques based of CAD parametrization: analyti
al mesh sensitivity with respe
t to geometry design variables  the NURBS knots  
an in prin
iple be 
omputed with the use of automati
 dierentiation of the CAD software but this is not possible without CAD sour
e 
ode and is unlikely to be pra
ti
al even then, given the size of 2 Distributions for whi
h pressure remains almost un
hanged over a signi
ant 
hordwise dis- tan
e. 14 the CAD 
ode base. This derivative 
an also be 
omputed using nite dieren
es, but the risk of poor a

ura
y still exists, and 
omputational 
osts are higher [76, 77℄. 2.1.4 Free form deformation (FFD) Computer graphi
s requires large graphi
al deformations su
h as stret
hing, twisting and other surfa
e morphing operations; soft obje
t animation (SOA) algorithms were developed to help with the geometry morphing required in graphi
s animation [86℄. In SOA algorithms, the obje
t surfa
e is treated as a pie
e of rubber and the desired deformation 
an be obtained by applying loads on it, so geometry morphing 
an be obtained without a 
hange in surfa
e topology. The surfa
e itself 
an be parametrized using Bezier or B-splines or even NURBS splines. A related approa
h, the free form deformation algorithm (FFD) treats the geometry as a void in a box-shaped pie
e of rubber. Deformation 
an be 
ontrolled by moving 
ontrol points pla
ed on the outer surfa
e of the box; the interior of the rubber box with its void is parametrized using a tensor produ
t of three spline representations (one in ea
h 
oordinate dire
tion). Sederberg and Parry developed an algorithm that uses the FFD 
on
ept with Bezier tri-variate volume representation [78℄. A disadvantage of the FFD method is that it requires large numbers of 
ontrol points to obtain lo
al deformation in the deformed geometry. However, Borrel and Rappoport presented a method to allow lo
al shape deformation via FFD by introdu
ing a set of 
ontrol points with 
onstrained lo
al B-splines that 
an be used to obtain deformation in a radius of inuen
e determined by the designer [10℄. 2.1.5 Multidis
iplinary aerodynami
/stru
tural shape optimization using deformations (MASSOUD) The MASSOUDmethod is an analogy to analyti
al methods that tries to parametrize the deformations in the geometry rather than the geometry itself. It also uti- lizes SOA algorithms and allows strong lo
al deformation 
ontrol. The MASSOUD parametrization requires few design variables be
ause it parametrizes the deforma- tion. Samareh has applied the MASSOUD method to parametrize a simple wing, a wing body blend, and a 
omplex air
raft 
onguration with su

ess [76℄. Nielsen and Anderson su

essfully adopted the MASSOUD parametrization s
heme for aerody- 15 nami
 optimization of turbulent ow using unstru
tured grids [64℄. 2.2 New Least Squares Parametrization Te
hnique A new least-squares based surfa
e parametrization method is presented in this se
- tion; the airfoil surfa
e is parametrized using pie
e-wise polynomials whose 
oe- 
ients are found by solving a least squares problem [6℄. This se
tion also des
ribes how to implement a thi
kness 
onstraint with the new parametrization and presents some validation test 
ases. 2.2.1 Airfoil geometry parametrization In the proposed te
hnique, the geometry is parametrized using pie
e-wise polyno- mials found by a least-squares t. The parametrization polynomials are 
ontrolled by a set of 
ontrol points and satisfy C2 
ontinuity at their meeting points. The airfoil upper and lower surfa
es are represented using two least square splines for ea
h surfa
e as shown in Fig. 2.1. The polynomials 3 used are P1(x) = a0 √ x+ a1x+ a2x 2 + a3x 3 0 < x < L1 P2(x) = b0 + b1x+ b2x 2 + b3x 3 L1 < x < L (2.4) where x is the normalized 
hord-wise position, L1 is 
hord-wise position that sep- arates the polynomial regions, and L = 1. These polynomials are suitable for an airfoil with a rounded leading edge due to the existen
e of the √ x term in P1 (x), whi
h gives an innite slope at x = 0. The x and y 
oordinates of the design 
ontrol points shown in Fig. 2.1 are used to nd the values of the polynomial 
oe
ients. These polynomials must satisfy 
ontinuity of value, slope, and 
urvature at their meeting point x = L1. These 
onditions 
an be written as ˆ Value 
ontinuity: P1(L1)− P2(L1) = 0 (2.5) 3 Te
hni
ally, P1 is not a polynomial be
ause of the presen
e of the √ x term used to give innite slope and nite radius of 
urvature at x = 0. However, the label is 
onvenient and not overly 
onfusing. 16 ˆ Slope 
ontinuity: P ′ 1(L1)− P ′ 2(L1) = 0 (2.6) ˆ Curvature 
ontinuity: P ′′ 1 (L1)− P ′′ 2 (L1) = 0 (2.7) An additional 
onstraint should be added on P2(L) to assure zero thi
kness at the trailing edge. The above hard 
onstraints should be stri
tly satised by the geom- etry parametrization polynomials; they 
an be written in matrix form as BP = 0 where B =  √ L1 L1 L 2 1 L 3 1 −1 −L1 −L21 −L31 1 2 √ L1 1 2L1 3L 2 1 0 −1 −2L1 −3L21 −1 4 √ L31 0 2 6L1 0 0 −2 −6L1 0 0 0 0 1 L L2 L3  , P =  a0 a1 a2 a3 b0 b1 b2 b3  The free parameters are 
hosen to best approximate the y 
oordinates of the airfoil shape 
ontrol points; the x 
oordinates of the 
ontrol points are xed. The resulting least square system with 
onstraints applied 
an be expressed as min ‖AP − c‖2 subje
t to BP = 0 (2.8) where A 
ontains powers of the x 
oordinates at the design points so that AP gives the y 
oordinates of the parametrized shape at the design points and c 
ontain the a
tual y 
oordinates of the design points. I have used a set of 
ontrol points that 
ontrols the airfoil shape fun
tions instead of using the 
oe
ients of the shape fun
tions to ease 
onstraints and boundary denitions. To give an example of how this least-squares system is 
onstru
ted, 
onsider the parametrized airfoil surfa
e 17 Figure 2.1: Least square surfa
e presentation of RAE 2822 airfoil using two polyno- mials P1 (x) & P2 (x) tted using nine 
ontrol points. shown in Fig. 2.1. Six 
ontrol points lies in the region of the polynomial P1 (x), while three points lies in the P2 (x) region. The 
orresponding least-squares system is  √ x1 x1 x 2 1 x 3 1 0 0 0 0√ x2 x2 x 2 2 x 3 2 0 0 0 0√ x3 x3 x 2 3 x 3 3 0 0 0 0√ x4 x4 x 2 4 x 3 4 0 0 0 0√ x5 x5 x 2 5 x 3 5 0 0 0 0√ x6 x6 x 2 6 x 3 6 0 0 0 0 0 0 0 0 1 x7 x 2 7 x 3 7 0 0 0 0 1 x8 x 2 8 x 3 8 0 0 0 0 1 x9 x 2 9 x 3 9   a0 a1 a2 a3 b0 b1 b2 b3  =  y1 y2 y3 y4 y5 y6 y7 y8 y9  (2.9) Be
ause the 
onstraint equation, Eq. 2.8, has a zero right hand side, the solution ve
tor P must lie in the null spa
e of the 
onstraint equations, i.e. it should be a linear 
ombination of the null spa
e basis of the 
onstraint equations. The matrix B is full row rank and to nd the null spa
e basis of it, QR fa
torization 
an be used: BT = QR 18 Q = [ −→q1 −→q2 −→q3 −→q4 | −→q5 −→q6 −→q7 −→q8 ] = [ Q1 | Q2 ] The ve
tors of the [Q2] matrix are unit ve
tors forming a basis of the null spa
e of matrix B. The solution ve
tor P must be a linear 
ombination of the ve
tors of the matrix Q2 P = z1 · −→q5 + z2 · −→q6 + z3 · −→q7 + z4 · −→q8 = Q2z (2.10) Substituting Eq. 2.10 into Eq. 2.8 AP = AQ2z = c (2.11) Solving the least-squares system by singular value de
omposition for numeri
al sta- bility, z = [AQ2] † c (2.12) where the re
tangular matrix [AQ2] † is the pseudo-inverse of [AQ2]. Finally, P = Q2 [AQ2] † c (2.13) Equation 2.13 gives the relationship between the polynomial 
oe
ients P and the y lo
ations of the 
ontrol points c. The sensitivity of the polynomial 
oe
ients with respe
t to the y lo
ation of the ith 
ontrol point, whi
h is needed to 
al
ulate the mesh sensitivity ∂M/∂D, is the ith 
olumn of the matrix Q2 [AQ2] † .The dependen
y of an airfoil surfa
e point −→rb = [ xa ya ] on a design variable Di 
an be found from ∂−→rb ∂Di = [ 0 ∂a0 ∂Di √ xa + ∂a1 ∂Di xa + ∂a2 ∂Di x2a + ∂a3 ∂Di x3a ] 0 < xa < L1 (2.14) ∂−→rb ∂Di = [ 0 ∂b0 ∂Di + ∂b1 ∂Di xa + ∂b2 ∂Di x2a + ∂b3 ∂Di x3a ] L1 < xa < L 19 where [ ∂a0 ∂Di ∂a1 ∂Di ∂a2 ∂Di ∂a3 ∂Di ∂b0 ∂Di ∂b1 ∂Di ∂b2 ∂Di ∂b3 ∂Di ]T is the ith ve
tor of the 
on- strained pseudo inverse matrix Q2 [AQ2] † . This pro
edure 
an be extended to parametrize airfoil surfa
e using any number of pie
e-wise polynomials, with an a
- 
ompanying in
rease in system size. Leading edge radius and trailing edge thi
kness 
onstraints 
an be added to the C2 
ontinuity requirements forming the 
onstraint system BP = d, and the 
ontrol point 
oordinates are used to 
onstru
t the system AP = {yc}. In pra
ti
e, all airfoil surfa
e points are used in the least squares system to nd the polynomials 
oe
ients. After nding the polynomials 
oe
ients, the designer 
an sele
t a set of 
ontrol points whi
h lies on the parametrized surfa
e; the least number of 
ontrol points is two per polynomial. The less 
ontrol points used, the less 
ontrol of airfoil geometry. I have sele
ted nine 
ontrol points at spe
i
 
hord- wise x-stations based on my engineering sense; it turns out that my sele
ted set of 
ontrol points were able to produ
e sho
k free optimal proles as will be shown in the next 
hapters, however, sele
tion of the 
ontrol points x-stations 
an be made by formulating a simple minimization problem. In this problem, a set of airfoil geometri
 data gathered and the obje
tive is to minimize the RMS error between the original airfoil surfa
es and the parametrized airfoil surfa
es, where the design variables are the 
ontrol points x-stations. 2.2.2 Thi
kness 
onstraint Some air
raft fuel tanks are pla
ed inside the wing, and the main landing gear of some air
raft are stored in the wings after being retra
ted. In addition, the wing must have su
ient bending rigidity from a stru
tural point of view; therefore, the wing thi
kness is a 
onstraint at some 
hord-wise stations. This subse
tion demonstrates how to add a thi
kness 
onstraint to the parametrization. Suppose the airfoil is parametrized using four polynomials, as shown in Fig. 2.2: two for the upper surfa
e P1(x) = a0 √ x+ a1x+ a2x 2 + a3x 3 0 < x < L1 P2(x) = b0 + b1x+ b2x 2 + b3x 3 L1 < x < L 20 Figure 2.2: Parametrized airfoil using four polynomials and two for the lower surfa
e P3(x) = c0 √ x+ c1x+ c2x 2 + c3x 3 0 < x < L1 P4(x) = d0 + d1x+ d2x 2 + d3x 3 L1 < x < L If the thi
kness tc need to be 
onstrained at a station xc where 0 < xc ≤ L1 , this 
onstraint will be Ct : P1(xc)− P3(xc) = tc (2.15) Ct : ( a0 √ xc + a1xc + a2x 2 c + a3x 3 c )− (c0√xc + c1xc + c2x2c + c3x3c) = tc Equation 2.15 provides a link between the upper and lower surfa
e polynomials; therefore, a 
oupled least-squares system needs to be 
onstru
ted and solved; the 21 hard 
onstraint equations are expressed as  B 00 B P1(xc) {0}1×4 −P3(xc) {0}1×4   a0 a1 . . . d2 d3  =  0 0 . . . 0 tc  (2.16) [Ba] {Pa} = {da} (2.17) The global least-squares system is[ A 0 0 A ] {Pa} = { cu cl } (2.18) [Aa] {Pa} = {ca} The above hard 
onstraint equations do not have a null right hand side due to the thi
kness 
onstraint; therefore, the solution (polynomials 
oe
ients) does not belong to the null spa
e of Ba. The solution pro
edure for the 
onstrained least- squares problem expressed by Eq. (2.16,2.18) is des
ribed by Masuda et al [52℄: ˆ Apply QR fa
torization to get BTa = QR. ˆ Let Q = [ Q1 Q2 ] , R = [ R1 0 ] where Q1
ontains the rst 7 
olumns of Q and Q2 
ontains the rest of the 
olumns, and R1is the rst 7× 7 sub matrix of R. ˆ The hard 
onstraints 
an now be written as BPa ≡ RTQTPa = da ˆ Let QTPa = [ y z ] . Then Pa = Q [ y z ] = Q1y +Q2z ∴ RTQTPa = [ RT1 0 ] · [ Q1 Q2 ]T Pa = da 7−→ RT1 y + 0 = da ∴ y = R−T1 da (2.19) 22 ˆ Be
ause R is upper triangular, Equation 2.19 
an be solved using forward substitution. ˆ The soft 
onstraints in equation 2.18 
an be rewritten as AaPa = Aa [ Q1y Q2z ] = ca (2.20) ∴ AaQ2z = ca −AaQ1y (2.21) ∴ z = [AaQ2] † [ca −AaQ1y] (2.22) ˆ Finally, the polynomial 
oe
ients 
an be found as Pa = Q1y +Q2z: Pa = [ Q1 −Q2 [AaQ2]†AaQ1 ] R−T1 da +Q2 [AaQ2] † ca (2.23) Using the last equation we noti
e that the sensitivity of the polynomial 
oe- 
ients with respe
t to the ith 
ontrol point y lo
ation, yci, is the i th ve
tor of the matrix [Q2] [AQ2] † . 2.2.3 Validation 
ases In this subse
tion, testing of the proposed least-squares surfa
e parametrization is 
arried out; parametrization of various types of airfoils is done to show that the proposed geometry parametrization s
heme has the required exibility to represent various types of airfoils. These airfoils in
lude NACA 4-, 5-, and 6-digit series, laminar ow, supersoni
, and super-
riti
al airfoil se
tions. 4 Figures 2.32.10 show the least-squares tting polynomials for various airfoils, while Table 2.1 shows the RMS error in the parametrized geometry for upper and lower surfa
es. The last parametrized airfoil, the laminar LV2 5 airfoil of German Aerospa
e Center, was parametrized using 10 polynomials, 5 for ea
h of its surfa
es. Although the RMS error is small (order of 10−5
hord for all parametrization method), as shown in Table 2.1, Fig. 2.11 shows that the dieren
e of the surfa
e pressure distribution, espe
ially at the peak velo
ity at the leading edge upper surfa
e, is still signi
ant. There are many x-stations at whi
h LV2 airfoil 
hanges its 
urvature, and also a 4 Readers interested in parti
ulars for the NACA airfoils are referred to Abbott and von Doen- ho's textbook [36℄. The seminal referen
e for the RAE airfoil is Cook et al [18℄. 5 LV2 geometry obtained by personal 
onta
t with German Aerospa
e Center (DLR) resear
hers. 23 Figure 2.3: Parametrized NACA 0011SC, 'O' are the original airfoil ordinates Figure 2.4: Parametrized NACA 0012, 'O' are the original airfoil ordinates 24 Figure 2.5: Parametrized NACA 6509, 'O' are the original airfoil ordinates Figure 2.6: Parametrized NACA 16006, 'O' are the original airfoil ordinates 25 Figure 2.7: Parametrized NACA 63412, 'O' are the original airfoil ordinates Figure 2.8: Parametrized NACA 644421, 'O' are the original airfoil ordinates 26 Figure 2.9: Parametrized RAE2822, 'O' are the original airfoil ordinates long interval of almost innite 
urvature value whi
h makes it hard for the geometry parametrization s
heme to a

urately present it and this 
auses the u
tuation in pressure resulted in the parametrized airfoil be
ause of la
k of a

urate presentation of 
urvature u
tuations. Better mat
hing 
an be obtained but this will in
rease the number of geometry design variables signi
antly. This 
ase illustrates 
learly the trade o between a

urately representing the geometry (whi
h will 
hange during optimization iteration) and 
hoosing a reasonable number of design variables; this 
hoi
e must be left for the designer. If the 
hange in pressure distribution be
omes una

eptable, or if the pressure distribution requires large number of airfoil 
ontrol points (design variables) to be a

urately represented, parametrization the geometry perturbation is the more attra
tive option. The values of RMS error in Table 2.1 are small 
ompared to the maximum airfoil thi
kness value: only the NACA 0011SC has an RMS error of more than 0.2% of maximum thi
kness. However, if the RMS error was large, in
reasing the number of surfa
e parametrization polynomials would redu
e the error. 27 Figure 2.10: Parametrized LV2 laminar airfoil, parametrized using 5 polynomials per surfa
e, 'O' are the original airfoil ordinates Figure 2.11: LV2 airfoil parametrized using dierent methods with 20 design vari- ables, presented with permission of Brezillon. 28 Airfoil Upper surfa
e error Lower surfa
e error Maximum Thi
kness NACA 0011SC 6.37 · 10−4 6.37 · 10−4 11 · 10−2 NACA 0012 2.35 · 10−4 2.35 · 10−4 12 · 10−2 NACA 6409 9.28 · 10−5 1.24 · 10−4 9 · 10−2 NACA 16-006 1.31 · 10−6 2.03 · 10−5 6 · 10−2 NACA 63-412 1.17 · 10−5 5.63 · 10−5 12 · 10−2 NACA 64-421 1.70 · 10−4 7.50 · 10−5 21 · 10−2 RAE 2822 1.20 · 10−5 3.20 · 10−5 slightly > 12 · 10−2 LV2 1.45 · 10−5 3.22 · 10−5 slightly > 12 · 10−2 Table 2.1: RMS error in dierent parametrized airfoil geometries 29 Chapter 3 MESH MORPHING AND MESH SENSITIVITY 3.1 Mesh Morphing The modi
ation of the aerodynami
 shape during optimization requires a 
hange of the mesh that presents the shape. This 
an be done by grid regeneration around the new geometry but this is time 
onsuming and will 
hange the dis
retization error [51℄. Another strategy is to adapt the old mesh to t the new shape of the airfoil using mesh movement. The tension spring analogy is one of the most widely used mesh deformation strategies for aerodynami
 optimization. The main idea is to repla
e the grid edges by springs with stiness inversely proportional to their length. The boundary points that lie on the airfoil surfa
e are moved with displa
ement spe
ied by the optimizer, far eld points are kept xed, and the interior point displa
ements are determined by equilibrium of the spring network [7℄. For large grid displa
ements, the linear spring analogy is not robust, and negative area 
ells 
an present after mesh morphing. Farhat et. al. improved the tension spring analogy by adding torsional springs at the grid nodes to prevent element ipping [25℄. Ea
h edge fa
es two angles as shown in Fig. 3.1; the edge stiness is modied to in
lude terms with the re
ipro
al of the sine of these angles. This allows the edge stiness to grow to innity if the angle tends to be zero and therefore prevent element ipping. Another strategy is to modify the mesh by solving a linear elasti
ity problem in whi
h the boundary displa
ements are known [14℄; the element modulus of elasti
ity 
an be the re
ipro
al of the distan
e from the wall, or it 
an be the re
ipro
al of the element size [82℄. The later strategy was applied by Stein et al. and the results showed that this method is robust, espe
ially for vis
ous 
al
ulations. In this 
ase, 30 Figure 3.1: S
hemati
 drawing of an edge pq and its fa
ing angles a and b. the elements of the boundary layer experien
ed small geometri
al 
hange while the elements away from the airfoil experien
ed larger 
hanges [80℄. However, the linear elasti
ity mesh movement s
heme is 
omputationally expensive 
ompared to the spring analogy method. The semi-torsional spring analogy method is adopted in this resear
h due to its simpli
ity and robustness. Consider the edge pq shown in Fig 3.1. The relationship between the for
es and the node displa
ements when treating this edge as a spring follows Hooke's Law as  Fpx Fpy Fqx Fqy  = 1lpq ( 1 + 1 sin (θa) + 1 sin (θb) ) −1 0 1 0 0 −1 0 1 1 0 −1 0 0 1 0 −1   upx upy uqx uqy  (3.1) where lpq is the length of edge pq, and θa and θb are the angles fa
ing the edge. After the assembly of the global stiness matrix, the system of equations that relates grid point displa
ement with nodal for
es 
an be written as [ Kii Kib Kbi Kbb ]{ Ui Ub } = { 0 Fb } (3.2) 31 where Ui and Ub are the interior and boundary mesh point displa
ements respe
tively. We do not need to know the values of boundary nodal for
es Fb be
ause the boundary points displa
ement ve
tor Ub is known expli
itly: it is the deformation required in the airfoil prole to minimize the obje
tive fun
tion. Therefore, Eq 3.2 
an be written as [ Kii Kib 0 I ]{ Ui Ub } = { 0 Ub } (3.3) Substituting Uki = −−→ rk+1i − −→ rki , Ub = −→rb −−→rboin Eq 3.3 we get[ Kii Kib 0 I ]k{ −−→ rk+1i−→rb } = [ Kii Kib 0 I ]k{ −→ rki−→rbo } + { 0 Ub } (3.4) where −→rbo is the initial position ve
tor of the boundary points before mesh morphing. The stiness matri
es [Kii] and [Kib] depend on the mesh fa
e lengths, whi
h 
hange during the mesh morphing stage; therefore, Eq. 3.4 is a non-linear equation and needs to be solved iteratively. 3.2 Testing Mesh Morphing Farhat et al [25℄ demonstrated the robustness of the semi-torsional mesh movement s
heme. In this se
tion, testing results are presented to demonstrate that the 
ur- rent implementation of this s
heme shows the same good behavior for several high- deformation 
ases. The rst test 
ase is an unstru
tured triangular mesh around a NACA 0012. The thi
kness of the airfoil is doubled, whi
h means that the airfoil surfa
e points will translate by several 
omputational 
ells in the y-dire
tion. Fig- ure 3.2 shows that even with this large displa
ement (multiple 
ells) of the boundary, no 
ells are inverted and no edges interse
ted with another. Displa
ement de
reases with the distan
e from the airfoil surfa
e, and there is almost no movement on the symmetry line. The se
ond 
ase tests mesh movement when the outer boundary is 
hanging. The original mesh is an unstru
tured triangular mesh (shown in Fig 3.3a). The outer boundaries are 
hanged, redu
ing the total area by almost 50% and turning the original right-angled 
orners nearly into 
usps. Figure 3.3b shows that the semi- 32 Figure 3.2: Mesh movement s
heme results of doubling the thi
kness of NACA 0012 torsional mesh movement was 
apable of adapting the mesh in the entire eld without element ipping. 3.3 Mesh Sensitivity For gradient-based optimization, as shown in the next 
hapter, the gradient 
ompu- tation requires 
al
ulation of the dependen
y of the residual on the design variables ∂R/∂D, whi
h in turn depends on the evaluation of mesh sensitivity. The mesh sensitivity tells how mesh points translate in the (x, y) plane with the perturbation of the geometri
 parametrization 
ontrol points (design variables). This translation, obviously, depends on the mesh movement s
heme. Truong et. al. 
ompared al- gebrai
 mesh movement to linear elasti
ity mesh movement s
hemes to study the impa
t of the adopted mesh movement s
heme on the nal optimized shape. There was a noti
eable dieren
e between the nal optimal airfoil shapes for a subsoni 
ase, although the dieren
e in the optimal obje
tive fun
tion value was of order of dis
retization error. For a transoni
 test 
ase, the dieren
e in the optimal shapes was almost negligible [83℄. The mesh sensitivity with respe
t to one of the design variables ∂M/∂Di (that is, 33 (a) Initial mesh (b) Final mesh Figure 3.3: Mesh movement results, large outer boundary deformation of a re
tan- gular domain the 
hange in mesh point lo
ations with a 
hange of a design variable) is 
al
ulated by dierentiating Eq. 3.3:[ Kii Kib 0 I ]{ ∂Ui ∂Di ∂Ub ∂Di } = { 0 ∂−→rb ∂Di } ∴ { ∂M ∂Di } = { ∂Ui ∂Di ∂Ub ∂Di } = [ Kii Kib 0 I ]−1{ 0 ∂−→rb ∂Di } (3.5) where ∂−→rb/∂Di is obtained using Eq. 2.14 and is related to the design variables via the pseudo inverse of the 
onstrained least-squares system solved in parametrizing the geometry. 34 Chapter 4 GRADIENT CALCULATION USING ADJOINT APPROACH Gradient 
al
ulation plays a key role in gradient-based optimization. The traditional nite dieren
e strategy is 
omputationally expensive as it requires at least as many CFD simulations as the number of design variables to 
ompute the gradient; ea
h of these is the solution to a large non-linear system of equations. The forward strategy to 
ompute the gradient is less expensive as it requires the 
al
ulation of the ow sensitivity with respe
t to the design variables, and uses the 
omputed ow sensitivity to 
ompute the gradient; the forward strategy requires solving a number of linear systems equal to the number of design variables to nd the ow property sensitivity with respe
t to all the design variables. The adjoint strategy is 
omputationally 
heaper; it requires the solution of one linear system whose right hand side is the dependen
y of the aerodynami
 obje
tive fun
tion on the ow eld properties. Due to its numeri
al e
ien
y and the 
orresponding redu
tion in 
omputational eort, the adjoint strategy is adopted in this resear
h. 4.1 Forward and Adjoint Formulations The obje
tive fun
tion, F , for aerodynami
 optimization is a fun
tion of the design variables, D, and the ow eld solution at the surfa
e points of the boundary 
ontrol volumes Us F = F (Us,D) (4.1) Us is expressed most 
onveniently in primitive variables: U = ( ρ u v P )T . Consider, for instan
e, the lift and drag for
es of a two dimensional airfoil, whi
h are perpendi
ular and parallel, respe
tively, to the in
oming ow, whi
h is in
lined 35 at an angle α to the airfoil 
hord. These 
an be evaluated as follows: FL = − {˛ Psnx ds } sinα+ {˛ Psny ds } cosα FD = {˛ Psnx ds } cosα+ {˛ Psny ds } sinα or in dis
rete form, FL = − ∑ wsPsnx sinα+ ∑ wsPsny cosα (4.2) FD = ∑ wsPsnx cosα+ ∑ wsPsny sinα (4.3) where Ps is the pressure at surfa
e integration point pressure, nx and ny are the unit normal 
omponents at the surfa
e integration point, α is the angle of atta
k, and ws is the ar
 length asso
iated with the surfa
e integration point. Note that this form uses the dimensional pressure and 
oordinates gives the dimensional lift and drag for
es. The lift and drag 
oe
ients, whi
h are the non-dimensional equivalents, are identi
al in form but use non-dimensional pressure and 
oordinates. These dis
rete integrals expressed as a fun
tion of geometri
 and ow properties of the 
ontrol volume su
h as the length of ea
h fa
e and the unit normal at ea
h Gauss point. The geometri
 properties depend on the design variables through the mesh sensitivity, while the ow properties at the Gauss points depend on the ow properties of the 
ontrol volume itself and its neighbors, whi
h in turns depend on the mesh and the boundary shape whi
h ultimately depends on the geometri
 design variables. The gradient of the obje
tive fun
tion 
an be obtained by using the 
hain rule dF dD = ∂F ∂Ubg ∂Ubg ∂U ∂U ∂M ∂M ∂D + ∂F ∂M ∂M ∂D (4.4) where Ubg is the boundary Gauss point ow properties, U is the 
ontrol volume averaged solution written in primitive variables, and M is the mesh point lo
ations. ∂M/∂D is the mesh dependen
y on the design variables whi
h 
omputations are presented in details in Chapter 3. The residual of the ow governing equations 
an be written as a fun
tion of the ow eld solution U and mesh geometri
 design variables D. If we apply the 
onstraint that the ow solution is 
onverged regardless 36 of variations in the design variables, we 
an write dR dD = ∂R ∂M ∂M ∂D + ∂R ∂U ∂U ∂M ∂M ∂D = 0 (4.5) Solving this for the solution sensitivity ∂U ∂D ≡ ∂U ∂M ∂M ∂D , we get ∂U ∂M ∂M ∂D = − [ ∂R ∂U ]−1 · ∂R ∂M ∂M ∂D = − [ ∂R ∂Q ∂Q ∂U ]−1 · { ∂R ∂M ∂M ∂D } (4.6) where the last equality expands the residual Ja
obian with respe
t to the primitive ow variables into the produ
t of the residual Ja
obian with respe
t to the 
onserved ow variables (whi
h is used in impli
it ow solvers) and a 
hange of variables for 
onserved to non-
onserved. Note that ∂Q ∂U is a blo
k diagonal matrix. Substituting Eq.4.6 in Eq. 4.4, we get the forward formulation for gradient 
omputation dF dD = − ∂F ∂Ubg ∂Ubg ∂U [ ∂R ∂Q ∂Q ∂U ]−1 · { ∂R ∂M ∂M ∂D } + ∂F ∂M ∂M ∂D (4.7) This form of the gradient requires solving as many linear systems as there are design variables in the optimization problem. Taking the transpose of Eq. 4.7, dF dD T = − { ∂R ∂M ∂M ∂D }T · [ ∂R ∂Q ∂Q ∂U ]−T { ∂F ∂Ubg ∂Ubg ∂U }T + { ∂F ∂M ∂M ∂D }T (4.8) where the residual sensitivity to mesh movement 
an be written using the 
hain rule as: ∂R ∂M = ∂R ∂AΩi ∂AΩi ∂M + ∑ fa
esΩi { ∂R ∂nx dnx dM + ∂R ∂ny dny dM + (4.9) ∂R ∂wi dwi dM + ∂R ∂U fa
e ∂U fa
e ∂M } We get the adjoint method, presented by A. Jameson [37, 41, 43, 38, 40℄. Now only one linear system solve is required. However, this linear system solve requires 37 expli
itly forming the global Ja
obian matrix ∂R ∂Q be
ause its transpose is required. Be
ause the CFD solver used in this thesis 
an form the global Ja
obian matrix expli
itly (author?) [57℄, the transpose of the Ja
obian 
an easily be formed as well. The Ja
obian from the last GMRES iteration is re-used, so the 
omputational eort of solving the adjoint problem is redu
ed to solving this linear system, a 
ost on the order of 1% of the CFD simulation 
omputational eort. To ease the programming eort when 
hanging the obje
tive fun
tion, we form three adjoint problems, one ea
h for the lift 
oe
ient CL, the drag 
oe
ient CD, and the moment 
oe
ient CM to nd ∂CL/∂D, ∂CD/∂D and ∂CM/∂D, respe
- tively. These aerodynami
 
oe
ient gradients 
an be used to evaluate any obje
tive fun
tion gradient for a fun
tion that depends on aerodynami
 for
es: dF dD = f ( dCL dD , dCD dD , dCM dD ) The solution pro
edure of the three adjoint problems 
an be summarized as follows, ˆ Using the steady state ow solution, we 
onstru
t the CFD simulation Ja
obian matrix ∂R/∂Q expressed in Eq. 1.6. ˆ We 
onstru
t ∂R/∂U as: ∂R ∂U = ∂R ∂Q ∂Q ∂U where ∂Q/∂U is the transformation matrix from 
onservative to primitive ow variables and is blo
k diagonal. ˆ We 
onstru
t ∂F ∂Ubg ∂Ubg ∂U , where ∂F/∂Ubg is the analyti
 dependen
y of the ob- je
tive fun
tion on the primitive ow properties at the airfoil surfa
e points, and ∂Ubg/∂U is the dependen
y of surfa
e point primitive ow properties on the 
ontrol volume average values of the primitive ow properties. The latter is known as a side ee
t of solution re
onstru
tion. ˆ We solve the three linear systems [ ∂R ∂U ]T ΨL,D,M = [ ∂CL,D,M ∂Ubg ∂Ubg ∂U ]T to get the adjoint ve
tors ΨL,D,M . ˆ We 
onstru
t the obje
tive fun
tion gradient ve
tor dCL,D,M dD T = { ∂F ∂M ∂M ∂D }T −{ ∂R ∂M ∂M ∂D }T ΨL,D,M , whi
h requires only a ve
tor dot produ
t for ea
h design 38 variable. As an example of how to use dCL,D,M dD to 
onstru
t the gradient of an aerodynami fun
tion, 
onsider the following obje
tive fun
tion whi
h represents a typi
al drag minimization fun
tion with a lift 
onstraint applied using a penalty term F = CD + k1 (CL − CLc)2 (4.10) The gradient of the above fun
tion 
an be written as dF dD = dCD dD + 2k1 (CL − CLc) dCL dD (4.11) also ∂F ∂D = ∂CD ∂D + 2k1 (CL − CLc) ∂CL ∂D (4.12) The dis
rete (invis
id) forms of CL and CD follow from the non-dimensionalization of Eqs. 4.2 and 4.3, where now the pressure and 
oordinates are non-dimensionalized: CL = − (∑ wsPsnx ) sinα+ (∑ wsPsny ) cosα (4.13) CD = (∑ wsPsnx ) cosα+ (∑ wsPsny ) sinα and their partial derivatives with respe
t to the geometry design variables are ∂CL ∂D = − (∑ ∂ws ∂D Psnx ) sinα+ (∑ ∂ws ∂D Psny ) cosα (4.14) − (∑ wsPs ∂nx ∂D ) sinα+ (∑ wsPs ∂ny ∂D ) cosα ∂CD ∂D = (∑ ∂ws ∂D Psnx ) cosα+ (∑ ∂ws ∂D Psny ) sinα+(∑ wsPs ∂nx ∂D ) cosα+ (∑ wsPs ∂ny ∂D ) sinα Using Eq. 4.9 to 
al
ulate ∂R/∂M , we need to 
ompute the terms ∂AΩi/∂M , ∂nx/∂M , ∂ny/∂M , and ∂wi/∂M , whi
h only depend on mesh points' spatial lo
a- tions, while the terms ∂R/∂nx and ∂R/∂ny are obtained by dire
t dierentiation of the Roe ux, and nally (∂R/∂Uface) · (∂Uface/∂M) is obtained by dierentiating 39 Figure 4.1: S
hemati
 drawing of an element fa
e κ, and illustration of its left and right sides the fa
e property re
onstru
tion s
heme whi
h depends on CFD solver te
hnology. The remainder of this se
tion will fo
us on the derivatives of the residual R, with the derivatives on the geometri
 terms in the following se
tion. Let us 
onsider the residual 
ontribution of fa
e κ of the 
ontrol volume Ωi shown in Fig. 4.1. The dis
rete form of this edge's residual 
ontribution to the 
ontrol volume Ωi is Rκ,Ωi = 1 AΩi J∑ j=1 { 1 2 ( ~FRj + ~FLj − |△F1|j − |△F4|j − |△F5|j ) wj } (4.15) where J = 1 when using one Gauss integration point at the middle of fa
eκ for a se
ond order a

urate s
heme. For the fourth order s
heme, we have two Gauss integration points, i.e. J = 2, on the fa
e κ. wj is the integration weight asso
iated with the Gauss point j. 40 To nd ∂Rκ,Ωi/∂AΩi , Eq. 4.15 
an be dire
tly dierentiated to get ∂Rκ,Ωi ∂AΩi = − 1 (AΩi) 2 ∑ j { 1 2 ( ~FRj + ~FLj − |△F1|j − |△F4|j − |△F5|j ) wj } = − 1 AΩi Rκ,Ωi (4.16) The term ∂Rκ,Ωi/∂nx,y 
an be written as ∂Rκ,Ωi ∂nx,y = 1 AΩi ∑ j { 1 2 ( ∂ ~FRj ∂nx,y + ∂ ~FLj ∂nx,y − ∂ |△F1|j ∂nx,y − ∂ |△F4|j ∂nx,y − ∂ |△F5|j ∂nx,y ) wj } (4.17) where 
omponent terms 
an be expanded using the denition of the Roe ux to give: ∂ ~FR,Lj ∂nx =  (ρu)R,L( ρu2 + p ) R,L (ρuv)R,L ({E + p}u)R,L  , ∂ ~FR,Lj ∂ny =  (ρv)R,L (ρuv)R,L( ρv2 + p ) R,L ({E + p} v)R,L  and ∂ |△F1|j ∂nx = uŨn√ Ũ2n  ( △ρ− △P ã2 ) 1 ũ ṽ ũ2+ṽ2 2 + ρ̃  0 △u− nx△U △v − ny△U ũ△u+ ṽ△v − Ũn△U  + ∣∣∣Ũn∣∣∣ ρ̃  0 −△U − nx△u −ny△v −ũ△U − Ũn△u   41 ∂ |△F1|j ∂ny = vŨn√ Ũ2n  ( △ρ− △P ã2 ) 1 ũ ṽ ũ2+ṽ2 2 + ρ̃  0 △u− nx△U △v − ny△U ũ△u+ ṽ△v − Ũn△U  + ∣∣∣Ũn∣∣∣ ρ̃  0 −nx△u −△U − ny△v −ṽ△U − Ũn△v   ∂ ∣∣∣△F̃4,5∣∣∣ ∂nx = ( Ũn ± ã ) ũ√( Ũn ± ã )2 (△P ± ρ̃ã△U 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã + ∣∣∣Ũn ± ã∣∣∣ (±ρ̃ã△u 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã + ∣∣∣Ũn ± ã∣∣∣ (△P ± ρ̃ã△U 2ã2 ) 0 ±ã 0 ±ũã  42 ∂ ∣∣∣△F̃4,5∣∣∣ ∂ny = ( Ũn ± ã ) ṽ√( Ũn ± ã )2 (△P ± ρ̃ã△U 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã + ∣∣∣Ũn ± ã∣∣∣ (±ρ̃ã△v 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã + ∣∣∣Ũn ± ã∣∣∣ (△P ± ρ̃ã△U 2ã2 ) 0 0 ±ã ±ṽã  To nd ∂Rκ,Ωi/∂wj , Eq. 4.15 is dierentiated; we get ∂Rκ,Ωi ∂wj = 1 (AΩi) { 1 2 ( ~FRj + ~FLj − |△F1|j − |△F4|j − |△F5|j )} (4.18) The dependen
y of the fa
e residual 
ontribution Rκ,Ωi on the re
onstru
ted ow properties at the fa
e U fa
e 
an found as the sum of the dependen
y on the right and left fa
e ow properties as ∂Rκ,Ωi ∂U fa
e = ∂Rκ,Ωi ∂UR + ∂Rκ,Ωi ∂UL (4.19) = 1 (AΩi) ∑ 1 2  ∂ ~FRj ∂UR + ∂ ~FLj ∂UL − ∂ { |△F1|j + |△F4|j + |△F5|j } ∂UR − ∂ { |△F1|j + |△F4|j + |△F5|j } ∂UL wj  The terms in Eq. 4.19 
an be obtained using symboli
 manipulator like Maple® or Matlab®. Another possible approa
h is to use an automati
 dierentiation pa
k- age to dierentiate the C, C++, or FORTRAN fun
tion that 
omputes the fa
e ux. 43 Figure 4.2: General triangular element with unit normal n̂ on one of its fa
es κ. Finite dieren
es 
an also be easily implemented to evaluate ∂RκΩi/∂U fa
e ∂Rκ,Ωi ∂U fa
e = Rκ,Ωi (UR + ǫ)−Rκ,Ωi (UR − ǫ) 2ǫ + Rκ,Ωi (UL + ǫ)−Rκ,Ωi (UL − ǫ) 2ǫ +o (ǫ)2 (4.20) where ǫ is 
hosen to be 10−8 in the above 
entral nite dieren
e formula so that the error will be on the order of ma
hine zero. The ANSLib CFD 
ode 
an use both approa
hes; they lead to nearly identi
al answers. 4.2 Element and Fa
e Geometri
 Properties Dependen
y on Mesh Coordinates Element and fa
e geometri
 properties depend dire
tly on the spatial 
oordinates of the verti
es. Fig. 4.2 shows a s
hemati
 drawing of a triangular element used for 
ell 
entered nite volume simulations, with its three verti
es and one of its three fa
es labeled for later referen
e. In the next subse
tions, the pro
edure for evaluating geometri
 properties like element area, fa
e length, and fa
e normals will be presented in addition to the mesh 
oordinate dependen
y of these properties. 44 4.2.1 Element area mesh dependen
y The element area of the triangular element shown in Fig 4.2
an be obtained as half of the 
ross produ
t of the two ve
tors −→r12,−→r13 A = σ (x2 − x1) (y3 − y1)− (y2 − y1) (x3 − x1) 2 (4.21) The verti
es' spatial lo
ations are related to the airfoil surfa
e points via the mesh movement s
heme, while the airfoil surfa
e points are related to the design vari- ables through the pseudo inverse matrix in the least-squares t during geometry parametrization. It worth mentioning that any general polygon 
an be split into several triangles and the area of ea
h triangle 
an be 
al
ulated using Eq. 4.21. The element area mesh dependen
y 
an be 
al
ulated by dire
t dierentiation of Eq. 4.21 dA dM = ∂A ∂x1 ∂x1 ∂M + ∂A ∂x2 ∂x2 ∂M + ∂A ∂x3 ∂x3 ∂M (4.22) + ∂A ∂y1 ∂y1 ∂M + ∂A ∂y2 ∂y2 ∂M + ∂A ∂y3 ∂y3 ∂M where ∂A ∂x1 = σ (y2 − y3) 2 ∂A ∂x2 = σ (y3 − y1) 2 ∂A ∂x3 = σ (y1 − y2) 2 ∂A ∂y1 = σ (x3 − x2) 2 ∂A ∂y2 = σ (x1 − x3) 2 ∂A ∂y3 = σ (x2 − x1) 2 45 σ = sgn ((x2 − x1) (y3 − y1)− (y2 − y1) (x3 − x1)) 4.2.2 Fa
e length mesh dependen
y The length of fa
e κ in Fig. 4.2, whi
h is a general straight fa
e in the mesh, 
an be 
al
ulated as L = √ (x3(M)− x2(M))2 + (y3(M)− y2(M))2 (4.23) The fa
e length mesh dependen
y 
an be obtained as dL dM = ∂L ∂x2 ∂x2 ∂M + ∂L ∂x3 ∂x3 ∂M + ∂L ∂y2 ∂y2 ∂M + ∂ L ∂y3 ∂y3 ∂M where ∂L ∂x2 = −x3 + x2√ (x3 − x2)2 + (y3 − y2)2 = − (x3 − x2) L ∂L ∂x3 = x3 − x2√ (x3 − x2)2 + (y3 − y2)2 = (x3 − x2) L ∂L ∂y2 = −y3 + y2√ (x3 − x2)2 + (y3 − y2)2 = − (y3 − y2) L ∂L ∂y3 = y3 − y2√ (x3 − x2)2 + (y3 − y2)2 = (y3 − y2) L The previous relations are for straight fa
es whi
h exist all over the domain interior, but for higher order s
hemes, the boundary fa
es are 
urved to enable higher order ux integration. A boundary element with one 
urved fa
e is shown in Fig. 4.3 The 
urved fa
e length is obtained using numeri
al integration, therefore to avoid dierentiating the numeri
al integration s
heme, ∂L/∂x2,3 and ∂L/∂y2,3 are evalu- 46 Figure 4.3: Boundary element with 
urved fa
e for high order integration s
heme ated using nite dieren
es as ∂L ∂x2,3 = L(x2,3 + ǫ)− L(x2,3 − ǫ) 2ǫ ∂L ∂y2,3 = L(y2,3 + ǫ)− L(y2,3 − ǫ) 2ǫ (4.24) For a se
ond order s
heme, only one Gauss integration point exists at the middle of the fa
e with wi = L, therefore dwi dM = dL dM For the fourth order s
heme, two Gauss integration points are required for ux integration. Figure 4.2.2 shows a s
hemati
 drawing of a 
urved edge with two Gauss points on it. Table 4.1 shows their integration weights and their parametri lo
ations ki on the fa
e starting from the point x2, y2. 47 Figure 4.4: S
hemati
 drawing of a 
urved fa
e with two Gauss points used Point 1 Point 2 wi L/2 L/2 ki ( 1 2 − 1√12)L ( 1 2 + 1√ 12 )L Table 4.1: Two point Gauss quadrature rule 48 4.2.3 Fa
e normal mesh dependen
y The unit normal ve
tor of fa
e κ that points outward from element Ωi (shown in Fig. 4.2) 
an be found using the three verti
es of the element as n̂ = [ nx ny ] = (~c3 − βr̂23) ‖~c3 − βr̂23‖ (4.25) where ~c3 = [ 2/3x3 − (1/3x1 + 1/3x2) 2/3 y3 − (1/3 y1 + 1/3 y2) ] r̂23 =  x3−x2√(|−x3+x2|)2+(|−y3+y2|)2 y3−y2√ (|−x3+x2|)2+(|−y3+y2|)2  β = (x3 − x2) (2/3x3 − (1/3x1 + 1/3x2))√ (|−x3 + x2|)2 + (|−y3 + y2|)2 + (y3 − y2) (2/3 y3 − (1/3 y1 + 1/3 y2))√ (|−x3 + x2|)2 + (|−y3 + y2|)2 The mesh dependen
y of the unit normal 
an be written as dn̂ dM = ∂n̂ ∂x1 ∂x1 ∂M + ∂n̂ ∂x2 ∂x2 ∂M + ∂n̂ ∂x3 ∂x3 ∂M + ∂n̂ ∂y1 ∂y1 ∂M + ∂n̂ ∂y2 ∂y2 ∂M + ∂n̂ ∂y3 ∂y3 ∂M (4.26) where ∂n̂ ∂xi = α1 − α2 ‖~c3 − βr̂23‖2 α1 = ‖~c3 − βr̂23‖ · ∂ ∂xi (~c3 − βr̂23) α2 = (~c3 − βr̂23) · ∂ ∂xi ‖~c3 − βr̂23‖ 49 ∂n̂ ∂yi = α3 − α4 ‖~c3 − βr̂23‖2 α3 = ‖~c3 − βr̂23‖ · ∂ ∂yi (~c3 − βr̂23) α4 = (~c3 − βr̂23) · ∂ ∂yi ‖~c3 − βr̂23‖ i = 1, 2, 3 The terms in the last equations 
an be found using automati
 dierentiation of the unit normal expression. 4.3 Fa
e Flow Properties Dependen
y on The Mesh The evaluation of the ∂U fa
e /∂M term in the mesh sensitivity of the residual, using Eq. 4.9, depends on the details of the CFD solver. The CFD solver used in this resear
h is the Advan
ed Numeri
al Simulation Library (ANSLib) whi
h is a multi- physi
s nite volume solver 
apable of 
ondu
ting CFD simulations up to fourth or- der a

ura
y; this solver has been built by Ollivier-Goo
h and 
o-workers [61, 56, 54℄. In ANSLib, the ow properties are assumed to 
hange within ea
h 
ontrol volume a

ording to a linear polynomial for se
ond order simulations, or a 
ubi
 polynomial for fourth order simulations. These polynomials are found using a least-squares re- 
onstru
tion; the least-squares system depends on the 
ontrol volume properties like moments of area whi
h eventually depend on mesh point lo
ations. In the next two subse
tions, the re
onstru
tion and mesh dependen
y of ow properties at fa
es will be presented. 4.3.1 Fa
e ow properties re
onstru
tion Solution re
onstru
tion is the key to determine solver a

ura
y. The solution is assumed to vary linearly within the 
ontrol volume for se
ond order a

ura
y; for higher order methods, ow properties are assumed to vary a

ording to higher order polynomial; variation a

ording to 
ubi
 polynomial is fourth order a

urate. The ow solution re
onstru
tion method presented in this subse
tion follows the method des
ribed in Ollivier-Goo
h and Van Altena [68℄ and Ollivier-Goo
h et al [69℄. 50 For se
ond order solution re
onstru
tion, the primitive variables U = [ ρ u v p ]T are re
onstru
ted at the fa
e Gauss points as UR2 (x, y) = Uref + (x− xref ) ∂U ∂x ∣∣∣∣ ref + (y − yref ) ∂U ∂y ∣∣∣∣ ref (4.27) For fourth order solution re
onstru
tion, the re
onstru
tion polynomial takes the form of a third degree Taylor series expansion in two variables: UR4 (x, y) = U R 2 (x, y) + (x− xref )2 2 ∂2U ∂x2 ∣∣∣∣ ref (4.28) + (x− xref) . (y − yref) ∂ 2U ∂x∂y ∣∣∣∣ ref + (y − yref )2 2 ∂2U ∂y2 ∣∣∣∣ ref + (x− xref )3 6 ∂3U ∂x3 ∣∣∣∣ ref + (x− xref)2 (y − yref) 2 ∂3U ∂x2∂y ∣∣∣∣ ref + (x− xref ) (y − yref )2 2 ∂3U ∂x∂y2 ∣∣∣∣ ref + (y − yref)3 6 ∂3U ∂y3 ∣∣∣∣ ref The referen
e point ~xref is 
hosen to be the triangle 
enter for 
ell-
entered nite volume s
heme and the vertex lo
ation for vertex-
entered s
heme. A 
onstrained least-squares system is 
onstru
ted to nd the re
onstru
tion poly- nomial 
oe
ient. The hard 
onstraint is that the re
onstru
tion polynomial should preserve the 
omputed 
ontrol volume average solution. The rest of the equations  whi
h are satised only in a least-squares sense  are obtained by requiring that the average of the re
onstru
tion polynomial approximates the 
ontrol volume average solution for neighboring 
ontrol volumes. Neighbors are 
hosen based on their topo- logi
al distan
e from the element, and entire layers of neighbors is added until the number of neighbor elements in the sten
il equals or ex
eeds 4 for the se
ond order s
heme and 16 for the fourth order s
heme. Figure 4.5 shows a s
hemati
 drawing of an element k and its neighbor layers. The mean 
onstraint equation is obtained by requiring 
onservation of the 
ontrol 51 Figure 4.5: S
hemati
 drawing of rst, se
ond and third neighbor layers of triangular element k. volume average by the re
onstru
tion polynomial and 
an be written as Uk = 1 Ak ˆ Ak URk dA ≡ Uk ref + xk ∂U ∂x ∣∣∣∣ k ref + yk ∂U ∂y ∣∣∣∣ k ref + (4.29) x2k 2 ∂2U ∂x2 ∣∣∣∣ k ref + xyk ∂2U ∂x∂y ∣∣∣∣ k ref + y2k 2 ∂2U ∂y2 ∣∣∣∣ k ref + ... where xnymk = 1 Ak ¨ Ak (x− xk ref )n (y − yk ref )m dA 52 Mat
hing the 
ontrol volume average in a neighbor j would require that U j = 1 Aj ˆ Aj URk dA ≡ Uk ref + 1 Aj ˆ Aj (x− xk ref ) ∂U ∂x ∣∣∣∣ ref dA+ (4.30) 1 Aj ˆ Aj (y − yk ref ) ∂U ∂y ∣∣∣∣ ref dA+ 1 Aj ˆ Aj (x− xk ref )2 2 ∂2U ∂x2 ∣∣∣∣ ref dA+ 1 Aj ˆ Aj (x− xk ref ) (y − yk ref ) ∂ 2U ∂x∂y ∣∣∣∣ ref dA+ 1 Aj ˆ Aj (y − yk ref )2 2 ∂2U ∂y2 ∣∣∣∣ ref dA+ ... Using (x− xk ref ) = (x− xj ref )− (xk ref − xj ref ) and (y − yk ref ) = (y − yj ref )− (yk ref − yj ref ) , we get U j = 1 Aj ˆ Aj URk dA ≡ Uk ref + x̂k,j ∂U ∂x ∣∣∣∣ k ref + ŷk,j ∂U ∂y ∣∣∣∣ k ref + (4.31) x̂2k 2 ∂2U ∂x2 ∣∣∣∣ k ref + x̂yk,j ∂2U ∂x∂y ∣∣∣∣ k ref + ŷ2k,j 2 ∂2U ∂y2 ∣∣∣∣ k ref + ... where x̂nymk,j = 1 Aj ˆ Aj ((x− xj ref )− (xk ref − xj ref ))n · ((y − yj ref )− (yk ref − yj ref ))m dA = m∑ r=0 n∑ q=0 m! r! (m− r) ! n! q! (n− q) ! (xj ref − xk ref ) q · (yj ref − yk ref )r xn−qj ref ym−rj ref for simpli
ity, the term x, yk ref will be written as x, yk, and x, yj ref will be repla
ed 53 with x, yj . The resulting 
onstrained least-squares system 
an be written as  1 xk yk x 2 k xyk y 2 k · · · 1 x̂k,1 ŷk,1 x̂2k,1 x̂yk,1 ŷ 2 k,1 · · · 1 x̂k,2 ŷk,2 x̂2k,2 x̂yk,2 ŷ 2 k,2 · · · 1 x̂k,3 ŷk,3 x̂2k,3 x̂yk,3 ŷ 2 k,3 · · · . . . . . . . . . . . . . . . . . . . . . 1 x̂k,N ŷk,N x̂2k,N x̂yk,N ŷ 2 k,N · · · · · ·   U ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . .  k =  Uk U1 U2 U3 . . . UN  (4.32) This 
onstrained least squares system 
an be rewritten as min ∥∥∥[ARec](P̃Rec)− (URec − Uk)∥∥∥ satisfying [BRec] (PRec) = ( Uk ) where [ARec] =  x̂k,1 − xk, ŷk,1 − yk, x̂2k,1 − x2k, x̂yk,1 − xyk, ŷ2k,1 − y2k, · · · x̂k,2 − xk, ŷk,2 − yk, x̂2k,2 − x2k, x̂yk,2 − xyk, ŷ2k,2 − y2k, · · · x̂k,3 − xk, ŷk,3 − yk, x̂2k,3 − x2k, x̂yk,3 − xyk, ŷ2k,3 − y2k, · · · . . . . . . . . . . . . . . . . . . x̂k,N − xk, ŷk,N − yk, x̂2k,N − x2k, x̂yk,N − xyk, ŷ2k,N − y2k, · · · · · ·  [BRec] = [ 1 xk yk x 2 k xyk y 2 k · · · ] and PRec =  U ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . .  k , P̃Rec =  ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . .  k , ( URec − Uk ) =  U1 − Uk U2 − Uk U3 − Uk . . . UN − Uk  54 The least-squares system is solved to obtain the re
onstru
tion polynomial 
oef- 
ients using singular value de
omposition method to nd the pseudo inverse of the least-squares 
oe
ient matrix as follows ARec = UΣV T A†Rec = V Σ †UT  ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . .  k =  x̂k,1 − xk, ŷk,1 − yk, x̂2k,1 − x2k, x̂yk,1 − xyk, ŷ2k,1 − y2k, · · · x̂k,2 − xk, ŷk,2 − yk, x̂2k,2 − x2k, x̂yk,2 − xyk, ŷ2k,2 − y2k, · · · x̂k,3 − xk, ŷk,3 − yk, x̂2k,3 − x2k, x̂yk,3 − xyk, ŷ2k,3 − y2k, · · · . . . . . . . . . . . . . . . . . . x̂k,N − xk, ŷk,N − yk, x̂2k,N − x2k, x̂yk,N − xyk, ŷ2k,N − y2k, · · · · · ·  † ·  U1 − Uk U2 − Uk U3 − Uk . . . UN − Uk  (4.33) Uk = Uk − [ xk yk x 2 k xyk y 2 k · · · ] ·  ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . .  (4.34) The unlimited fa
e ow properties 
an be found at fa
e Gauss points ( xfi , yff ) 55 as follows Ufi = Uk + (xfi − xk) ∂U ∂x ∣∣∣∣ k + (yfi − yk) ∂U ∂y ∣∣∣∣ k + (4.35) (xfi − xk)2 2 ∂2U ∂x2 ∣∣∣∣ k + (xfi − xk) (yfi − yk) ∂2U ∂x∂y ∣∣∣∣ k + (yfi − yk)2 2 ∂2U ∂y2 ∣∣∣∣ k + ... To ensure monotoni
ity and solution 
onvergen
e, a limiter is applied in order not to 
reate a new extremum at the 
ontrol volume fa
e; the limited ow properties at the fa
e 
an be written as Ufi = Uk + φk ( (xfi − xk) ∂U ∂x ∣∣∣∣ k + (yfi − yk) ∂U ∂y ∣∣∣∣ k + (4.36) (xfi − xk)2 2 ∂2U ∂x2 ∣∣∣∣ k + (xfi − xk) (yfi − yk) ∂2U ∂x∂y ∣∣∣∣ k + (yfi − yk)2 2 ∂2U ∂y2 ∣∣∣∣ k + ... ) where φk is the limiter value in the 
ontrol volume k. ANSLib uses two limiter fun
tions to 
al
ulate φ: Venkatakrishnan's limiter [84℄ and a higher-order limiter [56℄ designed not to degrade the a

ura
y of high-order s
hemes for smooth ows. We will derive the fa
e properties-mesh dependen
y for the limited 
ase; the unlimited 
ase 
an be obtained by simpli
ation of the limited 
ase by setting φk = 1. 56 4.3.2 Mesh dependen
e of the fa
e ow property re
onstru
tion The term ∂U fa
e /∂M for the limited 
ase 
an be obtained using dire
t dierentiation of Eq. 4.36 ∂Ufi ∂M = ∂Uk ∂M + ∂φk ∂M ( (xfi − xk) ∂U ∂x ∣∣∣∣ k + (yfi − yk) ∂U ∂y ∣∣∣∣ k + (4.37) (xfi − xk)2 2 ∂2U ∂x2 ∣∣∣∣ k + (xfi − xk) (yfi − yk) ∂2U ∂x∂y ∣∣∣∣ k + (yfi − yk)2 2 ∂2U ∂y2 ∣∣∣∣ k + ... ) + φk ( ∂ (xfi − xk) ∂M ∂U ∂x ∣∣∣∣ k + ∂ (yfi − yk) ∂M ∂U ∂y ∣∣∣∣ k + ∂ (xfi−xk) 2 2 ∂M ∂2U ∂x2 ∣∣∣∣ k + ∂ (xfi − xk) (yfi − yk) ∂M ∂2U ∂x∂y ∣∣∣∣ k + ∂ (yfi−yk) 2 2 ∂M ∂2U ∂y2 ∣∣∣∣ k + ... + φk (xfi − xk) ∂ ∂U∂x ∣∣ k ∂M + (yfi − yk) ∂ ∂U ∂y ∣∣∣ k ∂M + (xfi − xk)2 2 ∂ ∂ 2U ∂x2 ∣∣∣ k ∂M + (xfi − xk) (yfi − yk) ∂ ∂ 2U ∂x∂y ∣∣∣ k ∂M + (yfi − yk)2 2 ∂ ∂ 2U ∂y2 ∣∣∣ k ∂M + ...  The unlimited fa
e ow property-mesh dependen
y 
an be obtained by setting φk = 1 and ∂φk/∂M = 0 in Eq. 4.37. In Eq. 4.37, the term ∂φk/∂M 
an be obtained by dierentiating the limiter 57 expression, and the terms ∂ ∂ nU ∂xm∂yn−m ∣∣∣ k /∂M 
an be evaluated as follows: [ARec] { P̃Rec } = { URec } ∂ [ARec] ∂M { P̃Rec } + [ARec] ∂ { P̃Rec } ∂M = ∂ { URec } ∂M = {0} ∴ ∂ { P̃Rec } ∂M = − [ARec]† { ∂ [ARec] ∂M { P̃Rec }} ∂ { P̃Rec } ∂M = − [ARec]† ∂ [ARec] ∂M [ARec] † {URec}(4.38) and ∵ Uk =Uk − [ xk yk x 2 k xyk y 2 k · · · ] · { P̃Rec } ∴ ∂Uk ∂M =− ∂ ∂M [ xk yk x 2 k xyk y 2 k · · · ] · { P̃Rec } − (4.39) [ xk yk x 2 k xyk y 2 k · · · ] · ∂ { P̃Rec } ∂M where ∂ [ARec] ∂M =  ∂x̂k,1 ∂M ∂ŷk,1 ∂M ∂x̂2k,1 ∂M ∂x̂yk,1 ∂M ∂ŷ2k,1 ∂M · · · ∂x̂k,2 ∂M ∂ŷk,2 ∂M ∂x̂2k,2 ∂M ∂x̂yk,2 ∂M ∂ŷ2k,2 ∂M · · · ∂x̂k,3 ∂M ∂ŷk,3 ∂M ∂x̂2k,3 ∂M ∂x̂yk,3 ∂M ∂ŷ2k,3 ∂M · · · . . . . . . . . . . . . . . . . . . ∂x̂k,N ∂M ∂ŷk,N ∂M ∂x̂2k,N ∂M ∂x̂yk,N ∂M ∂ŷ2k,N ∂M · · · · · ·  ∂ { P̃Rec } ∂M = ∂ ∂M  ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . .  58 and ∂ ∂M {(xfi − xk)m (yfi − yk)n} = m (xfi − xk)m−1 (yfi − yk)n ( ∂xfi ∂M − ∂xk ∂M ) + n (xfi − xk)m (yfi − yk)n−1 ( ∂yfi ∂M − ∂yk ∂M ) ∂x̂nymk,j ∂M = m∑ r=0 n∑ q=0 { m! r! (m− r) ! n! q! (n− q) ! q (xj ref − xk ref ) q−1 · (yj ref − yk ref ) · xn−qj ref y m−r j ref · ( ∂xj ref ∂M − ∂xk ref ∂M ) + m! r! (m− r) ! n! q! (n− q) ! r (xj ref − xk ref ) q (yj ref − yk ref )r−1 · xn−qj ref y m−r j ref · ( ∂yj ref ∂M − ∂yk ref ∂M ) 59 Chapter 5 GRADIENT VALIDATION In this se
tion, the 
omparison of the gradient a

ura
y evaluated using se
ond and fourth order s
hemes is 
arried out. Both se
ond and fourth order gradients are evaluated using the adjoint approa
h and 
ompared to their 
orresponding nite dieren
e gradient. We present three test 
ases: subsoni
, non-limited transoni
 and limited transoni
 test 
ases. In all test 
ases we use 18 design points to parametrize the airfoil geometry. 5.1 Subsoni
 Test Case In this test 
ase, the evaluation of the lift 
oe
ient gradient with respe
t to the airfoil geometri
 design 
ontrol points is presented for both se
ond and fourth order s
hemes, in
luding a 
omparison with a nite dieren
e gradient 
al
ulated using the same order of a

ura
y in the ow solver. The airfoil used in this test 
ase is NACA 0012 at subsoni
 
onditions M = 0.5 and α = 2o. Figure 5.1 shows a representative 
omparison between the eld of pressure sensitivity with respe
t to one of the geometry design points 
omputed using nite dieren
es and the solution sensitivity 
al
ulation of Eq. 4.6 for both se
ond and fourth order a

urate 
ompu- tations. Agreement is ex
ellent in all 
ases; this is also true for the other design 
ontrol points. The ex
ellent mat
hing of the pressure sensitivity when 
omparing the se
ond order and the fourth order results indi
ates that the two s
hemes will give similar gradient ve
tors and similar optimization des
ent dire
tions for subsoni ow optimization. Table 5.1 shows quantitatively an ex
ellent mat
hing between the obje
tive fun
- tion gradient magnitude (less than half a per
ent dieren
e) and dire
tion (less than a half a degree dieren
e) when 
omparing nite dieren
es and the adjoint ap- proa
h of Eq. 4.8 for both se
ond and fourth order s
hemes. Figure 5.1 shows good 60 (a) Se
ond order sensitivity 
al
ulation (b) Se
ond order nite dieren
e (
) Fourth order sensitivity 
al
ulation (d) Fourth order nite dieren
e Figure 5.1: The pressure sensitivity with respe
t to one of the design 
ontrol points 
omputed for subsoni
 ow over NACA 0012, 
omparing the sensitivity 
al
ulation of Eq. 4.6 with nite dieren
e results. agreement between the obje
tive fun
tion gradient 
omponents ∂CL/∂yi, where yi are the y lo
ations of the design 
ontrol points, 
omputed using adjoint and nite dieren
e approa
hes for both se
ond and fourth order 
omputations. The maxi- mum error, normalized by the gradient magnitude, is only 0.005 whi
h gives a high level of a

ura
y for the 
omputed gradient. Taken together, these results imply that the order of dis
retization error has little ee
t on the 
omputed gradient ve
tor for subsoni
 ow. 5.2 Transoni
 Test Case with No Limiter In this test 
ase, the sensitivity analysis is 
arried out for NACA 0012 airfoil with M = 0.8 and α = 2◦. The drag 
oe
ient gradient is evaluated without using limiters in the CFD simulation, so overshoot/undershoot at the sho
k lo
ation is expe
ted. Figure 5.3 shows good agreement of the pressure sensitivity 
omputed using Eq. 4.6 61 Se
ond order Fourth order Adjoint FD Adjoint FD Gradient ve
tor magnitude 18.605 18.561 18.575 18.655 Angle with se
ond order adjoint 0 ◦ 0.247 ◦ 1.117 ◦ 1.330 ◦ Angle with fourth order adjoint 1.117 ◦ 0.917 ◦ 0 ◦ 0.249 ◦ Table 5.1: The magnitude of se
ond and fourth order CL gradients and angles be- tween the evaluated gradients for NACA 0012 in subsoni
 ow. Figure 5.2: CL gradient error in se
ond and fourth order s
hemes with respe
t to the design points normalized by gradient magnitude. 62 (a) Se
ond order sensitivity 
al
ulation (b) Se
ond order nite dieren
e (
) Fourth order sensitivity 
al
ulation (d) Fourth order nite dieren
e Figure 5.3: The pressure sensitivity with respe
t to one of the design 
ontrol points 
omputed for an unlimited transoni
 ow around NACA 0012, 
omparing the sensi- tivity 
al
ulation of Eq. 4.6 with nite dieren
e results. and nite dieren
e; the gure also shows that for transoni
 ows, the pressure sensitivity 
omputed using se
ond order and fourth order a

urate adjoint s
heme are dierent, espe
ially near the sho
k lo
ation. Table 5.2 shows again the ex
ellent mat
hing between the 
omputed adjoint and nite dieren
e gradient with a dire
tion dieren
e less than a degree and nearly identi
al magnitude. Figure 5.4 shows that for unlimited transoni
 ow, both se
- ond and fourth order adjoint gradients are an ex
ellent mat
h to the 
orresponding nite dieren
e gradient, with the se
ond order s
hemes mat
hing more 
losely than the fourth order s
heme. Comparison of the se
ond and fourth order gradient ve
tors shows little dieren
e between them (about 4% in magnitude and 2◦ in dire
tion). Again, we expe
t that the se
ond and fourth order s
hemes will give similar opti- mization sear
h dire
tions. 63 Se
ond order Fourth order Adjoint FD Adjoint FD Gradient ve
tor magnitude 1.895 1.895 1.814 1.813 Angle with se
ond order adjoint 0 ◦ 0.252 ◦ 2.236 ◦ 1.997 ◦ Angle with fourth order adjoint 2.236 ◦ 2.354 ◦ 0 ◦ 0.518 ◦ Table 5.2: The magnitude of se
ond and fourth order CD gradients and angles between the evaluated gradients for NACA 0012 in an unlimited transoni
 ow. Figure 5.4: CD gradient error in se
ond and fourth order s
hemes with respe
t to the design points, normalized by gradient magnitude. 64 (a) Se
ond order sensitivity 
al
ulation (b) Se
ond order nite dieren
e (
) Fourth order sensitivity 
al
ulation (d) Fourth order nite dieren
e Figure 5.5: The pressure sensitivity with respe
t to one of the design 
ontrol points 
omputed for an unlimited transoni
 ow around NACA 0012. 5.3 Transoni
 Test Case with Limiter In this test 
ase, the impa
t of using a limiter in the CFD simulation on the a

ura
y of the 
omputed gradient for se
ond order and fourth order s
hemes is studied. Two dierent limiters are used, the Venkatakrishnan limiter [84℄ and the higher order limiter of Mi
halak and Ollivier-Goo
h [56℄. Figure 5.5 shows very good mat
hing of the se
ond order pressure sensitivity 
omputed using adjoint and nite dieren
e te
hniques with the use of Venkatakrishnan limiter. Mat
hing is less good between the pressure sensitivity and nite dieren
e results for the fourth order a

urate s
heme; this lower level of pressure sensitivity mat
hing will lead to less a

urate gradient values when using the limited fourth order s
heme. Table 5.3 shows that with the use of the Venkatakrishnan limiter, the se
ond order gradient magnitude is a very good mat
h with the 
orresponding nite dieren
e gradient; the larger error in gradient value observed for the fourth order s
heme is 
omparable to the dieren
e in magnitude between se
ond and fourth order s
hemes 65 Se
ond order Fourth order Adjoint FD Adjoint FD Gradient ve
tor magnitude 1.897 1.895 1.887 1.812 Angle with se
ond order adjoint 0 ◦ 3.612 ◦ 2.292 ◦ 3.334 ◦ Angle with fourth order adjoint 2.292 ◦ 5.275 ◦ 0 ◦ 4.873 ◦ Table 5.3: Magnitudes of se
ond and fourth order CD gradients and angles between the evaluated gradients for NACA 0012 in Venkatakrishnan limited transoni
 ow Se
ond order Fourth order Adjoint FD Adjoint FD Gradient ve
tor magnitude 1.676 1.753 1.576 1.672 Angle with se
ond order adjoint 0 ◦ 2.845 ◦ 30.1 ◦ 18.8 ◦ Angle with fourth order adjoint 30.1 ◦ 21.05 ◦ 0 ◦ 18.9 ◦ Table 5.4: The magnitude of se
ond and fourth order CD gradients and angles be- tween the evaluated gradients for NACA 0012 using higher order limiter in transoni ow. for the unlimited transoni
 
ase. Also, the dieren
e in nite dieren
e and adjoint gradient dire
tion grows to several degrees with the use of Venkatakrishnan limiter. Table 5.4 shows the same behavior with the higher order limiter of C. Mi
halak and C. Ollivier-Goo
h. It shows also that the error in the gradient magnitude is larger 
ompared to the Venkatakrishnan results. The mat
hing in adjoint and nite dieren
e gradient dire
tions is good for the se
ond order s
heme but it is not for the fourth order s
heme with the use of higher order limiter. It is well known that the use of limiters 
auses 
onvergen
e problem due to the non-dierentiability of the limiting pro
edure [55℄. The limiter ae
ts both the right and left hand sides of Eq. 1.6. We spe
ulate that the poor mat
hing between nite dieren
e and adjoint gradients for the high-order limiter is related to the non-dierentiability of the limiting pro
edure and that this ee
t will vary in strength from one limiting pro
edure to another. To redu
e the inuen
e of ∂Ri/∂D in 
ontrol volume i when 
al
ulating the adjoint gradient, the fourth order Ja
obian was modied numeri
ally by making the non zero stru
ture of the fourth order Ja
obian matrix the same as the non zero stru
ture of the se
ond order Ja
obian, and dropping the rest of values in the fourth order Ja
obian matrix. The right hand side is still 
onstru
ted with fourth order 66 Mod 4th order Limiter Venkat. HO Gradient s
heme Adj. FD Adj. FD Gradient ve
tor magnitude 1.735 1.812 1.966 1.672 Angle with mod 4th order Venkat. adjoint 0 ◦ 2.845 ◦ 6.247 ◦ 3.338 ◦ Angle with mod 4th order HO. adjoint 6.247 ◦ 7.256 ◦ 0 ◦ 7.207 ◦ Table 5.5: The magnitudes and angles between the evaluated modied fourth order CD gradients using adjoint, and nite dieren
e for NACA 0012 using Venkatakr- ishnan and higher order limiters in transoni
 ow. a

ura
y. The above modi
ation doesn't ae
t the a

ura
y of the CFD simulation be
ause the right hand side remains fourth order a

urate. The 
omputed gradient using this approa
h is presented in Table 5.5 and shows a redu
tion in the error in the fourth order 
omputed adjoint gradient, espe
ially in the dire
tion of the gradient. The limiter used from now on is the Venkatakrishnan limiter as it produ
es an adjoint 
omputed gradient with better mat
hing to the nite dieren
e gradient 
ompared with the high order limiter. 5.4 Sensitivity of nite dieren
e gradient to perturbation magnitude In the previous se
tions, I have used perturbation amplitude ǫ = 10−6 and using the 
entral dieren
e formula to evaluate the nite dieren
e gradient. In this se
tion, the sensitivity of the evaluated gradient with respe
t to ǫ is presented for the hardest 
ase, transoni
 ow with limiter. The limiter used in the transoni
 simulation is Venkatakrishnan's limiter. Table 5.6 shows the norm of the drag gradient for dierent values of perturbation amplitude ǫ; it shows also that de
reasing ǫ to be less than 10−6 will not pra
ti
ally 
hange the value of the 
omputed gradient norm. Therefore, ǫ = 10−6 is 
hosen. 67 Figure 5.6: The normalized CD gradient error in se
ond, fourth, and modied fourth order s
hemes with respe
t to the design points in a limited transoni
 ow (Venkatakrishnan limiter) Figure 5.7: The normalized CD gradient error in se
ond, fourth, and modied fourth order s
hemes with respe
t to the design points in a limited transoni
 ow (higher order limiter). 68 ǫ Se
ond order ∥∥∥∂CD∂D ∥∥∥ Fourth order ∥∥∥∂CD∂D ∥∥∥ 10−3 1.90237 1.86814 10−6 1.90275 1.87105 10−9 1.09278 1.87114 Table 5.6: Finite dieren
e drag gradient sensitivity with respe
t to perturbation amplitude ǫ 69 Chapter 6 GRADIENT BASED OPTIMIZATION TEST CASES In this se
tion we present four optimization test 
ases; in all of them we 
ompare the optimal shape resulting from using se
ond and fourth order s
hemes. The rst two 
ases are inverse design problems, one subsoni
 and the other transoni
. In both test 
ases, a target pressure distribution is obtained using CFD simulation of a parametrized NACA 2412 airfoil and the starting geometry is a NACA 0012 airfoil. The optimizer will try to nd the geometry whose surfa
e pressure distribution mat
hes the target pressure distribution. Those two test 
ases are important, as the design spa
e has only one solution point at whi
h the resulting pressure distribution from the optimal geometry will mat
h the target pressure distribution. The third test 
ase is a transoni
 drag minimization with no lift 
onstraint starting from the RAE 2822 airfoil. The obje
tive of this test 
ase is to minimize CD at M = 0.73 and angle of atta
k = 2◦. In this test 
ase a strong sho
k wave is formed near mid 
hord of the initial airfoil geometry and we are seeking a sho
k free geometry, or at least a geometry that produ
es a mu
h weaker sho
k wave. Geometri
 
onstraints are applied so the airfoil thi
kness will be positive all the way along the airfoil se
tion. The fourth test 
ase repeats test 
ase three but adds the lift 
oe
ient as a 
onstraint; in this 
ase, we 
ompare the resulting optimal shape with the optimization results of Brezillon and Gauger [11℄ . 6.1 Subsoni
 Inverse Design In this test 
ase, the target pressure distribution is obtained for the parametrized NACA 2412 at a subsoni
 
ondition, M = 0.5 and α = 2◦ , using se
ond and fourth order CFD simulations. The starting geometry is the parametrized NACA 0012. 70 The obje
tive fun
tion to be minimized is F = ˛ (PT − Pi)2 dS (6.1) The above obje
tive fun
tion and its gradient 
an be expressed in dis
rete form as F = ∑ (PT − Pi)2 ws (6.2) dF dxd = ∑ 2 (PT − Pi) (−∂Pi ∂xd ) ws + (PT − Pi)2 ∂ws ∂xd (6.3) where ws is the ar
 length asso
iated with the surfa
e Gauss point. Figure 6.1 shows the target pressure distribution of the NACA 2412, the initial pressure distribution of NACA 0012, and the optimized airfoil pressure distribution obtained by the se
ond order and fourth order s
hemes; both s
hemes su

essfully rea
hed the target pressure distribution. Figures 6.2 and 6.3 show the 
onvergen
e history and gradient norm history for both s
hemes. Both s
hemes took about the same number of optimization iterations (28 for se
ond order and 32 for fourth order) to drop the obje
tive fun
tion value by eight orders of magnitude. Figures 6.4 and 6.5 show how 
lose the optimized NACA 0012 is to the NACA 2412 using both se
ond order and fourth order s
hemes. The error between the target geometry and the optimized prole is larger in the fourth order s
heme due to ina

ura
y of the pressure interpolation s
heme used to evaluate the obje
tive fun
tion; nevertheless the resulting geometry is an ex
ellent mat
h with the NACA 2412, with an error of less 0.1% of the surfa
e movement. 6.2 Transoni
 Inverse Design In this test 
ase, the target pressure is obtained using CFD simulation of the ow over a NACA 2412 airfoil at transoni
 
onditions, M = 0.73, α = 2◦. The obje
- tive fun
tion to be minimized is again the integration of the square of the pressure dieren
e between the target pressure and the optimized pressure as expressed in Eq. 6.1. The fourth order optimization is based on the modied fourth order gra- dient evaluation strategy. Figure 6.6 shows the initial, target, and the optimized 71 (a) Se
ond order (b) Fourth order Figure 6.1: Subsoni
 NACA 2412 inverse design pressure distributions for the initial, target, and optimized airfoil proles 0 5 10 15 20 25 30 35 10−12 10−10 10−8 10−6 10−4 10−2 Iteration F   Fourth order Second order Figure 6.2: Se
ond and fourth order optimization 
onvergen
e history. 72 0 5 10 15 20 25 30 35 10−6 10−5 10−4 10−3 10−2 10−1 Iteration Fi rs t− or de r o pt im al ity   Fourth order Second order Figure 6.3: Gradient norm history Figure 6.4: Subsoni
 inverse design optimal airfoil shapes 73 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 x 10−5 Upper surface Dy x/c D y   4th order 2nd order 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 x 10−5 Lower surface Dy x/c D y   4th order 2nd order (a) Upper surfa
e (b) Lower surfa
e Figure 6.5: The dieren
e between the target prole and the optimized proles, se
ond and fourth order pressure distribution. The 
onvergen
e history and gradient norm history are shown in Figs. 6.7 and 6.8; the fourth order s
heme is slower to 
onverge 
ompared to the se
ond order s
heme, due to the use of larger number of airfoil surfa
e Gauss quadrature points used in the fourth order 
omputations (double the number used for the se
ond order s
heme) whi
h makes the minimizer slower to 
onverge. The obje
tive fun
tion dropped only three order of magnitudes before 
onvergen
e stall. This stall is due to the high non-linearity in the target pressure distribution be
ause of the existen
e of a strong sho
k wave in it. The noise in the pressure sensitivity generated by the sho
k wave in the target pressure distribution doesn't allow further 
onvergen
e; however, the gradient magnitude dropped four order of magnitudes from its initial value. Figures 6.9 and 6.10 show that the optimal shapes for the two s
hemes dier by less than 10−4 of the 
hord length on the lower surfa
e and less than 10−3 
hord length on the upper surfa
e (where the strong sho
k wave exists); both optimal shapes are in good agreement with the NACA 2412: the maximum deviation is about 5% of maximum surfa
e movement. 6.3 Drag Minimization without Lift Constraint In this test 
ase, minimization of the drag 
oe
ient will be 
arried out with no lift 
onstraint applied. The airfoil to be optimized is an RAE 2822 at transoni 74 (a) Se
ond order (b) Fourth order Figure 6.6: Subsoni
 NACA 2412 inverse design pressure distributions for the initial, target, and optimized airfoil proles Figure 6.7: Se
ond and fourth order optimization 
onvergen
e history. 75 0 5 10 15 20 25 30 35 40 45 10−4 10−3 10−2 10−1 100 G ra de  N or m iteration   2nd order 4th order Figure 6.8: Gradient norm history (a) Upper surfa
e (b) Lower surfa
e Figure 6.9: The dieren
e between the target prole and the optimized proles, se
ond and fourth order 76 Figure 6.10: Transoni
 inverse design optimal airfoil shapes. CL CL optimized airfoil CD CD optimized airfoil Se
ond order 0.865 0.765 0.0081 0.0046 6 Fourth order 0.849 0.759 0.0099 0.0047 Table 6.1: Aerodynami
 
oe
ients of original and optimized RAE 2822 airfoil at transoni
 
onditions. 
onditions: M = 0.73 and α = 2◦. Fig. 6.11 shows the initial solution with a strong sho
k wave standing at 70% 
hord. Geometri
 
onstraints are added to insure that there is no interse
tion between the airfoil upper and lower surfa
es along the airfoil. Figure 6.12 shows the optimal shapes and the optimized pressure elds resulting from using se
ond and fourth order s
hemes. Figures 6.12 and 6.13 show that the dieren
e between the se
ond and fourth order optimal proles is notable on both upper and lower surfa
es. The surfa
e pressure distribution of the original RAE 2822 airfoil and of the optimized airfoil for se
ond and fourth order s
hemes is shown in Fig. 6.14. Both s
hemes su

essfully produ
ed a similar sho
k free pressure distribution but with dierent proles. The original and optimized CD are shown in Table 6.1. A drag redu
tion of about 50% is a
hieved by both s
hemes. The value of the CD of the se
ond order optimized airfoil is obtained from fourth order a

urate CFD simulation over the optimized se
ond order airfoil prole to eliminate the dieren
es in dis
retization error in the nal results. Figure 6.16 shows the 
onvergen
e history of the optimization. Both s
hemes rea
hed their optimal solution in about the same number of optimization iterations. 77 Figure 6.11: Pressure 
ontours of RAE 2822 at Ma
h 0.73 and angle of atta
k 2 Figure 6.12: Optimized RAE 2822. 78 (a) Upper surfa
e (b) Lower surfa
e Figure 6.13: Optimization surfa
e displa
ements of the original RAE2822 surfa
es. Figure 6.14: Pressure distribution 
omparison for the original and the optimized geometries. 79 0 5 10 15 20 25 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 −grade     −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−−    +grad C D BFGS optimal profile 0 5 10 15 20 25 0 0.002 0.004 0.006 0.008 0.01 0.012 −grade     −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−−    +grad C D BFGS optimal profile (a) Se
ond order response surfa
e (b) Fourth order response surfa
e Figure 6.15: Obje
tive fun
tion response surfa
e along positive and negative gradient dire
tions 
entered at the optimal prole of RAE 2822 transoni
 drag minimization without lift 
onstraint. The number of CFD simulations for the se
ond order s
heme is 50 
ompared to 45 for the fourth order s
heme, yet the overall 
omputational 
ost of the se
ond order s
heme is 40% less than the fourth order 
omputations. The gradient of the obje
tive fun
tion is not zero at the optimal point; to understand the reasons for that, I have plotted the obje
tive fun
tion values along positive and negative gradient dire
tion at the optimal solution point as shown in Fig 6.15. The se
ond order response surfa
e shows that the obtained optimal solution is a possible lo
al optimum and the shape of the obje
tive fun
tion shows that the gradient 
omputed using nite dieren
e will not be zero at the optimal point. For the fourth order 
omputations, the gure shows that the optimizer rea
hed a near optimal solution. While a slightly lower value is available,the small dieren
e between the obje
tive fun
tion value at the obtained optimal point and the minimum value along the negative gradient dire
tion (on the order of the dis
retization error) may 
ause a line sear
h failure due to insu
ient de
rease in the obje
tive fun
tion. 6.4 Drag Minimization with Lift Constraint This test 
ase represents a typi
al optimization task required in aerospa
e industry, as it is required to minimize the drag while the lift is un
hanged. Therefore the lift 
oe
ient CL will be an aerodynami
 
onstraint. The original RAE 2822 geometry 80 Figure 6.16: Se
ond and fourth order optimization 
onvergen
e history. will be used as a starting shape; the obje
tive fun
tion to be minimized is F = CD + 10 · (CL − CLC )2 (6.4) where CLC is the original lift 
oe
ient of RAE 2822 atM = 0.73 and α = 2 ◦ . In this test 
ase both geometri
 and aerodynami
 
onstraints are applied to avoid getting a non-feasible airfoil geometry, while not ae
ting airfoil lift 
oe
ient. The nonlinear 
onstraint, the lift 
oe
ient, is added as a penalty term in the obje
tive fun
tion as shown in Eq. 6.4. Figure 6.17 show the optimized airfoil's pressure distribution for both se
ond and fourth order s
hemes. Both se
ond and fourth order s
hemes travel 99% of the way towards their optimal obje
tive fun
tion after 7 iterations; the se
ond order s
heme then be
omes very slow to rea
h its optimal value. The se
ond order s
heme 
osts 112 CFD simulations to rea
h its optimum 
ompared to 49 CFD simulations for the fourth order. Figure 6.20 shows a noti
eable dieren
e between the optimized proles, in
luding a notably larger nose radius and less aft 
hamber for the fourth order s
heme. A sho
k free geometry is obtained with weak 
ompression waves on the airfoil upper surfa
e. Drag is redu
ed by about 50%, as shown in Table 6.2, while the lift 
oe
ient is about 1% higher than its original value. 81 CL CL optimized airfoil CD CD optimized airfoil Se
ond order 0.865 0.871 0.0081 0.0047 7 Fourth order 0.849 0.853 0.0099 0.0051 Table 6.2: Aerodynami
 
oe
ients of original and optimized RAE 2822 airfoil at transoni
 
onditions. Figure 6.17: Surfa
e pressure distribution, of the initial and optimized geometries Figure 6.18: Se
ond and fourth order optimization 
onvergen
e history. 82 Figure 6.19: Optimal shapes 
omparison: se
ond order, fourth order, and optimized prole by Brezillon and Gauger 
ompared with the original RAE2822. (a) Upper surfa
e (b) Lower surfa
e Figure 6.20: Se
ond and fourth order optimized prole surfa
e displa
ements from the initial shape. 83 0 5 10 15 20 25 0 1 2 3 4 5 6 7 8 −grade     −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−−    +grad F BFGS Optimal profile 0 5 10 15 20 25 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 −grade     −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−−    +grad F BFGS Optimal Profile (a) Se
ond order response surfa
e (b) Fourth order response surfa
e Figure 6.21: Obje
tive fun
tion response surfa
e along positive and negative gradient dire
tions 
entered at the optimal prole of RAE 2822 transoni
 drag minimization without lift 
onstraint. At the optimization breakdown point (optimal solution found), the gradient is not zero, possibly due to the existen
e of the sho
k wave whi
h 
auses high non-linearity in the obje
tive fun
tion and lift 
onstraint behavior. Figure 6.21_shows plots of the obje
tive fun
tion response along negative and positive gradient dire
tion. It shows also that the obtained optimal solutions are lo
al minima along the gradient dire
tion, and shows also that the gradient 
omputed by nite dieren
e is not zero. It is also worth noting that the gradient of the obje
tive fun
tion at the breakdown point is dominated by the gradient of the lift 
onstraint. Figures 6.19 and 6.20 show a 
omparison between the optimized airfoil proles based on se
ond and fourth order 
omputations. Dieren
es in shape are espe
ially noti
eable near the leading edge, and in the reexed part of the lower surfa
e near the trailing edge. The se
ond order optimized prole Cd is four drag 
ounts less than the fourth order prole drag 
oe
ient, nearly 10%. The same problem was studied by Brezillon and Gauger [11℄; they used the se
ond order solver MEGAFLOW of the German Aerospa
e Center and obtained their optimal prole by 
ontrolling the 
amber line of the RAE 2822 via 20 
ontrol points. Comparing the optimal pressure distribution of our se
ond order 
omputations shown in Fig 6.17 and that of Brezil- lon, Fig 6.22, shows that both solvers tend to a

elerate the ow near the leading edge upper surfa
e more than the original airfoil shape, followed by a gradual de
el- eration till the trailing edge; on the other hand, the surfa
e pressure on the lower 84 Figure 6.22: initial and optimal pressure distribution obtained by Brezillon and Gauger [11℄ (presented with permission) airfoil surfa
e 
hanges only slightly from the original surfa
e pressure distribution of the RAE2822. This similarity in the optimized pressure distribution suggests that whatever geometry parametrization te
hnique is used, the surfa
e pressure 
hanges will be qualitatively similar. The lift penalty fa
tor used in this test 
ase is sele
ted su
h that violating the lift 
onstraint by one lift 
ount (10−2) 
auses 10 drag 
ounts (10−3) in
rease in the obje
tive fun
tion value. Table 6.3 shows the ee
t of the lift 
onstraint penalty weighting fa
tor on the optimal obje
tive fun
tion value, as well as the optimal drag and lift 
oe
ients; these results are based on se
ond order 
omputations. 85 Lift penalty weight F opt CD opt CL opt Lift penalty 5 0.0052 0.0047 0.873 3.2 × 10−4 10 0.0051 0.0047 0.871 3.6 × 10−4 20 0.0082 0.0080 0.868 1.8 × 10−4 Table 6.3: Lift penalty weight ee
t on Drag minimization of RAE 2822 with lift 
onstraint 6.5 Mesh Renement Study of a Drag Minimization with Lift Constraint In the last test 
ase, Table 6.2 shows that the optimal value of CD is dierent when 
omparing se
ond and fourth order s
heme results. This subse
tion examines the impa
t of mesh renement on the optimal value of CD. Three mesh grid sizes are used, with 5000, 11000, and 15000 triangles, respe
tively. In all test 
ases, CLC = 0.84 is used as a lift 
onstraint value for both se
ond and fourth order 
omputations. Figure 6.23 
ompares the dieren
e in the optimal CD value for se
ond and fourth order s
hemes with mesh renement; the plotted value of se
ond order CD optimal are obtained using fourth order CFD simulation on the optimal shape obtained using se
ond order s
heme. 86 Figure 6.23: Optimal CD value with mesh renement The gure shows that the two s
hemes tend to rea
h the same value of the optimal CD with mesh renement with only one drag 
ount dieren
e (whi
h is on order of the dis
retization error). This behavior is expe
ted as transoni
 CFD simulations using se
ond and fourth order s
hemes for the same geometry produ
e almost the same surfa
e pressure distribution and hen
e obje
tive fun
tion value. The root mean square of dieren
e between the se
ond and fourth order optimized airfoil proles of drops from 3 · 10−5 on the 
oarsest mesh to 2 · 10−5 on the nest grid. 87 Chapter 7 PARTICLE SWARM OPTIMIZATION AND A NEW HYBRID OPTIMIZATION METHOD 7.1 Introdu
tion to Gradient Free Optimization Gradient based optimization results depend on the optimization starting point. To nd a global optimal solution, gradient independent methods are an obvious 
an- didate. Perhaps the best-known non-gradient based optimization te
hnique is the geneti
 algorithm [30℄, a biologi
ally-inspired approa
h whi
h 
an nd the global optimal point of an obje
tive fun
tion with multiple lo
al optima. In geneti
 al- gorithm optimization, a randomly generated population of 
andidate solutions are generated; their tness is determined by the value of the obje
tive fun
tion [29℄. So- lution points with the best tness are sele
ted to be parents of the next generation. A new generation of ospring is generated by 
rossover and mutation of the parents' geneti
 
hara
teristi
s (that is, the parents' solution point lo
ations in the design spa
e). This pro
edure is repeated until the globally optimal solution is obtained. The number of individuals in ea
h generation should not be less than the number of design variables, otherwise poor optimization 
onvergen
e will result. Figure 7.1 shows a ow
hart of a simple geneti
 algorithm optimizer. The dependen
y of the population size on the number of design variables and the large number of genera- tions required to 
onverge to the global optimum 
ombine to make geneti
 algorithm impra
ti
ally expensive for large aerodynami
 optimization problems. Geneti
 algo- 88 rithm optimization has been used for aerodynami
 optimization problems by many resear
hers. A geneti
 algorithm was used for transoni
 wing drag minimization with a lift 
onstraint by Gage and Kroo [29℄. Anderson used GA for wing aerody- nami
 shape optimization with stru
tural 
onstraints [1℄. Jang and Lee used GA to maximize lift to drag ratio of an airfoil using the Euler ow model with the NACA 0012 as a starting geometry [44℄. Oyama et al applied a geneti
 algorithm with a Navier-Stokes solver for transoni
 wing optimization [70℄. They also explored the use of fra
tal analysis in GA aerodynami
 optimization [71℄. Parti
le swarm optimization is a random sear
h te
hnique rst introdu
ed in 1995 as a te
hnique that simulate the behavior of a herd of predators hunting for food [46℄. This te
hnique has proved to be faster than geneti
 algorithm to rea
h the global optimal by many resear
hers (for example, [34, 59℄); swarm intelligen
e has the advantage of not being highly sensitive to problem dimensions [79℄, whi
h is not true for evolutionary optimization like geneti
 algorithm. The parti
le swarm optimization method (des
ribed in Se
tion 7.2) 
an be made more e
ient by hy- bridizing it with sequential quadrati
 programming, as des
ribed in Se
tion 7.3. The ee
tiveness of this hybrid approa
h for un
onstrained and 
onstrained aerodynami optimization is demonstrated in Se
tion 7.4. 7.2 Swarm Intelligen
e Parti
le swarm optimization is a random sear
h method that sear
hes for the global optimal solution of an obje
tive fun
tion. It uses a number of swarm parti
les that s
ans the design spa
e looking for the global optimal solution. Ea
h swarm parti- 
le represents a point in the design spa
e and it moves through the design spa
e a

ording to a simple formula. A parti
le's velo
ity depends on its inertia, its own experien
e, and so
ial experien
e gained by 
ommuni
ating with other swarm par- ti
les. The parti
le swarm te
hnique keeps a re
ord of lo
ation of the best obje
tive fun
tion value for ea
h parti
le ( −→pi ) as well the global best lo
ation (−→g ), whi
h is the best lo
ation among the −→pi 's. Algorithm 7.1 is the general pseudo 
ode of the parti
le swarm te
hnique. The velo
ity of the parti
le in an optimization iteration depends on ˆ Its velo
ity in the previous iteration (momentum part) 89 Figure 7.1: Geneti
 Algorithm optimization 90 Algorithm 7.1 Parti
le swarm pseudo 
ode For ea
h parti
le i in the swarm Initialize parti
le position −→xi Set its best position to be its initial position: −→pi = −→xi Initialize its velo
ity −→ Vi = ~0 End Set the global best position: −→g = −→pi where F (−→pi ) < F (−→pj ) for all i 6= j Do For ea
h parti
le i Cal
ulate its velo
ity −→ Vi Update its position −→xi Cal
ulate its obje
tive fun
tion value (tness) F (−→xi) If the 
urrent tness value is better than its best tness (F (−→xi) < F (−→pi )) Set the 
urrent parti
le position to be its best position: −→pi = −→xi End if End for Sele
t the best parti
le obje
tive fun
tion and position to be the global best if it is better than the stored value While max iteration 
ount not rea
hed or 
onvergen
e 
riterion not met. 91 ˆ Its distan
e from its known best value −→pi (
ognitive part) ˆ Its distan
e from the swarm global best value −→g (so
ial part) Figure 7.2 shows s
hemati
ally how the above three fa
tors ae
t the parti
le ve- lo
ity, the velo
ity of ea
h parti
le and its position update are obtained using the following formulas −→ Vi(k + 1) = w −→ Vi(k) + c1 [ −→r1 ] · {−→pi −−→xi(k)} + c2 [−→r2 ] · {−→g −−→xi(k)} (7.1) −→xi(k + 1) = −→xi(k) +−→Vi(k + 1), i = 1..N (7.2) where w = 0.73 is the inertial weight fa
tor, c1 = c2 = 1.49 are the 
ognitive and so
ial a

eleration fa
tors respe
tively; these 
onstants are sele
ted a

ording to Cler
's 
onstri
tion models [16℄. −→r1 and −→r2 are two random ve
tors of size N whose elements are independent and uniformly randomly sele
ted from the uniform distribution on the interval [0, 1]. The  · operation is a 
omponent by 
omponent multipli
ation not a dot produ
t and this is true from this point through the rest of the 
hapter. Parti
le swarm optimization suers from solution stagnation on
e the parti
les have 
onverged to a high quality lo
ally optimal solution whi
h is not a true globally optimal point [23℄. To avoid this problem, a restart strategy is suggested when premature 
onvergen
e is dete
ted [9℄. At ea
h restart, swarm parti
les start from new positions in the design spa
e and 
onverge to an optimal solution; restarting is 
arried out m times, and the optimal solution is the best of all global bests. This te
hnique is 
alled the multi-start parti
le swarm optimization (MPSO). A drawba
k of this me
hanism is that restarting may 
ause repetition of the sear
h 
omputations. Evers and Ben Ghalia suggested a more e
ient me
hanism to es
ape from high quality lo
al wells. They re-s
atter the parti
les around the 
onvergen
e point in a region with radius smaller than the design spa
e norm, and large enough some at least some parti
les to be outside the basin of attra
tion of the lo
al optimum [24, 23℄; their s
heme is 
alled the regrouped parti
le swarm optimization (RegPSO). The main three elements in their s
heme are ˆ Prematurity dete
tion. 92 Figure 7.2: Velo
ity 
omponents and position update of a parti
le ˆ Redening upper and lower boundaries for de
ision (design) variables. ˆ Regrouping the parti
les by s
attering them on a sear
h domain dened by the new upper and lower boundaries. In the following subse
tions, the above three elements of the RegPSO strategy will be presented; interested readers are en
ouraged to 
onsult Evers and Ben Ghalia [23℄ for more detail. All 
onstants presented in the next subse
tions were found by Evers and Ben Ghalia using numeri
al optimization experiments on standard ben
hmark optimization problems. 7.2.1 Premature 
onvergen
e dete
tion Convergen
e is dete
ted when all parti
les 
onverge to an optimal solution and start to lose their momentum; this loss of momentum prevents es
aping from the neigh- borhood of a high quality lo
al optimal solution. Premature 
onvergen
e 
an be dete
ted when the greatest distan
e between a swarm parti
le and the global best point −→g falls below a 
ertain threshold. This was suggested by Van den Bergh and adopted in Ever's work with su

ess [9, 23℄. 93 For multidimensional optimization problem of size n andm parti
les, the de
ision variables −→x = [ x1 x2 x3 · · · xn ]T have initial upper and lower limits −→ x0U =[ x1U x2U x3U · · · xnU ]T and −→ x0L = [ x1L x2L x3L · · · xnL ]T . First, the range is dened as −−−→ range (Ω) = −→xU −−→xL (7.3) Also the normalized swarm radius 
an be found a

ording to δ norm = maxi=1...m ‖−→xi (k)−−→g (k)‖ ‖−−−→range (Ω)‖ (7.4) Convergen
e is dete
ted when the normalized swarm radius δ norm drops below ε = 1.1 · 10−4. 7.2.2 Redening de
ision variables range On
e 
onvergen
e is dete
ted, the upper and lower limits of the de
ision variables for the regrouped swarm are 
omputed. The initial range of parti
le j 
an be found using its upper limit x0jU and lower limit x 0 jLas rangej ( Ω0 ) = x0jU − x0jL (7.5) Upper and lower limits of the de
ision variables are 
hanged when 
onvergen
e is dete
ted in order to s
atter swarm parti
les in a ball surrounding the 
onverged solution. Changing the range of ea
h de
ision variable j in a regrouping stage r is done a

ording to the following formula [23℄: rangej (Ω r) = max ( rangej ( Ω0 ) , 1.2 × 104 max i=1...m ∣∣∣∣−−→xr−1i,j −−−→gr−1j ∣∣∣∣) (7.6) When the 
urrent range of design variables rangej is su
iently small, the optimiza- tion problem is 
onsidered 
onverged. 7.2.3 Regrouping and position 
lipping Regrouping parti
les in a region surrounding the best position −→g not only moves them away from the lo
al minimum but allows them to re
over their momentum, 94 whi
h had been lost due to 
onvergen
e. Regrouping is done a

ording to the fol- lowing formula −→xi = −→g + [−→ri ]T · −−−−−−−→ range (Ωr)− 1 2 −−−−−−−→ range (Ωr) (7.7) where −−−−−−−→ range (Ωr) = [range1 (Ω r) , range2 (Ω r) , . . . , rangen (Ω r)]T . Equation 7.7 requires position 
lipping to keep the parti
les within the 
on- strained design spa
e, i.e xrj,L ≦ xi,j ≦ x r j,U . The upper and lower limits of the 
omponents of parti
le's de
ision ve
tor 
an be found as xrj,L = max ( x0j,L, g r−1 j − 1 2 rangej (Ω r) ) (7.8) xrj,U = min ( x0j,U , g r−1 j + 1 2 rangej (Ω r) ) 7.3 Hybrid SQP-RegPSO Te
hnique Swarm intelligen
e starts with a set of randomly generated position points in the design spa
e and denes its initial best lo
ation, −→g , as the point of best tness among the swarm parti
les; this best point a
t as a gravity 
enter in the design spa
e and attra
ts other parti
les to it; if the initial best point is of high quality, this redu
es the total 
omputational eort to nd the global optimal solution. The te
hnique we are proposing here for hybridization of sequential quadrati
 programing (SQP) and parti
le swarm optimization is to use SQP to nd the initial high quality −→g before using RegPSO te
hnique. The pseudo-
ode of the hybrid s
heme is summarized in Algorithm 7.2. 7.4 Optimization Test Cases This se
tion will examine the e
ien
y of the proposed hybrid s
heme and 
ompare it with the original RegPSO te
hnique. To do so, a drag optimization that seeks a sho
k free prole with CL = 0.4 at M = 0.8 and angle of atta
k α = 1.5 ◦ is 
arried out. Two distin
t starting points were sele
ted to explore the existen
e of a unique 95 Algorithm 7.2 Hybrid SQP-RegPSO pseudo 
ode For ea
h parti
le in the swarm Initialize parti
le position −→xi Set its best position to be its initial position: −→pi = −→xi Initialize its velo
ity −→ Vi = 0 End Set the global best position: −→g = −→pi where F (−→pi ) < F (−→pj ) for all i 6= j, or the SQP optimal solution, whi
hever is better Set δ norm = 1 Set t = 1 Do For ea
h parti
le i Cal
ulate its velo
ity −→ Vi Update its position −→xi Clip position su
h that xL,j ≤ xi,j ≤ xU,j Cal
ulate its obje
tive fun
tion value (tness) F (−→xi) If the 
urrent tness value is better than its best tness (F (−→xi) < F (−→pi )) Set the 
urrent parti
le position to be its best position: −→pi = −→xi End if End for Sele
t the best parti
le obje
tive fun
tion and position to be the global best if it is better than the stored value t = t+ 1 Compute δ norm using Eq 7.4 If δ norm ≤ ε Compute the range of ea
h de
ision variable range ( Ωrj ) S
atter the parti
les around −→g a

ording to Eq 7.7 Apply position 
lipping with limits dened a

ording to Eq 7.8 End if While (δ norm > ε and t ≤ max iterations). 96 global best 8 ; the rst is to start from NACA 0012, the se
ond one is to start from NACA 00083 airfoil. 9 These 
ases were 
hosen be
ause experiments showed that SQP optimization with these two starting points 
onverged to distin
t optimal solutions with signi
antly dierent obje
tive fun
tion values. All the 
omputations here are se
ond order as we aiming to nd a global optimization method with reasonable 
omputational 
ost, not to 
ompare se
ond and fourth order optimal results. 7.4.1 Constrained drag optimization of the NACA 0012 In this test 
ase, drag optimization of the NACA 0012 at a lift 
oe
ient of CL = 0.4, M = 0.8 and angle of atta
k α = 2◦ is presented. The obje
tive fun
tion 
an be written as F = CD + 10 · (CL − 0.4)2 (7.9) SQP optimization is used rst to nd a lo
al optimum and thus provide an initial high quality −→g as a starting point for the RegPSO s
heme. The optimal prole found by SQP optimization with BFGS Hessian approximation was not sho
k free. RegPSO is used after SQP and su

essfully rea
hed a sho
k free prole. In this test 
ase, the swarm size was 10 parti
les and the number of de
ision variables is 18. The rst phase of the hybrid s
heme (SQP optimization) 
osts 92 CFD simulations while the se
ond phase 
osts 490 CFD simulations whi
h makes the total 
ost 582 CFD simulations to rea
h sho
k free prole. Using RegPSO alone as an optimization te
hnique with a maximum of 1000 CFD simulations does not result in a sho
k free prole, as the optimization has not yet 
onverged to a global optimum. Fig 7.4.1 shows the pressure distribution of the original NACA 0012 at transoni
 
onditions, the SQP optimal prole, the RegPSO optimal shape, and the hybrid sho
k free prole. Drag has been redu
ed from 150 drag 
ounts at the initial lo
al optimum found by SQP optimization to 8.89 drag 
ounts at the nal global optimum. Figure 7.4 shows the surfa
e pressure distribution of the original NACA 0012 air- foil, and the optimized prole using SQP with BFGS approximate Hessian, RegPSO, 8 In the 
ase of a nite region in the design spa
e for whi
h the ow is sho
k free, there may be a region of global best solutions diering only in the amount of dis
retization error present. 9 That is, an airfoil from the 4-digit NACA series with 8.3% maximum thi
kness. 97 NACA 0012 Original BFGS Optimal RegPSO 1000 CFD simulations Hybrid 492 Total CFD simulations Figure 7.3: Pressure Field of the original NACA 0012 and the optimized proles Figure 7.4: Surfa
e pressure distribution of original NACA 0012 and the optimized proles 98 Figure 7.5: Airfoil shapes of original NACA 0012 and the optimized proles and hybrid BFGS-RegPSO te
hniques. Figure 7.5 shows the dierent optimized air- foil proles resulted from the three optimization te
hniques, note that the initial global best point of the hybrid s
heme is the SQP optimal shape. To study the fun
tion variation between the initial best known shape, whi
h results from SQP op- timization, and the nal best shape of the hybrid s
heme, the ve
tor that 
onne
ts these two points in design spa
e was divided into 100 steps and the obje
tive fun
- tion value at ea
h of these points 
omputed. Figure 7.7 shows that there is a des
ent dire
tion (the ve
tor that 
onne
ts the initial and nal −→g ) at the SQP optimal point but the SQP optimizer 
an not nd it. Possible reasons for that in
lude: 1. An error in the 
omputed adjoint gradient. 2. The SQP gradient and the ve
tor that 
onne
ts the initial and nal −→g are perpendi
ular due to the inuen
e of numeri
al noise and of the behavior of the Euler ow model physi
s in the SQP gradient as the des
ent dire
tion resulting from adjoint 
omputations tends to minimize the drag by shifting the sho
k wave forward so it o

urs at lower Ma
h number to redu
e its strength, while this me
hanism to redu
e the drag is not stri
tly followed by RegPSO s
heme . 99 Figure 7.6: Fun
tion minimization iterations in the se
ond phase of the hybrid s
heme (RegPSO) with a NACA 0012 airfoil as the starting point. To determine the real reason, we 
omputed the adjoint gradient and the nite dif- feren
e gradient at the SQP lo
al optimum, then 
omputed the angles between the ˆ Adjoint 
omputed gradient ve
tor. ˆ Finite dieren
e 
omputed gradient ve
tor. ˆ Ve
tor that 
onne
ts the SQP (lo
al) optimum and the hybrid (global) opti- mum. Table 7.1 shows the angles between those three ve
tors. The adjoint and nite dieren
e gradient dire
tions are in ex
ellent agreement but are almost perpendi
ular to the third ve
tor that 
onne
ts lo
al and global optima. This observation supports the hypothesis that the inuen
e of numeri
al noise and Euler ow model physi
s prevents nding that des
ent dire
tion; however, this needs to be veried by another test 
ase. 100 0 20 40 60 80 100 120 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 BFGS Optimal to Hybrid optimal steps F BFGS optimal Hybrid optimal Figure 7.7: Obje
tive fun
tion values between SQP and hybrid optimal points (start- ing geometry: NACA 0012). Angle in degrees Adjoint and nite dieren
e gradients 1.78 Adjoint gradient and SQP-hybrid ve
tor 86.6 Finite dieren
e gradient and SQP-hybrid ve
tor 86.7 Table 7.1: Angles between adjoint, FD., and SQP-hybrid ve
tors 101 7.4.2 Constrained drag minimization of NACA 00083 In this test 
ase, the same drag optimization problem is studied at the same transoni 
onditions, but starting from a thinner airfoil. Again, SQP 
ould not eliminate the sho
k wave on the airfoil's upper surfa
e, but a sho
k free optimal shape is obtained using RegPSO in the se
ond phase as shown in Fig 7.8. Figure 7.9 shows a 
omparison of the surfa
e pressure of the original NACA 00083, SQP with BFGS approximate Hessian (lo
al) optimum, and hybrid BFGS- RegPSO (global) optimum results. The BFGS optimal prole weakens the sho
k strength to redu
e the drag, but the hybrid s
heme eliminates the sho
k wave 
om- pletely. Figure 7.10 shows a prole 
omparison of the original NACA 00083, BFGS optimal, and the hybrid optimal prole. The drag 
oe
ient is redu
ed from 32 drag 
ounts at the BFGS optimal solution to 8.89 drag 
ounts at the hybrid optimal solution. The optimal value of the drag 
oe
ient found by the hybrid s
heme in this test 
ase is the same as the value obtained in the previous 
ase and in both 
ases, it is of order of the dis
retization error. The total number of CFD simulations is 1080, 80 for the SQP phase and 1000 for the RegPSO phase. Again, to study the obje
tive fun
tion variation between the BFGS (lo
al) opti- mum, and the nal (global) optimum of the hybrid s
heme, we took 100 steps on the ve
tor that 
onne
ts these points and 
ompute the obje
tive fun
tion values. Figure 7.12 shows the variation in the obje
tive fun
tion value along that ve
tor, whi
h is not a des
ent ve
tor at the BFGS optimal point. Also, there is numeri
al noise in the obje
tive fun
tion values at some stations along that ve
tor, in addition to an ap- pearan
e then vanishing of a sho
k wave. As in the previous test 
ase, we 
omputed the obje
tive fun
tion gradient at the initial gBest using adjoint and nite dieren
e strategy and 
al
ulated its angle with the ve
tor that 
onne
ts the initial and the nal best position ve
tors −→g obtained using the hybrid s
heme. Table7.2 shows that the ve
tor 
onne
t the initial and nal best solutions is almost perpendi
ular to the adjoint and nite dieren
e gradients; it shows also that adjoint and nite dieren
e gradients are a very good mat
h. The orthogonality of these ve
tors was previously noted in the last test 
ase as well. This strongly suggests that numeri
al noise in the gradient and the inuen
e of the ow model physi
s prevents nding a des
ent 102 NACA 00083 original BFGS Optimal Hybrid Optimal Figure 7.8: Pressure surfa
e of NACA 00083 and the optimized proles Angle in degree Adjoint and Finite dieren
e 0.27 Adjoint and SQP-hybrid ve
tor 86.93 Finite dieren
e and SQP-hybrid ve
tor 86.87 Table 7.2: Angles between adjoint, FD, and BFGS-hybrid best solutions ve
tors dire
tion that leads from the SQP lo
al optimum to the sho
k free global optimum. 7.4.3 Aerodynami
 and thi
kness 
onstraint drag minimization of NACA 0012 In this test 
ase thi
kness is 
onstrained to be tc = 0.1 at xtc = 0.35, in addition to the lift 
onstraint CLc = 0.4 at transoni
 
onditions, M = 0.8 and α = 1.5 ◦ . SQP with BFGS Hessian approximation optimization te
hnique is rst used, the obje
tive 103 Figure 7.9: Surfa
e pressure of NACA 00083 and the optimized proles Figure 7.10: Prole 
omparison of NACA 00083, BFGS optimal, and hybrid optimal airfoils 104 Figure 7.11: Fun
tion minimization iterations in the se
ond phase of the hybrid s
heme (RegPSO) of NACA 00083 start 0 20 40 60 80 100 120 0 0.005 0.01 0.015 0.02 0.025 BFGS optimal to hybrid optimal steps F BFGS optimal Hybrid optimal Figure 7.12: Obje
tive fun
tion values between BFGS and hybrid optimal points, starting geometry NACA 00083 105 fun
tion to be minimized is F = CD + 10 (CL − 0.4)2 (7.10) subje
t tot (xtc) = tc The thi
kness 
onstraint 
an be 
ast as an equality 
onstraint equation using the surfa
e parametrization matrixQ2 [AQ1] † ; for this test 
ase, xc < L1, so the thi
kness 
onstraint equation 
an be written as t (xtc) = P1(xtc)− P3(xtc) = tc (7.11) ∴ ( a0 √ xc + a1xc + a2x 2 c + a3x 3 c )− (c0√xc + c1xc + c2x2c + c3x3c) = tc Re
all that the relationship of polynomial 
oe
ients and the design variables 
an be found using Q2 [AQ1] † as ai = ∑( Q2 [AQ1] † ) i,j yj (7.12) Equation (7.12) 
an be used with Eq (7.11) to derive an algebrai
 linear equation that relates the thi
kness 
onstraint tc to design variables dire
tly; this last equation 
an be used as a 
onstraint equation added to the optimization problem and this 
onstraint will be satised using a Lagrange multiplier. For the se
ond phase of the hybrid s
heme, the thi
kness 
onstraint is satised using a penalty term in the obje
tive fun
tion; the obje
tive fun
tion for RegPSO phase 
an be written as F = CD + 10 (CL − 0.4)2 + (t (xtc)− tc)2 (7.13) The penalty method is used to approximately satisfy the thi
kness 
onstraint term due to the la
k of Lagrange multiplier in non-gradient based optimization te
hniques. Drag 
oe
ient have been redu
ed from 230 to 120 drag 
ounts in the rst phase, then redu
ed to 9 drag 
ounts using the hybrid s
heme, a sho
k free prole is obtained as shown in Figs (7.13, 7.14) ; this sho
k free prole is obtained with 900 CFD simulations in total. 106 NACA 0012 original Hybrid Optimal Figure 7.13: Pressure surfa
e of NACA 0012 and the optimized proles with thi
kness 
onstraint Figure 7.14: Surfa
e pressure of NACA 0012 and the optimized prole with thi
kness 
onstraint 107 NACA 00083 original Hybrid Optimal Figure 7.15: Pressure surfa
e of NACA 00083 and the optimized proles with thi
k- ness 
onstraint 7.4.4 Aerodynami
 and thi
kness 
onstraint drag minimization of NACA 00083 This test 
ase is a repeat of the last one, but the starting geometry is dierent. The starting geometry is NACA 00083 whi
h does not satisfy the thi
kness 
onstraint. The SQP optimizer is used rst;it satises the thi
kness 
onstraint after the rst few iterations, then it starts to minimize the drag. Drag is redu
ed from 60 drag 
ounts to 40 drag 
ounts in the rst phase. The se
ond (RegPSO) phase redu
ed the drag to 9 drag 
ounts and produ
ed a sho
k free optimal shape, as shown in Figs. 7.15 and 7.16. Comparison of the optimized proles of the last two test 
ases are shown in Fig. 7.17. The two proles are very 
lose in upper surfa
e shape, whi
h is tightly 
onstrained by the need to eliminate the sho
k wave; there is more dieren
e in the lower surfa
e on whi
h the ow is subsoni
. The last observation suggests that for invis
id ow simulations, there is a region of global optimal solution rather than one point global optimum. More pre
isely, on
e a sho
k-free solution is obtained, the obje
tive fun
tion is too at to distinguish reliably between the drag of sho
k-free solutions; this is not surprising, 
onsidering that the drag for sho
k-free invis
id ows is due solely to dis
retization error. In all test 
ases, the new hybrid s
heme was able to nd a sho
k free airfoil prole regardless of the starting points. Only ten parti
les were using in the swarm to s
an the design spa
e, whi
h is roughly half of the number of design variables. This is lower than the re
ommended population size required by GA to 
onverge well. The 108 Figure 7.16: Surfa
e pressure of NACA 00083 and the optimized prole with thi
k- ness 
onstraint optimization 
omputational 
ost is four to nine times the 
ost of SQP optimization. Using SQP optimization with ten dierent starting geometries might lead to a sho
k free prole, but there is no guarantee for that. As the results above show, there is a high 
han
e of getting trapped in lo
al minima or suering a breakdown in the optimization pro
ess be
ause of solution noise. Table 7.3 shows the ee
t of the thi
kness penalty weight on the optimization results using se
ond order a

urate 
omputations. It shows that a thi
kness penalty weight of 10 is su
ient to satisfy the thi
kness geometri
 
onstraint, but it will in
rease the value of CD opt and may prevent obtaining a sho
k free optimal prole. Figure Shows a 
omparison of the optimal proles using two dierent thi
kness penalty weights; the maximum thi
kness of the optimized airfoils were 0.077 and 0.095 respe
tively for thi
kness penalty of 1 and 10. Therefore it is re
ommended to use a thi
kness penalty of order of 10 in order to satisfy a maximum thi
kness 
onstraint. Figure 7.18 shows a 
omparison of the optimal pressure distribution of the two optimal proles; a sho
k free pressure distribution is obtained using thi
kness penalty weight of 1. The drag 
oe
ient for this prole is 9 drag 
ounts whi
h is 109 Figure 7.17: Prole 
omparison of the hybrid s
heme optimal proles, starting ge- ometries NACA 0012 and NACA 00083 of order of dis
retization error. It shows also that with thi
kness penalty weight 10, there is a weak 
ompression wave on the airfoil upper surfa
e whi
h 
auses an in
rease in the drag 
oe
ient to 17 drag 
ounts. Thi
k. penalty F opt CD opt CL opt t (xtc) Lift penalty Thi
k. penalty weight 1 0.00134 0.00088 0.3997 0.077 0.000001 0.00048 10 0.00208 0.0017 0.4005 0.095 0.0000026 0.00025 Table 7.3: Thi
kness penalty weight impa
t on the optimization results of NACA 00083 using hybrid s
heme 110 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1.5 −1 −0.5 0 0.5 1 1.5 x/c − Cp   Penalty weight 10 Penalty weight 1 Figure 7.18: Thi
kness penalty weight impa
t on the Optimal pressure distribution of NACA 00083 using hybrid optimization s
heme. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 0.04 x/c y/ c   Thickness penalty 1 Thickness renalty 10 Figure 7.19: Thi
kness penalty weight impa
t on the Optimal prole with NACA 00083 starting geometry using hybrid optimization s
heme. 111 Algorithm 7.3 Pseudo 
ode of RegPSO-SQP hybrid T  optimization s
heme 1. Use the initial geometry as the initial global best −→g . 2. Randomly generate swarm. 3. Add the initial geometry as a swarm member. 4. Use RegPSO strategy with stopping 
riterion of max 900 CFD simulations to nd an estimated global best −→g . 5. Use SQP with starting point −→g to nd better global best solution. 7.5 Comparing SQP-RegPSO with RegPSO-SQP optimization strategies I have presented results of the proposed hybrid optimization s
heme in whi
h SQP optimization is used to obtain a reasonably good upper bound on the global optimum. Based on the results of previous test 
ases, this s
heme seems to be a promising global optimization s
heme. Another hybrid s
heme in whi
h RegPSO pre
edes SQP optimization is possible. In this se
tion I will 
ompare the optimization results of the last ve optimization problems using the proposed hybrid s
heme with the results of RegPSO-SQP s
heme; I will refer to the latter s
heme as hybrid T . Algorithm 7.3 shows a pseudo 
ode of hybrid T optimization. Both s
hemes has two phases, SQP phase and RegPSO phase, but they dier in the order of the phases. Tables 7.4 and 7.5 summarize the optimization results for both s
hemes. Figures 7.207.24 show a 
omparison of the obtained pressure distribution of the optimal shape obtained using these s
hemes. These 
omparisons show that the hybrid s
heme is better than the hybrid T for all 
ases; in Case II, the s
hemes are nearly identi
al in drag 
oe
ient though dierent in nal shape. This was expe
ted be
ause in hybrid T if the rst (RegPSO) phase 
ould not rea
h the basin of attra
tion of the global minimum, the se
ond (SQP) phase will not be able to nd it. Tables 7.4, 7.5 show a 
omparison of Phase I and Phase II optimal obje
tive fun
tion value for both hybrid and hybrid T for the ve test 
ases in this 
hapter. As a hybridization quality measure, the ratio of the optimal obje
tive fun
tion phase 112 Case Phase I Opt (SQP) Phase II Opt (RegPSO) Phase II/I I 10 0.0150 0.0009 0.0600 II 11 0.0032 0.0009 0.2813 III 12 0.0120 0.0013 0.1083 IV 13 0.0040 0.0013 0.3250 V 14 0.0040 0.0021 0.525 Table 7.4: Optimization results of the hybrid s
heme in Phase I, II Case Phase I Opt (RegPSO) Phase II Opt (SQP) Phase II/I I 0.0137 0.0049 0.358 II 0.0019 0.0011 0.579 III 0.0125 0.0125 1.0 IV 0.0017 0.0017 1.0 V 0.0052 0.0052 1.0 Table 7.5: Optimization results of the hybrid T s
heme in Phase I, II II to Phase I optimal is used and the smaller this value the more su

essful the hybridization sequen
e is. The tables show also that the sequen
e used in hybrid is better than hybrid T . 113 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 x/c − Cp   RegPSO−SQP SQP−RegPSO Figure 7.20: Hybrid and hybrid T optimal pressure distribution for NACA 0012 (Case I). 114 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 x/c − Cp   RegPSO−SQP SQP−RegPSO Figure 7.21: Hybrid and hybrid T optimal pressure distribution for NACA 00083 (Case II). 115 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 x/c − Cp   RegPSO−SQP SQP−RegPSO Figure 7.22: Hybrid and hybrid T optimal pressure distribution for NACA 0012 with 10% thi
kness 
onstraint and unit thi
kness penalty weight (Case III). 116 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 x/c − Cp   RegPSO−SQP SQP−RegPSO Figure 7.23: Hybrid and hybrid T optimal pressure distribution for NACA 00083 with 10% thi
kness 
onstraint and unit thi
kness penalty weight (Case IV). 117 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 x/c − Cp   RegPSO−SQP SQP−RegPSO Figure 7.24: Hybrid and hybrid T optimal pressure distribution for NACA 00083 with 10% thi
kness 
onstraint and ten thi
kness penalty weight (Case V). 118 Chapter 8 CONCLUSIONS AND RECOMMENDATIONS 8.1 Contributions and Con
lusions This thesis presents the rst study of the use of high order nite volume CFD methods in aerodynami
 optimization. I have developed the rst nite volume based optimization 
ode that uses high order Ja
obian in gradient 
omputations using the adjoint te
hnique. A new geometry parametrization te
hnique has been developed and presented; this new te
hnique is based on least squares surfa
e tting. To develop an e
ient global optimizer, a new hybrid SQP-RegPSO s
heme has been developed to nd the global optimal solution in the feasible design spa
e. I presented both 
onstrained and un
onstrained aerodynami
 optimization using sequential quadrati
 programming with the use of adjoint method to 
ompute the obje
tive fun
tion gradient. The 
omputations were based on higher order nite volume method on unstru
tured grids. I took advantage of evaluating the exa
t Ja
obian matrix to 
ompute the ow adjoint eld and ow solution sensitivity. The 
omputed gradient based on the adjoint method is an ex
ellent mat
h with the 
orresponding nite dieren
e gradient in subsoni
 ow for both se
ond and fourth order s
hemes. For transoni
 ow, the se
ond order Ja
obian mat
hes well with the 
orresponding order of a

ura
y nite dieren
e gradient; on the other hand, the fourth order Ja
obian does not mat
h the nite dieren
e gradient well when using a limiter. This error in the fourth order gradient is attributed to the high sensitivity of the fourth order Ja
obian matrix to the use of the limiter in transoni ow whi
h tends to ae
t the diagonal dominan
e of the Ja
obian matrix. To ta
kle this problem and to enhan
e the diagonal dominan
e of the fourth order Ja
obian matrix, we used the same non-zero stru
ture of the se
ond order Ja
obian matrix 119 and dropped the rest of the non-zero elements; this tends to in
rease the diagonal dominan
e of the modied fourth order Ja
obian matrix and improves its 
ondition number. I have developed a new approa
h for smooth parametrization of airfoil shapes, based on a least-squares t of the parameters of two polynomials to a set of 
ontrol points. Compared with using all mesh points as design variables, this approa
h redu
es the size of the design spa
e and eliminates os
illations in the shape. We use a semi-torsional spring analogy to deform the mesh grid in the entire ow eld when the surfa
e shape evolves during optimization. It is used also to 
al
ulate the mesh sensitivity terms in the adjoint gradient. To test the developed gradient based optimizer; I have used the optimizer to ta
kle inverse design problems in whi
h only one optimal solution point in the de- sign spa
e exists. Additional transoni
 drag minimization test 
ases have been pre- sented with and without lift 
onstraint. As a summary of test 
ases observations: A subsoni
 inverse design test 
ase shows that both se
ond and fourth order s
hemes rea
hed their target geometry after almost the same number of iterations. The same behavior is observed for transoni
 inverse design test 
ase. This indi
ates that the 
onvergen
e is independent of the order of a

ura
y of the spatial dis
retization. In a drag minimization test 
ase without lift 
onstraint, both se
ond and fourth order s
heme rea
hed their optimum shape after almost the same number of optimiza- tion iterations. The dieren
e of the resulting optimal airfoil shape is small. Both s
hemes redu
ed the drag by almost 50% of its original value but the lift also went down. In a drag minimization test 
ase with lift 
onstraint, the fourth order s
heme was faster to rea
h its optimal shape with 8 optimization iterations (whi
h 
ost 49 CFD simulations), while the se
ond order s
heme took 13 iterations (112 CFD simulations) to rea
h its optimum. Both s
hemes rea
hed their optimum with almost the same wall 
lo
k run time and found nearly identi
al airfoil optimal shapes. Based on all the test 
ases we presented, we 
on
lude that the spatial dis
retization error produ
es a systemati
 error in the obje
tive fun
tion and has a little ee
t on the nal optimal shape. A mesh renement study shows that both se
ond and fourth order s
hemes tend to give the same optimal value for the obje
tive fun
tion as we rene the mesh. 120 Based on the results of this resear
h, I re
ommend the use of a high order method for obje
tive fun
tion value 
omputations and for a

urate predi
tion of for
es and sho
k wave 
apture in transoni
 ows, and the use of a se
ond order method (but based on the obtained forth order solution) to 
ompute the obje
tive fun
tion gradient needed for gradient based optimizer. If the mesh resolution makes se
ond and fourth order simulation results similar, using the se
ond order s
heme for optimization will be e
onomi
al. I have developed also a new hybrid gradient/non-gradient based optimization te
hnique that uses sequential quadrati
 programming with BFGS Hessian approx- imation te
hnique in an initial gradient based optimization phase, followed by the regrouped parti
le swarm optimization method in the non-gradient random sear
h phase. The SQP optimization phase leads to a high quality initial best-solution point that redu
es the total 
omputational eort of the RegPSO to rea
h a sho
k free prole. Unlike with the SQP s
heme, the transoni
 drag optimization with lift 
onstraint of the NACA 0012 and NACA 00083 ta
kled with the hybrid s
heme lead to sho
k free optimal shapes; in both 
ases, the obje
tive fun
tion gradient ve
tor is almost perpendi
ular to the ve
tor that 
onne
ts the SQP optimum and the hybrid opti- mum. This explains why the gradient based optimization 
an not rea
h the global optimal, and also shows that its optimal solution depends on the starting geometry. The obje
tive fun
tion behavior along the ve
tor that 
onne
ts the SQP and the hybrid optima shows that numeri
al noise may have a harmful impa
t on gradient based optimization as well. Transoni
 drag optimization of NACA 0012 and NACA 00083 have been ta
k- led with both aerodynami
 and geometri
 
onstraints. For RegPSO, the geometri 
onstraint is satised using penalty method; in both 
ases, sho
k free proles have been obtained. The last test 
ase (drag redu
tion of NACA 00083 with a thi
kness 
onstraint) shows that the hybrid s
heme was able to nd a sho
k free prole despite in
reasing the thi
kness during optimization. The 
omputational expense of the hy- brid s
heme is four to nine times more than SQP optimization, but in all 
ases, a global optimal (sho
k free) prole obtained, and in aerospa
e industry, this is worth the in
rease in 
omputational eort. 121 8.2 Future Work We intend to apply non gradient-based optimization te
hniques like parti
le swarm algorithm to study the ee
t of CFD s
heme order of a

ura
y on the nal optimal shape in the near future. Another obvious extension of the 
urrent work is aero- dynami
 optimization using RANS with turbulen
e models. The ow solver should solve the turbulen
e equations 
oupled with the mean ow equations; otherwise, the ow sensitivity values will be less ina

urate. Thus, for a three dimensional ow simulation with k − ǫ turbulen
e model, the blo
k size of the Ja
obian matrix is 7× 7. Fully 
oupled solvers have also been shown to 
onverge more e
iently. For three-dimensional problems, the geometry parametrization should be re- pla
ed with a perturbation parametrization (parametrization of the 
hange in the shape). This will redu
e the required number of design variables whi
h is very important for three dimensional optimization problems. Importantly, perturbation parametrization also prevents the 
hange in the initial pressure distribution when parametrizing the geometry itself. Exploring the ee
t of the geometry parametrization te
hnique used on SQP opti- mization results is needed. We spe
ulate that using the right geometry parametriza- tion te
hnique may lead to a design spa
e in whi
h aerodynami
 
onstraints are simply 
onne
ted and an SQP optimizer 
an nd a globally optimal solution. 122 Bibliography [1℄ Anderson, M.B. The Potential of Geneti
 Algorithms for Subsoni
 Wing De- sign. In AIAA Paper 95-3925, presented at the 1st AIAA Air
raft Engineering, Te
hnology, and Operations Congress, Los Angeles, CA, 1995. [2℄ Anderson, Murray B. Geneti
 Algorithms In Aerospa
e Design: Substantial Progress, Tremendous Potential. In NATO/Von-Karmen-Institute Workshop on Intelligent System, 2002. [3℄ Anderson, W. Kyle, and Bonhaus, Daryl L. Aerodynami
 Design on Unstru
- tured Grids for Turbulent Flows. Te
hni
al report, NASA Langley Resear
h Center, Hampton, Virginia, 1997. [4℄ Avriel, M. Nonlinear Programming Analysis and Methods. Dover Publishing. ISBN 0-486-43227-0., 2003. [5℄ Azab, M. B., and Ollivier-Goo
h, C. F. Higher Order Two Dimensional Aerody- nami
 Optimization Using Unstru
tured Grids and Adjoint Sensitivity Compu- tations. In AIAA-2010-1430, 48th AIAA Aerospa
e s
ien
es meeting, Orlando, FL, USA, 4-7Jan 2010, Orlando FL, 2010. [6℄ Azab, M. B., and Ollivier-Goo
h, Carl. Constrained and Un
onstrained Aerody- nami
 Quadrati
 Programming Optimization Using High Order Finite Volume Method and Adjoint Sensitivity Computations. In 49th Aerospa
e s
ien
e meet- ing, Orlando FL, AIAA-2011-0183, 2011. [7℄ Batina, J. T. Unsteady Euler Airfoil Solutions Using Unstru
tured Dynami Meshes. AIAA Journal ,, 28 No 8.:13811388, 1990. 123 [8℄ Beall, M. W., Walsh, J., and Shephard, Mark S. A

essing CAD Geometry for Mesh Generation. In Pro
eedings of the 12th International Meshing Roundtable, Sandia National Laboratories, pages 20033030, 2003. [9℄ Bergh, F., and Engelbre
ht, A. A New Lo
ally Convergent Parti
le Swarm Optimiser. In Conferen
e on Systems, Man and Cyberneti
s, 2002. [10℄ Borrel, P., and Rappoport, A. Simple Constrained Deformations for Geometri Modeling and Intera
tive Design. ACM Transa
tions on Graphi
s, 13:137155, 1994. [11℄ Brezillon, J. and Gauger, N. R. 2D and 3D Aerodynami
 Shape Optimisation Using the Adjoint Approa
h. Aerospa
e S
ien
e and Te
hnology, 8:715727, 2004. [12℄ Broyden, C. G. The Convergen
e of a Class of Double-rank Minimization Algo- rithms. Journal of the Institute of Mathemati
s and Its Appli
ations, 6:7690, 1970. [13℄ Byrd,R. H., Lu, P., and No
edal, J. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM Journal on S
ienti
 and Statisti
al Comput- ing, 16, 5:11901208., 1995. [14℄ Cavallo, P. A., Hosangadi, A., Lee, R. A., and Dash, S. M. Dynami
 Unstru
- tured Grid Methodology with Appli
ation to Aero/Propulsive Flow Field. In AIAA Paper 1997-2310, 2002. [15℄ Chandrashekarappa, P., and Duvigneau, R. Metamodel-Assisted Parti
le Swarm Optimization and Appli
ation to Aerodynami
 Shape Optimization. Re- sear
h Report RR-6397, INRIA, 2007. [16℄ Cler
, M., and Kennedy, J. The Parti
le Swarm - Explosion, Stability, and Convergen
e in a Multidimensional Complex Spa
e. Evolutionary Computation, IEEE Transa
tions on, 6(1):5873, 2002. [17℄ Consentino, G., and Holst, T. Numeri
al Optimization Design of Advan
ed Transoni
 Wing Congurations. In AIAA 85-0424, 1985. 124 [18℄ Cook, P.H., M
Donald M.A., and Firmin, M.C.P. . Aerofoil RAE 2822: Pres- sure Distributions, and Boundary Layer and Wake Measurements. In Experi- mental Data Base for Computer Program Assessment, AGARD Report AR 138. AGARD, 1979. [19℄ Dadone, Andrea, Mohammadi, Bijan, and Petruzzelli, Ni
ola. In
omplete Sen- sitivities and BFGS Methods for 3D Aerodynami
 Shape Design. Resear
h Report RR-3633, INRIA, 1999. Projet M3N. [20℄ Duvigneau, R. Aerodynami
 Shape Optimization with Un
ertain Operating Conditions using Metamodels. Resear
h Report RR-6143, INRIA, 2007. [21℄ Eberhart, R. C., and Kennedy, J. A New Optimizer Using Parti
le Swarm The- ory. In Sixth International Symposium on Mi
ro Ma
hine and Human S
ien
e, Nagoya, Japan, pp. 39-43., 1995. [22℄ Engelbre
ht, A. P. Fundamentals of Computational Swarm Intelligen
e. John Wiley and Sons, 2006. [23℄ Evers, G. I., and Ben Ghalia, M. Regrouping Parti
le Swarm Optimization a New Global Optimization Algorithm with Improved Performan
e Consisten
y A
ross Ben
hmarks. In SMC IEEE International Conferen
e, pages 39013908, o
t. 2009. [24℄ Evers, George. An Automati
 Regrouping Me
hanism to Deal with Stagna- tion in Parti
le Swarm Optimization. Master's thesis, University of Texas-Pan Ameri
an, 2009. [25℄ Farhat, C., Degand, C., Koobus, B., and Lesoinne, M. Torsional Springs for Two-Dimensional Dynami
 Unstru
tured Fluid Meshes. Computer methods in applied me
hani
s and engineering, 163(1-4):231245, 1998. [26℄ Frink, N. T. Three-Dimensional Upwind S
heme for Solving The Euler Equa- tions on Unstru
tured Tetrahedral Grids. PhD thesis, Virginia Polyte
hni
 In- stitute and State University, 1991. 125 [27℄ Frink, N. T., Parikh, P., and Pirzadeh, S. A Fast Upwind Solver for The Euler Equations on Three-Dimensional Unstru
tured Meshes. In AIAA paper 91-0102, 1991. [28℄ Gage, P., and Kroo, I. Development of the Quasi-Pro
edural Method for Use in Air
raft Conguration Optimization. In AIAA Paper No. 92-4693, Presented at the Fourth AIAA/USAF/NASA/OAI Symposium on Multidis
iplinary Analysis and Optimization, 1992. [29℄ Gage, P., and Kroo, I. A Role of Geneti
 Algorithms in a Preliminary Design Environment. In AIAA Paper No. 93-3933, 1993. [30℄ Goldberg, David E. Geneti
 Algorithms in Sear
h, Optimization, and Ma
hine Learning. Addison-Wesley Publishing Company, 1989. [31℄ Goldfarb, D. Fa
torized Variable Metri
 Methods for Un
onstrained Optimiza- tion. Mathemati
s of Computation, 30(136):796811, 1976. [32℄ Gregg, R.D., and Misegades, K.P. Transoni
 Wing Optimization Using Evolu- tion Theory. In presented at the AIAA 25th Aerospa
e S
ien
es Meeting, AIAA 87-0520, 1987. [33℄ Hart, Willian E. . Adaptive Global Optimization with Lo
al Sear
h. PhD thesis, University of California, San Diego, 1994. [34℄ Hassan, R., Cohanim, A. , and de We
k,O. A Comparison of Parti
le Swarm Optimization and The Geneti
 Algorithm. In AIAA-2005-1897, page 1897, 2005. [35℄ Hi
ks, R. M. , and Henne, P. A. . Wing Design by Numeri
al Optimization. Journal of Air
raft, 15:407412., 1978. [36℄ Ira H. Abbott and Albert Edward Von Doenho. Theory of Wing Se
tions. Dover, 1959. [37℄ Jameson, A. Aerodynami
 Design via Control Theory. Journal of S
ienti Computing, 3:233260, 1988. 10.1007/BF01061285. 126 [38℄ Jameson, A. A Perspe
tive on Computational Algorithms for Aerodynami Analysis and Design. Progress in Aerospa
e S
ien
es, 37(2):197243, 2001. doi: DOI: 10.1016/S0376-0421(01)00004-5. [39℄ Jameson, A. Aerodynami
 Shape Optimization Using The Adjoint Method. VKI Le
tur seriese, 2003. [40℄ Jameson, A. Optimum Transoni
 Wing Design Using Control Theory. In IU- TAM Symposium held in Guttingen, Germany 2-6 September 2002, page 253. Springer Netherlands, 2003. IUTAM Symposium Transsoni
um IV: pro
eedings of the IUTAM Symposium held in Guttingen, Germany 2-6 September 2002. [41℄ Jameson, A. , and Reuther, J. Control theory based airfoil design using Eu- ler equations. In AIAA/USAF/NASA/ISSMO Symposium on Multidis
iplinary Analysis and Optimization, 1994. [42℄ Jameson, A., and Kim, S. Redu
tion of The Adjoint Gradient Formula in The Continuous Limit. AIAA paper, 40:2003, 2003. [43℄ Jameson, A., Pier
e,N. A., and Martinelli, L. Optimum Aerodynami
 Design Using The Navier-Stokes Equations. In 35th AIAA Aerospa
e S
ien
es Meeting & Exhibit, Reno, NV, 1997. [44℄ Jang, M., and Lee, J. Geneti
 Algorithm Based Design of Transoni
 Air- foils Using Euler Equations. In AIAA Paper 2000-1584, Presented at the 41st AIAA/ASME/ASCE/AHS/ASC Stru
tures, Stru
tural Dynami
s, and Materi- als, 2000. [45℄ Karthik, M. Appli
ation Of The Dis
rete Adjoint Method To Coupled Multidis
i- plinary Unsteady Flow Problems For Error Estimation And Optimization. PhD thesis, Department of Me
hani
al Engineering, The University of Wyoming, 2009. [46℄ Kennedy, J., and Eberhart, R. Parti
le Swarm Optimization. In Pro
eedings of IEEE International Conferen
e on Neural Networks. IV. pp. 1942-1948, 1995. [47℄ Kroo, I. Drag Due to Lift: Con
epts for Predi
tion and Redu
tion. Annual Review of Fluid Me
hani
, 33:587617, 2001. 127 [48℄ Kulfan, B. M. Re
ent Extensions and Appli
ations of the CST Universal Para- metri
 Geometry Representation Method. In 7th AIAA Aviation Te
hnology, Integration and Operations Conferen
e, 2007. [49℄ Kulfan, B. M., Bussoletti, J. E. . Fundamental Parametri
 Geometry Represen- tations for Air
raft Component Shapes. In 11th AIAA/ISSMO Multidis
iplinary Analysis and Optimization Conferen
e, 2006. [50℄ Li, W. , Huyse, L. and Padula, S. Robust airfoil optimization to a
hieve drag redu
tion over a range of ma
h numbers. Stru
tural and Multidis
iplinary Op- timization, 24:3850, 2002. 10.1007/s00158-002-0212-4. [51℄ Martineau, D. G., Georgala, J. M. A Mesh Movement Algorithm for High Quality Generalised Meshes. In 42nd AIAA Fluid Dynami
s Conferen
e and Exhibit, AIAA Paper 2004-0614, Reno, NV, 2004. [52℄ Masuda, H., Yoshioka, Y., and Furukawa, Y. Intera
tive Mesh Deformation Us- ing Equality-Constrained Least Squares. Computers and Graphi
s, 30(6):936 946, 2006. [53℄ Mengistu, Temesgen, and Ghaly, Wahid. Aerodynami
 Optimization of Tur- boma
hinery Blades Using Evolutionary Methods and ANN-Based Surrogate Models. Optimization and Engineering, 9:239255, 2008. 10.1007/s11081-007- 9031-1. [54℄ Mi
halak, C. E
ient High-Order A

urate Unstru
tured Finite-Volume Algo- rithms for Vis
ous and Invis
id Compressible Flows. PhD thesis, The University of British Columbia, 2009. [55℄ Mi
halak, C., and Ollivier-Goo
h, C. Dierentiability of Slope Limiters on Unstru
tured Grids. In Fourteenth Annual Conferen
e of the Computational Fluid Dynami
s So
iety of Canada, 2006. [56℄ Mi
halak, C., and Ollivier-Goo
h, C. A

ura
y Preserving Limiter for The High- Order A

urate Solution of The Euler Equations. Journal of Computational Physi
s, 228(23):86938711, 2009. 128 [57℄ Mi
halak, Christopher, and Ollivier-Goo
h, Carl. Globalized Matrix-Expli
it Newton-GMRES for The High-Order A

urate Solution of The Euler Equations. Computers and Fluids, 39:11561167, 2010. [58℄ Mousavi, A., Castonguay, P., and Nadarajah, S. K. Survey of Shape Parame- terization Te
hniques and its Ee
t on Three-Dimensional Aerodynami
 Shape Optimization. In Pro
eedings of the AIAA 18th Computational Fluid Dynami
s Conferen
e. AIAA 2007-3837, 25 - 28 June 2007 2007. [59℄ Mouser, C. R., and Dunn, S. A. Comparing Geneti
 Algorithms and Parti- 
le Swarm optimisation for an Inverse Problem Exer
ise. In Rob May and A. J. Roberts, editors, Pro
. of 12th Computational Te
hniques and Appli- 
ations Conferen
e CTAC-2004, volume 46, pages C89C101, Mar
h 2005. http://anziamj.austms.org.au/V46/CTAC2004/Mous [Mar
h 21, 2005℄. [60℄ Nadarajah, S., and Jameson, A. A Comparison of The Continuous and Dis- 
rete Adjoint Approa
h to Automati
 Aerodynami
 Optimization. AIAA paper, 667:2000, 2000. [61℄ Nejat, A., and Ollivier-Goo
h, C. A High-Order A

urate Unstru
tured Finite Volume Newton-Krylov Algorithm for Invis
id Compressible Flows. Journal of Computational Physi
s, 227(4):25822609, 2008. [62℄ Neme
, M., and Aftosmis, M. J. Adjoint Sensitivity Computations for an Embedded-Boundary Cartesian Mesh Method. Journal of Computational Physi
s, 227(4):27242742, 2008. [63℄ Neme
, M., and Zingg, D. W. . Newton-Krylov Algorithm for Aerodynami Design Using the Navier-Stokes Equations. AIAA Journal, 40:11461154, 2002. [64℄ Nielsen, E., and Anderson, K. Aerodynami
 Design Optimization on Unstru
- tured Meshes Using the Navier-Stokes Equations. Te
hni
al report, NASA Lan- gley Te
hni
al Report Server, 1998. [65℄ Nielsen, Eri
 J., and Kleb, Bil. E
ient Constru
tion of Dis
rete Adjoint Op- erators on Unstru
tured Grids by Using Complex Variables. In 43rd AIAA Aerospa
e S
ien
es Meeting and Exhibit Reno, Nevada, 2005. 129 [66℄ No
edal, J., and Wright Stephen, J. Numeri
al Optimization. Springer New York, 2006. [67℄ Obayashi, S. Geneti
 Algorithm for Aerodynami
 Inverse Optimization Prob- lems. In Geneti
 Algorithms in Engineering Systems: Innovations and Appli
a- tions, 1995. [68℄ Ollivier-Goo
h, C. and Van Altena, M. A High-Order-A

urate Unstru
tured Mesh Finite-Volume S
heme for The Adve
tion-Diusion Equation. Journal of Computational Physi
s, 181(2):729752, 2002. [69℄ Ollivier-Goo
h, C., Nejat, A., and Mi
halak, C. On Obtaining and Verifying High-Order Finite-Volume Solutions to the Euler Equations on Unstru
tured Meshes. AIAA Journal, 47(9):21052120, 2009. [70℄ Oyama, A., Obayashi, S., and Nakahashi, K. Transoni
 Wing Optimization Using Geneti
 Algorithm. In AIAA Paper 97-1854, 13th Computational Fluid Dynami
s Conferen
e, 1997. [71℄ Oyama, A., Obayashi, S., and Nakahashi, K. Fra
tional fa
torial design of geneti
 
oding for aerodynami
 optimization. In AIAA Paper 99-3298, 3rd AIAA Weakly Ionized Gases Workshop, 1999. [72℄ Piegl, L., and Tiller, W. The NURBS Book. Springer, Berlin, 1995. [73℄ Reuther, J., Jameson, A., Farmer, J., Martinelli, L., and Saunders, D. Aerody- nami
 Shape Optimization of Complex Air
raft Congurations via an Adjoint Formulation. In 34rd AIAA Aerospa
e S
ien
es Meeting and Exhibit, no. AIAA 1996-0094, 1996. [74℄ Roe, P. L. Approximate Riemann Solvers, Parameter Ve
tors, and Dieren
e S
hemes. Journal of Computational Physi
s, 43(2):357372, 1981. doi: DOI: 10.1016/0021-9991(81)90128-5. [75℄ Saad, You
ef, and S
hultz, Martin H. GMRES: a Generalized Minimal Resid- ual Algorithm for Solving Nonsymmetri
 Linear Systems. SIAM J. S
i. Stat. Comput., 7:856869, July 1986. 130 [76℄ Samareh, Jamshid A. A Survey of Shape Parameterization Te
hniques. In CEAS/AIAA/ICASE/NASA Langley International Forum on Aeroelasti
ity and Stru
tural Dynami
s Williamsburg, VA, Also, NASA/CP-1999-209136, June 1999, pp. 333-343., 1999. [77℄ Samareh, Jamshid A. Geometry and Grid/Mesh Generation Issues for CFD and CSM Shape Optimization. Optimization and Engineering, 6:2132, 2005. 10.1023/B:OPTE.0000048535.08259.a8. [78℄ Sederberg, Thomas W., and Parry, S
ott R. Free-Form Deformation of Solid Geometri
 Models. SIGGRAPH Comput. Graph., 20:151160, August 1986. [79℄ Shi, Y., and Eberhart, R. C. Empiri
al Study of Parti
le Swarm Optimization. In Pro
eedings of Congress onEvolutionary Computation, CEC 99., volume 3, page 1950 Vol. 3, 1999. [80℄ Stein, K., Tezduyar, T., and Benney, R. Mesh Moving Te
hniques for Fluid- Stru
ture Intera
tions with Large Displa
ements. J. Appl. Me
h., 70:5863, 2003. [81℄ Tautges,Timothy J. The Common Geometry Module (CGM): A Generi
, Ex- tensible Geometry Interfa
e. In Pro
eedings of the 9th International Meshing Roundtable, pages 337348, 2000. [82℄ Tezduyar, T. E., Behr, M., Mittal, S., and Johnson, A. A. Computation of Unsteady In
ompressible Flows with The Stabilized Finite Element Methods: Spa
e-Time Formulations, Iterative Strategies and Massively Parallel Imple- mentations. New Methods in Transient Analysis, 246:0, 1992. [83℄ Truong, Anh H., Oldeld, Chad A., and Zingg, David W. Mesh Movement for a Dis
rete-Adjoint Newton Krylov Algorithm for Aerodynami
 Optimization. AIAA Journal, 46, issue 7:16951704, 2008. [84℄ Venkatakrishnan, V. Convergen
e to Steady State Solutions of The Euler Equa- tions on Unstru
tured Grids with Limiters. Journal of Computational Physi
s, 118(1):120130, 1995. 131 [85℄ Venter, G., and Sobiesz
zanski-Sobieski, J. Multidis
iplinary Optimization of a Transport Air
raft Wing Using Parti
le Swarm Optimization. Stru
t Multidis- 
iplinary Optimization, 26:pp 121131, 2004. [86℄ Watt, A., and Watt, M. Advan
ed Animation and Rendering Te
hniques. Addison-Wesley Publishing Company, New York, 1992. [87℄ Zhu, C., Byrd, R. H., Lu, P., and No
edal, J. Algorithm 778: L-BFGS-B: Fortran Subroutines for Large-S
ale Bound-Constrained Optimization. ACM Transa
tions on Mathemati
al Software (TOMS), 23(4):550560, 1997. [88℄ Zymaris,A.S., Papadimitriou, D.I., Giannakoglou, K.C., and Othmer,C. Feasi- bility Study of Constant Eddy-Vis
osity Assumption in Gradient-Based Design Optimization. Journal of Computational Physi
s, 229:52285245, 2010. 132

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0072307/manifest

Comment

Related Items