Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

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

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

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

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

Full Text

Aerodynami  Optimization Using  High-Order Finite-Volume CFD Simulations by Mohammad Baher Azab  B.S ., Aerospa e Engineering Cairo University (Egypt), 1999 M.S ., Aerospa e Engineering Cairo University (Egypt), 2002  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Fa ulty of Graduate Studies (Me hani al Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Van ouver) O tober 2011  ©  Mohammad Baher Azab 2011  Abstra t The growth of  omputer power and storage  apa ity allowed engineers to ta kle  engineering design as an optimization problem. mization is  For transport air raft, drag mini-  riti al to in rease range and redu e operating  osts. Lift and geometri  onstraints are added to the optimization problem to meet payload and rigidity  on-  straints of the air raft. Higher order methods in CFD simulations have proved to be a valuable tool and are expe ted to repla e  urrent se ond order CFD methods  in the near future; therefore, exploring the use of higher order CFD methods in aerodynami  optimization is of great resear h interest and is one goal of this thesis.  Gradient-based optimization te hniques are well known for fast  onvergen e, but  they are only lo al minimizers; therefore their results depend on the starting point in the design spa e.  The gradient-independent optimization te hniques  the global minimum of an obje tive fun tion but require vast therefore, for global optimization with reasonable  an nd  omputational eort;  omputational  ost, a hybrid op-  timization strategy is needed. A new least-squares based geometry parametrization is used to des ribe airfoil shapes and a semi-torsional spring analogy mesh morphing tool updates the grid everywhere when the airfoil geometry  hanges during shape optimization.  For the gradient based optimization s heme, both se ond and fourth order simulations have been used to well known for its low  ompute the obje tive fun tion; the adjoint approa h,  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  minimization without and with lift  inverse design problems and for drag  onstraint to validate the developed optimizer.  The optimization s heme used is Sequential Quadrati  Programming (SQP) with the  BFGS approximation of the Hessian matrix. A mesh renement study is presented for an aerodynami ally  onstrained drag minimization problem to show how se ond  ii  and fourth order optimal results behave with mesh renement. A hybrid parti le swarm / BFGS s heme has been developed for use as a global optimizer. It has been tested on a drag minimization problem with a lift  onstraint;  the hybrid s heme obtained a sho k free proles, while gradient-based optimization ould not in general.  iii  Prefa e The resear h ideas and methods explored in the three thesis are the fruits of a  o-authored manus ripts of this  lose working relationship between Dr Carl Ollivier-Goo h  and Mohammad Azab. The implementation of the methods, the data analysis, and the manus ript preparation were done by Mohammad Azab with invaluable guidan e from Carl Ollivier-Goo h throughout the pro ess.  iv  Table of Contents Abstra t Table of Contents List of Tables List of Figures List of Algorithms Nomen lature 1 INTRODUCTION  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  v  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  viii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  x  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xv  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xvi  . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  1.1  Finite Volume Flow Solver  1.2  Numeri al Aerodynami  1.3  . . . . . . . . . . . . . . . . . . . . . . .  Optimization  Gradient-based aerodynami  optimization  . . . . . . . . . .  6  1.2.2  Gradient free optimization  . . . . . . . . . . . . . . . . . . .  8  . . . . . . . . . . . . . . . . . . . . . . .  10  Contributions of the Thesis  Survey of Geometri  . . . . . . . . . . . . . . . . .  Parametrization Te hniques  12  . . . . . . . . . . .  12  . . . . . . . . . . . . . . . . . . .  13  2.1.1  Analyti al parametrization  2.1.2  Pie e-wise spline parametrization  2.1.3  CAD parametrization  2.1.4  Free form deformation (FFD)  2.1.5  Multidis iplinary aerodynami /stru tural shape optimization  . . . . . . . . . . . . . . .  14  . . . . . . . . . . . . . . . . . . . . . .  14  . . . . . . . . . . . . . . . . .  using deformations (MASSOUD) 2.2  5  1.2.1  2 GEOMETRY PARAMETRIZATION 2.1  . . . . . . . . . . . . . . . .  2  . . . . . . . . . . . . . . . .  New Least Squares Parametrization Te hnique  v  . . . . . . . . . . . .  15  15 16  2.2.1  Airfoil geometry parametrization  . . . . . . . . . . . . . . . .  16  2.2.2  Thi kness  . . . . . . . . . . . . . . . . . . . . . . .  20  2.2.3  Validation  . . . . . . . . . . . . . . . . . . . . . . . . .  23  onstraint ases  3 MESH MORPHING AND MESH SENSITIVITY  . . . . . . . . .  30  . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  30  3.1  Mesh Morphing  3.2  Testing Mesh Morphing  3.3  Mesh Sensitivity  . . . . . . . . . . . . . . . . . . . . . . . .  32  . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  33  4 GRADIENT CALCULATION USING ADJOINT APPROACH 4.1  Forward and Adjoint Formulations  4.2  Element and Fa e Geometri dinates  4.3  35  . . . . . . . . . . . . . . . . . .  35  Properties Dependen y on Mesh Coor-  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4.2.1  Element area mesh dependen y  4.2.2 4.2.3  44  . . . . . . . . . . . . . . . .  45  Fa e length mesh dependen y  . . . . . . . . . . . . . . . . .  46  Fa e normal mesh dependen y  . . . . . . . . . . . . . . . . .  49  Fa e Flow Properties Dependen y on The Mesh  . . . . . . . . . . .  50  . . . . . . . . . . . . . .  50  4.3.1  Fa e ow properties re onstru tion  4.3.2  Mesh dependen e of the fa e ow property re onstru tion  . .  57  . . . . . . . . . . . . . . . . . . . . . . .  60  . . . . . . . . . . . . . . . . . . . . . . . . . . .  60  5 GRADIENT VALIDATION 5.1  Subsoni  Test Case  5.2  Transoni  Test Case with No Limiter  5.3  Transoni  Test Case with Limiter  5.4  Sensitivity of nite dieren e gradient to perturbation magnitude  . . . . . . . . . . . . . . . . .  61  . . . . . . . . . . . . . . . . . . .  65  .  67  . . . . . .  70  . . . . . . . . . . . . . . . . . . . . . . . . .  70  . . . . . . . . . . . . . . . . . . . . . . . .  71  6 GRADIENT BASED OPTIMIZATION TEST CASES 6.1  Subsoni  Inverse Design  6.2  Transoni  6.3  Drag Minimization without Lift Constraint  6.4  Drag Minimization with Lift Constraint  6.5  Inverse Design  . . . . . . . . . . . . . .  74  . . . . . . . . . . . . . . . .  80  Mesh Renement Study of a Drag Minimization with Lift Constraint  86  vi  7 PARTICLE SWARM OPTIMIZATION AND A NEW HYBRID OPTIMIZATION METHOD . . . . . . . . . . . . . . . . . . . . . .  7.1  Introdu tion to Gradient Free Optimization  7.2  Swarm Intelligen e  . . . . . . . . . . . . . .  88  . . . . . . . . . . . . . . . . . . . . . . . . . . .  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  . . . . . . . . . . . . . . .  94  . . . . . . . . . . . . . . . . . . . .  95  . . . . . . . . . . . . . . . . . . . . . . . .  95  7.3  Hybrid SQP-RegPSO Te hnique  7.4  Optimization Test Cases  lipping  7.4.1  Constrained drag optimization of the NACA 0012  7.4.2  Constrained drag minimization of NACA 00083  7.4.3  Aerodynami NACA 0012  7.4.4  Aerodynami NACA 00083  7.5  88  and thi kness  . . . . . .  97  . . . . . . .  102  onstraint drag minimization of  . . . . . . . . . . . . . . . . . . . . . . . . . . . and thi kness  103  onstraint drag minimization of  . . . . . . . . . . . . . . . . . . . . . . . . . .  108  Comparing SQP-RegPSO with RegPSO-SQP optimization strategies  112  8 CONCLUSIONS AND RECOMMENDATIONS  . . . . . . . . . .  119  . . . . . . . . . . . . . . . . . . . . .  119  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  122  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  123  8.1  Contributions and Con lusions  8.2  Future Work  Bibliography  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 5.2  The magnitude of se ond and fourth order  CD  ow. . .  62  gradients and angles  between the evaluated gradients for NACA 0012 in an unlimited transoni 5.3  ow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  Magnitudes of se ond and fourth order  CD  64  gradients and angles be-  tween the evaluated gradients for NACA 0012 in Venkatakrishnan limited transoni 5.4  ow  . . . . . . . . . . . . . . . . . . . . . . . . . .  The magnitude of se ond and fourth order  CD  66  gradients and angles  between the evaluated gradients for NACA 0012 using higher order limiter in transoni 5.5  ow. . . . . . . . . . . . . . . . . . . . . . . . . .  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 5.6  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  Aerodynami  oe ients of original and optimized RAE 2822 airfoil onditions.  Aerodynami at transoni  6.3  .  ǫ  at transoni 6.2  ow.  67  Finite dieren e drag gradient sensitivity with respe t to perturbation amplitude  6.1  66  . . . . . . . . . . . . . . . . . . . . . . . . .  69  77  oe ients of original and optimized RAE 2822 airfoil onditions.  . . . . . . . . . . . . . . . . . . . . . . . . .  82  Lift penalty weight ee t on Drag minimization of RAE 2822 with lift onstraint  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  viii  86  7.1  Angles between adjoint, FD., and SQP-hybrid ve tors  . . . . . . . .  101  7.2  Angles between adjoint, FD, and BFGS-hybrid best solutions ve tors  103  7.3  Thi kness penalty weight impa t on the optimization results of NACA 00083 using hybrid s heme . . . . . . . . . . . . . . . . . . . . . . . .  110  7.4  Optimization results of the hybrid s heme in Phase I, II  . . . . . . .  113  7.5  T Optimization results of the hybrid s heme in Phase I, II  . . . . . .  113  ix  List of Figures 1.1  Airfoil aerodynami tion  2.1  optimization  y le using gradient-based optimiza-  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  7  Least square surfa e presentation of RAE 2822 airfoil using two polynomials  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 variables, presented with permission of Brezillon.  pq  a  b.  3.1  S hemati  3.2  Mesh movement s heme results of doubling the thi kness of NACA 0012  3.3  drawing of an edge  . . . . . . . . . . . . .  and its fa ing angles  and  . . .  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  S hemati right sides  4.2  31  33  Mesh movement results, large outer boundary deformation of a re tangular domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4.1  28  drawing of an element fa e  κ, and illustration of its left and  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  General triangular element with unit normal  x  34  n̂  on one of its fa es  κ.  40 44  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  5.1  . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  The pressure sensitivity with respe t to one of the design points  omputed for subsoni  sensitivity 5.2  k. .  CL  ow over NACA 0012,  ontrol  omparing the  al ulation of Eq. 4.6 with nite dieren e results.  . . . .  . . . . . . . . .  The pressure sensitivity with respe t to one of the design points  omputed for an unlimited transoni  omparing the sensitivity  CD  ow around NACA 0012,  al ulation of Eq. 4.6 with nite dieren e  The pressure sensitivity with respe t to one of the design points  5.6  63  gradient error in se ond and fourth order s hemes with respe t to  the design points, normalized by gradient magnitude. . . . . . . . . . 5.5  62  ontrol  results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4  61  gradient error in se ond and fourth order s hemes with respe t to  the design points normalized by gradient magnitude. 5.3  52  omputed for an unlimited transoni  The normalized  CD  64  ontrol  ow around NACA 0012.  65  gradient error in se ond, fourth, and modied  fourth order s hemes with respe t to the design points in a limited transoni 5.7  ow (Venkatakrishnan limiter) . . . . . . . . . . . . . . . .  The normalized  CD  68  gradient error in se ond, fourth, and modied  fourth order s hemes with respe t to the design points in a limited transoni  6.1  Subsoni  ow (higher order limiter).  . . . . . . . . . . . . . . . . . .  68  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  6.5  The dieren e between the target prole and the optimized proles,  inverse design optimal airfoil shapes  se ond and fourth order  . . . . . . . . . . . . .  . . . . . . . . . . . . . . . . . . . . . . . . .  xi  73  74  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  . . . . . . . . . . . . . . . . . . . . . . . . .  6.10 Transoni  inverse design optimal airfoil shapes.  6.11 Pressure  ontours of RAE 2822 at Ma h 0.73 and angle of atta k 2  6.12 Optimized RAE 2822.  . . . . . . . . . . . .  geometries.  77  .  78  . . . . . . . . . . . . . . . . . . . . . . . . . .  78  6.13 Optimization surfa e displa ements of the original RAE2822 surfa es. 6.14 Pressure distribution  76  79  omparison for the original and the optimized  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  79  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.  6.16 Se ond and fourth order optimization  . . . . . . . . . . . . . . .  80  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  82  6.19 Optimal shapes  omparison:  onvergen e history. . . . . . .  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 gradient 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) 6.23 Optimal  . . . . . . . . . . . . . . . .  85  value with mesh renement . . . . . . . . . . . . . . . .  87  Algorithm optimization . . . . . . . . . . . . . . . . . . . . .  90  CD  7.1  Geneti  7.2  Velo ity  7.3  Pressure Field of the original NACA 0012 and the optimized proles  omponents and position update of a parti le  xii  . . . . . . . .  93 98  7.4  Surfa e pressure distribution of original NACA 0012 and the optimized 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.  7.7  .  100  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 kness  onstraint  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  107  7.14 Surfa e pressure of NACA 0012 and the optimized prole with thi kness  onstraint  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  107  7.15 Pressure surfa e of NACA 00083 and the optimized proles with thi kness  onstraint  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  108  7.16 Surfa e pressure of NACA 00083 and the optimized prole with thi kness  onstraint  7.17 Prole  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  109  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. 7.20 Hybrid and hybrid  . . . . .  111  T optimal pressure distribution for NACA 0012  (Case I). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  114  T optimal pressure distribution for NACA 00083 7.21 Hybrid and hybrid (Case II).  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xiii  115  7.22 Hybrid and hybrid with 10% thi kness III).  T optimal pressure distribution for NACA 0012 onstraint and unit thi kness penalty weight (Case  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  116  T optimal pressure distribution for NACA 00083 7.23 Hybrid and hybrid with 10% thi kness IV).  onstraint and unit thi kness penalty weight (Case  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  117  T optimal pressure distribution for NACA 00083 7.24 Hybrid and hybrid with 10% thi kness  onstraint and ten thi kness penalty weight (Case  V). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xiv  118  List of Algorithms 7.1  Parti le swarm pseudo  7.2  Hybrid SQP-RegPSO pseudo  7.3  Pseudo  ode  . . . . . . . . . . . . . . . . . . . . . . . ode  . . . . . . . . . . . . . . . . . . .  T ode of RegPSO-SQP hybrid  optimization s heme  xv  . . . .  91 96 112  Nomen lature α  Angle of atta k  △Q  Change in  △t  Time step  δnorm  Normalized swarm radius  ∂R ∂Q  ja obian matrix  ∂U ∂M  Dependen y of premitive ow variables solution on mesh.  onservative ow properties  dCL,D,M Gradient of lift, drag. and moment dD  oe ient  dF dD  Obje tive fun tion gradient  dR dD  Residual dependen y of design variable  γ  Spe i  n̂  Fa e unit normal  Ωi  Control volume i  Uk  Control volume avaraged primitive ow properties  − → F  Flux ve tor a ross  − → 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  heats ratio  ontrol volume boundary  xvi  − → Vi  Parti le i velo ity  − → xi (k)  Parti le position in iteration k  − → xi (k + 1) ∂M/∂D ∂R/∂D → ∂− rb /∂Di ΨL,D,M  Parti le position in iteration k+1  Mesh sensitivity Residual senstivity with respe t to design variables Dependen y of boundary point position  − → rb ve  Adjoint ve tor of of lift, drag. and moment  ρ  Density  ρ̃  Roe avaraged density  à  Roe averaged ux ja obian matrix  ã  Roe avaraged speed of sound  h̃  Roe avaraged enthalpy  ũ  Roe avaraged x-velo ity  ṽ  Roe avaraged y-velo ity  F~ (QR ) , F~ (QL )  tor on design variable D  oe ient  Right and left fa e uxes  CD  Drag  CL  Lift  CM  Moment  Cp , Cv  Spe i  Di  Desigen variable i  et  Spe i  Kib  Stiness matrix intries relates interior nodes with boundary points  oe ient oe ient oe ient heats at  onstant pressure and volume  total energy  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  ws  Gaussian integration weight  yc  Geometry  a  Speed of sound  E  Energy per unit volume  e  Spe i  F  Aerodynami  ontrol volume boundary  ontrol point variable y ordinate (design veriable)  internal energy 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 OllivierGoo h. You are su h a great resear h supervisor and gentleman. I would like also to a knowledge my supervision  ommittee, Dr Mi hael Friedlander, Dr Philip Loewen,  and Dr James Olson for their valuable remarks that helped me to justify my resear h strategy. I would like also to thank my fellows in the resear h and development department of the Egyptian Air For e espe ially Col. Dr Eng. Mohamed Ibrahim Mustafa and Major General Dr Abdalla El-Ramsisy for their great support during my resear h in Egypt whi h has had a great impa t on my PhD resear h progress. Great thanks goes also to my dear professor Dr Chen Grief for being there when I needed his advi e. Finally, I would like to thank my wife and the light of my life, Deena, for her support and bringing happiness to my life.  Mohammad Baher Azab  xix  Dedi ation To my home  ountry, the land of Pharaohs  Egypt To the land of ho key  Canada To the souls of my parents  Baher and Laila  xx  Chapter 1  INTRODUCTION Drag redu tion of transport air raft is of great importan e be ause it redu es air raft fuel  onsumption and thereby redu es operating  osts and environmental impa t in  the form of pollution and global warming. Drag, lift and other aerodynami  for es  an be predi ted using CFD simulations, whi h have be ome an essential tool for aerodynami  analysis and design. CFD simulations  grids give a  urate aerodynami  for e predi tions, and unstru tured grids have the  advantage of easily representing any travel at transoni  arried out using unstru tured  omplex shape.  As most transport air raft  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 the optimization problem is to minimize an aerodynami hanging geometri to some geometri  obje tive (often drag) by  design variables, given an initial aerodynami and aerodynami  assessment of the aerodynami  onstraints.  design;  shape and subje t  This pro ess requires a  urate  hara teristi s of a given geometry. The ow model  used in this resear h is Euler's ow model; although this model negle ts vis ous ow ee ts, redu ing the drag using this invis id ow model will ultimately redu e the drag on real  ongurations as the invis id drag is about 40% of the total drag [47℄.  A higher order nite volume CFD solver developed by C. Ollivier-Goo h and  o-  workers [68, 61, 56℄ is used in this resear h. Optimization s hemes  an be divided into two main  ategories: gradient-based  and non-gradient-based s hemes. Gradient-based s hemes require  omputing the ob-  je tive fun tion value and its gradient with respe t to the design variables, while the non-gradient-based te hniques require only Gradient-based s hemes are fast to spa e  omputing the obje tive fun tion value.  onverge to a lo al optimal point in the design  ompared to the non-gradient-based methods, but the optimal solution found  by gradient-based optimization depends on the starting point [33℄.  1  Non-gradient-  based s hemes  an nd global minimum solution regardless of the starting point,  but with larger  omputational  ost  ompared to gradient-based s heme.  1.1 Finite Volume Flow Solver The two dimensional integral form of Euler's equation volume  Ωi  an be written for a  ontrol  as  ∂ ∂t  ¨  ˛  QdV +  Ωi  dΩi  − → F · n̂dl = 0  n̂ = nx i + ny j is the outward pointing normal → − Q and F are respe tively the onserved variable ve where  boundary of  ontrol volume  Ωi   boundaries. These  ρ    (1.1)  to the  ontrol volume fa es;  tor and the ux a ross the  an be expressed as    ρUn      ρu  →  ρUn u + nx p −   Q =   ρv  , F · n̂ =  ρU v + n p n y    E (E + p) Un "  u2 + v 2 p = (γ − 1) E − ρ 2  #         Un = n x u + n y v where is the uid density, pressure, and  E  normal to the spe i heats  γ  heats at  u  and  v  are the Cartesian velo ity  is the energy per unit volume. ontrol volume boundary. onstant pressure  (Cp ) and  Cp , Cν  Cp =  onstant volume  2  is the  Cν =  relations like the  (Cv ), the ratio of spe  onstant for air,  γR , γ−1  p  is the velo ity in the dire tion  Other thermodynami  an be expressed in terms of the gas  γ =  Un  omponents,  R:  R γ −1  i  For a thermally and  alori ally perfe t gas, thermodynami  properties  an be related  by  P = ρRT e = Cv T et = e + E = ρet  where  T  energy,  h  is the temperature, is the spe i  e   1 2 u + v2 2  P h=e+ ρ p a = γRT is the spe i  enthalpy, and  a  internal energy,  et  is the spe i  total  is the speed of sound.  For dis rete solution of the Euler equations, ow properties are normalized by some referen e values in order to redu e the round-o errors in the dis retized linear system; for external aerodynami s, the normalized ow properties are as follows  ρ̄ = ū = x̄ =  a P , P̄ = a∞ ρ∞ a2∞  v P̄ 1 v̄ = , Ē = + ρ̄ ū2 + v̄ 2 a∞ γ−1 2 y ta∞ ȳ = , t̄ = L L  ρ , ρ∞ u , a∞ x , L  ā =  With this normalization, the normalized Euler equations are identi al to their dimensional 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 evaluated at ea h  ell fa e  κ  using the following formula  i h ~ = 1 F~ (QR ) + F~ (QL ) − à (QR − QL ) F 2 κ where  QR , QL are the  onserved ow properties at right and left of the  the Roe averaged matrix  à  is the ux Ja obian  3  ∂ F~ /∂Q  (1.2)  ell fa e  κ and  at Roe-averaged quantities  as follows  ρ̃ = ũ =  h̃ =  √  ρL ρR  q uL + uR ρρRL q , 1 + ρρRL q hL + hR ρρRL q , 1 + ρρRL  q vL + vR ρρRL q ṽ = 1 + ρρRL  ũ2 + ṽ 2 ã2 = (γ − 1) h̃ − 2  The Roe-averaged Ja obian matrix tin t, but the eigensystem is  à  !  (1.3)  has four eigenvalues. Three of these are dis-  omplete. The Roe dissipation matrix  an be written  in terms the eigenvalues and eigenve tors of the Ja obian as       à = X̃     where  X̃  term  à (QR − QL )  Ũn  0  0  0  0  Ũn  0  0  0  0  Ũn + ã  0  0  0  0  Ũn − ã       −1  X̃     (1.4)  is the eigenve tor matrix evaluated using Roe-averaged ow quantities; the an be written a  à (QR − QL ) = =  ording to Frink as follows [27, 26℄  à △Q △F̃1 + △F̃4 + △F̃5  4  (1.5)  where  △F̃1  △F̃4,5  △P △ρ − 2  ã      Ũn ± ã  =        1        ũ ṽ ũ2 +ṽ2 2           + ρ̃       1   △P ± ρ̃ã△U  ũ ± nx ã  ṽ ± n ã 2ã2 y  ˜ ho ± Ũn ã        △u − nx △U   △v − ny △U    ˜ ũ△u + ṽ△v − Un △U   0        △ρ = ρR − ρL ,  △u = uR − uL ,  △P  △U = nx △u + ny △v  Higher-order a variables  Ũn  =         U =  h  = PR − PL ,  △v = vR − vL ,  ura y is obtained by least-squares re onstru tion of the non- onserved  ρ u v p  iT  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  where  △Q =  h  an be written as,    I ∂R + △t ∂Q  △ρ △ρu △ρv △E  n  iT  △ Qn = R n ,  (1.6)  ∂R ∂Q is the global Ja obian matrix, and  is the residual. The steady state solution is obtained iteratively when  △Q → 0.  R In  pra ti e, we use a quasi-Newton generalization of Eq. 1.6 that in ludes residual-based lo al time stepping [57℄ and solve the system using GMRES [75℄.  1.2 Numeri al Aerodynami Optimization Aerodynami  design used to rely on CFD simulations in  onjun tion with experimen-  tal testing and engineering intuition of the designer. With the growth of high speed omputers, integrating numeri al optimization s hemes with CFD simulations has  5  be ome possible and is now used for aerodynami  design and optimization. Gradient-  based optimization te hniques are widely used be ause they rea h an optimized shape after a reasonable exe ution time; however, the nal optimal shape is the lo al minimum lo ated downhill from the optimization starting point. methods like geneti optimum but  Non-gradient-based  algorithm (GA) or parti le swarm (PS) are slower to nd an  an nd the global minimum regardless of the starting point; their  drawba k is the large number of iterations required to rea h this minimum  om-  pared to gradient-based s hemes.  1.2.1 Gradient-based aerodynami optimization Gradient-based optimization depends on evaluating the gradient of the obje tive fun tion with respe t to the design variables and using the gradient in a linear model (steepest des ent) or a quadrati  model (Newton or quasi-Newton model)  to nd a sear h dire tion; this sear h dire tion is the dire tion in whi h the design variables should  hange their values to minimize the obje tive fun tion [66℄.  Gradient-based te hniques have been widely used in aerodynami their fast  onvergen e to an optimal solution.  The obtained optimal shape is bi-  ased by the optimization starting point (initial aerodynami no guarantee that gradient-based methods  optimization due  shape), and there is  an nd the global best optimal shape  in the design spa e [2℄. Hi ks and Henne where among the rst to apply gradientbased optimization te hniques in aerodynami  design in the late 1970's [35℄. Sin e  then many resear hers have investigated the use of gradient-based optimization te hniques like steepest des ent and quadrati tion [17, 31, 28℄. gradient  programming in aerodynami  optimiza-  The most expensive part of gradient-based optimization is the  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 whi h is  ompute the gradient (using a  entral nite dieren e rule),  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 dient redu ed the  omputational  ost of gradient  ompute the gra-  al ulation to the  ost of one ow  simulation regardless of the number of design variables [41, 43, 60, 38, 39, 42, 11, 62℄.  6  Figure 1.1: Airfoil aerodynami  optimization  y le using gradient-based optimization  Numerous resear hers have applied gradient-based optimization using the adjoint approa h in aerodynami  optimization sin e the early 1990's. Reuther and his  o-workers applied the adjoint approa h for aerodynami  optimization of air raft  onguration using Euler's ow model [73℄. Jameson developed an adjoint formulation for the Navier-Stokes equations and applied it to transoni  wing optimization  [43℄; in this work, Jameson suggested that an optimal pressure distribution rst be obtained using Euler's ow model then used as a target pressure for an inverse design optimization problem using a Navier-Stokes adjoint optimizer to redu e the overall omputational  ost. However, Jameson treated the eddy vis osity as  onstant whi h  was later shown to be a bad assumption. Anderson and Bonhaus examined the ee t of the strength of  oupling of the turbulen e model to the ow equations. They  om-  pared the adjoint gradient of the ow equations, with the eddy vis osity frozen, with the nite dieren e gradient of the found that freezing the eddy vis osity  ombined ow and turbulen e equations, and an lead to signi ant error in the  gradient. Therefore, they developed a ow solver that  ouples the Spalart-Allmaras  one-equation turbulen e model with the ow equations; their  5×5  blo ks for two dimensional ows and  They found that this  6×6  omputed  oupled solver used  blo ks for three dimensional  ases.  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 Zymaris et al developed a  k−ε  hemi ally rea ting ows [65℄.  ontinuous adjoint optimizer for turbulent ow using the  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 point, the sequential quadrati  onverge to minimal  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. Hessian is expensive to  The exa t  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. aerodynami  optimization for transoni  Dadone et.  and supersoni  wings and  al used BFGS in 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 approximation in gradient by Neme  ura y  aused by  omputations [19℄. BFGS optimization has also been used  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 heuristi  optimization methods, are optimization te hniques that do not require obje tive  fun tion gradient problems. They  omputation and therefore an be  an be applied to non-dierentiable  ategorized as evolutionary s hemes (in luding geneti  rithms) and random sear h s hemes (like the parti le swarm te hnique).  8  algo-  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 mutating these  hromosomes and optimization is  arried out by  rossing and  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 next optimization iteration. Transoni  hromosomes of the ospring of the  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 a geneti  algorithm in subsoni  he added the geometri  onstraint.  Somewhat later, Anderson applied  wing optimization with stru tural  and aerodynami  onstraints [1℄;  onstraints to the obje tive fun tion as  penalty terms. Jang and Lee applied a geneti  algorithm in subsoni  and transoni  invis id airfoil optimization; their obje tive was to maximize the lift-to-drag ratio  1 [44℄. Oyama et al applied a geneti  of an airfoil, beginning from the NACA 0012  algorithm with a Navier-Stokes solver for transoni  wing optimization [70℄.  also explored the use of fra tal analysis in GA aerodynami The parti le swarm method (PSOpt) is a sto hasti  They  optimization [71℄.  optimization sear h method  developed by Eberhart and Kennedy in 1995, inspired by the so ial behavior of bird o ks [46, 21℄. The general idea of the parti le swarm optimization is to randomly generate a swarm of parti les in the design spa e. For ea h parti le a tness value is  al ulated based on CFD simulation.  Then parti les y in the design spa e  a  ording to a simple formula that takes into a  ount that parti le's own best t-  ness position and the swarm's overall best tness position [22℄. PSOpt is known to suer from premature  onvergen e prior to dis overing the true global minimizer;  Evers suggests an automati  regrouping PSO (RegPSO) that automati ally triggers  swarm regrouping when premature  onvergen e is dete ted. The suggested regroup-  ing strategy aims to liberate parti les from sub-optimal solutions and enables nding the global minimum [24℄. Although the PSO algorithm has been applied to a wide 1  The National Advisory Committee for Aeronauti s (NACA) airfoil family geometri  oordinates  an be found at University of Illinois Urbana-Champaign website http://www.ae.illinois.edu/mselig/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 optimization te hnique in the multidis iplinary optimization of a wing; the obje tive was to maximize the air raft range by maximizing lift-to-drag ratio and redu ing the wing weight subje t to geometri  onstraints [85℄. Chandrashekarappa and Du-  vigneau [15℄ used parti le swarm optimization s heme to aerodynami ally optimize wings in supersoni to aerodynami  onditions; Duvigneau also applied parti le swarm optimization  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 fun tion at lower  an  ompute an a  omputational  used. The rst major  urate value of an aerodynami  obje tive  ost than required when a se ond order method is  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 gradientbased 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 rea h a true global optimum (for invis id transoni subje t to aerodynami  and geometri  an  aerodynami s, a sho k free airfoil  onstraints) after a reasonable number of CFD  simulations. This s heme is a hybrid (BFGS + regrouped parti le swarm) s heme that takes the advantages of gradient-based and gradient-free optimization s hemes. To a hieve these goals, the following pre-existing high-order a  omponents were needed, in addition to the  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. Chapter 2 des ribes the requirements on su h a parametrization and presents a new least-squares spline parametrization method developed for this work.  Robust mesh movement.  As the airfoil shape  hanges during optimization, the  omputational mesh must be updated. Be ause mesh regeneration will intro-  10  du e una  eptably large  hanges in the dis retization error, mesh movement is  strongly preferred. Chapter 3 des ribes the semi-torsional spring mesh movement s heme used in this thesis.  Mesh sensitivity al ulation. a  Cal ulating the gradient of the obje tive fun tion  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  bines aspe ts of the geometri  om-  parametrization and mesh movement s hemes.  Obje tive fun tion gradient al ulation.  E ient gradient  al ulation requires  the solution of the adjoint to the governing ow equations. Chapter 4 des ribes how this is done.  A key innovation of the thesis is e ient solution of the  dis rete adjoint problem for a high-order a  urate ow solution s heme.  addition to showing the formulation of the gradient  al ulation, this  ompares nite dieren e and adjoint gradients for subsoni  are presented both with and without limiting of the  hapter  and transoni  for se ond and fourth order nite volume s hemes. For transoni  In  ows  ow, results  omputational solution to  prevent overshoots.  Optimization drivers.  The gradient-based optimizer using in this thesis is the  BFGS-based quasi-Newton solver in Matlab's optimization toolbox. Gradient based optimization test  ases for inverse design problems and for drag mini-  mization with and without a lift  onstraint is shown in Chapter 6. Also, the  impa t of mesh renement on the se ond and fourth order optimal results is studied. The regrouped parti le swarm optimization  ode was written by the author.  This s heme and its hybridization with the BFGS s heme is presented in Chapter 7. Examples are given to show the ee tiveness of the hybrid s heme in nding global optima when the gradient-based s heme is unable to.  11  Chapter 2  GEOMETRY PARAMETRIZATION The geometry of engineered obje ts is dened mathemati ally in  omputer-aided  design (CAD) software, then exported as a group of points or polygons whi h approximates the original geometry. This dis rete data provides input to mesh generation software that  mesh ),  reates a dis rete representation of the  omputational domain (a  whi h in turn is used as input for CFD analysis. Aerodynami  optimization modies the aerodynami variables. An obvious approa h  shape by  design and  hanging a set of geometri  design  hoi e is to use all the surfa e grid points of a wing, but this  auses two problems. First, it makes the design spa e very large and this  may lead to a highly expensive optimization. Se ond, it may lead to a non-smooth geometry due to the displa ement independen e of a surfa e mesh point; this problem  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 geometri  parametrization. Geometri  parametrization te hniques  an be  lassied as  analyti al; pie e-wise spline tting; CAD based; and free form deformation (FFD) approa hes. In Se tion 2.1, a review of various geometry parametrization te hniques is presented. Se tion 2.2 des ribes in detail the geometry parametrization te hnique used in this resear h, a novel pie e-wise least-squares tting te hnique [5, 6℄.  2.1 Survey of Geometri Parametrization Te hniques In this review, various te hniques applied to aerodynami  optimization are presented;  the basi s of these te hniques and their appli ation limitations are dis ussed, in luding 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 optimization [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 optimization design variables. The sinusoidal displa ement shape fun tion is expressed as  h(x) = where  a  and  b  are      sin πx  onstants to  ln(0.5) ln(b)  0<x<1  (2.1)  ontrol the peak lo ation and the width of the  sinusoidal displa ement shape fun tion.  b  a  a=4  is re ommended for most  ases, while  must be between 0 and 1 [45℄.  The Class-Shape fun tion-Transformation method (CST) is another analyti al geometry parametrization te hnique presented by Kulfan [49, 48℄. CST parametrizes the airfoil geometry using the following formula  y = x0.5 (1 − x)S (x) + x△ZT E where  S (x)  is the shape fun tion and  △ZT E  (2.2)  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 X i=0  where the weights proles, it is not  bi  bi  n! xi (1 − x)n−i i! (n − i)!  (2.3)  are used as the design variables. Although CST gives non-wavy  apable of representing  omplex geometries. Mousavi and Nadiragah  ompared the impa t of using dierent geometry parametrization te hniques on the optimal wing shape; using CST parametrization gave drag  oe ients higher by  about 15% than the optimal geometry from a B-spline parametrization [58℄ for a three-dimensional lift- onstrained transoni  13  drag minimization problem.  2.1.2 Pie e-wise spline parametrization Bezier  urves  an be used for airfoil shape parametrization. Obayashi used Bezier  urve for aerodynami the Bezier  optimization using geneti  algorithm [67℄; he noti ed that  urve representation fails to represent geometries that gives rooftop pres-  2 be ause Bezier  sure distributions ization of Bezier  urves are always  onvex.  urves, were found to be more suitable.  B-splines, a general-  A  resentation is a very good geometry parametrization te hnique. the  ubi  spline representation  ubi  B-spline rep-  Better  ontrol of  an be obtained by in reasing the number of splines  that represent the airfoil; Li et al. optimized the NACA 0012 at single point and multi-point operating  onditions with a lift  onstraint using spline representation as  geometry parametrization te hnique [50℄. CAD systems typi ally use non-uniform rational B-spline (NURBS) representation for geometry modeling, allowing them to represent any in  omplex geometry; a detailed dis ussion of NURBS  The NURBS Book  [72℄.  Mengistu and Ghaly applied su  parametrization s heme to turboma hine blade aerodynami  an be found  essfully a NURBS  optimization using a  gradient-free method [53℄. While pie e-wise spline parametrization is well-suited for two-dimensional shapes and simple three-dimensional geometries, require a large number of  omplex shapes  ontrol points, redu ing the ee tiveness of gradient-based  optimization te hniques [76℄.  2.1.3 CAD parametrization Computer aided design pa kages have evolved to implement NURBS for geometry representation due the ex ellent properties of NURBS. Linking the CAD and grid generation software  an be done using an API that allows a  internal interfa e [81, 8℄. However, imposing geometri le. Mesh sensitivity  ess to the CAD system's  onstraints is still an obsta-  al ulation is another obsta le for gradient-based optimization  te hniques based of CAD parametrization: analyti al mesh sensitivity with respe t to geometry design variables  the NURBS knots  with the use of automati without CAD sour e 2  an in prin iple be  omputed  dierentiation of the CAD software but this is not possible  ode and is unlikely to be pra ti al even then, given the size of  Distributions for whi h pressure remains almost un hanged over a signi ant  tan e.  14  hordwise dis-  the CAD  ode base. This derivative  but the risk of poor a  an also be  ura y still exists, and  omputed using nite dieren es,  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  obtained without a hange in surfa e topology. The surfa e itself  an be  an be parametrized  using Bezier or B-splines or even NURBS splines. A related approa h, the free form deformation algorithm (FFD) treats the geometry as a void in a box-shaped pie e of rubber. Deformation  an be  ontrolled by moving  ontrol points pla ed on the outer  surfa e of the box; the interior of the rubber box with its void is parametrized using a tensor produ t of three spline representations (one in ea h  oordinate dire tion).  Sederberg and Parry developed an algorithm that uses the FFD  on ept with Bezier  tri-variate volume representation [78℄. A disadvantage of the FFD method is that it requires large numbers of  ontrol points to obtain lo al deformation in the deformed  geometry. However, Borrel and Rappoport presented a method to allow lo al shape deformation via FFD by introdu ing a set of B-splines that  ontrol points with  onstrained lo al  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 MASSOUD method is an analogy to analyti al methods that tries to parametrize the deformations in the geometry rather than the geometry itself. lizes SOA algorithms and allows strong lo al deformation  It also uti-  ontrol. The MASSOUD  parametrization requires few design variables be ause it parametrizes the deformation. Samareh has applied the MASSOUD method to parametrize a simple wing, a wing body blend, and a Anderson su  omplex air raft  onguration with su  ess [76℄. Nielsen and  essfully adopted the MASSOUD parametrization s heme for aerody-  15  nami  optimization of turbulent ow using unstru tured grids [64℄.  2.2 New Least Squares Parametrization Te hnique A new least-squares based surfa e parametrization method is presented in this se tion; the airfoil surfa e is parametrized using pie e-wise polynomials whose  oe-  ients are found by solving a least squares problem [6℄. This se tion also des ribes how to implement a thi kness some validation test  onstraint with the new parametrization and presents  ases.  2.2.1 Airfoil geometry parametrization In the proposed te hnique, the geometry is parametrized using pie e-wise polynomials found by a least-squares t. The parametrization polynomials are by a set of  2 ontrol points and satisfy C  ontrolled  ontinuity at their meeting points.  The  airfoil upper and lower surfa es are represented using two least square splines for  3 used are  ea h surfa e as shown in Fig. 2.1. The polynomials  √ P1 (x) = a0 x + a1 x + a2 x2 + a3 x3 2  P2 (x) = b0 + b1 x + b2 x + b3 x where  x  is the normalized  hord-wise position,  arates the polynomial regions, and  L = 1.  3  L1  0 < x < L1 L1 < x < L  is  hord-wise position that sep-  These polynomials are suitable for an  x = 0.  The  x and y  √  term in  P1 (x),  oordinates of the design  ontrol  airfoil with a rounded leading edge due to the existen e of the whi h gives an innite slope at  (2.4)  x  points shown in Fig. 2.1 are used to nd the values of the polynomial These polynomials must satisfy meeting point  ˆ  Value  x = L1 .  These  ontinuity of value, slope, and  onditions  urvature at their  an be written as  ontinuity:  P1 (L1 ) − P2 (L1 ) = 0 3  oe ients.  Te hni ally,  P1  is not a polynomial be ause of the presen e of the  slope and nite radius of  urvature at  x = 0.  (2.5)  √  However, the label is  onfusing.  16  x  term used to give innite  onvenient and not overly  ˆ  Slope  ontinuity: ′  ′  ′′  ′′  P1 (L1 ) − P2 (L1 ) = 0 ˆ  Curvature  (2.6)  ontinuity:  P1 (L1 ) − P2 (L1 ) = 0 An additional  P2 (L)  onstraint should be added on  trailing edge. The above hard  (2.7)  to assure zero thi kness at the  onstraints should be stri tly satised by the geom-  etry parametrization polynomials; they  an be written in matrix form as  BP = 0 where   √ L1  1  √  2 L B =  −1 1  √ 3  4 L1 0  L1 1  L31  L21  −1 −L1  2L1 3L21  0  0  2  6L1  0  0  0  0  1  −1 0  L  −L21  −L31   −2L1 −3L21   , −2 −6L1   L2 L3  The free parameters are  hosen to best approximate the  shape  x  ontrol points; the  least square system with  oordinates of the  onstraints applied  min kAP − ck2 where the  y  a tual  A  ontains powers of the  x    y  P =  3  b0      b1      b2     b  3  ontrol points are xed. The resulting an be expressed as  subje t to  BP = 0  (2.8)  oordinates at the design points so that  oordinates of the design points. I have used a set of  ontrols the airfoil shape fun tions instead of using the fun tions to ease                    a0     a1      a2     a   oordinates of the airfoil  oordinates of the parametrized shape at the design points and  y                   c  AP  gives  ontain the  ontrol points that  oe ients of the shape  onstraints and boundary denitions. To give an example of how  this least-squares system is  onstru ted,  17  onsider the parametrized airfoil surfa e  Figure 2.1: Least square surfa e presentation of RAE 2822 airfoil using two polynomials  P1 (x)  &  P2 (x)  shown in Fig. 2.1.  tted using nine  Six  ontrol points lies in the region of the polynomial  while three points lies in the is   √                   P2 (x)  region. The  x1 x1 x21 x31 0  0  0  0  x2 x22 x32 0  0  0  0  √ x2 √ x3 √ x4 √ x5 √ x6  Be ause the  ontrol points.  x33 x34 x35 x36  0  0  0  0  0  0  0  0  0  0  0  0  x6  x23 x24 x25 x26  0  0  0  0  0  0  0  0  1 x7 x27 x37  0  0  0  0  1 x8 x28 x38  0  0  0  0  1 x9 x29 x39  x3 x4 x5  P1 (x),  orresponding least-squares system                                          a0           a1           a2          a3 =  b0            b1           b2          b3    y1     y2      y3      y4    y5    y6      y7      y8     y9  (2.9)  onstraint equation, Eq. 2.8, has a zero right hand side, the solution  ve tor  P  linear  ombination of the null spa e basis of the  must lie in the null spa e of the  onstraint equations, i.e. it should be a onstraint equations. The matrix  is full row rank and to nd the null spa e basis of it, QR fa torization  B T = QR 18  B  an be used:  Q=  h  − → → → → → → → → q1 − q2 − q3 − q4 | − q5 − q6 − q7 − q8  The ve tors of the matrix  B.  matrix  Q2  [Q2 ]  i  =  h  Q1 | Q2  i  matrix are unit ve tors forming a basis of the null spa e of  The solution ve tor  P  must be a linear  ombination of the ve tors of the  → → → → P = z1 · − q5 + z2 · − q6 + z3 · − q7 + z4 · − q 8 = Q2 z  (2.10)  Substituting Eq. 2.10 into Eq. 2.8  AP = AQ2 z = c  (2.11)  Solving the least-squares system by singular value de omposition for numeri al stability,  z = [AQ2 ]† c where the re tangular matrix  [AQ2 ]†  (2.12)  is the pseudo-inverse of  [AQ2 ].  Finally,  P = Q2 [AQ2 ]† c  (2.13)  Equation 2.13 gives the relationship between the polynomial  y  lo ations of the  the  i  th the i "  − → rb =  oe ients with respe t to the  ontrol point, whi h is needed to olumn of the matrix  xa ya  #  Q2 [AQ2 ]  on a design variable  → ∂− rb = ∂Di  "  → ∂− rb = ∂Di  "  ∂a0 √ ∂Di xa  Di  †  ∂a1 ∂Di xa  +  +  ∂b1 ∂Di xa  +  and the  lo ation of  ∂M/∂D ,  is  .The dependen y of an airfoil surfa e point  an be found from  ∂a2 2 ∂Di xa  +  ∂a3 3 ∂Di xa  0 ∂b0 ∂Di  y  al ulate the mesh sensitivity  0 +  P  c.  ontrol points  The sensitivity of the polynomial th  oe ients  ∂b2 2 ∂Di xa  +  ∂b3 3 ∂Di xa  19  #  #  0 < xa < L 1  L 1 < xa < L  (2.14)  where  h  ∂a0 ∂Di  ∂a1 ∂Di  ∂a2 ∂Di  ∂a3 ∂Di  ∂b0 ∂Di  strained pseudo inverse matrix  ∂b1 ∂Di  ∂b2 ∂Di †  Q2 [AQ2 ]  .  ∂b3 ∂Di  iT  is the  This pro edure  ith  ve tor of the  on-  an be extended to  parametrize airfoil surfa e using any number of pie e-wise polynomials, with an a ompanying in rease in system size. Leading edge radius and trailing edge thi kness onstraints system  an be added to the  BP = d,  and the  C2  ontinuity requirements forming the  ontrol point  oordinates are used to  onstraint  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  an sele t a set of number of less  oe ients, the designer  ontrol points whi h lies on the parametrized surfa e; the least  ontrol points is two per polynomial. The less  ontrol of airfoil geometry. I have sele ted nine  ontrol points used, the  ontrol points at spe i  hord-  wise x-stations based on my engineering sense; it turns out that my sele ted set of ontrol points were able to produ e sho k free optimal proles as will be shown in the next  hapters, however, sele tion of the  ontrol points x-stations  by formulating a simple minimization problem. geometri  an be made  In this problem, a set of airfoil  data gathered and the obje tive is to minimize the RMS error between  the original airfoil surfa es and the parametrized airfoil surfa es, where the design variables are the  ontrol points x-stations.  2.2.2 Thi kness onstraint Some air raft fuel tanks are pla ed inside the wing, and the main landing gear of some air raft are stored in the wings after being retra ted. In addition, the wing must have su ient bending rigidity from a stru tural point of view; therefore, the wing thi kness is a  onstraint at some  how to add a thi kness  hord-wise stations. This subse tion demonstrates  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 + a1 x + a2 x2 + a3 x3 P2 (x) = b0 + b1 x + b2 x2 + b3 x3 20  0 < x < L1 L1 < x < L  Figure 2.2: Parametrized airfoil using four polynomials  and two for the lower surfa e  √ P3 (x) = c0 x + c1 x + c2 x2 + c3 x3 P4 (x) = d0 + d1 x + d2 x2 + d3 x3  If the thi kness this  tc  need to be  C t : a0 x c + a1 x c +  a2 x2c  +  L1 < x < L  onstrained at a station  onstraint will be  √  0 < x < L1  a3 x3c    − c0  √  xc  where  0 < xc ≤ L 1  Ct : P1 (xc ) − P3 (xc ) = tc  xc + c1 xc + c2 x2c + c3 x3c = tc  ,  (2.15)  Equation 2.15 provides a link between the upper and lower surfa e polynomials; therefore, a  oupled least-squares system needs to be  21  onstru ted and solved; the  hard  onstraint equations are expressed as      B  0  0  B  P1 (xc ) {0}1×4 −P3 (xc ) {0}1×4                       a0      a1    . . .     d2     d3   =      0           0     . .  (2.16)  .         0        t  c  [Ba ] {Pa } = {da }  (2.17)  The global least-squares system is  "  A  0  0  A  #  (  {Pa } =  cu cl  )  (2.18)  [Aa ] {Pa } = {ca } The above hard the thi kness  onstraint equations do not have a null right hand side due to  onstraint; therefore, the solution (polynomials  belong to the null spa e of  Ba .  oe ients) does not  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  ˆ  Let  ˆ ˆ  Q=  Q  and  of  R.  h  Q2  The hard  Let  i  Q1 Q2  ,  R=  BaT = QR. # " R1  ontains the rest of the  onstraints  QT Pa =  "  y z  #  .  ∴ RT QT Pa =  where  0  olumns, and  an now be written as  Then  "  R1T 0  Pa = Q #  ·  h  Q1  "  y z  #  Q1 Q2  R1 is the  rst  7×7  olumns of sub matrix  BPa ≡ RT QT Pa = da  = Q1 y + Q2 z iT  Pa = da 7−→ R1T y + 0 = da  ∴ y = R1−T da 22  ontains the rst 7  (2.19)  ˆ  Be ause  R  is upper triangular, Equation 2.19  an be solved using forward  substitution.  ˆ  The soft  onstraints in equation 2.18  Aa Pa = Aa  an be rewritten as  h  Q1 y Q2 z  ∴ Aa Q2 z = ca − Aa Q1 y  i  = ca  (2.20) (2.21)  ∴ z = [Aa Q2 ]† [ca − Aa Q1 y]  ˆ  Finally, the polynomial  oe ients  an be found as  (2.22)  Pa = Q1 y + Q2 z :  i h Pa = Q1 − Q2 [Aa Q2 ]† Aa Q1 R1−T da + Q2 [Aa Q2 ]† ca  (2.23)  Using the last equation we noti e that the sensitivity of the polynomial ients with respe t to the the matrix  [Q2 ] [AQ2 ]  †  ith  ontrol point  y  lo ation,  yci ,  is the  ith  oe-  ve tor of  .  2.2.3 Validation ases In this subse tion, testing of the proposed least-squares surfa e parametrization is arried out; parametrization of various types of airfoils is done to show that the proposed geometry parametrization s heme has the required exibility to represent various types of airfoils.  These airfoils in lude NACA 4-, 5-, and 6-digit series,  4 Figures 2.32.10 show  laminar ow, supersoni , and super- riti al airfoil se tions.  the least-squares tting polynomials for various airfoils, while Table 2.1 shows the RMS error in the parametrized geometry for upper and lower surfa es. The last parametrized airfoil, the laminar LV2  5 airfoil of German Aerospa e  Center, was parametrized using 10 polynomials, 5 for ea h of its surfa es. Although the RMS error is small (order of  10−5  hord for all parametrization method), as shown  in Table 2.1, Fig. 2.11 shows that the dieren e of the surfa e pressure distribution, espe ially at the peak velo ity at the leading edge upper surfa e, is still signi ant. There are many x-stations at whi h LV2 airfoil 4  hanges its  urvature, and also a  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 parametrization s heme to a  urvature value whi h makes it hard for the geometry urately present it and this  auses the u tuation in  pressure resulted in the parametrized airfoil be ause of la k of a of  urvature u tuations. Better mat hing  an be obtained but this will in rease the  number of geometry design variables signi antly. This trade o between a  ase illustrates  urately representing the geometry (whi h will  optimization iteration) and  learly the  hange during  hoosing a reasonable number of design variables; this  hoi e must be left for the designer. If the una  urate presentation  hange in pressure distribution be omes  eptable, or if the pressure distribution requires large number of airfoil  points (design variables) to be a  ontrol  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 variables, presented with permission of Brezillon.  28  Airfoil NACA 0011SC NACA 0012 NACA 6409 NACA 16-006 NACA 63-412 NACA 64-421 RAE 2822 LV2  Upper surfa e error  Lower surfa e error  · 10−4  · 10−4  6.37 2.35 · 10−4 9.28 · 10−5 1.31 · 10−6 1.17 · 10−5 1.70 · 10−4 1.20 · 10−5 1.45 · 10−5  6.37 2.35 · 10−4 1.24 · 10−4 2.03 · 10−5 5.63 · 10−5 7.50 · 10−5 3.20 · 10−5 3.22 · 10−5  Maximum Thi kness  11 · 10−2 12 · 10−2 9 · 10−2 6 · 10−2 12 · 10−2 21 · 10−2 −2 slightly > 12 · 10 −2 slightly > 12 · 10  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  of the mesh that presents the shape. This the new geometry but this is time error [51℄.  hange  an be done by grid regeneration around  onsuming and will  hange the dis retization  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  30  al ulations. In this  ase,  Figure 3.1: S hemati  drawing of an edge  pq  and its fa ing angles  the elements of the boundary layer experien ed small geometri al elements away from the airfoil experien ed larger elasti ity mesh movement s heme is spring analogy method.  a  and  b.  hange while the  hanges [80℄. However, the linear  omputationally expensive  ompared to the  The semi-torsional spring analogy method is adopted in  this resear h due to its simpli ity and robustness. Consider the edge  pq  shown in Fig 3.1. The relationship between the for es and  the node displa ements when treating this edge as a spring follows Hooke's Law as    Fpx    Fpy   F  qx Fqy      −1      0  1 1  = 1 1+ +  lpq sin (θa ) sin (θb )   1  0  where lpq is the length of edge  pq, and θa  and  θb        −1 0 1  0 −1 0      1 0 −1  0  1  0   upx     u  py  uqx     u   (3.1)  qy  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  "  Kii  Kib  Kbi Kbb  an be written as  #(  Ui Ub  31  )  =  (  0 Fb  )  (3.2)  where  Ui and Ub are the interior and boundary mesh point displa  We do not need to know the values of boundary nodal for es points displa ement ve tor  Ub  ements respe tively.  Fb be  ause the boundary  is known expli itly: it is the deformation required in  the airfoil prole to minimize the obje tive fun tion. Therefore, Eq 3.2  "  as  where  0  I  #(  )  Ui Ub  −−→ − → → Uik = rik+1 − rik , Ub = − rb − − r→ bo in  Substituting  "  Kii Kib  Kii Kib 0  I  0 Ub  )  (3.3)  Eq 3.3 we get  #k ( #k ( −−→ ) " k+1 K K ri ii ib = − → 0 I rb  − r→ bo is the initial position ve  The stiness matri es  =  (  an be written  ) − ) ( → 0 rik + − Ub r→ bo  (3.4)  tor of the boundary points before mesh morphing.  [Kii ] and [Kib ] depend on the mesh fa  e lengths, whi h  hange  during the mesh morphing stage; therefore, Eq. 3.4 is a non-linear equation and needs to be solved iteratively.  3.2 Testing Mesh Morphing Farhat et al [25℄ demonstrated the robustness of the semi-torsional mesh movement s heme. In this se tion, testing results are presented to demonstrate that the  ur-  rent implementation of this s heme shows the same good behavior for several highdeformation NACA 0012.  ases. The rst test  ase is an unstru tured triangular mesh around a  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  ure 3.2 shows that even with this large displa ement (multiple no  y -dire  tion. Fig-  ells) of the boundary,  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). outer boundaries are  The  hanged, redu ing the total area by almost 50% and turning  the original right-angled  orners nearly into  32  usps. Figure 3.3b shows that the semi-  Figure 3.2: Mesh movement s heme results of doubling the thi kness of NACA 0012  torsional mesh movement was  apable of adapting the mesh in the entire eld without  element ipping.  3.3 Mesh Sensitivity For gradient-based optimization, as shown in the next tation requires  ∂R/∂D ,  whi h in turn depends on the evaluation of mesh sensitivity.  of the geometri  ompu-  al ulation of the dependen y of the residual on the design variables  sensitivity tells how mesh points translate in the parametrization  (x, y)  The mesh  plane with the perturbation  ontrol points (design variables). This translation,  obviously, depends on the mesh movement s heme. gebrai  hapter, the gradient  Truong et.  al.  ompared al-  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  33  ∂M/∂Di  (that is,  (a) Initial mesh  (b) Final mesh  Figure 3.3: Mesh movement results, large outer boundary deformation of a re tangular domain  the  hange in mesh point lo ations with a  hange of a design variable) is  al ulated  by dierentiating Eq. 3.3:  "  ∴ where  → ∂− rb /∂Di    ∂M ∂Di  Kii Kib    0  =  I (  ∂Ui ∂Di ∂Ub ∂Di  #(  ∂Ui ∂Di ∂Ub ∂Di  )  "  =  )  =  (  Kii Kib 0  I  0  → ∂− rb ∂Di  )  #−1 (  0  → ∂− rb ∂Di  )  (3.5)  is obtained using Eq. 2.14 and is related to the design variables via  the pseudo inverse of the  onstrained least-squares system solved in parametrizing  the geometry.  34  Chapter 4  GRADIENT CALCULATION USING ADJOINT APPROACH Gradient  al ulation plays a key role in gradient-based optimization. The traditional  nite dieren e strategy is  omputationally expensive as it requires at least as many  CFD simulations as the number of design variables to  ompute the gradient; ea h of  these is the solution to a large non-linear system of equations. The to  ompute the gradient is less expensive as it requires the  sensitivity with respe t to the design variables, and uses the to  forward strategy  al ulation of the ow  omputed ow sensitivity  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, variables, volumes  F,  for aerodynami  optimization is a fun tion of the design  D , and the ow eld solution at the surfa  e points of the boundary  ontrol  Us F = F (Us , D)  Us  is expressed most  onveniently in primitive variables:  (4.1)  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:    ˛    ˛  Ps ny ds cos α Ps nx ds sin α + FL = − ˛  ˛  FD = Ps nx ds cos α + Ps ny ds sin α or in dis rete form,  where  Ps  normal  X X FL = − ws Ps nx sin α + ws Ps ny cos α X X FD = ws Ps nx cos α + ws Ps ny sin α is the pressure at surfa e integration point pressure,  omponents at the surfa e integration point,  is the ar  α is the  nx  (4.3)  and  ny  are the unit  angle of atta k, and  ws  length asso iated with the surfa e integration point. Note that this form  uses the dimensional pressure and for es. The lift and drag  oordinates gives the dimensional lift and drag  oe ients, whi h are the non-dimensional equivalents, are  identi al in form but use non-dimensional pressure and  oordinates.  These dis rete integrals expressed as a fun tion of geometri of the  (4.2)  and ow properties  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 variables. The gradient of the obje tive fun tion rule  where  an be obtained by using the  ∂F ∂Ubg ∂U ∂M ∂F ∂M dF = + dD ∂Ubg ∂U ∂M ∂D ∂M ∂D Ubg  is the boundary Gauss point ow properties,  averaged solution written in primitive variables, and  ∂M/∂D  design  M  U  hain  (4.4) is the  ontrol volume  is the mesh point lo ations.  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 be written as a fun tion of the ow eld solution variables  D.  If we apply the  U  and mesh geometri  onstraint that the ow solution is  36  an  design  onverged regardless  of variations in the design variables, we  an write  ∂R ∂M ∂R ∂U ∂M dR = + =0 dD ∂M ∂D ∂U ∂M ∂D Solving this for the solution sensitivity  ∂U ∂D  ≡  (4.5)  ∂U ∂M ∂M ∂D , we get        ∂U ∂M ∂R −1 ∂R ∂M ∂R ∂Q −1 ∂R ∂M =− =− · · ∂M ∂D ∂U ∂M ∂D ∂Q ∂U ∂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 ow variables (whi h is used in impli it ow solvers) and a  onserved  hange of variables for  ∂Q onserved to non- onserved. Note that ∂U is a blo k diagonal matrix. Substituting Eq.4.6 in Eq.  4.4, we get the forward formulation for gradient  omputation      ∂R ∂M ∂F ∂Ubg ∂R ∂Q −1 ∂F ∂M dF =− · + dD ∂Ubg ∂U ∂Q ∂U ∂M ∂D ∂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 T =− dD    ∂R ∂M ∂M ∂D   T      ∂R ∂Q −T ∂F ∂M T ∂F ∂Ubg T · + ∂Q ∂U ∂Ubg ∂U ∂M ∂D  where the residual sensitivity to mesh movement  an be written using the  (4.8)  hain rule  as:  ∂R ∂M  =  X  ∂R dnx ∂R dny + + ∂nx dM ∂ny dM fa es Ωi  ∂R ∂Ufa e ∂R dwi + ∂wi dM ∂Ufa e ∂M  ∂R ∂AΩi + ∂AΩi ∂M  We get the adjoint method, presented by A. Jameson [37, 41, 43, 38, 40℄.  (4.9)  Now  only one linear system solve is required. However, this linear system solve requires  37  expli itly forming the global Ja obian matrix  ∂R ∂Q be ause its transpose is required.  Be ause the CFD solver used in this thesis  an form the global Ja obian matrix  expli itly  (author?)  [57℄, the transpose of the Ja obian  an easily be formed as  well. The Ja obian from the last GMRES iteration is re-used, so the  omputational  eort of solving the adjoint problem is redu ed to solving this linear system, a on the order of 1% of the CFD simulation To ease the programming eort when  omputational eort. hanging the obje tive fun tion, we form  three adjoint problems, one ea h for the lift and the moment  oe ient  tively. These aerodynami  CM  to nd  oe ient  CL ,  the drag  ∂CL /∂D , ∂CD /∂D  oe ient gradients    dCL dCD dCM , , dD dD dD  The solution pro edure of the three adjoint problems  ˆ  Using the steady state ow solution, we matrix  ˆ  We  ∂R/∂Q  onstru t  and  oe ient  ∂CM /∂D ,  CD ,  respe -  an be used to evaluate any obje tive  fun tion gradient for a fun tion that depends on aerodynami  dF =f dD  ost  for es:   an be summarized as follows,  onstru t the CFD simulation Ja obian  expressed in Eq. 1.6.  ∂R/∂U  as:  ∂R ∂Q ∂R = ∂U ∂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  the  ontrol volume average values of the primitive ow properties. The latter  is the dependen y of surfa e point primitive ow properties on  is known as a side ee t of solution re onstru tion.  ˆ  We solve the three linear systems adjoint ve tors  ˆ  We   ∂R  ΨL,D,M .   ∂R T ∂U  ΨL,D,M =  h  onstru t the obje tive fun tion gradient ve tor  ∂M ∂M ∂D  T  ΨL,D,M ,  i ∂CL,D,M ∂Ubg T to get the ∂Ubg ∂U dCL,D,M T dD  =   ∂F  ∂M ∂M ∂D  T  −  whi h requires only a ve tor dot produ t for ea h design  38  variable.  As an example of how to use fun tion,  dCL,D,M to dD  onstru t the gradient of an aerodynami  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 The gradient of the above fun tion  also  (4.10)  an be written as  dCD dCL dF = + 2k1 (CL − CLc ) dD dD dD  (4.11)  ∂CD ∂CL ∂F = + 2k1 (CL − CLc ) ∂D ∂D ∂D  (4.12)  The dis rete (invis id) forms of  CL and CD follow from the non-dimensionalization  of Eqs. 4.2 and 4.3, where now the pressure and  oordinates are non-dimensionalized:    X X ws Ps ny cos α CL = − ws Ps nx sin α +   X X ws Ps ny sin α CD = ws Ps nx cos α +  (4.13)  and their partial derivatives with respe t to the geometry design variables are  ∂CL ∂D  ∂CD ∂D  X   X  ∂ws ∂ws = − Ps nx sin α + Ps ny cos α ∂D ∂D X  X  ∂nx ∂ny ws Ps − ws Ps sin α + cos α ∂D ∂D X  X  ∂ws ∂ws = Ps nx cos α + Ps ny sin α + ∂D ∂D X  X  ∂ny ∂nx ws Ps ws Ps cos α + sin α ∂D ∂D  Using Eq. 4.9 to  ∂nx /∂M , ∂ny /∂M ,  al ulate and  tions, while the terms  ∂wi /∂M ,  ∂R/∂nx  the Roe ux, and nally  ∂R/∂M , and  we need to  ompute the terms  (4.14)  ∂AΩi /∂M ,  whi h only depend on mesh points' spatial lo a-  ∂R/∂ny  are obtained by dire t dierentiation of  (∂R/∂Uf ace ) · (∂Uf ace /∂M ) 39  is obtained by dierentiating  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 the derivatives on the geometri  Let us  onsider the residual  R,  with  terms in the following se tion.  ontribution of fa e  κ of the  ontrol volume  Ωi  shown  in Fig. 4.1.  The dis rete form of this edge's residual  Rκ,Ωi  ontribution to the  ontrol volume  J    1 X 1 ~ FRj + F~Lj − |△F1 |j − |△F4 |j − |△F5 |j wj = AΩi 2  Ωi  is  (4.15)  j=1  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. with the Gauss point  J = 2,  on the fa e  κ. wj  j. 40  is the integration weight asso iated  To nd  ∂Rκ, Ωi /∂AΩi ,  Eq. 4.15  an be dire tly dierentiated to get     ∂Rκ, Ωi 1 1 X 1 ~ ~ FRj + FLj − |△F1 |j − |△F4 |j − |△F5 |j wj = − =− Rκ,Ωi 2 ∂AΩi 2 AΩi (AΩi ) j  (4.16)  The term  ∂Rκ, Ωi /∂nx,y  ∂Rκ,Ωi 1 X = ∂nx,y AΩi j  where  (  1 2    ∂ F~R,Lj =  ∂nx   ∂ |△F1 |j ∂ |△F4 |j ∂ |△F5 |j ∂ F~Lj ∂ F~Rj + − − − ∂nx,y ∂nx,y ∂nx,y ∂nx,y ∂nx,y  !  wj  )  (4.17)  omponent terms    an be written as  an be expanded using the denition of the Roe ux to give:   (ρu)R,L   ρu2 + p R,L  ,  (ρuv)R,L  ({E + p} u)R,L      ∂ F~R,Lj =  ∂ny   (ρv)R,L (ρuv)R,L  ρv 2 + p R,L  ({E + p} v)R,L         and  ∂ |△F1 |j ∂nx  =    1      0       △u − nx △U ũ △P   + ρ̃  △ρ − 2      ã △v − ny △U    ṽ     ũ2 +ṽ2 ũ△u + ṽ△v − U˜n △U 2      0           −△U − n △u  x  ρ̃     −ny △v         ˜ −ũ△U − U △u   uŨ q n Ũn2  Ũn         n  41         +       ∂ |△F1 |j ∂ny  =           1           ũ v Ũn △P   + ρ̃  q △ρ − 2      ã  ṽ   Ũn2     ũ2 +ṽ2 2      0             −n △u x   Ũn ρ̃       −△U − ny △v       ˜ −ṽ△U − U △v         △u − nx △U  +  △v − ny △U    ˜ ũ△u + ṽ△v − Un △U  0  n  ∂ △F̃4,5 ∂nx  =   1   Ũn ± ã ũ  △P ± ρ̃ã△U    ũ ± nx ã  r  ṽ ± n ã 2 2 2ã y  Ũn ± ã h˜o ± Ũn ã   1     ũ ± n ã ±ρ̃ã△u  x  + Ũn ± ã   2ã2  ṽ ± ny ã  h˜o ± Ũn ã   0    ±ã  △P ± ρ̃ã△U    Ũn ± ã  0  2ã2   ±ũã 42      +    ∂ △F̃4,5 ∂ny  To nd   1   Ũn ± ã ṽ  △P ± ρ̃ã△U    ũ ± nx ã  r  ṽ ± n ã 2 2 2ã y  Ũn ± ã ˜ ho ± Ũn ã   1     ±ρ̃ã△v  ũ ± nx ã   Ũn ± ã  ṽ ± n ã  + 2ã2 y   ˜ ho ± Ũn ã   0    0  △P ± ρ̃ã△U    Ũn ± ã  ±ã  2ã2   ±ṽã  =  ∂Rκ, Ωi /∂wj ,  ∂Rκ, Ωi 1 = ∂wj (AΩi )    +    Eq. 4.15 is dierentiated; we get     1 ~ ~ FRj + FLj − |△F1 |j − |△F4 |j − |△F5 |j 2  The dependen y of the fa e residual properties at the fa e    Ufa  e  ontribution  Rκ, Ωi  (4.18)  on the re onstru ted ow  an found as the sum of the dependen y on the right  and left fa e ow properties as  ∂Rκ, Ωi ∂Rκ, Ωi ∂Rκ, Ωi = (4.19) + ∂Ufa e ∂UR ∂UL o n   ∂ |△F | + |△F | + |△F | ~ ~ X 1 4 5 ∂ FLj j j j 1  ∂ FRj 1 + − =  (AΩi ) 2 ∂UR ∂UL ∂UR o  n  ∂ |△F1 |j + |△F4 |j + |△F5 |j  wj −  ∂UL  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 age to dierentiate the C, C++, or FORTRAN fun tion that  43  dierentiation pa k-  omputes the fa e ux.  Figure 4.2: General triangular element with unit normal  Finite dieren es  n̂  an also be easily implemented to evaluate  on one of its fa es  ∂Rκ Ωi /∂Ufa  κ.  e  Rκ, Ωi (UR + ǫ) − Rκ, Ωi (UR − ǫ) Rκ, Ωi (UL + ǫ) − Rκ, Ωi (UL − ǫ) ∂Rκ, Ωi = + +o (ǫ)2 ∂Ufa e 2ǫ 2ǫ (4.20)  where  ǫ  is  −8 in the above hosen to be 10  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 of the verti es. for  ell  properties depend dire tly on the spatial  Fig. 4.2 shows a s hemati  oordinates  drawing of a triangular element used  entered nite volume simulations, with its three verti es and one of its  three fa es labeled for later referen e. evaluating geometri  In the next subse tions, the pro edure for  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  A=σ  − − → r→ 12 , r13  (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 variables through the pseudo inverse matrix in the least-squares t during geometry parametrization.  It worth mentioning that any general polygon  several triangles and the area of ea h triangle The element area mesh dependen y  an be  an be  an be split into  al ulated using Eq. 4.21.  al ulated by dire t dierentiation of  Eq. 4.21  dA dM  =  ∂A ∂x1 ∂A ∂x2 ∂A ∂x3 + + ∂x1 ∂M ∂x2 ∂M ∂x3 ∂M ∂A ∂y1 ∂A ∂y2 ∂A ∂y3 + + + ∂y1 ∂M ∂y2 ∂M ∂y3 ∂M  where  ∂A ∂x1 ∂A ∂x2 ∂A ∂x3  ∂A ∂y1 ∂A ∂y2 ∂A ∂y3  = = =  = = =  σ (y2 − y3 ) 2 σ (y3 − y1 ) 2 σ (y1 − y2 ) 2  σ (x3 − x2 ) 2 σ (x1 − x3 ) 2 σ (x2 − x1 ) 2 45  (4.22)  σ = 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,  al ulated as  L=  q  (x3 (M ) − x2 (M ))2 + (y3 (M ) − y2 (M ))2  The fa e length mesh dependen y  dL dM  =  an be  (4.23)  an be obtained as  ∂L ∂x2 ∂L ∂x3 ∂L ∂y2 ∂ L ∂y3 + + + ∂x2 ∂M ∂x3 ∂M ∂y2 ∂M ∂y3 ∂M  where  ∂L ∂x2  =  ∂L ∂x3  =  ∂L ∂y2  =  ∂L ∂y3  =  q q q q  −x3 + x2  =  − (x3 − x2 ) L  x3 − x2  =  (x3 − x2 ) L  −y3 + y2  =  − (y3 − y2 ) L  y3 − y2  =  (y3 − y2 ) L  (x3 − x2 )2 + (y3 − y2 )2 (x3 − x2 )2 + (y3 − y2 )2 (x3 − x2 )2 + (y3 − y2 )2 (x3 − x2 )2 + (y3 − y2 )2  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 higher order ux integration. A boundary element with one  urved to enable  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,  46  ∂L/∂x2,3  and  ∂L/∂y2,3  are evalu-  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 ∂y2,3  = =  L(x2,3 + ǫ) − L(x2,3 − ǫ) 2ǫ 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  dL dwi = dM 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  47  x2 , y 2 .  Figure 4.4: S hemati  drawing of a  urved fa e with two Gauss points used  Point 1  wi ki  ( 12  L/2 − √112 )L  Point 2  ( 12  L/2 + √112 )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 Fig. 4.2)  that points outward from element  Ωi  (shown in  an be found using the three verti es of the element as  n̂ =  "  nx ny  #  =  (~c3 − βr̂23 ) k~c3 − βr̂23 k  (4.25)  where  ~c3 =  r̂23  "  2/3 x3 − (1/3 x1 + 1/3 x2 )  2/3 y3 − (1/3 y1 + 1/3 y2 )  x3 −x2 √ 2 2 +(|−y3 +y2 |)  =  (|−x3 +x2y|)3 −y 2 √ 2 2   β =  (|−x3 +x2 |) +(|−y3 +y2 |)  (x3 − x2 ) (2/3 x3 − (1/3 x1 + 1/3 x2 )) q (|−x3 + x2 |)2 + (|−y3 + y2 |)2  +  (y3 − y2 ) (2/3 y3 − (1/3 y1 + 1/3 y2 )) q (|−x3 + x2 |)2 + (|−y3 + y2 |)2  The mesh dependen y of the unit normal  dn̂ dM  #  =  an be written as  ∂ n̂ ∂x1 ∂ n̂ ∂x2 ∂ n̂ ∂x3 + + ∂x1 ∂M ∂x2 ∂M ∂x3 ∂M ∂ n̂ ∂y1 ∂ n̂ ∂y2 ∂ n̂ ∂y3 + + + ∂y1 ∂M ∂y2 ∂M ∂y3 ∂M  where  ∂ n̂ ∂xi  =  α1 − α2 k~c3 − βr̂23 k2  ∂ (~c3 − βr̂23 ) ∂xi ∂ k~c3 − βr̂23 k = (~c3 − βr̂23 ) · ∂xi  α1 = k~c3 − βr̂23 k · α2  49  (4.26)  ∂ n̂ ∂yi  =  α3 − α4 k~c3 − βr̂23 k2  ∂ (~c3 − βr̂23 ) ∂yi ∂ k~c3 − βr̂23 k = (~c3 − βr̂23 ) · ∂yi i = 1, 2, 3  α3 = k~c3 − βr̂23 k · α4  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  ∂Ufa 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 multiphysi s nite volume solver der a  apable of  ondu ting CFD simulations up to fourth or-  ura y; this solver has been built by Ollivier-Goo h and  In ANSLib, the ow properties are assumed to a  o-workers [61, 56, 54℄.  hange within ea h  ording to a linear polynomial for se ond order simulations, or a  ontrol volume  ubi  polynomial  for fourth order simulations. These polynomials are found using a least-squares reonstru tion; the least-squares system depends on the  ontrol volume properties like  moments of area whi h eventually depend on mesh point lo ations. In the next two subse tions, the re onstru tion and mesh dependen y of ow properties at fa es will be presented.  4.3.1 Fa e ow properties re onstru tion Solution re onstru tion is the key to determine solver a assumed to vary linearly within the  ording to  The solution is  ontrol volume for se ond order a  higher order methods, ow properties are assumed to vary a polynomial; variation a  ura y.  ubi  ura y; for  ording to higher order  polynomial is fourth order a  urate. The  ow solution re onstru tion method presented in this subse tion follows the method des ribed in Ollivier-Goo h and Van Altena [68℄ and Ollivier-Goo h et al [69℄.  50  h  U=  For se ond order solution re onstru tion, the primitive variables are re onstru ted at the fa e Gauss points as  U2R (x, y) = Uref + (x − xref )  ∂U ∂x  ref  + (y − yref )  ∂U ∂y  ρ u v p  (4.27)  ref  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:  U4R (x, y) = U2R (x, y) +  (x − xref )2 ∂ 2 U 2 ∂x2  ∂2U + (x − xref ) . (y − yref ) ∂x∂y +  (x − xref ) 6  3  + ref 2  +  ref  ~xref  is  ref  ∂3U  (x − xref ) (y − yref ) 2 ∂x2 ∂y ∂3U  (x − xref ) (y − yref ) 2 ∂x∂y 2  The referen e point  (y − yref )2 ∂ 2 U + 2 ∂y 2 2  ∂3U ∂x3  (4.28)  ref  + ref  (y − yref ) 6  hosen to be the triangle  enter for  3  ref  ∂3U ∂y 3  ref  ell- entered nite  volume s heme and the vertex lo ation for vertex- entered s heme.  A nomial  onstrained least-squares system is oe ient. The hard  preserve the  omputed  onstru ted to nd the re onstru tion poly-  onstraint is that the re onstru tion polynomial should  ontrol volume average solution. The rest of the equations   whi h are satised only in a least-squares sense  are obtained by requiring that the average of the re onstru tion polynomial approximates the solution for neighboring  ontrol volumes. Neighbors are  ontrol volume average  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 of an element  The mean  k  drawing  and its neighbor layers.  onstraint equation is obtained by requiring  51  onservation of the  ontrol  iT  Figure 4.5: S hemati element  drawing of rst, se ond and third neighbor layers of triangular  k.  volume average by the re onstru tion polynomial and  Uk =  1 Ak  ˆ  Ak  UkR dA ≡ Uk ref + xk  x2 k ∂ 2 U 2 ∂x2 where  xn y m  k  + xy k k ref  1 = Ak  ¨  Ak  ∂2U ∂x∂y  ∂U ∂x  + yk k ref  + k ref  an be written as  ∂U ∂y  y2k ∂ 2 U 2 ∂y 2  + ... k ref  (x − xk ref )n (y − yk ref )m dA 52  + k ref  (4.29)  Mat hing the  Uj =  j  ontrol volume average in a neighbor  1 Aj  ˆ  UkR dA  1 Aj  ˆ  (y − yk ref )  1 Aj  ˆ  (x − xk ref ) (y − yk ref )  1 Aj  ˆ  (y − yk ref )2 ∂ 2 U 2 ∂y 2  Aj  Aj  Aj  Aj  ≡ Uk ref ∂U ∂y  1 + Aj  ˆ  Aj  dA + ref  would require that  ∂U ∂x  (x − xk ref ) 1 Aj  ˆ  ∂2U ∂x∂y  (x −  dA +  ref xk ref )2  Aj  ∂2U ∂x2  2  (4.30)  dA + ref  dA + ref  dA + ... ref  Using  (x − xk ref ) = (x − xj ref ) − (xk ref − xj ref ) and  (y − yk ref ) = (y − yj ref ) − (yk ref − yj ref ) , we get  Uj  =  1 Aj  ˆ  Aj  UkR dA ≡ Uk ref + x bk,j  c2 k ∂ 2 U x 2 ∂x2  k,j  =  1 Aj  k ref  + ybk,j  ∂U ∂y  yb2 k,j ∂ 2 U 2 ∂y 2  k ref  +x cy k,j  ˆ  ((x − xj ref ) − (xk ref − xj ref ))n ·  where  \ xn y m  ∂2U ∂x∂y  ∂U ∂x  Aj  + k ref  +  (4.31)  k ref  + ... k ref  ((y − yj ref ) − (yk ref − yj ref ))m dA m X n X n! m! (xj ref − xk ref )q · = r! (m − r) ! q! (n − q) ! r=0 q=0  m−r (yj ref − yk ref )r xjn−q ref yj ref  for simpli ity, the term  x, yk ref  will be written as  53  x, yk ,  and  x, yj ref  will be repla ed  with  x, yj .    1  The resulting  xk  This  x2k c x2 k,1 c x2 k,2 c x2 k,3  yk    1 x bk,1    1 x bk,2   1 x bk,3   . .  .. . .  1 x bk,N  onstrained least-squares system  ybk,1  ybk,2  ybk,3 . . .  xyk x cyk,1  x cyk,2  x cyk,3  . . .  c x2 k,N  ybk,N  yk2 yb2 k,1 yb2 k,2 yb2 k,3  . . .  ···  ···  ···  ···  . . .  ..  yb2 k,N  x cy k,N  onstrained least squares system    .  ······                  an be written as  U ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2  satisfying [BRec ] (PRec ) = U k  where  x bk,1 − xk ,  ybk,1 − y k ,              =            k  Uk     U1   U2    U3   .  . .  UN (4.32)  an be rewritten as     min [ARec ] PeRec − U Rec − U k    . . .    c x2 k,1 − x2k , c x2 k,2 − x2 ,    x cyk,1 − xyk ,  yb2 k,1 − yk2 , yb2 − y 2 ,  ···    x x cyk,2 − xyk , ··· k,2 k k  bk,2 − xk , ybk,2 − y k ,  x 2 2 c b 2 2 cyk,3 − xyk , y k,3 − yk , ··· [ARec ] =  bk,3 − xk , ybk,3 − y k , x k,3 − xk , x  . . . . . ..  . . . . . . . . . . .  2 2 c b 2 2 x bk,N − xk , ybk,N − y k , x k,N − xk , x cyk,N − xyk , y k,N − yk , · · · · · · h i [BRec ] = 1 xk y k x2k xyk yk2 · · · and    PRec        =        U ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2 . . .           ,       k    PeRec       =       ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2 . . .          ,      k  54    U1 − Uk    U2 − Uk    U Rec − U k =  U 3 − U k  .  . .  UN − Uk                      The least-squares system is solved to obtain the re onstru tion polynomial  oef-   ients using singular value de omposition method to nd the pseudo inverse of the least-squares  oe ient matrix as follows  ARec = U ΣV T A†Rec = V Σ† U T               ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2 . . .      x2 k,1 − x2k , x cy k,1 − xyk , yb2 k,1 − yk2 , ··· x bk,1 − xk , ybk,1 − y k , c    2 2 c b 2 2  x cy k,2 − xyk , y k,2 − yk , ···   bk,2 − xk , ybk,2 − y k , x k,2 − xk , x   2 2 c b 2 2  = x cy k,3 − xyk , y k,3 − yk , ··· bk,3 − xk , ybk,3 − y k , x k,3 − xk , x    . . . . . ..  . . . . .  . . . . . .    2 2 c2 k,N − x , x x bk,N − xk , ybk,N − y k , x cyk,N − xyk , yb2 k,N − yk , · · · · · · k k    U1 − Uk    U2 − Uk   U −U  3 k  .  . .  UN − Uk  Uk = U k −  h            (4.33)    xk y k x2k xyk yk2  The unlimited fa e ow properties     i   ··· ·        ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2 . . .               an be found at fa e Gauss points  55  (4.34)  xf i , y f f    †       ·     as follows  Ufi = Uk + (xfi − xk )  ∂U ∂x  (xfi − xk )2 ∂ 2 U 2 ∂x2  k  (yfi − yk )2 ∂ 2 U 2 ∂y 2  k  k  reate a new extremum at the  at the fa e  (4.35)  ∂2U ∂x∂y  onvergen e, a limiter is applied in order  an be written as  (xfi − xk )2 ∂ 2 U 2 ∂x2  k  (yfi − yk )2 ∂ 2 U 2 ∂y 2  k  is the limiter value in the  fun tions to  al ulate  φ:  ∂U ∂x  k  + (yfi − yk )  ∂U ∂y  +  (4.36)  k ∂2U  + (xfi − xk ) (yfi − yk ) ∂x∂y !  + k  + ...  ontrol volume  k.  ANSLib uses two limiter  Venkatakrishnan's limiter [84℄ and a higher-order limiter [56℄  designed not to degrade the a  ura y of high-order s hemes for smooth ows. We  will derive the fa e properties-mesh dependen y for the limited ase  + k  ontrol volume fa e; the limited ow properties  Ufi = Uk + φk (xfi − xk )  φk  + k  + ...    where  ∂U ∂y  + (xfi − xk ) (yfi − yk )  To ensure monotoni ity and solution not to  + (yfi − yk )  an be obtained by simpli ation of the limited  56  ase; the unlimited  ase by setting  φk = 1.  4.3.2 Mesh dependen e of the fa e ow property re onstru tion  The term  ∂Ufa e /∂M  for the limited  ase  an be obtained using dire t dierentiation  of Eq. 4.36  ∂Ufi ∂M  =  ∂Uk ∂φk + ∂M ∂M    (xfi − xk )  (xfi − xk )2 ∂ 2 U 2 ∂x2  ∂U ∂x  k  + (yfi − yk )  ∂U ∂y  +  ∂2U + (xfi − xk ) (yfi − yk ) ∂x∂y k !  (yfi − yk )2 ∂ 2 U + ... + 2 ∂y 2 k  ∂ (yfi − yk ) ∂U ∂ (xfi − xk ) ∂U + φk ∂M ∂x k ∂M ∂y  (4.37)  k  + k  + k  2  ∂  (xfi −xk ) 2  ∂M  ∂2U ∂x2  k  ∂2U ∂y 2  k  +  2  ∂  (yfi −yk ) 2  ∂M   φk (xfi − xk )  ∂ (xfi − xk ) (yfi − yk ) ∂ 2 U ∂M ∂x∂y   + ... +  ∂U ∂x k  ∂  ∂M  ∂  + (yfi − yk )  ∂U ∂y k  ∂M  ∂2U  (xfi − xk )2 ∂ ∂x2 2 ∂M  k  ∂2U  (yfi − yk )2 ∂ ∂y2 2 ∂M  k  ∂φk /∂M = 0  ∂2U ∂x∂y k  ∂M  +  + ...  The unlimited fa e ow property-mesh dependen y and  + ∂  + (xfi − xk ) (yfi − yk )   + k  an be obtained by setting  φk = 1  in Eq. 4.37.  In Eq. 4.37, the term  ∂φk /∂M  an be obtained by dierentiating the limiter  57  expression, and the terms  and  ∂  ∂nU ∂xm ∂y n−m k  /∂M  an be evaluated as follows:  o n  = U Rec [ARec ] PeRec n o  eRec n o ∂ P ∂ U Rec ∂ [ARec ] e PRec + [ARec ] = = {0} ∂M ∂M n∂M o  o n ∂ PeRec † ∂ [ARec ] PeRec = − [ARec ] ∴ ∂M ∂M n o ∂ PeRec  ∂ [ARec ] = − [ARec ]† [ARec ]† U Rec (4.38) ∂M ∂M i n o xk y k x2k xyk yk2 · · · · PeRec i n o ∂ h ∂Uk =− ∴ xk yk x2k xyk yk2 · · · · PeRec − ∂M ∂M n o i ∂ PeRec h xk y k x2k xyk yk2 · · · · ∂M ∵ Uk =U k −  where       ∂ [ARec ]   =  ∂M     h  ∂ ybk,1 ∂M ∂ ybk,2 ∂M ∂ ybk,3 ∂M  ∂b xk,1 ∂M ∂b xk,2 ∂M ∂b xk,3 ∂M . . .  . . .  ∂ ybk,N ∂M  ∂b xk,N ∂M  n  ∂ PeRec ∂M  o  c2 ∂x k,1 ∂M c2 ∂x k,2 ∂M c2 ∂x k,3  ∂M . . .  . . .  c2 ∂x k,N ∂M        ∂   = ∂M      58  c 2 ∂y k,1 ∂M c 2 ∂y k,2 ∂M c 2 ∂ y k,3 ∂M  ∂x c y k,1 ∂M ∂x c y k,2 ∂M ∂x c y k,3 ∂M  . . .  c 2 ∂y k,N ∂M  ∂x c y k,N ∂M  ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2 . . .               ··· ··· ··· ..  .  ······  (4.39)              and    ∂ ∂xk m n m−1 n ∂xfi {(xfi − xk ) (yfi − yk ) } = m (xfi − xk ) (yfi − yk ) − + ∂M ∂M ∂M   ∂yk m n−1 ∂yfi − n (xfi − xk ) (yfi − yk ) ∂M ∂M  nym \ ∂x k,j ∂M  =  m X n  X  n! m! q (xj ref − xk ref )q−1 · (yj ref − yk ref ) · r! (m − r) ! q! (n − q) ! r=0 q=0   ∂xk ref ∂xj ref n−q m−r − + xj ref yj ref · ∂M ∂M m! n! r (xj ref − xk ref )q (yj ref − yk ref )r−1 · r! (m − r) ! q! (n − q) !   ∂yj ref ∂yk ref n−q m−r xj ref yj ref · − ∂M ∂M  59  Chapter 5  GRADIENT VALIDATION In this se tion, the  omparison of the gradient a  fourth order s hemes is  arried out.  Both se ond and fourth order gradients are  evaluated using the adjoint approa h and dieren e gradient. We present three test limited transoni  test  ases. In all test  ura y evaluated using se ond and  ompared to their  orresponding nite  ases: subsoni , non-limited transoni  and  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  airfoil geometri  design  oe ient gradient with respe t to the  ontrol points is presented for both se ond and fourth order  s hemes, in luding a  omparison with a nite dieren e gradient  the same order of a  ura y in the ow solver.  is NACA 0012 at subsoni representative  tations.  The airfoil used in this test  onditions M = 0.5 and  α = 2o .  ase  Figure 5.1 shows a  omparison between the eld of pressure sensitivity with respe t to  one of the geometry design points sensitivity  al ulated using  omputed using nite dieren es and the solution  al ulation of Eq. 4.6 for both se ond and fourth order a  Agreement is ex ellent in all  urate  ompu-  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 omputed for subsoni  ow over NACA 0012,  ontrol points  omparing the sensitivity  al ulation  of Eq. 4.6 with nite dieren e results.  agreement between the obje tive fun tion gradient are the y lo ations of the design  ontrol points,  omponents  ∂CL /∂yi ,  where  yi  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 subsoni  omputed gradient ve tor for  ow.  5.2 Transoni Test Case with No Limiter In this test = 0.8 and  ase, the sensitivity analysis is  α=  2◦ . The drag  arried out for NACA 0012 airfoil with M  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  61  omputed using Eq. 4.6  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  tween the evaluated gradients for NACA 0012 in subsoni  Figure 5.2:  CL  gradients and angles beow.  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 omputed for an unlimited transoni tivity  ow around NACA 0012,  ontrol points  omparing the sensi-  al ulation of Eq. 4.6 with nite dieren e results.  and nite dieren e; the gure also shows that for transoni sensitivity  omputed using se ond order and fourth order a  ows, the pressure 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 optimization 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  Figure 5.4:  CD  ow.  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 omputed for an unlimited transoni  ontrol points  ow around NACA 0012.  5.3 Transoni Test Case with Limiter In this test of the  ase, the impa t of using a limiter in the CFD simulation on the a  ura y  omputed gradient for se ond order and fourth order s hemes is studied. Two  dierent limiters are used, the Venkatakrishnan limiter [84℄ and the higher order limiter of Mi halak and Ollivier-Goo h [56℄. Figure 5.5 shows very good mat hing of the se ond order pressure sensitivity  omputed using adjoint and nite dieren e  te hniques with the use of Venkatakrishnan limiter. Mat hing is less good between the pressure sensitivity and nite dieren e results for the fourth order a will lead to less a  urate s heme; this lower level of pressure sensitivity mat hing  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  Se ond order Adjoint Gradient ve tor magnitude  ◦  0  Angle with fourth order adjoint  30.1  ◦  Fourth order  FD  1.676  Angle with se ond order adjoint  Adjoint  FD  1.576  1.672  1.753 2.845  ◦  30.1  21.05  ◦  0  Table 5.4: The magnitude of se ond and fourth order  ow  CD  ◦  ◦  18.8  ◦  18.9  ◦  gradients and angles be-  tween the evaluated gradients for NACA 0012 using higher order limiter in transoni ow.  for the unlimited transoni  ase. Also, the dieren e in nite dieren e and adjoint  gradient dire tion grows to several degrees with the use of Venkatakrishnan limiter. Table 5.4 shows the same behavior with the higher order limiter of C. Mi halak and C. Ollivier-Goo h. larger  It shows also that the error in the gradient magnitude is  ompared to the Venkatakrishnan results. The mat hing in adjoint and nite  dieren e gradient dire tions is good for the se ond order s heme but it is not for the fourth order s heme with the use of higher order limiter. It is well known that the use of limiters  auses  onvergen e problem due to the non-dierentiability of the limiting  pro edure [55℄. The limiter ae ts both the right and left hand sides of Eq. 1.6. We spe ulate that the poor mat hing between nite dieren e and adjoint gradients for the high-order limiter is related to the non-dierentiability of the limiting pro edure and that this ee t will vary in strength from one limiting pro edure to another. To redu e the inuen e of  ∂Ri /∂D  in  ontrol volume  i  when  al ulating the  adjoint gradient, the fourth order Ja obian was modied numeri ally by making the non zero stru ture of the fourth order Ja obian matrix the same as the non zero stru ture of the se ond order Ja obian, and dropping the rest of values in the fourth order Ja obian matrix.  The right hand side is still  66  onstru ted with fourth order  Mod 4th order Limiter  Venkat.  Gradient s heme  Adj.  Gradient ve tor magnitude  1.735  Angle with mod 4th order Venkat. adjoint Angle with mod 4th order HO. adjoint  0  ◦  6.247  ◦  HO  FD  Adj.  1.812  1.966  2.845  ◦  6.247  7.256  ◦  0  ◦  ◦  FD 1.672 3.338  ◦  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  a  ow.  ura y. The above modi ation doesn't ae t the a  be ause the right hand side remains fourth order a  ura y of the CFD simulation  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 ase, transoni  ow with limiter.  ǫ is presented for the hardest  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 order limiter).  68  ow (higher  ǫ 10−3 10−6 10−9  Se ond order  ∂CD ∂D  Fourth order  1.90237  1.86814  1.90275  1.87105  1.09278  1.87114  ∂CD ∂D  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  both test  and the other transoni .  In  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  RAE 2822 airfoil. The obje tive of this test  ◦ angle of atta k = 2 . In this test  onstraint starting from the  ase is to minimize  CD  at M = 0.73 and  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 in this  ase, we  ase repeats test  ase three but adds the lift  oe ient as a  onstraint;  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 order CFD simulations.  ondition,  M = 0.5  and  α = 2◦  , using se ond and fourth  The starting geometry is the parametrized NACA 0012.  70  The obje tive fun tion to be minimized is  F =  ˛  (PT − Pi )2 dS  The above obje tive fun tion and its gradient  F =  where  ws  X  X dF = 2 (PT − Pi ) dxd  is the ar  (6.1)  an be expressed in dis rete form as  (PT − Pi )2 ws    −∂Pi ∂xd    (6.2)  ws + (PT − Pi )2  ∂ws ∂xd  (6.3)  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 rea hed the target pressure distribution. Figures 6.2 and 6.3 show the history and gradient norm history for both s hemes.  essfully  onvergen e  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 gradient evaluation strategy.  Figure 6.6 shows the initial, target, and the optimized  71  (a) Se ond order Figure 6.1: Subsoni  (b) Fourth order  NACA 2412 inverse design pressure distributions for the initial,  target, and optimized airfoil proles  −2  10  Fourth order Second order −4  10  −6  F  10  −8  10  −10  10  −12  10  0  5  10  15  20  25  30  35  Iteration  Figure 6.2: Se ond and fourth order optimization  72  onvergen e history.  −1  10  Fourth order Second order −2  First−order optimality  10  −3  10  −4  10  −5  10  −6  10  0  5  10  15  20  25  30  Iteration  Figure 6.3: Gradient norm history  Figure 6.4: Subsoni  inverse design optimal airfoil shapes  73  35  −5  1.5  −5  Upper surface Dy  x 10  1  Lower surface Dy  x 10  4th order 2nd order  4th order 2nd order 1  0.5  0.5  Dy  Dy  0 0  −0.5 −0.5 −1  −1  −1.5  0  0.2  0.4  0.6  0.8  −1.5  1  0  0.2  0.4  0.6  x/c  (a) Upper surfa e Figure 6.5:  0.8  1  x/c  (b) Lower surfa e  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 obje tive fun tion dropped only three order of magnitudes before  onverge. The  onvergen e stall.  This stall is due to the high non-linearity in the target pressure distribution be ause of the existen e of a strong sho k wave in it. The noise in the pressure sensitivity generated by the sho k wave in the target pressure distribution doesn't allow further onvergen e; however, the gradient magnitude dropped four order of magnitudes from its initial value. Figures 6.9 and 6.10 show that the optimal shapes for the two s hemes dier by less than  −3 than 10  10−4  of the  hord length on the lower surfa e and less  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 lift  ase, minimization of the drag  onstraint applied.  oe ient will be  arried out with no  The airfoil to be optimized is an RAE 2822 at transoni  74  (a) Se ond order Figure 6.6: Subsoni  (b) Fourth order  NACA 2412 inverse design pressure distributions for the initial,  target, and optimized airfoil proles  Figure 6.7: Se ond and fourth order optimization  75  onvergen e history.  0  10  2nd order 4th order  −1  Grade Norm  10  −2  10  −3  10  −4  10  0  5  10  15  20  25  30  35  40  45  iteration  Figure 6.8: Gradient norm history  (a) Upper surfa e Figure 6.9:  (b) Lower surfa e  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  CD  CL  optimized airfoil  Se ond order  0.865  0.765  0.0081  Fourth order  0.849  0.759  0.0099  Table 6.1: Aerodynami transoni  onditions:  CD  optimized airfoil 0.0046  6  0.0047  oe ients of original and optimized RAE 2822 airfoil at  onditions.  M = 0.73  and  α = 2◦ .  sho k wave standing at 70%  Fig. 6.11 shows the initial solution with a strong  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  airfoil is obtained from fourth order a  CD  of the se ond order optimized  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.035  0.012  0.03  0.01  0.025 0.008  D  C  C  D  0.02 BFGS optimal profile  0.006  0.015 0.004 0.01 BFGS optimal profile 0.002  0.005  0  0 −grade  5 10 15 20 −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−−  0  25 +grad  0 −grade  (a) Se ond order response surfa e  5 10 15 20 −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−−  25 +grad  (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  without lift  drag minimization  onstraint.  The number of CFD simulations for the se ond order s heme is 50 for the fourth order s heme, yet the overall s heme is 40% less than the fourth order  omputational  ompared to 45  ost of the se ond 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 where test  C LC  is the original lift  ase both geometri  oe ient of RAE 2822 at  and aerodynami  M = 0.73 and α = 2◦ .  In this  onstraints are applied to avoid getting a  non-feasible airfoil geometry, while not ae ting airfoil lift onstraint, the lift  (6.4)  oe ient. The nonlinear  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. se ond order s heme  osts 112 CFD simulations to rea h its optimum  The  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  Se ond order  0.865  0.871  0.0081  Fourth order  0.849  0.853  0.0099  Table 6.2: Aerodynami transoni  CD  optimized airfoil 0.0047  7  0.0051  oe ients of original and optimized RAE 2822 airfoil at  onditions.  Figure 6.17: Surfa e pressure distribution, of the initial and optimized geometries  Figure 6.18: Se ond and fourth order optimization  82  onvergen e history.  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  8  0.7  7  0.6  6 0.5  5  F  F  0.4  4  0.3  3 0.2  2  BFGS Optimal Profile  1 0  0.1  BFGS Optimal profile  0 −grade  5 10 15 20 −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−−  0  25  0 −grade  +grad  (a) Se ond order response surfa e  5 10 15 20 −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−−  25 +grad  (b) Fourth order response surfa e  Figure 6.21: Obje tive fun tion response surfa e along positive and negative gradient dire tions without lift  entered at the optimal prole of RAE 2822 transoni  drag minimization  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 in the obje tive fun tion and lift  onstraint behavior.  auses high non-linearity  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 the fourth order prole drag  Cd  is four drag  ounts less than  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 amber line of the RAE 2822 via 20 distribution of our se ond order  ontrolling the  ontrol points. Comparing the optimal pressure  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 eleration 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 lift  onstraint by one lift  obje tive fun tion value.  −2 ) ount (10  ase is sele ted su h that violating the  auses 10 drag  −3 ) in rease in the  ounts (10  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  85  omputations.  Lift penalty weight  Fopt  CD  5  0.0052  0.0047  0.873  10  0.0051  0.0047  0.871  20  0.0082  0.0080  0.868  opt  CL  opt  Lift penalty  3.2 × 10−4 3.6 × 10−4 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 when  ase, Table 6.2 shows that the optimal value of  CD  is dierent  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.  CLC = 0.84  is used as a lift  omputations.  Figure 6.23  In all test  ases,  onstraint value for both se ond and fourth order  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  on order of the dis retization error).  ount dieren e (whi h is  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 grid.  3 · 10−5  on the  87  oarsest mesh to  2 · 10−5  on the nest  Chapter 7  PARTICLE SWARM OPTIMIZATION AND A NEW HYBRID OPTIMIZATION METHOD 7.1 Introdu tion to Gradient Free Optimization Gradient based optimization results depend on the optimization starting point. To nd a global optimal solution, gradient independent methods are an obvious  an-  didate.  Perhaps the best-known non-gradient based optimization te hnique is the  geneti  algorithm [30℄, a biologi ally-inspired approa h whi h  an nd the global  optimal point of an obje tive fun tion with multiple lo al optima. gorithm optimization, a randomly generated population of  In geneti  al-  andidate solutions are  generated; their tness is determined by the value of the obje tive fun tion [29℄. Solution points with the best tness are sele ted to be parents of the next generation. A new generation of ospring is generated by geneti  rossover and mutation of the parents'  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 shows a ow hart of a simple geneti  onvergen e will result. Figure 7.1  algorithm optimizer. The dependen y of the  population size on the number of design variables and the large number of generations required to  onverge to the global optimum  impra ti ally expensive for large aerodynami  88  ombine to make geneti  algorithm  optimization problems. Geneti  algo-  rithm optimization has been used for aerodynami resear hers. with a lift nami  A geneti  optimization problems by many  algorithm was used for transoni  wing drag minimization  onstraint by Gage and Kroo [29℄. Anderson used GA for wing aerody-  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 Navier-Stokes solver for transoni  wing optimization [70℄.  use of fra tal analysis in GA aerodynami  algorithm with a  They also explored the  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 optimization method (des ribed in Se tion 7.2) bridizing it with sequential quadrati  algorithm. The parti le swarm  an be made more e ient by hy-  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 partile 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  → −  − → g ),  fun tion value for ea h parti le ( pi ) as well the global best lo ation ( the best lo ation among the  − → pi 's.  Algorithm 7.1 is the general pseudo  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  whi h is  ode of the  Figure 7.1: Geneti  Algorithm optimization  90  Algorithm 7.1  For ea h parti le  Parti le swarm pseudo  i  ode  in the swarm  − → xi  Initialize parti le position  Set its best position to be its initial position: Initialize its velo ity  − ~ → Vi = 0  End Set the global best position:  − → → g =− pi  where  − → → pi = − xi  − − F (→ pi ) < F (→ pj )  for all  i 6= j  Do For ea h parti le  i  − → → − x  Cal ulate its velo ity Vi Update its position  i  Cal ulate its obje tive fun tion value (tness) If the  − F (→ xi )  − − (F (→ xi ) < F ( → pi )) → − → − position: pi = xi  urrent tness value is better than its best tness  Set the  urrent parti le position to be its best  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  91  onvergen e  riterion not met.  − → pi  ˆ  Its distan e from its known best value  ˆ  Its distan e from the swarm global best value  ( ognitive part)  − → g  (so ial part)  Figure 7.2 shows s hemati ally how the above three fa tors ae t the parti le velo 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)} → − − → → xi (k + 1) = − xi (k) + Vi (k + 1), where  w = 0.73  and so ial a to Cler 's  is the inertial weight fa tor,  onstri tion models [16℄.  − → r1  and  − → r2  i = 1..N  c1 = c2 = 1.49  eleration fa tors respe tively; these  (7.1)  (7.2) are the  ognitive  onstants are sele ted a  ording  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℄. premature  To avoid this problem, a restart strategy is suggested when  onvergen e is dete ted [9℄. At ea h restart, swarm parti les start from  new positions in the design spa e and arried out te hnique is  m  onverge to an optimal solution; restarting is  times, and the optimal solution is the best of all global bests. This  alled the multi-start parti le swarm optimization (MPSO). A drawba k  of this me hanism is that restarting may  ause repetition of the sear h  omputations.  Evers and Ben Ghalia suggested a more e ient me hanism to es ape from high quality lo al wells. They re-s atter the parti les around the  onvergen e point in a  region with radius smaller than the design spa e norm, and large enough some at least some parti les to be outside the basin of attra tion of the lo al optimum [24, 23℄; their s heme is  alled the regrouped parti le swarm optimization (RegPSO). The  main three elements in their s heme are  ˆ  Prematurity dete tion.  92  Figure 7.2: Velo ity  omponents and position update of a parti le  ˆ  Redening upper and lower boundaries for de ision (design) variables.  ˆ  Regrouping the parti les by s attering them on a sear h domain dened by the new upper and lower boundaries.  In the following subse tions, the above three elements of the RegPSO strategy will be presented; interested readers are en ouraged to for more detail. All  onsult Evers and Ben Ghalia [23℄  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 neighborhood 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 variables  h  x1U  iT  h  n and m parti  les, the de ision  − → 0 have initial upper and lower limits xU = x1 x2 x3 · · · xn iT iT h − → 0 and xL = . First, x3U · · · xnU x1L x2L x3L · · · xnL  − → x =  x2U  the range is dened as  − → − → − −−→ range (Ω) = x − x U  Also the normalized swarm radius  δnorm =  (7.3)  L  an be found a  ording to  → → maxi=1...m k− xi (k) − − g (k)k −−→ k− range (Ω)k  Convergen e is dete ted when the normalized swarm radius  ε = 1.1  (7.4)  δnorm  drops below  · 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 using its upper limit  x0jU  j  omputed. The initial range of parti le  and lower limit  rangej  an be found  x0jL as   Ω0 = x0jU − x0jL  Upper and lower limits of the de ision variables are  (7.5)  hanged when  onvergen e is  dete ted in order to s atter swarm parti les in a ball surrounding the solution. Changing the range of ea h de ision variable done a  j  onverged  in a regrouping stage  r  is  ording to the following formula [23℄:  rangej  When the  r  (Ω ) = max    rangej  Ω  0    −− → −− → r−1 , 1.2 × 10 max xi,j − gjr−1 4  i=1...m    (7.6)  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 lowing formula  onvergen e. Regrouping is done a  ording to the fol-  −−−−−−→ 1 −−−−−−−→ T − − → → → xi = − g + [− ri ] · range (Ωr ) − range (Ωr ) 2  (7.7)  −−−−−−−→ r ) = [range (Ωr ) , range (Ωr ) , . . . , range (Ωr )]T . 1 2 n  where range (Ω  Equation 7.7 requires position  r strained design spa e, i.e xj,L  lipping to keep the parti les within the  ≦ xi,j ≦  omponents of parti le's de ision ve tor  xrj,U .  on-  The upper and lower limits of the  an be found as    1 xrj,L = max x0j,L , gjr−1 − rangej (Ωr ) 2   1 r 0 r−1 r xj,U = min xj,U , gj + rangej (Ω ) 2  (7.8)  7.3 Hybrid SQP-RegPSO Te hnique Swarm intelligen e starts with a set of randomly generated position points in the design spa e and denes its initial best lo ation,  − → g , as the point of best tness among  the swarm parti les; this best point a t as a gravity  enter in the design spa e and  attra ts other parti les to it; if the initial best point is of high quality, this redu es the total  omputational eort to nd the global optimal solution. The te hnique we  are proposing here for hybridization of sequential quadrati  programing (SQP) and  parti le swarm optimization is to use SQP to nd the initial high quality  − → g  before  using RegPSO te hnique. The pseudo- ode of the hybrid s heme is summarized in Algorithm 7.2.  7.4 Optimization Test Cases This se tion will examine the e ien y of the proposed hybrid s heme and  ompare  it with the original RegPSO te hnique. To do so, a drag optimization that seeks a sho k free prole with  CL = 0.4  at  M = 0.8  and angle of atta k  α = 1.5◦  is  arried  out. Two distin t starting points were sele ted to explore the existen e of a unique  95  Algorithm 7.2  Hybrid SQP-RegPSO pseudo  ode  For ea h parti le in the swarm Initialize parti le position  − → xi  Set its best position to be its initial position: Initialize its velo ity  − → Vi = 0  End Set the global best position:  − → → pi = − xi  − − − → → pi ) < F (→ pj ) for all i 6= j , g =− pi where F (→  or the SQP  optimal solution, whi hever is better Set Set  δnorm = 1 t=1  Do For ea h parti le  i  − → → − x  Cal ulate its velo ity Vi Update its position  i  Clip position su h that  xL,j ≤ xi,j ≤ xU,j  Cal ulate its obje tive fun tion value (tness) If the  − F (→ xi )  − − (F (→ xi ) < F ( → pi )) → − → − position: pi = xi  urrent tness value is better than its best tness  Set the  urrent parti le position to be its best  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 δnorm δnorm ≤ ε  Compute If  using Eq 7.4  Compute the range of ea h de ision variable range S atter the parti les around  Apply position  − → g  a  ording to Eq 7.7  lipping with limits dened a  End if While (δnorm  >ε  and  t ≤ max    Ωrj  iterations).  96  ording to Eq 7.8  8  global best ; the rst is to start from NACA 0012, the se ond one is to start from  9 These ases were hosen be ause experiments showed that SQP  NACA 00083 airfoil.  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  M = 0.8  oe ient of CL = 0.4, ◦ 2 is presented. The obje tive fun tion an be  ase, drag optimization of the NACA 0012 at a lift  and angle of atta k  α =  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) the se ond phase  osts 92 CFD simulations while  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. original NACA 0012 at transoni  Fig 7.4.1 shows the pressure distribution of the 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 airfoil, and the optimized prole using SQP with BFGS approximate Hessian, RegPSO, 8  In the  ase of a nite region in the design spa e for whi h the ow is sho k free, there may be  a region of global best solutions diering only in the amount of dis retization error present.  9  That is, an airfoil from the 4-digit NACA series with 8.3% maximum thi kness.  97  NACA 0012 Original  BFGS Optimal  RegPSO 1000 CFD simulations  Hybrid 492 Total CFD simulations  Figure 7.3: Pressure Field of the original NACA 0012 and the optimized proles  Figure 7.4: Surfa e pressure distribution of original NACA 0012 and the optimized proles  98  Figure 7.5: Airfoil shapes of original NACA 0012 and the optimized proles  and hybrid BFGS-RegPSO te hniques. Figure 7.5 shows the dierent optimized airfoil 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 optimization, 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 dire tion (the ve tor that but the SQP optimizer 1. An error in the  omputed. Figure 7.7 shows that there is a des ent  onne ts the initial and nal  − → g ) at the SQP optimal point  an not nd it. Possible reasons for that in lude:  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  ˆ  Adjoint  ˆ  Finite dieren e  ˆ  Ve tor that  omputed the angles between the  omputed gradient ve tor.  omputed gradient ve tor.  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.02  0.018  BFGS optimal  0.016  0.014  F  0.012  0.01  0.008  0.006  Hybrid optimal  0.004  0.002  0  0  20  40  60 BFGS Optimal to Hybrid optimal steps  80  100  120  Figure 7.7: Obje tive fun tion values between SQP and hybrid optimal points (starting geometry: NACA 0012).  Angle in degrees Adjoint and nite dieren e gradients  1.78  Adjoint gradient and SQP-hybrid ve tor  86.6  Finite dieren e gradient and SQP-hybrid ve tor  86.7  Table 7.1: Angles between adjoint, FD., and SQP-hybrid ve tors  101  7.4.2 Constrained drag minimization of NACA 00083 In this test  ase, the same drag optimization problem is studied at the same transoni  onditions, but starting from a thinner airfoil. Again, SQP  ould not eliminate the  sho k wave on the airfoil's upper surfa e, but a sho k free optimal shape is obtained using RegPSO in the se ond phase as shown in Fig 7.8. Figure 7.9 shows a  omparison of the surfa e pressure of the original NACA  00083, SQP with BFGS approximate Hessian (lo al) optimum, and hybrid BFGSRegPSO (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 pletely. Figure 7.10 shows a prole  om-  omparison of the original NACA 00083, BFGS  optimal, and the hybrid optimal prole. The drag to 8.89 drag  oe ient is redu ed from 32 drag  ounts at the BFGS optimal solution  ounts at the hybrid optimal solution. The optimal value of the drag  oe ient found by the hybrid s heme in this test obtained in the previous  ase and in both  ase is the same as the value  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) optimum, 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 appearan 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  nal best position ve the ve tor  − → tors g obtained using the hybrid s  onne ts the initial and the heme. Table7.2 shows that  onne t the initial and nal best solutions is almost perpendi ular to the  adjoint and nite dieren e gradients; it shows also that adjoint and nite dieren e gradients are a very good mat h. The orthogonality of these ve tors was previously noted in the last test  ase as well. This strongly suggests that numeri al noise in  the gradient and the inuen e of the ow model physi s prevents nding a des ent  102  NACA 00083 original  BFGS Optimal  Hybrid Optimal Figure 7.8: Pressure surfa e of NACA 00083 and the optimized proles  Angle in degree Adjoint and Finite dieren e  0.27  Adjoint and SQP-hybrid ve tor  86.93  Finite dieren e and SQP-hybrid ve tor  86.87  Table 7.2: Angles between adjoint, FD, and BFGS-hybrid best solutions ve tors  dire tion that leads from the SQP lo al optimum to the sho k free global optimum.  7.4.3 Aerodynami and thi kness onstraint drag minimization of NACA 0012 In this test the lift  ase thi kness is  onstraint  CLc = 0.4  onstrained to be at transoni  tc = 0.1  onditions,  at  xtc = 0.35,  M = 0.8  and  in addition to  α = 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.025  0.02  F  0.015  0.01  BFGS optimal 0.005 Hybrid optimal  0  0  20  40  60 BFGS optimal to hybrid optimal steps  80  100  120  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 subje t tot (xtc )  The thi kness  onstraint  an be  surfa e parametrization matrix onstraint equation  = tc  ast as an equality  Q2 [AQ1 ]  †  (7.10)  ; for this test  onstraint equation using the ase,  xc < L1 , so the thi  an be written as  t (xtc ) = P1 (xtc ) − P3 (xtc ) = tc   √ √ ∴ a0 xc + a1 xc + a2 x2c + a3 x3c − c0 xc + c1 xc + c2 x2c + c3 x3c = tc  Re all that the relationship of polynomial be found using  Q2 [AQ1 ]  †  an be used as a  oe ients and the design variables  X  Q2 [AQ1 ]†    i,j  yj  an  (7.12)  an be used with Eq (7.11) to derive an algebrai  that relates the thi kness  (7.11)  as  ai = Equation (7.12)  kness  onstraint  tc  linear equation  to design variables dire tly; this last equation  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 The penalty method is used to approximately satisfy the thi kness  (7.13)  onstraint term  due to the la k of Lagrange multiplier in non-gradient based optimization te hniques. Drag  oe ient have been redu ed from 230 to 120 drag  then redu ed to 9 drag  ounts in the rst phase,  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 kness  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 The SQP optimizer is used rst;it satises the thi kness  onstraint.  onstraint after the rst few  iterations, then it starts to minimize the drag. Drag is redu ed from 60 drag to 40 drag to 9 drag  ounts  ounts in the rst phase. The se ond (RegPSO) phase redu ed the 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 Fig. 7.17. The two proles are very  ases are shown in  lose in upper surfa e shape, whi h is tightly  onstrained by the need to eliminate the sho k wave; there is more dieren e in the lower surfa e on whi h the ow is subsoni . The last observation suggests that for invis id ow simulations, there is a region of global optimal solution rather than one point global optimum.  More pre isely, on e a sho k-free solution is obtained, the  obje tive fun tion is too at to distinguish reliably between the drag of sho k-free solutions; this is not surprising,  onsidering that the drag for sho k-free invis id ows  is due solely to dis retization error. In all test  ases, the new hybrid s heme was able to nd a sho k free airfoil prole  regardless of the starting points. Only ten parti les were using in the swarm to s an the design spa e, whi h is roughly half of the number of design variables. This is lower than the re ommended population size required by GA to  108  onverge well. The  Figure 7.16: Surfa e pressure of NACA 00083 and the optimized prole with thi kness  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 in rease the value of Figure Shows a  onstraint, but it will  CD  and may prevent obtaining a sho k free optimal prole. opt 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  109  ounts whi h is  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  in rease in the drag  oe ient to 17 drag  CL  ounts.  Fopt  CDopt  1  0.00134  0.00088  0.3997  10  0.00208  0.0017  0.4005  Thi k. penalty  auses an  opt  t (xtc )  Lift penalty  Thi k. penalty  0.077  0.000001  0.00048  0.095  0.0000026  0.00025  weight  Table 7.3: Thi kness penalty weight impa t on the optimization results of NACA 00083 using hybrid s heme  110  1.5 Penalty weight 10 Penalty weight 1 1  −Cp  0.5  0  −0.5  −1  −1.5  0  0.1  0.2  0.3  0.4  0.5 x/c  0.6  0.7  0.8  0.9  1  Figure 7.18: Thi kness penalty weight impa t on the Optimal pressure distribution of NACA 00083 using hybrid optimization s heme.  0.04 Thickness penalty 1 Thickness renalty 10 0.03  0.02  0.01  y/c  0  −0.01  −0.02  −0.03  −0.04  −0.05  −0.06  0  0.1  0.2  0.3  0.4  0.5 x/c  0.6  0.7  0.8  0.9  1  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  T  ode of RegPSO-SQP hybrid  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 nd an estimated global best 5. Use SQP with starting point  − → g.  − → g  riterion of max 900 CFD simulations to  to nd better global best solution.  7.5 Comparing SQP-RegPSO with RegPSO-SQP optimization strategies I have presented results of the proposed hybrid optimization s heme in whi h SQP optimization is used to obtain a reasonably good upper bound on the global optimum. Based on the results of previous test  ases, this s heme seems to be a promising  global optimization s heme. Another hybrid s heme in whi h RegPSO pre edes SQP optimization is possible. In this se tion I will  ompare the optimization results of the  last ve optimization problems using the proposed hybrid s heme with the results  T  of RegPSO-SQP s heme; I will refer to the latter s heme as hybrid . 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  T show that the hybrid s heme is better than the hybrid for all ases; in Case II, the s hemes are nearly identi al in drag was expe ted be ause in hybrid  oe ient though dierent in nal shape. This  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  10  0.0150  0.0009  0.0600  0.0032  0.0009  0.2813  0.0120  0.0013  0.1083  0.0040  0.0013  0.3250  0.0040  0.0021  0.525  I  11  II  12 III 13 IV 14  V  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  T s heme in Phase I, II  Table 7.5: Optimization results of the hybrid  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  T  better than hybrid .  113  1.5 RegPSO−SQP SQP−RegPSO 1  −Cp  0.5  0  −0.5  −1  −1.5  0  0.2  0.4  0.6  0.8  1  x/c T optimal pressure distribution for NACA 0012 (Case  Figure 7.20: Hybrid and hybrid I).  114  1 RegPSO−SQP SQP−RegPSO 0.5  −Cp  0  −0.5  −1  −1.5  0  0.2  0.4  0.6  0.8  1  x/c  Figure 7.21:  T optimal pressure distribution for NACA 00083  Hybrid and hybrid  (Case II).  115  1.5 RegPSO−SQP SQP−RegPSO 1  −Cp  0.5  0  −0.5  −1  −1.5  0  0.2  0.4  0.6  0.8  1  x/c  Figure 7.22: Hybrid and hybrid 10% thi kness  T optimal pressure distribution for NACA 0012 with  onstraint and unit thi kness penalty weight (Case III).  116  1 RegPSO−SQP SQP−RegPSO 0.5  −Cp  0  −0.5  −1  −1.5  0  0.2  0.4  0.6  0.8  1  x/c T optimal pressure distribution for NACA 00083 with  Figure 7.23: Hybrid and hybrid 10% thi kness  onstraint and unit thi kness penalty weight (Case IV).  117  1.5 RegPSO−SQP SQP−RegPSO 1  −Cp  0.5  0  −0.5  −1  −1.5  0  0.2  0.4  0.6  0.8  1  x/c T optimal pressure distribution for NACA 00083 with  Figure 7.24: Hybrid and hybrid 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  optimization. I have developed the rst nite volume based  ode that uses high order Ja obian in gradient  omputations using the  adjoint te hnique. A new geometry parametrization te hnique has been developed and presented; this new te hnique is based on least squares surfa e tting. To develop an e ient global optimizer, a new hybrid SQP-RegPSO s heme has been developed to nd the global optimal solution in the feasible design spa e. I presented both sequential quadrati  onstrained and un onstrained aerodynami  programming with the use of adjoint method to  obje tive fun tion gradient.  The  ompute the  omputations were based on higher order nite  volume method on unstru tured grids. Ja obian matrix to  optimization using  I took advantage of evaluating the exa t  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 order s hemes. the  ow for both se ond and fourth  For transoni  ow, the se ond order Ja obian mat hes well with  orresponding order of a  ura y nite dieren e gradient; on the other hand,  the fourth order Ja obian does not mat h the nite dieren e gradient well when using a limiter.  This error in the fourth order gradient is attributed to the high  sensitivity of the fourth order Ja obian matrix to the use of the limiter in transoni ow whi h tends to ae t the diagonal dominan e of the Ja obian matrix. To ta kle this problem and to enhan e the diagonal dominan e of the fourth order Ja obian matrix, we used the same non-zero stru ture of the se ond order Ja obian matrix  119  and dropped the rest of the non-zero elements; this tends to in rease the diagonal dominan e of the modied fourth order Ja obian matrix and improves its  ondition  number. I have developed a new approa h for smooth parametrization of airfoil shapes, based on a least-squares t of the parameters of two polynomials to a set of points.  ontrol  Compared with using all mesh points as design variables, this approa h  redu es the size of the design spa e and eliminates os illations in the shape. We use a semi-torsional spring analogy to deform the mesh grid in the entire ow eld when the surfa e shape evolves during optimization. It is used also to  al ulate the mesh  sensitivity terms in the adjoint gradient. To test the developed gradient based optimizer; I have used the optimizer to ta kle inverse design problems in whi h only one optimal solution point in the design spa e exists. Additional transoni  drag minimization test  ases have been pre-  sented with and without lift  onstraint. As a summary of test  subsoni  ase shows that both se ond and fourth order s hemes  inverse design test  ases observations: A  rea hed their target geometry after almost the same number of iterations. The same behavior is observed for transoni  inverse design test  onvergen e is independent of the order of a a drag minimization test  ase without lift  ase. This indi ates that the  ura y of the spatial dis retization. In  onstraint, both se ond and fourth order  s heme rea hed their optimum shape after almost the same number of optimization 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 on all the test  lo k run time and found nearly identi al airfoil optimal shapes. Based ases we presented, we  produ es a systemati  on lude that the spatial dis retization error  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 wave  apture in transoni  omputations and for a  urate predi tion of for es and sho k  ows, and the use of a se ond order method (but based on  the obtained forth order solution) to  ompute the obje tive fun tion gradient needed  for gradient based optimizer. If the mesh resolution makes se ond and fourth order simulation results similar, using the se ond order s heme for optimization will be e onomi al. I have developed also a new hybrid gradient/non-gradient based optimization te hnique that uses sequential quadrati  programming with BFGS Hessian approx-  imation te hnique in an initial gradient based optimization phase, followed by the regrouped parti le swarm optimization method in the non-gradient random sear h phase.  The SQP optimization phase leads to a high quality initial best-solution  point that redu es the total  omputational eort of the RegPSO to rea h a sho k  free prole. Unlike with the SQP s heme, the transoni  drag optimization with lift  onstraint  of the NACA 0012 and NACA 00083 ta kled with the hybrid s heme lead to sho k free optimal shapes; in both  ases, the obje tive fun tion gradient ve tor is almost  perpendi ular to the ve tor that  onne ts the SQP optimum and the hybrid opti-  mum. This explains why the gradient based optimization  an not rea h the global  optimal, and also shows that its optimal solution depends on the starting geometry. The obje tive fun tion behavior along the ve tor that  onne ts the SQP and the  hybrid optima shows that numeri al noise may have a harmful impa t on gradient based optimization as well. Transoni  drag optimization of NACA 0012 and NACA 00083 have been ta k-  led with both aerodynami  and geometri  onstraints. For RegPSO, the geometri  onstraint is satised using penalty method; in both been obtained. The last test  ases, sho k free proles have  ase (drag redu tion of NACA 00083 with a thi kness  onstraint) shows that the hybrid s heme was able to nd a sho k free prole despite  in reasing  the thi kness during optimization. The  omputational expense of the hy-  brid s heme is four to nine times more than SQP optimization, but in all  ases, a  global optimal (sho k free) prole obtained, and in aerospa e industry, this is worth the in rease in  omputational eort.  121  8.2 Future Work We intend to apply non gradient-based optimization te hniques like parti le swarm algorithm to study the ee t of CFD s heme order of a shape in the near future. dynami  Another obvious extension of the  urrent work is aero-  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 simulation with  7 × 7.  ura y on the nal optimal  Fully  k−ǫ  urate.  Thus, for a three dimensional ow  turbulen e model, the blo k size of the Ja obian matrix is  oupled solvers have also been shown to  onverge more e iently.  For three-dimensional problems, the geometry parametrization should be re-  pla ed with a perturbation parametrization (parametrization of the shape).  hange in the  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 optimization results is needed. We spe ulate that using the right geometry parametrization te hnique may lead to a design spa e in whi h aerodynami simply  onne ted and an SQP optimizer  onstraints are  an nd a globally optimal solution.  122  Bibliography [1℄ Anderson, M.B. The Potential of Geneti  Algorithms for Subsoni  Wing De-  AIAA Paper 95-3925, presented at the 1st AIAA Air raft Engineering, Te hnology, and Operations Congress, Los Angeles, CA, 1995. sign. In  [2℄ Anderson, Murray B.  Geneti  Algorithms In Aerospa e Design:  Progress, Tremendous Potential.  on Intelligent System, 2002.  In  NATO/Von-Karmen-Institute Workshop  [3℄ Anderson, W. Kyle, and Bonhaus, Daryl L. Aerodynami tured Grids for Turbulent Flows.  Substantial  Design on Unstru -  Te hni al report, NASA Langley Resear h  Center, Hampton, Virginia, 1997.  [4℄ Avriel, M.  Nonlinear Programming Analysis and Methods.  Dover Publishing.  ISBN 0-486-43227-0., 2003.  [5℄ Azab, M. B., and Ollivier-Goo h, C. F. Higher Order Two Dimensional Aerodynami  Optimization Using Unstru tured Grids and Adjoint Sensitivity Compu-  AIAA-2010-1430, 48th AIAA Aerospa e s ien es meeting, Orlando, FL, USA, 4-7Jan 2010, Orlando FL, 2010.  tations. In  [6℄ Azab, M. B., and Ollivier-Goo h, Carl. Constrained and Un onstrained Aerodynami  Quadrati  Programming Optimization Using High Order Finite Volume  Method and Adjoint Sensitivity Computations. In  ing, Orlando FL, AIAA-2011-0183, 2011. [7℄ Batina, J. T. Meshes.  49th Aerospa e s ien e meet-  Unsteady Euler Airfoil Solutions Using Unstru tured Dynami  AIAA Journal ,, 28 No 8.:13811388, 1990. 123  [8℄ Beall, M. W., Walsh, J., and Shephard, Mark S. A  essing CAD Geometry for  Pro eedings of the 12th International Meshing Roundtable, Sandia National Laboratories, pages 20033030, 2003.  Mesh Generation. In  [9℄ Bergh, F., and Engelbre ht, A. Optimiser. In  A New Lo ally Convergent Parti le Swarm  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  Aerospa e S ien e and Te hnology,  Using the Adjoint Approa h.  8:715727,  2004. [12℄ Broyden, C. G. The Convergen e of a Class of Double-rank Minimization Algorithms.  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.  ing, 16, 5:11901208., 1995.  SIAM Journal on S ienti and Statisti al Comput-  [14℄ Cavallo, P. A., Hosangadi, A., Lee, R. A., and Dash, S. M. Dynami  Unstru -  tured Grid Methodology with Appli ation to Aero/Propulsive Flow Field. In  AIAA Paper 1997-2310, 2002. [15℄ Chandrashekarappa,  P.,  and Duvigneau,  R.  Metamodel-Assisted Parti le  Swarm Optimization and Appli ation to Aerodynami  Shape Optimization. Re-  sear h Report RR-6397, INRIA, 2007. [16℄ Cler , M., and Kennedy, J.  The Parti le Swarm - Explosion, Stability, and  Convergen e in a Multidimensional Complex Spa e.  IEEE Transa tions on, 6(1):5873, 2002. [17℄ Consentino, G., and Holst, T. Transoni  Evolutionary Computation,  Numeri al Optimization Design of Advan ed  Wing Congurations. In  AIAA 85-0424, 1985. 124  [18℄ Cook, P.H., M Donald M.A., and Firmin, M.C.P. . Aerofoil RAE 2822: Pres-  Experimental Data Base for Computer Program Assessment, AGARD Report AR 138. sure Distributions, and Boundary Layer and Wake Measurements. In  AGARD, 1979.  [19℄ Dadone, Andrea, Mohammadi, Bijan, and Petruzzelli, Ni ola. In omplete Sensitivities and BFGS Methods for 3D Aerodynami  Shape Design.  Resear h  Report RR-3633, INRIA, 1999. Projet M3N.  [20℄ Duvigneau, R.  Aerodynami  Shape Optimization with Un ertain Operating  Conditions using Metamodels. Resear h Report RR-6143, INRIA, 2007.  [21℄ Eberhart, R. C., and Kennedy, J. A New Optimizer Using Parti le Swarm The-  Sixth International Symposium on Mi ro Ma hine and Human S ien e, Nagoya, Japan, pp. 39-43., 1995. ory. In  [22℄ Engelbre ht, A. P.  Fundamentals of Computational Swarm Intelligen e.  John  Wiley and Sons, 2006.  [23℄ Evers, G. I., and Ben Ghalia, M. Regrouping Parti le Swarm Optimization a New Global Optimization Algorithm with Improved Performan e Consisten y A ross Ben hmarks. In  SMC IEEE International Conferen e, pages 39013908,  o t. 2009.  [24℄ Evers, George.  An Automati  Regrouping Me hanism to Deal with Stagna-  tion in Parti le Swarm Optimization. Master's thesis, University of Texas-Pan Ameri an, 2009.  [25℄ Farhat, C., Degand, C., Koobus, B., and Lesoinne, M. Two-Dimensional Dynami  Unstru tured Fluid Meshes.  applied me hani s and engineering, 163(1-4):231245,  Torsional Springs for  Computer methods in  1998.  Three-Dimensional Upwind S heme for Solving The Euler Equations on Unstru tured Tetrahedral Grids. PhD thesis, Virginia Polyte hni In-  [26℄ Frink, N. T.  stitute and State University, 1991.  125  [27℄ Frink, N. T., Parikh, P., and Pirzadeh, S. A Fast Upwind Solver for The Euler Equations on Three-Dimensional Unstru tured Meshes. In  AIAA paper 91-0102,  1991. [28℄ Gage, P., and Kroo, I. Development of the Quasi-Pro edural Method for Use in  AIAA Paper No. 92-4693, Presented at the Fourth AIAA/USAF/NASA/OAI Symposium on Multidis iplinary Analysis and Optimization, 1992. Air raft Conguration Optimization. In  [29℄ Gage, P., and Kroo, I. A Role of Geneti Environment. In  AIAA Paper No. 93-3933, 1993.  [30℄ Goldberg, David E.  Learning.  Geneti Algorithms in Sear h, Optimization, and Ma hine  Addison-Wesley Publishing Company, 1989.  [31℄ Goldfarb, D. Fa torized Variable Metri tion.  Algorithms in a Preliminary Design  Methods for Un onstrained Optimiza-  Mathemati s of Computation, 30(136):796811,  [32℄ Gregg, R.D., and Misegades, K.P. Transoni tion Theory. In  87-0520, 1987.  1976.  Wing Optimization Using Evolu-  presented at the AIAA 25th Aerospa e S ien es Meeting, AIAA  [33℄ Hart, Willian E. .  Adaptive Global Optimization with Lo al Sear h.  PhD thesis,  University of California, San Diego, 1994. [34℄ Hassan, R., Cohanim, A. , and de We k,O. A Comparison of Parti le Swarm Optimization and The Geneti  Algorithm.  In  AIAA-2005-1897,  page 1897,  2005. [35℄ Hi ks, R. M. , and Henne, P. A. .  Wing Design by Numeri al Optimization.  Journal of Air raft, 15:407412., 1978. [36℄ Ira H. Abbott and Albert Edward Von Doenho.  Theory of Wing Se tions.  Dover, 1959. [37℄ Jameson, A.  Aerodynami  Design via Control Theory.  Computing, 3:233260, 1988.  10.1007/BF01061285.  126  Journal of S ienti  [38℄ Jameson, A.  A Perspe tive on Computational Algorithms for Aerodynami  Analysis and Design.  Progress in Aerospa e S ien es, 37(2):197243, 2001.  doi:  DOI: 10.1016/S0376-0421(01)00004-5. [39℄ Jameson, A.  Aerodynami Shape Optimization Using The Adjoint Method.  VKI  Le tur seriese, 2003. [40℄ Jameson, A. Optimum Transoni  Wing Design Using Control Theory. In  TAM Symposium held in Guttingen, Germany 2-6 September 2002,  IU-  page 253.  Springer Netherlands, 2003. IUTAM Symposium Transsoni um IV: pro eedings of the IUTAM Symposium held in Guttingen, Germany 2-6 September 2002. [41℄ Jameson, A. , and Reuther, J.  Control theory based airfoil design using Eu-  AIAA/USAF/NASA/ISSMO Symposium on Multidis iplinary Analysis and Optimization, 1994.  ler equations. In  [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 Using The Navier-Stokes Equations. In  & Exhibit, Reno, NV, 1997. [44℄ Jang, M., and Lee, J.  Geneti  Design  35th AIAA Aerospa e S ien es Meeting  Algorithm Based Design of Transoni  Air-  AIAA Paper 2000-1584, Presented at the 41st AIAA/ASME/ASCE/AHS/ASC Stru tures, Stru tural Dynami s, and Materials, 2000.  foils Using Euler Equations. In  Appli ation Of The Dis rete Adjoint Method To Coupled Multidis iplinary Unsteady Flow Problems For Error Estimation And Optimization. PhD  [45℄ Karthik, M.  thesis, Department of Me hani al Engineering, The University of Wyoming, 2009.  Pro eedings of IEEE International Conferen e on Neural Networks. IV. pp. 1942-1948, 1995.  [46℄ Kennedy, J., and Eberhart, R. Parti le Swarm Optimization. In  [47℄ Kroo, I.  Drag Due to Lift: Con epts for Predi tion and Redu tion.  Review of Fluid Me hani  , 33:587617, 2001.  127  Annual  [48℄ Kulfan, B. M. Re ent Extensions and Appli ations of the CST Universal Parametri  Geometry Representation Method. In  Integration and Operations Conferen e, 2007.  7th AIAA Aviation Te hnology,  [49℄ Kulfan, B. M., Bussoletti, J. E. . Fundamental Parametri tations for Air raft Component Shapes. In  Geometry Represen-  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.  timization, 24:3850, 2002.  Stru tural and Multidis iplinary Op-  10.1007/s00158-002-0212-4.  [51℄ Martineau, D. G., Georgala, J. M.  A Mesh Movement Algorithm for High  42nd AIAA Fluid Dynami s Conferen e and Exhibit, AIAA Paper 2004-0614, Reno, NV, 2004.  Quality Generalised Meshes.  In  [52℄ Masuda, H., Yoshioka, Y., and Furukawa, Y. Intera tive Mesh Deformation Us-  Computers and Graphi s,  ing Equality-Constrained Least Squares.  30(6):936  946, 2006. [53℄ Mengistu, Temesgen, and Ghaly, Wahid.  Aerodynami  Optimization of Tur-  boma hinery Blades Using Evolutionary Methods and ANN-Based Surrogate Models.  Optimization and Engineering,  9:239255, 2008. 10.1007/s11081-007-  9031-1.  E ient High-Order A urate Unstru tured Finite-Volume Algorithms for Vis ous and Invis id Compressible Flows. PhD thesis, The University  [54℄ Mi halak, C.  of British Columbia, 2009. [55℄ Mi halak, C., and Ollivier-Goo h, C.  Dierentiability of Slope Limiters on  Fourteenth Annual Conferen e of the Computational Fluid Dynami s So iety of Canada, 2006.  Unstru tured Grids.  In  [56℄ Mi halak, C., and Ollivier-Goo h, C. A Order A  ura y Preserving Limiter for The High-  urate Solution of The Euler Equations.  Physi s, 228(23):86938711,  2009.  128  Journal of Computational  [57℄ Mi halak, Christopher, and Ollivier-Goo h, Carl. Newton-GMRES for The High-Order A  Computers and Fluids, 39:11561167,  Globalized Matrix-Expli it  urate Solution of The Euler Equations.  2010.  [58℄ Mousavi, A., Castonguay, P., and Nadarajah, S. K. Survey of Shape Parameterization Te hniques and its Ee t on Three-Dimensional Aerodynami Optimization. In  Shape  Pro eedings of the AIAA 18th Computational Fluid Dynami s  Conferen e. AIAA 2007-3837, 25 - 28 June 2007 2007. [59℄ Mouser, C. R., and Dunn, S. A.  Comparing Geneti  Algorithms and Parti-  le Swarm optimisation for an Inverse Problem Exer ise.  In Rob May and  Pro . of 12th Computational Te hniques and Appliations Conferen e CTAC-2004, volume 46, pages C89C101, Mar h 2005.  A. J. Roberts, editors,  http://anziamj.austms.org.au/V46/CTAC2004/Mous [60℄ Nadarajah, S., and Jameson, A.  [Mar h 21, 2005℄.  A Comparison of The Continuous and Dis-  rete Adjoint Approa h to Automati  Aerodynami  Optimization.  AIAA paper,  667:2000, 2000. [61℄ Nejat, A., and Ollivier-Goo h, C. A High-Order A  urate Unstru tured Finite  Volume Newton-Krylov Algorithm for Invis id Compressible Flows.  Computational Physi s, 227(4):25822609, [62℄ Neme , M., and Aftosmis, M. J. Embedded-Boundary  2008.  Adjoint Sensitivity Computations for an  Cartesian Mesh Method.  Physi s, 227(4):27242742,  Journal of  Journal of Computational  2008.  [63℄ Neme , M., and Zingg, D. W. .  Newton-Krylov Algorithm for Aerodynami  Design Using the Navier-Stokes Equations.  AIAA Journal, 40:11461154, 2002.  [64℄ Nielsen, E., and Anderson, K. Aerodynami  Design Optimization on Unstru -  tured Meshes Using the Navier-Stokes Equations. Te hni al report, NASA Langley 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.  Aerospa e S ien es Meeting and Exhibit Reno, Nevada, 2005. 129  In  43rd AIAA  [66℄ No edal, J., and Wright Stephen, J.  Numeri al Optimization.  Springer New  York, 2006. [67℄ Obayashi, S. Geneti lems. In  Algorithm for Aerodynami  Inverse Optimization Prob-  Geneti Algorithms in Engineering Systems: Innovations and Appli a-  tions, 1995.  [68℄ Ollivier-Goo h, C. and Van Altena, M. A High-Order-A  urate Unstru tured  Mesh Finite-Volume S heme for The Adve tion-Diusion Equation.  Computational Physi s, 181(2):729752,  Journal of  2002.  [69℄ Ollivier-Goo h, C., Nejat, A., and Mi halak, C.  On Obtaining and Verifying  High-Order Finite-Volume Solutions to the Euler Equations on Unstru tured Meshes.  AIAA Journal, 47(9):21052120,  2009.  [70℄ Oyama, A., Obayashi, S., and Nakahashi, K. Using Geneti  Algorithm. In  Dynami s Conferen e, 1997.  oding for aerodynami  optimization.  AIAA Weakly Ionized Gases Workshop, 1999. [72℄ Piegl, L., and Tiller, W.  Wing Optimization  AIAA Paper 97-1854, 13th Computational Fluid  [71℄ Oyama, A., Obayashi, S., and Nakahashi, K. geneti  Transoni  The NURBS Book.  Fra tional fa torial design of In  AIAA Paper 99-3298, 3rd  Springer, Berlin, 1995.  [73℄ Reuther, J., Jameson, A., Farmer, J., Martinelli, L., and Saunders, D. Aerodynami  Shape Optimization of Complex Air raft Congurations via an Adjoint  Formulation. In  34rd AIAA Aerospa e S ien es Meeting and Exhibit, no. AIAA  1996-0094, 1996.  [74℄ Roe, P. L. Approximate Riemann Solvers, Parameter Ve tors, and Dieren e S hemes.  Journal of Computational Physi s,  43(2):357372, 1981.  doi: DOI:  10.1016/0021-9991(81)90128-5. [75℄ Saad, You ef, and S hultz, Martin H. GMRES: a Generalized Minimal Residual Algorithm for Solving Nonsymmetri  Comput., 7:856869, July 1986. 130  Linear Systems.  SIAM J. S i. Stat.  [76℄ Samareh, Jamshid A.  A Survey of Shape Parameterization Te hniques.  In  CEAS/AIAA/ICASE/NASA Langley International Forum on Aeroelasti ity and Stru tural Dynami s Williamsburg, VA, Also, NASA/CP-1999-209136, June 1999, pp. 333-343., 1999. [77℄ Samareh, Jamshid A.  Geometry and Grid/Mesh Generation Issues for CFD  and CSM Shape Optimization.  Optimization and Engineering,  6:2132, 2005.  10.1023/B:OPTE.0000048535.08259.a8. [78℄ Sederberg, Thomas W., and Parry, S ott R. Free-Form Deformation of Solid Geometri  Models.  SIGGRAPH Comput. Graph., 20:151160, August 1986.  [79℄ Shi, Y., and Eberhart, R. C. Empiri al Study of Parti le Swarm Optimization. In  Pro eedings of Congress onEvolutionary Computation, CEC 99.,  volume 3,  page 1950 Vol. 3, 1999. [80℄ Stein, K., Tezduyar, T., and Benney, R. Mesh Moving Te hniques for FluidStru 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 , Extensible Geometry Interfa e.  In  Roundtable, pages 337348, 2000.  Pro eedings of the 9th International Meshing  [82℄ Tezduyar, T. E., Behr, M., Mittal, S., and Johnson, A. A.  Computation of  Unsteady In ompressible Flows with The Stabilized Finite Element Methods: Spa e-Time Formulations, Iterative Strategies and Massively Parallel Implementations.  New Methods in Transient Analysis, 246:0, 1992.  [83℄ Truong, Anh H., Oldeld, Chad A., and Zingg, David W. Mesh Movement for a Dis rete-Adjoint Newton Krylov Algorithm for Aerodynami  Optimization.  AIAA Journal, 46, issue 7:16951704, 2008. [84℄ Venkatakrishnan, V. Convergen e to Steady State Solutions of The Euler Equations on Unstru tured Grids with Limiters. 118(1):120130, 1995.  131  Journal of Computational Physi s,  [85℄ Venter, G., and Sobiesz zanski-Sobieski, J. Multidis iplinary Optimization of a Transport Air raft Wing Using Parti le Swarm Optimization.  iplinary Optimization, 26:pp 121131, 2004. [86℄ Watt, A., and Watt, M.  Stru t Multidis-  Advan ed Animation and Rendering Te hniques.  Addison-Wesley Publishing Company, New York, 1992. [87℄ Zhu, C., Byrd, R. H., Lu, P., and No edal, J.  Algorithm 778:  L-BFGS-B:  Fortran Subroutines for Large-S ale Bound-Constrained Optimization.  Transa tions on Mathemati al Software (TOMS), 23(4):550560,  ACM  1997.  [88℄ Zymaris,A.S., Papadimitriou, D.I., Giannakoglou, K.C., and Othmer,C. Feasibility Study of Constant Eddy-Vis osity Assumption in Gradient-Based Design Optimization.  Journal of Computational Physi s, 229:52285245, 2010.  132  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items