Aerodynami
Optimization Using HighOrder FiniteVolume 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. Gradientbased 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 gradientindependent 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 leastsquares based geometry parametrization is used to des
ribe airfoil shapes and a semitorsional 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 gradientbased optimization
ould not in general. iii Prefa
e The resear
h ideas and methods explored in the three
oauthored manus
ripts of this thesis are the fruits of a
lose working relationship between Dr Carl OllivierGoo
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 OllivierGoo
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 Gradientbased 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
ewise 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 SQPRegPSO 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 SQPRegPSO with RegPSOSQP 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 SQPhybrid ve
tors . . . . . . . . 101 7.2 Angles between adjoint, FD, and BFGShybrid 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 gradientbased 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 SQPRegPSO pseudo
ode . . . . . . . . . . . . . . . . . . . 96 7.3 Pseudo
ode of RegPSOSQP 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 xvelo
ity ṽ Roe avaraged yvelo
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 ElRamsisy 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. OllivierGoo
h and
o workers [68, 61, 56℄ is used in this resear
h. Optimization s
hemes
an be divided into two main
ategories: gradientbased and nongradientbased s
hemes. Gradientbased s
hemes require
omputing the ob je
tive fun
tion value and its gradient with respe
t to the design variables, while the nongradientbased te
hniques require only
omputing the obje
tive fun
tion value. Gradientbased s
hemes are fast to
onverge to a lo
al optimal point in the design spa
e
ompared to the nongradientbased methods, but the optimal solution found by gradientbased optimization depends on the starting point [33℄. Nongradient 1 based s
hemes
an nd global minimum solution regardless of the starting point, but with larger
omputational
ost
ompared to gradientbased 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 roundo 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 Roeaveraged 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 Roeaveraged 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 Roeaveraged 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 Higherorder a
ura
y is obtained by leastsquares 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 quasiNewton generalization of Eq. 1.6 that in
ludes residualbased 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. Nongradientbased 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 gradientbased s
hemes. 1.2.1 Gradientbased aerodynami
optimization Gradientbased 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 quasiNewton 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℄. Gradientbased 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 gradientbased 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 gradientbased optimization te
h niques like steepest des
ent and quadrati
programming in aerodynami
optimiza tion [17, 31, 28℄. The most expensive part of gradientbased 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 gradientbased optimization Numerous resear
hers have applied gradientbased optimization using the ad joint approa
h in aerodynami
optimization sin
e the early 1990's. Reuther and his
oworkers 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 NavierStokes 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 NavierStokes 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 SpalartAllmaras oneequation 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 nondierentiable 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 lifttodrag ratio of an airfoil, beginning from the NACA 0012 1 [44℄. Oyama et al applied a geneti algorithm with a NavierStokes 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 suboptimal 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 UrbanaChampaign 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 lifttodrag 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 gradientbased and gradientfree optimization s
hemes. To a
hieve these goals, the following
omponents were needed, in addition to the preexisting highorder 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 leastsquares 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 semitorsional 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 highorder 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 gradientbased optimizer using in this thesis is the BFGSbased quasiNewton 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 gradientbased s
heme is unable to. 11 Chapter 2 GEOMETRY PARAMETRIZATION The geometry of engineered obje
ts is dened mathemati
ally in
omputeraided 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 nonsmooth 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
ewise 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
ewise leastsquares 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 ClassShape fun
tionTransformation 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 nonwavy 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 Bspline parametrization [58℄ for a threedimensional lift
onstrained transoni
drag minimization problem. 13 2.1.2 Pie
ewise 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. Bsplines, a general ization of Bezier
urves, were found to be more suitable. A
ubi
Bspline 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 multipoint operating
onditions with a lift
onstraint using spline representation as geometry parametrization te
hnique [50℄. CAD systems typi
ally use nonuniform rational Bspline (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 gradientfree method [53℄. While pie
ewise spline parametrization is wellsuited for twodimensional shapes and simple threedimensional geometries,
omplex shapes require a large number of
ontrol points, redu
ing the ee
tiveness of gradientbased 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 gradientbased 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 Bsplines or even NURBS splines. A related approa
h, the free form deformation algorithm (FFD) treats the geometry as a void in a boxshaped 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 trivariate 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 Bsplines 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 leastsquares based surfa
e parametrization method is presented in this se
 tion; the airfoil surfa
e is parametrized using pie
ewise 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
ewise polyno mials found by a leastsquares 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
hordwise position, L1 is
hordwise 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 leastsquares 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 leastsquares 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 leastsquares system by singular value de
omposition for numeri
al sta bility, z = [AQ2] † c (2.12) where the re
tangular matrix [AQ2] † is the pseudoinverse 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
ewise 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 xstations 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 xstations
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 xstations. 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
hordwise 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 leastsquares 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 leastsquares 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 leastsquares 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 6digit series, laminar ow, supersoni
, and super
riti
al airfoil se
tions. 4 Figures 2.32.10 show the leastsquares 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 xstations 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 16006 1.31 · 10−6 2.03 · 10−5 6 · 10−2 NACA 63412 1.17 · 10−5 5.63 · 10−5 12 · 10−2 NACA 64421 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 semitorsional 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 nonlinear equation and needs to be solved iteratively. 3.2 Testing Mesh Morphing Farhat et al [25℄ demonstrated the robustness of the semitorsional 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 ydire
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 rightangled
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 gradientbased 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 leastsquares system solved in parametrizing the geometry. 34 Chapter 4 GRADIENT CALCULATION USING ADJOINT APPROACH Gradient
al
ulation plays a key role in gradientbased 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 nonlinear 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 nondimensional equivalents, are identi
al in form but use nondimensional 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 reused, 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 nondimensionalization of Eqs. 4.2 and 4.3, where now the pressure and
oordinates are nondimensionalized: 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 − △F1j − △F4j − △F5j ) 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 − △F1j − △F4j − △F5j ) 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 − ∂ △F1j ∂nx,y − ∂ △F4j ∂nx,y − ∂ △F5j ∂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 ∂ △F1j ∂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 ∂ △F1j ∂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 − △F1j − △F4j − △F5j )} (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 − ∂ { △F1j + △F4j + △F5j } ∂UR − ∂ { △F1j + △F4j + △F5j } ∂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 leastsquares 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 OllivierGoo
h and
oworkers [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 leastsquares re
onstru
tion; the leastsquares 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 OllivierGoo
h and Van Altena [68℄ and OllivierGoo
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 leastsquares 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 leastsquares 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 leastsquares 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 leastsquares 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 leastsquares
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 higherorder limiter [56℄ designed not to degrade the a
ura
y of highorder s
hemes for smooth ows. We will derive the fa
e propertiesmesh 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 propertymesh 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
, nonlimited 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 OllivierGoo
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. OllivierGoo
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 nondierentiability 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 highorder limiter is related to the nondierentiability 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 nonlinearity 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 nonfeasible 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 nonlinearity 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 bestknown nongradient based optimization te
hnique is the geneti
algorithm [30℄, a biologi
allyinspired 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 NavierStokes 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 multistart 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 res
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 SQPRegPSO 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 SQPRegPSO 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 4digit 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 BFGSRegPSO 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 SQPhybrid ve
tor 86.6 Finite dieren
e gradient and SQPhybrid ve
tor 86.7 Table 7.1: Angles between adjoint, FD., and SQPhybrid 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 SQPhybrid ve
tor 86.93 Finite dieren
e and SQPhybrid ve
tor 86.87 Table 7.2: Angles between adjoint, FD, and BFGShybrid 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 nongradient 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
kfree solution is obtained, the obje
tive fun
tion is too at to distinguish reliably between the drag of sho
kfree solutions; this is not surprising,
onsidering that the drag for sho
kfree 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 RegPSOSQP 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 SQPRegPSO with RegPSOSQP 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 RegPSOSQP 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 SQPRegPSO 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 nonzero stru
ture of the se
ond order Ja
obian matrix 119 and dropped the rest of the nonzero 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 leastsquares 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 semitorsional 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/nongradient 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 nongradient random sear
h phase. The SQP optimization phase leads to a high quality initial bestsolution 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 gradientbased 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 threedimensional 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 953925, 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/VonKarmenInstitute 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 0486432270., 2003. [5℄ Azab, M. B., and OllivierGoo
h, C. F. Higher Order Two Dimensional Aerody nami
Optimization Using Unstru
tured Grids and Adjoint Sensitivity Compu tations. In AIAA20101430, 48th AIAA Aerospa
e s
ien
es meeting, Orlando, FL, USA, 47Jan 2010, Orlando FL, 2010. [6℄ Azab, M. B., and OllivierGoo
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, AIAA20110183, 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 Doublerank 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 19972310, 2002. [15℄ Chandrashekarappa, P., and Duvigneau, R. MetamodelAssisted Parti
le Swarm Optimization and Appli
ation to Aerodynami
Shape Optimization. Re sear
h Report RR6397, 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 850424, 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 RR3633, INRIA, 1999. Projet M3N. [20℄ Duvigneau, R. Aerodynami
Shape Optimization with Un
ertain Operating Conditions using Metamodels. Resear
h Report RR6143, 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. 3943., 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 TexasPan Ameri
an, 2009. [25℄ Farhat, C., Degand, C., Koobus, B., and Lesoinne, M. Torsional Springs for TwoDimensional Dynami
Unstru
tured Fluid Meshes. Computer methods in applied me
hani
s and engineering, 163(14):231245, 1998. [26℄ Frink, N. T. ThreeDimensional 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 ThreeDimensional Unstru
tured Meshes. In AIAA paper 910102, 1991. [28℄ Gage, P., and Kroo, I. Development of the QuasiPro
edural Method for Use in Air
raft Conguration Optimization. In AIAA Paper No. 924693, 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. 933933, 1993. [30℄ Goldberg, David E. Geneti
Algorithms in Sear
h, Optimization, and Ma
hine Learning. AddisonWesley 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 870520, 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 AIAA20051897, 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/S03760421(01)000045. [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 26 September 2002, page 253. Springer Netherlands, 2003. IUTAM Symposium Transsoni
um IV: pro
eedings of the IUTAM Symposium held in Guttingen, Germany 26 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 NavierStokes 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 20001584, 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. 19421948, 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/s0015800202124. [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 20040614, Reno, NV, 2004. [52℄ Masuda, H., Yoshioka, Y., and Furukawa, Y. Intera
tive Mesh Deformation Us ing EqualityConstrained 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 ANNBased Surrogate Models. Optimization and Engineering, 9:239255, 2008. 10.1007/s11081007 90311. [54℄ Mi
halak, C. E
ient HighOrder A
urate Unstru
tured FiniteVolume Algo rithms for Vis
ous and Invis
id Compressible Flows. PhD thesis, The University of British Columbia, 2009. [55℄ Mi
halak, C., and OllivierGoo
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 OllivierGoo
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 OllivierGoo
h, Carl. Globalized MatrixExpli
it NewtonGMRES for The HighOrder 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 ThreeDimensional Aerodynami
Shape Optimization. In Pro
eedings of the AIAA 18th Computational Fluid Dynami
s Conferen
e. AIAA 20073837, 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 CTAC2004, 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 OllivierGoo
h, C. A HighOrder A
urate Unstru
tured Finite Volume NewtonKrylov 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 EmbeddedBoundary Cartesian Mesh Method. Journal of Computational Physi
s, 227(4):27242742, 2008. [63℄ Neme
, M., and Zingg, D. W. . NewtonKrylov Algorithm for Aerodynami Design Using the NavierStokes Equations. AIAA Journal, 40:11461154, 2002. [64℄ Nielsen, E., and Anderson, K. Aerodynami
Design Optimization on Unstru
 tured Meshes Using the NavierStokes 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℄ OllivierGoo
h, C. and Van Altena, M. A HighOrderA
urate Unstru
tured Mesh FiniteVolume S
heme for The Adve
tionDiusion Equation. Journal of Computational Physi
s, 181(2):729752, 2002. [69℄ OllivierGoo
h, C., Nejat, A., and Mi
halak, C. On Obtaining and Verifying HighOrder FiniteVolume 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 971854, 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 993298, 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 19960094, 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/00219991(81)901285. [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/CP1999209136, June 1999, pp. 333343., 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. FreeForm 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
eTime 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
reteAdjoint 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
zanskiSobieski, 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. AddisonWesley Publishing Company, New York, 1992. [87℄ Zhu, C., Byrd, R. H., Lu, P., and No
edal, J. Algorithm 778: LBFGSB: Fortran Subroutines for LargeS
ale BoundConstrained 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 EddyVis
osity Assumption in GradientBased Design Optimization. Journal of Computational Physi
s, 229:52285245, 2010. 132