Aerodynami Optimization Using High-Order Finite-Volume CFD Simulations by Mohammad Baher Azab B.S ., Aerospa e Engineering Cairo University (Egypt), 1999 M.S ., Aerospa e Engineering Cairo University (Egypt), 2002 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Fa ulty of Graduate Studies (Me hani al Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Van ouver) O tober 2011 © Mohammad Baher Azab 2011 Abstra t The growth of omputer power and storage apa ity allowed engineers to ta kle engineering design as an optimization problem. For transport air raft, drag mini- mization is riti al to in rease range and redu e operating osts. Lift and geometri onstraints are added to the optimization problem to meet payload and rigidity on- straints of the air raft. Higher order methods in CFD simulations have proved to be a valuable tool and are expe ted to repla e urrent se ond order CFD methods in the near future; therefore, exploring the use of higher order CFD methods in aerodynami optimization is of great resear h interest and is one goal of this thesis. Gradient-based optimization te hniques are well known for fast onvergen e, but they are only lo al minimizers; therefore their results depend on the starting point in the design spa e. The gradient-independent optimization te hniques an nd the global minimum of an obje tive fun tion but require vast omputational eort; therefore, for global optimization with reasonable omputational ost, a hybrid op- timization strategy is needed. A new least-squares based geometry parametrization is used to des ribe airfoil shapes and a semi-torsional spring analogy mesh morphing tool updates the grid everywhere when the airfoil geometry hanges during shape optimization. For the gradient based optimization s heme, both se ond and fourth order sim- ulations have been used to ompute the obje tive fun tion; the adjoint approa h, well known for its low omputational ost, has been used for gradient omputa- tion and mat hes well with nite dieren e gradient. The gradient based optimizer have been tested for subsoni and transoni inverse design problems and for drag minimization without and with lift onstraint to validate the developed optimizer. The optimization s heme used is Sequential Quadrati Programming (SQP) with the BFGS approximation of the Hessian matrix. A mesh renement study is presented for an aerodynami ally onstrained drag minimization problem to show how se ond ii and fourth order optimal results behave with mesh renement. A hybrid parti le swarm / BFGS s heme has been developed for use as a global optimizer. It has been tested on a drag minimization problem with a lift onstraint; the hybrid s heme obtained a sho k free proles, while gradient-based optimization ould not in general. iii Prefa e The resear h ideas and methods explored in the three o-authored manus ripts of this thesis are the fruits of a lose working relationship between Dr Carl Ollivier-Goo h and Mohammad Azab. The implementation of the methods, the data analysis, and the manus ript preparation were done by Mohammad Azab with invaluable guidan e from Carl Ollivier-Goo h throughout the pro ess. iv Table of Contents Abstra t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv Nomen lature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Finite Volume Flow Solver . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Numeri al Aerodynami Optimization . . . . . . . . . . . . . . . . 5 1.2.1 Gradient-based aerodynami optimization . . . . . . . . . . 6 1.2.2 Gradient free optimization . . . . . . . . . . . . . . . . . . . 8 1.3 Contributions of the Thesis . . . . . . . . . . . . . . . . . . . . . . . 10 2 GEOMETRY PARAMETRIZATION . . . . . . . . . . . . . . . . . 12 2.1 Survey of Geometri Parametrization Te hniques . . . . . . . . . . . 12 2.1.1 Analyti al parametrization . . . . . . . . . . . . . . . . . . . 13 2.1.2 Pie e-wise spline parametrization . . . . . . . . . . . . . . . 14 2.1.3 CAD parametrization . . . . . . . . . . . . . . . . . . . . . . 14 2.1.4 Free form deformation (FFD) . . . . . . . . . . . . . . . . . 15 2.1.5 Multidis iplinary aerodynami /stru tural shape optimization using deformations (MASSOUD) . . . . . . . . . . . . . . . . 15 2.2 New Least Squares Parametrization Te hnique . . . . . . . . . . . . 16 v 2.2.1 Airfoil geometry parametrization . . . . . . . . . . . . . . . . 16 2.2.2 Thi kness onstraint . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.3 Validation ases . . . . . . . . . . . . . . . . . . . . . . . . . 23 3 MESH MORPHING AND MESH SENSITIVITY . . . . . . . . . 30 3.1 Mesh Morphing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2 Testing Mesh Morphing . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3 Mesh Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4 GRADIENT CALCULATION USING ADJOINT APPROACH 35 4.1 Forward and Adjoint Formulations . . . . . . . . . . . . . . . . . . 35 4.2 Element and Fa e Geometri Properties Dependen y on Mesh Coor- dinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.2.1 Element area mesh dependen y . . . . . . . . . . . . . . . . 45 4.2.2 Fa e length mesh dependen y . . . . . . . . . . . . . . . . . 46 4.2.3 Fa e normal mesh dependen y . . . . . . . . . . . . . . . . . 49 4.3 Fa e Flow Properties Dependen y on The Mesh . . . . . . . . . . . 50 4.3.1 Fa e ow properties re onstru tion . . . . . . . . . . . . . . 50 4.3.2 Mesh dependen e of the fa e ow property re onstru tion . . 57 5 GRADIENT VALIDATION . . . . . . . . . . . . . . . . . . . . . . . 60 5.1 Subsoni Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.2 Transoni Test Case with No Limiter . . . . . . . . . . . . . . . . . 61 5.3 Transoni Test Case with Limiter . . . . . . . . . . . . . . . . . . . 65 5.4 Sensitivity of nite dieren e gradient to perturbation magnitude . 67 6 GRADIENT BASED OPTIMIZATION TEST CASES . . . . . . 70 6.1 Subsoni Inverse Design . . . . . . . . . . . . . . . . . . . . . . . . . 70 6.2 Transoni Inverse Design . . . . . . . . . . . . . . . . . . . . . . . . 71 6.3 Drag Minimization without Lift Constraint . . . . . . . . . . . . . . 74 6.4 Drag Minimization with Lift Constraint . . . . . . . . . . . . . . . . 80 6.5 Mesh Renement Study of a Drag Minimization with Lift Constraint 86 vi 7 PARTICLE SWARM OPTIMIZATION AND A NEW HYBRID OPTIMIZATION METHOD . . . . . . . . . . . . . . . . . . . . . . 88 7.1 Introdu tion to Gradient Free Optimization . . . . . . . . . . . . . . 88 7.2 Swarm Intelligen e . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 7.2.1 Premature onvergen e dete tion . . . . . . . . . . . . . . . . 93 7.2.2 Redening de ision variables range . . . . . . . . . . . . . . . 94 7.2.3 Regrouping and position lipping . . . . . . . . . . . . . . . 94 7.3 Hybrid SQP-RegPSO Te hnique . . . . . . . . . . . . . . . . . . . . 95 7.4 Optimization Test Cases . . . . . . . . . . . . . . . . . . . . . . . . 95 7.4.1 Constrained drag optimization of the NACA 0012 . . . . . . 97 7.4.2 Constrained drag minimization of NACA 00083 . . . . . . . 102 7.4.3 Aerodynami and thi kness onstraint drag minimization of NACA 0012 . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.4.4 Aerodynami and thi kness onstraint drag minimization of NACA 00083 . . . . . . . . . . . . . . . . . . . . . . . . . . 108 7.5 Comparing SQP-RegPSO with RegPSO-SQP optimization strategies 112 8 CONCLUSIONS AND RECOMMENDATIONS . . . . . . . . . . 119 8.1 Contributions and Con lusions . . . . . . . . . . . . . . . . . . . . . 119 8.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 vii List of Tables 2.1 RMS error in dierent parametrized airfoil geometries . . . . . . . . 29 4.1 Two point Gauss quadrature rule . . . . . . . . . . . . . . . . . . . . 48 5.1 The magnitude of se ond and fourth order CL gradients and angles between the evaluated gradients for NACA 0012 in subsoni ow. . . 62 5.2 The magnitude of se ond and fourth order CD gradients and angles between the evaluated gradients for NACA 0012 in an unlimited tran- soni ow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.3 Magnitudes of se ond and fourth order CD gradients and angles be- tween the evaluated gradients for NACA 0012 in Venkatakrishnan limited transoni ow . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4 The magnitude of se ond and fourth order CD gradients and angles between the evaluated gradients for NACA 0012 using higher order limiter in transoni ow. . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.5 The magnitudes and angles between the evaluated modied fourth order CD gradients using adjoint, and nite dieren e for NACA 0012 using Venkatakrishnan and higher order limiters in transoni ow. . 67 5.6 Finite dieren e drag gradient sensitivity with respe t to perturbation amplitude ǫ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.1 Aerodynami oe ients of original and optimized RAE 2822 airfoil at transoni onditions. . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.2 Aerodynami oe ients of original and optimized RAE 2822 airfoil at transoni onditions. . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.3 Lift penalty weight ee t on Drag minimization of RAE 2822 with lift onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 viii 7.1 Angles between adjoint, FD., and SQP-hybrid ve tors . . . . . . . . 101 7.2 Angles between adjoint, FD, and BFGS-hybrid best solutions ve tors 103 7.3 Thi kness penalty weight impa t on the optimization results of NACA 00083 using hybrid s heme . . . . . . . . . . . . . . . . . . . . . . . . 110 7.4 Optimization results of the hybrid s heme in Phase I, II . . . . . . . 113 7.5 Optimization results of the hybrid T s heme in Phase I, II . . . . . . 113 ix List of Figures 1.1 Airfoil aerodynami optimization y le using gradient-based optimiza- tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Least square surfa e presentation of RAE 2822 airfoil using two poly- nomials P1 (x) & P2 (x) tted using nine ontrol points. . . . . . . . 18 2.2 Parametrized airfoil using four polynomials . . . . . . . . . . . . . . 21 2.3 Parametrized NACA 0011SC, 'O' are the original airfoil ordinates . . 24 2.4 Parametrized NACA 0012, 'O' are the original airfoil ordinates . . . 24 2.5 Parametrized NACA 6509, 'O' are the original airfoil ordinates . . . 25 2.6 Parametrized NACA 16006, 'O' are the original airfoil ordinates . . . 25 2.7 Parametrized NACA 63412, 'O' are the original airfoil ordinates . . . 26 2.8 Parametrized NACA 644421, 'O' are the original airfoil ordinates . . 26 2.9 Parametrized RAE2822, 'O' are the original airfoil ordinates . . . . . 27 2.10 Parametrized LV2 laminar airfoil, parametrized using 5 polynomials per surfa e, 'O' are the original airfoil ordinates . . . . . . . . . . . . 28 2.11 LV2 airfoil parametrized using dierent methods with 20 design vari- ables, presented with permission of Brezillon. . . . . . . . . . . . . . 28 3.1 S hemati drawing of an edge pq and its fa ing angles a and b. . . . 31 3.2 Mesh movement s heme results of doubling the thi kness of NACA 0012 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Mesh movement results, large outer boundary deformation of a re t- angular domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.1 S hemati drawing of an element fa e κ, and illustration of its left and right sides . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2 General triangular element with unit normal n̂ on one of its fa es κ. 44 x 4.3 Boundary element with urved fa e for high order integration s heme 47 4.4 S hemati drawing of a urved fa e with two Gauss points used . . . 48 4.5 S hemati drawing of rst, se ond and third neighbor layers of trian- gular element k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.1 The pressure sensitivity with respe t to one of the design ontrol points omputed for subsoni ow over NACA 0012, omparing the sensitivity al ulation of Eq. 4.6 with nite dieren e results. . . . . 61 5.2 CL gradient error in se ond and fourth order s hemes with respe t to the design points normalized by gradient magnitude. . . . . . . . . . 62 5.3 The pressure sensitivity with respe t to one of the design ontrol points omputed for an unlimited transoni ow around NACA 0012, omparing the sensitivity al ulation of Eq. 4.6 with nite dieren e results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.4 CD gradient error in se ond and fourth order s hemes with respe t to the design points, normalized by gradient magnitude. . . . . . . . . . 64 5.5 The pressure sensitivity with respe t to one of the design ontrol points omputed for an unlimited transoni ow around NACA 0012. 65 5.6 The normalized CD gradient error in se ond, fourth, and modied fourth order s hemes with respe t to the design points in a limited transoni ow (Venkatakrishnan limiter) . . . . . . . . . . . . . . . . 68 5.7 The normalized CD gradient error in se ond, fourth, and modied fourth order s hemes with respe t to the design points in a limited transoni ow (higher order limiter). . . . . . . . . . . . . . . . . . . 68 6.1 Subsoni NACA 2412 inverse design pressure distributions for the ini- tial, target, and optimized airfoil proles . . . . . . . . . . . . . . . . 72 6.2 Se ond and fourth order optimization onvergen e history. . . . . . . 72 6.3 Gradient norm history . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.4 Subsoni inverse design optimal airfoil shapes . . . . . . . . . . . . . 73 6.5 The dieren e between the target prole and the optimized proles, se ond and fourth order . . . . . . . . . . . . . . . . . . . . . . . . . 74 xi 6.6 Subsoni NACA 2412 inverse design pressure distributions for the ini- tial, target, and optimized airfoil proles . . . . . . . . . . . . . . . . 75 6.7 Se ond and fourth order optimization onvergen e history. . . . . . . 75 6.8 Gradient norm history . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.9 The dieren e between the target prole and the optimized proles, se ond and fourth order . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.10 Transoni inverse design optimal airfoil shapes. . . . . . . . . . . . . 77 6.11 Pressure ontours of RAE 2822 at Ma h 0.73 and angle of atta k 2 . 78 6.12 Optimized RAE 2822. . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6.13 Optimization surfa e displa ements of the original RAE2822 surfa es. 79 6.14 Pressure distribution omparison for the original and the optimized geometries. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 6.15 Obje tive fun tion response surfa e along positive and negative gra- dient dire tions entered at the optimal prole of RAE 2822 transoni drag minimization without lift onstraint. . . . . . . . . . . . . . . . 80 6.16 Se ond and fourth order optimization onvergen e history. . . . . . . 81 6.17 Surfa e pressure distribution, of the initial and optimized geometries 82 6.18 Se ond and fourth order optimization onvergen e history. . . . . . . 82 6.19 Optimal shapes omparison: se ond order, fourth order, and opti- mized prole by Brezillon and Gauger ompared with the original RAE2822. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.20 Se ond and fourth order optimized prole surfa e displa ements from the initial shape. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.21 Obje tive fun tion response surfa e along positive and negative gra- dient dire tions entered at the optimal prole of RAE 2822 transoni drag minimization without lift onstraint. . . . . . . . . . . . . . . . 84 6.22 initial and optimal pressure distribution obtained by Brezillon and Gauger [11℄ (presented with permission) . . . . . . . . . . . . . . . . 85 6.23 Optimal CD value with mesh renement . . . . . . . . . . . . . . . . 87 7.1 Geneti Algorithm optimization . . . . . . . . . . . . . . . . . . . . . 90 7.2 Velo ity omponents and position update of a parti le . . . . . . . . 93 7.3 Pressure Field of the original NACA 0012 and the optimized proles 98 xii 7.4 Surfa e pressure distribution of original NACA 0012 and the opti- mized proles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 7.5 Airfoil shapes of original NACA 0012 and the optimized proles . . . 99 7.6 Fun tion minimization iterations in the se ond phase of the hybrid s heme (RegPSO) with a NACA 0012 airfoil as the starting point. . 100 7.7 Obje tive fun tion values between SQP and hybrid optimal points (starting geometry: NACA 0012). . . . . . . . . . . . . . . . . . . . . 101 7.8 Pressure surfa e of NACA 00083 and the optimized proles . . . . . 103 7.9 Surfa e pressure of NACA 00083 and the optimized proles . . . . . 104 7.10 Prole omparison of NACA 00083, BFGS optimal, and hybrid opti- mal airfoils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.11 Fun tion minimization iterations in the se ond phase of the hybrid s heme (RegPSO) of NACA 00083 start . . . . . . . . . . . . . . . . 105 7.12 Obje tive fun tion values between BFGS and hybrid optimal points, starting geometry NACA 00083 . . . . . . . . . . . . . . . . . . . . . 105 7.13 Pressure surfa e of NACA 0012 and the optimized proles with thi k- ness onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.14 Surfa e pressure of NACA 0012 and the optimized prole with thi k- ness onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.15 Pressure surfa e of NACA 00083 and the optimized proles with thi k- ness onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 7.16 Surfa e pressure of NACA 00083 and the optimized prole with thi k- ness onstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7.17 Prole omparison of the hybrid s heme optimal proles, starting ge- ometries NACA 0012 and NACA 00083 . . . . . . . . . . . . . . . . 110 7.18 Thi kness penalty weight impa t on the Optimal pressure distribution of NACA 00083 using hybrid optimization s heme. . . . . . . . . . . 111 7.19 Thi kness penalty weight impa t on the Optimal prole with NACA 00083 starting geometry using hybrid optimization s heme. . . . . . 111 7.20 Hybrid and hybrid T optimal pressure distribution for NACA 0012 (Case I). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 7.21 Hybrid and hybrid T optimal pressure distribution for NACA 00083 (Case II). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 xiii 7.22 Hybrid and hybrid T optimal pressure distribution for NACA 0012 with 10% thi kness onstraint and unit thi kness penalty weight (Case III). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 7.23 Hybrid and hybrid T optimal pressure distribution for NACA 00083 with 10% thi kness onstraint and unit thi kness penalty weight (Case IV). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.24 Hybrid and hybrid T optimal pressure distribution for NACA 00083 with 10% thi kness onstraint and ten thi kness penalty weight (Case V). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 xiv List of Algorithms 7.1 Parti le swarm pseudo ode . . . . . . . . . . . . . . . . . . . . . . . 91 7.2 Hybrid SQP-RegPSO pseudo ode . . . . . . . . . . . . . . . . . . . 96 7.3 Pseudo ode of RegPSO-SQP hybrid T optimization s heme . . . . 112 xv Nomen lature α Angle of atta k △Q Change in onservative ow properties △t Time step δ norm Normalized swarm radius ∂R ∂Q ja obian matrix ∂U ∂M Dependen y of premitive ow variables solution on mesh. dCL,D,M dD Gradient of lift, drag. and moment oe ient dF dD Obje tive fun tion gradient dR dD Residual dependen y of design variable γ Spe i heats ratio n̂ Fa e unit normal Ωi Control volume i Uk Control volume avaraged primitive ow properties −→ F Flux ve tor a ross ontrol volume boundary −→g Swarm best known position in the design spa e −→pi Parti le i best known position in the design spa e −→rb Position ve tor of point b xvi −→ Vi Parti le i velo ity −→xi(k) Parti le position in iteration k −→xi(k + 1) Parti le position in iteration k+1 ∂M/∂D Mesh sensitivity ∂R/∂D Residual senstivity with respe t to design variables ∂−→rb/∂Di Dependen y of boundary point position −→rbve tor on design variable D ΨL,D,M Adjoint ve tor of of lift, drag. and moment oe ient ρ Density ρ̃ Roe avaraged density Ã Roe averaged ux ja obian matrix ã Roe avaraged speed of sound h̃ Roe avaraged enthalpy ũ Roe avaraged x-velo ity ṽ Roe avaraged y-velo ity ~F (QR) , ~F (QL) Right and left fa e uxes CD Drag oe ient CL Lift oe ient CM Moment oe ient Cp, Cv Spe i heats at onstant pressure and volume Di Desigen variable i et Spe i total energy Kib Stiness matrix intries relates interior nodes with boundary points xvii Kii Stiness matrix intries relates interior points PRec Flow property re onstru tion polynomial Ps Surfa e pressure Q Conservative ow variables ve tor QR, QL onserved ow properties at right and left of the ell fa e Ubg Boundary Gauss point ow premitive properties Ufi Re onstru ted primitive ow properties at fa e Gaussian integration point i Ui, Ub Interior and boundary mesh point displa ements respe tively Un Velo ity in the dire tion normal to the ontrol volume boundary ws Gaussian integration weight yc Geometry ontrol point variable y ordinate (design veriable) a Speed of sound E Energy per unit volume e Spe i internal energy F Aerodynami obje tive fun tion GMRES Generalized minimal residual h Spe i enthalpy M Ma h number P Pressure R Residual RMS Root mean square T Temperature xviii A knowledgments I would like to express my deepest gratitude to my supervisor Dr Carl Ollivier- Goo h. You are su h a great resear h supervisor and gentleman. I would like also to a knowledge my supervision ommittee, Dr Mi hael Friedlander, Dr Philip Loewen, and Dr James Olson for their valuable remarks that helped me to justify my resear h strategy. I would like also to thank my fellows in the resear h and development department of the Egyptian Air For e espe ially Col. Dr Eng. Mohamed Ibrahim Mustafa and Major General Dr Abdalla El-Ramsisy for their great support during my resear h in Egypt whi h has had a great impa t on my PhD resear h progress. Great thanks goes also to my dear professor Dr Chen Grief for being there when I needed his advi e. Finally, I would like to thank my wife and the light of my life, Deena, for her support and bringing happiness to my life. Mohammad Baher Azab xix Dedi ation To my home ountry, the land of Pharaohs Egypt To the land of ho key Canada To the souls of my parents Baher and Laila xx Chapter 1 INTRODUCTION Drag redu tion of transport air raft is of great importan e be ause it redu es air raft fuel onsumption and thereby redu es operating osts and environmental impa t in the form of pollution and global warming. Drag, lift and other aerodynami for es an be predi ted using CFD simulations, whi h have be ome an essential tool for aerodynami analysis and design. CFD simulations arried out using unstru tured grids give a urate aerodynami for e predi tions, and unstru tured grids have the advantage of easily representing any omplex shape. As most transport air raft travel at transoni speed, it is of great importan e to redu e the wave drag by weakening or entirely eliminating sho k waves on the wing. Aerospa e engineers use CFD simulations with numeri al optimization te hniques for aerodynami design; the optimization problem is to minimize an aerodynami obje tive (often drag) by hanging geometri design variables, given an initial aerodynami shape and subje t to some geometri and aerodynami onstraints. This pro ess requires a urate assessment of the aerodynami hara teristi s of a given geometry. The ow model used in this resear h is Euler's ow model; although this model negle ts vis ous ow ee ts, redu ing the drag using this invis id ow model will ultimately redu e the drag on real ongurations as the invis id drag is about 40% of the total drag [47℄. A higher order nite volume CFD solver developed by C. Ollivier-Goo h and o- workers [68, 61, 56℄ is used in this resear h. Optimization s hemes an be divided into two main ategories: gradient-based and non-gradient-based s hemes. Gradient-based s hemes require omputing the ob- je tive fun tion value and its gradient with respe t to the design variables, while the non-gradient-based te hniques require only omputing the obje tive fun tion value. Gradient-based s hemes are fast to onverge to a lo al optimal point in the design spa e ompared to the non-gradient-based methods, but the optimal solution found by gradient-based optimization depends on the starting point [33℄. Non-gradient- 1 based s hemes an nd global minimum solution regardless of the starting point, but with larger omputational ost ompared to gradient-based s heme. 1.1 Finite Volume Flow Solver The two dimensional integral form of Euler's equation an be written for a ontrol volume Ωi as ∂ ∂t ¨ Ωi QdV + ˛ dΩi −→ F · n̂dl = 0 (1.1) where n̂ = nxi + nyj is the outward pointing normal to the ontrol volume fa es; Q and −→ F are respe tively the onserved variable ve tor and the ux a ross the boundary of ontrol volume Ωi boundaries. These an be expressed as Q = ρ ρu ρv E , −→F · n̂ = ρUn ρUnu+ nxp ρUnv + nyp (E + p)Un p = (γ − 1) [ E − ρ ( u2 + v2 ) 2 ] Un = nxu+ nyv where is the uid density, u and v are the Cartesian velo ity omponents, p is the pressure, and E is the energy per unit volume. Un is the velo ity in the dire tion normal to the ontrol volume boundary. Other thermodynami relations like the spe i heats at onstant pressure (Cp) and onstant volume (Cv), the ratio of spe i heats γ an be expressed in terms of the gas onstant for air, R: γ = Cp Cν , Cp = γR γ − 1 , Cν = R γ − 1 2 For a thermally and alori ally perfe t gas, thermodynami properties an be related by P = ρRT e = CvT et = e+ 1 2 ( u2 + v2 ) E = ρet h = e+ P ρ a = √ γRT where T is the temperature, e is the spe i internal energy, et is the spe i total energy, h is the spe i enthalpy, and a is the speed of sound. For dis rete solution of the Euler equations, ow properties are normalized by some referen e values in order to redu e the round-o errors in the dis retized linear system; for external aerodynami s, the normalized ow properties are as follows ρ̄ = ρ ρ∞ , ā = a a∞ , P̄ = P ρ∞a2∞ ū = u a∞ , v̄ = v a∞ , Ē = P̄ γ − 1 + 1 2 ρ̄ ( ū2 + v̄2 ) x̄ = x L , ȳ = y L , t̄ = ta∞ L With this normalization, the normalized Euler equations are identi al to their di- mensional form with the addition of [̄·] to every variable. From this point on, the normalized ow properties are used and therefore [̄·] will be dropped. The ux is evaluated using Roe ux dieren e splitting te hnique [74℄ and eval- uated at ea h ell fa e κ using the following formula ~F = 1 2 [ ~F (QR) + ~F (QL)− ∣∣∣Ã∣∣∣ (QR −QL)] κ (1.2) where QR, QL are the onserved ow properties at right and left of the ell fa e κ and the Roe averaged matrix Ã is the ux Ja obian ∂ ~F/∂Q at Roe-averaged quantities 3 as follows ρ̃ = √ ρLρR ũ = uL + uR √ ρR ρL 1 + √ ρR ρL , ṽ = vL + vR √ ρR ρL 1 + √ ρR ρL h̃ = hL + hR √ ρR ρL 1 + √ ρR ρL , ã2 = (γ − 1) ( h̃− ( ũ2 + ṽ2 ) 2 ) (1.3) The Roe-averaged Ja obian matrix Ã has four eigenvalues. Three of these are dis- tin t, but the eigensystem is omplete. The Roe dissipation matrix an be written in terms the eigenvalues and eigenve tors of the Ja obian as ∣∣∣Ã∣∣∣ = X̃ ∣∣∣Ũn∣∣∣ 0 0 0 0 ∣∣∣Ũn∣∣∣ 0 0 0 0 ∣∣∣Ũn + ã∣∣∣ 0 0 0 0 ∣∣∣Ũn − ã∣∣∣ X̃ −1 (1.4) where X̃ is the eigenve tor matrix evaluated using Roe-averaged ow quantities; the term ∣∣∣Ã∣∣∣ (QR −QL) an be written a ording to Frink as follows [27, 26℄∣∣∣Ã∣∣∣ (QR −QL) = ∣∣∣Ã∣∣∣△Q = ∣∣∣△F̃1∣∣∣+ ∣∣∣△F̃4∣∣∣+ ∣∣∣△F̃5∣∣∣ (1.5) 4 where ∣∣∣△F̃1∣∣∣ = ∣∣∣Ũn∣∣∣ ( △ρ− △P ã2 ) 1 ũ ṽ ũ2+ṽ2 2 + ρ̃ 0 △u− nx△U △v − ny△U ũ△u+ ṽ△v − Ũn△U ∣∣∣△F̃4,5∣∣∣ = ∣∣∣Ũn ± ã∣∣∣ (△P ± ρ̃ã△U 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã △ρ = ρR − ρL, △u = uR − uL, △v = vR − vL, △P = PR − PL, △U = nx△u+ ny△v Higher-order a ura y is obtained by least-squares re onstru tion of the non- onserved variables U = [ ρ u v p ]T and Gauss quadrature in ux integration [68, 61℄. After integrating uxes around ea h ontrol volume in the mesh, an impli it time dis retization leads to a sparse system of linear equations, whi h for the simplest ase of a global time step an be written as,[ I △t + ∂R ∂Q ]n △Qn = Rn (1.6) where △Q = [ △ρ △ρu △ρv △E ]T , ∂R ∂Q is the global Ja obian matrix, and R is the residual. The steady state solution is obtained iteratively when △Q → 0. In pra ti e, we use a quasi-Newton generalization of Eq. 1.6 that in ludes residual-based lo al time stepping [57℄ and solve the system using GMRES [75℄. 1.2 Numeri al Aerodynami Optimization Aerodynami design used to rely on CFD simulations in onjun tion with experimen- tal testing and engineering intuition of the designer. With the growth of high speed omputers, integrating numeri al optimization s hemes with CFD simulations has 5 be ome possible and is now used for aerodynami design and optimization. Gradient- based optimization te hniques are widely used be ause they rea h an optimized shape after a reasonable exe ution time; however, the nal optimal shape is the lo al min- imum lo ated downhill from the optimization starting point. Non-gradient-based methods like geneti algorithm (GA) or parti le swarm (PS) are slower to nd an optimum but an nd the global minimum regardless of the starting point; their drawba k is the large number of iterations required to rea h this minimum om- pared to gradient-based s hemes. 1.2.1 Gradient-based aerodynami optimization Gradient-based optimization depends on evaluating the gradient of the obje tive fun tion with respe t to the design variables and using the gradient in a linear model (steepest des ent) or a quadrati model (Newton or quasi-Newton model) to nd a sear h dire tion; this sear h dire tion is the dire tion in whi h the de- sign variables should hange their values to minimize the obje tive fun tion [66℄. Gradient-based te hniques have been widely used in aerodynami optimization due their fast onvergen e to an optimal solution. The obtained optimal shape is bi- ased by the optimization starting point (initial aerodynami shape), and there is no guarantee that gradient-based methods an nd the global best optimal shape in the design spa e [2℄. Hi ks and Henne where among the rst to apply gradient- based optimization te hniques in aerodynami design in the late 1970's [35℄. Sin e then many resear hers have investigated the use of gradient-based optimization te h- niques like steepest des ent and quadrati programming in aerodynami optimiza- tion [17, 31, 28℄. The most expensive part of gradient-based optimization is the gradient al ulation. Hi ks and Henne used nite dieren e rules to al ulate the obje tive fun tion gradient. This means that two CFD simulations were required for ea h design variable to ompute the gradient (using a entral nite dieren e rule), whi h is omputationally expensive. The same strategy was applied by Consentino and Holst to optimize transoni wings [17℄. The use of the adjoint method, whi h was originally applied to aerodynami s problems by Jameson, to ompute the gra- dient redu ed the omputational ost of gradient al ulation to the ost of one ow simulation regardless of the number of design variables [41, 43, 60, 38, 39, 42, 11, 62℄. 6 Figure 1.1: Airfoil aerodynami optimization y le using gradient-based optimization Numerous resear hers have applied gradient-based optimization using the ad- joint approa h in aerodynami optimization sin e the early 1990's. Reuther and his o-workers applied the adjoint approa h for aerodynami optimization of air raft onguration using Euler's ow model [73℄. Jameson developed an adjoint formula- tion for the Navier-Stokes equations and applied it to transoni wing optimization [43℄; in this work, Jameson suggested that an optimal pressure distribution rst be obtained using Euler's ow model then used as a target pressure for an inverse design optimization problem using a Navier-Stokes adjoint optimizer to redu e the overall omputational ost. However, Jameson treated the eddy vis osity as onstant whi h was later shown to be a bad assumption. Anderson and Bonhaus examined the ee t of the strength of oupling of the turbulen e model to the ow equations. They om- pared the adjoint gradient of the ow equations, with the eddy vis osity frozen, with the nite dieren e gradient of the ombined ow and turbulen e equations, and found that freezing the eddy vis osity an lead to signi ant error in the omputed gradient. Therefore, they developed a ow solver that ouples the Spalart-Allmaras one-equation turbulen e model with the ow equations; their oupled solver used 5× 5 blo ks for two dimensional ows and 6× 6 blo ks for three dimensional ases. They found that this oupling of the turbulen e equations with the ow equations led 7 to an a urate adjoint gradient [3℄. The same observation was veried by Nielsen and Kleb who extended their adjoint solver to deal with hemi ally rea ting ows [65℄. Zymaris et al developed a ontinuous adjoint optimizer for turbulent ow using the k − ε turbulen e model and applied it to du t optimization; they showed also that the assumption of onstant eddy vis osity leads to great ina ura y in the omputed gradient and this leads to a poor sear h dire tion [88℄. Regarding the optimization te hnique used, early resear hers used the steepest des ent s heme, in whi h the design variables are updated on a sear h dire tion exa tly opposite to the gradient of the obje tive fun tion [66, 4℄. This s heme has been implemented by Jameson and other resear hers [37, 41, 73, 3℄, but as the steepest des ent s heme requires a large number of iterations to onverge to minimal point, the sequential quadrati programming (SQP) s heme seems to be an attra tive andidate as an optimization te hnique. The use of SQP requires omputing the Hessian of the obje tive fun tion with respe t to the design variables. The exa t Hessian is expensive to ompute and may not be positive denite; therefore, the BroydenFlet herGoldfarbShanno (BFGS) approximate formula is often adopted to approximate the Hessian using obje tive fun tion gradient history [12, 13, 87℄. The BFGS approximation always gives a positive denite approximate Hessian, therefore a real optimization sear h dire tion is guaranteed. Dadone et. al used BFGS in aerodynami optimization for transoni and supersoni wings and ompared BFGS and steepest des ent; there results showed that the BFGS method is more e ient in nding the optimal solution and is less sensitive to any ina ura y aused by approximation in gradient omputations [19℄. BFGS optimization has also been used by Neme and Zingg in subsoni and transoni turbulent aerodynami optimization of two dimensional airfoil [63℄. 1.2.2 Gradient free optimization Gradient free or gradient independent optimization methods, also known as heuris- ti optimization methods, are optimization te hniques that do not require obje tive fun tion gradient omputation and therefore an be applied to non-dierentiable problems. They an be ategorized as evolutionary s hemes (in luding geneti algo- rithms) and random sear h s hemes (like the parti le swarm te hnique). 8 Geneti algorithm optimization is an evolutionary optimization algorithm that is inspired by Darwin's theory of evolution and natural sele tion [30℄. The design variables are treated as hromosomes and optimization is arried out by rossing and mutating these hromosome to nd a better solution that minimizes the obje tive fun tion. The initial population is randomly generated, the obje tive fun tion value is omputed using CFD for ea h population member and a tness value is omputed based on that; some members are sele ted based on their tness to be the parents of the next generation and are used to generate the hromosomes of the ospring of the next optimization iteration. Transoni wing optimization using a geneti algorithm was explored by Gregg and Misegades [32℄ and by Gage and Kroo [29℄ in the late 1980's to minimize drag with lift onstraint. Somewhat later, Anderson applied a geneti algorithm in subsoni wing optimization with stru tural onstraints [1℄; he added the geometri and aerodynami onstraints to the obje tive fun tion as penalty terms. Jang and Lee applied a geneti algorithm in subsoni and transoni invis id airfoil optimization; their obje tive was to maximize the lift-to-drag ratio of an airfoil, beginning from the NACA 0012 1 [44℄. Oyama et al applied a geneti algorithm with a Navier-Stokes solver for transoni wing optimization [70℄. They also explored the use of fra tal analysis in GA aerodynami optimization [71℄. The parti le swarm method (PSOpt) is a sto hasti optimization sear h method developed by Eberhart and Kennedy in 1995, inspired by the so ial behavior of bird o ks [46, 21℄. The general idea of the parti le swarm optimization is to randomly generate a swarm of parti les in the design spa e. For ea h parti le a tness value is al ulated based on CFD simulation. Then parti les y in the design spa e a ording to a simple formula that takes into a ount that parti le's own best t- ness position and the swarm's overall best tness position [22℄. PSOpt is known to suer from premature onvergen e prior to dis overing the true global minimizer; Evers suggests an automati regrouping PSO (RegPSO) that automati ally triggers swarm regrouping when premature onvergen e is dete ted. The suggested regroup- ing strategy aims to liberate parti les from sub-optimal solutions and enables nding the global minimum [24℄. Although the PSO algorithm has been applied to a wide 1 The National Advisory Committee for Aeronauti s (NACA) airfoil family geometri oordinates an be found at University of Illinois Urbana-Champaign website http://www.ae.illinois.edu/m- selig/ads/ oord_database.html. An explanation of the meaning of the digits in the NACA airfoil naming s heme an be found in [36℄. 9 range of engineering problems in the literature, very few aerodynami optimization appli ations are known. Venter and Sobiesz zanski applied the parti le swarm op- timization te hnique in the multidis iplinary optimization of a wing; the obje tive was to maximize the air raft range by maximizing lift-to-drag ratio and redu ing the wing weight subje t to geometri onstraints [85℄. Chandrashekarappa and Du- vigneau [15℄ used parti le swarm optimization s heme to aerodynami ally optimize wings in supersoni onditions; Duvigneau also applied parti le swarm optimization to aerodynami optimization of wings at transoni speed with free stream Ma h number un ertainty [20℄. 1.3 Contributions of the Thesis High order CFD methods an ompute an a urate value of an aerodynami obje tive fun tion at lower omputational ost than required when a se ond order method is used. The rst major ontribution of this thesis is a study of whether the e ien y of high order methods for CFD analysis translates into improved e ien y for gradient- based aerodynami optimization. Gradient based optimization te hniques are known to be lo al minimizers, with results depending on the starting airfoil geometry. The se ond major ontribution of this thesis is the development and study of an optimization s heme that an rea h a true global optimum (for invis id transoni aerodynami s, a sho k free airfoil subje t to aerodynami and geometri onstraints) after a reasonable number of CFD simulations. This s heme is a hybrid (BFGS + regrouped parti le swarm) s heme that takes the advantages of gradient-based and gradient-free optimization s hemes. To a hieve these goals, the following omponents were needed, in addition to the pre-existing high-order a urate ow solver: An e ient geometry parametrization method. For optimization purposes, the airfoil shape must be represented by a nite number of design variables. Chap- ter 2 des ribes the requirements on su h a parametrization and presents a new least-squares spline parametrization method developed for this work. Robust mesh movement. As the airfoil shape hanges during optimization, the omputational mesh must be updated. Be ause mesh regeneration will intro- 10 du e una eptably large hanges in the dis retization error, mesh movement is strongly preferred. Chapter 3 des ribes the semi-torsional spring mesh move- ment s heme used in this thesis. Mesh sensitivity al ulation. Cal ulating the gradient of the obje tive fun tion a urately and e iently requires information about the movement of the mesh with hanges in design variables. Se tion 3.3 des ribes this pro ess, whi h om- bines aspe ts of the geometri parametrization and mesh movement s hemes. Obje tive fun tion gradient al ulation. E ient gradient al ulation requires the solution of the adjoint to the governing ow equations. Chapter 4 des ribes how this is done. A key innovation of the thesis is e ient solution of the dis rete adjoint problem for a high-order a urate ow solution s heme. In addition to showing the formulation of the gradient al ulation, this hapter ompares nite dieren e and adjoint gradients for subsoni and transoni ows for se ond and fourth order nite volume s hemes. For transoni ow, results are presented both with and without limiting of the omputational solution to prevent overshoots. Optimization drivers. The gradient-based optimizer using in this thesis is the BFGS-based quasi-Newton solver in Matlab's optimization toolbox. Gradient based optimization test ases for inverse design problems and for drag mini- mization with and without a lift onstraint is shown in Chapter 6. Also, the impa t of mesh renement on the se ond and fourth order optimal results is studied. The regrouped parti le swarm optimization ode was written by the author. This s heme and its hybridization with the BFGS s heme is presented in Chap- ter 7. Examples are given to show the ee tiveness of the hybrid s heme in nding global optima when the gradient-based s heme is unable to. 11 Chapter 2 GEOMETRY PARAMETRIZATION The geometry of engineered obje ts is dened mathemati ally in omputer-aided design (CAD) software, then exported as a group of points or polygons whi h ap- proximates the original geometry. This dis rete data provides input to mesh gener- ation software that reates a dis rete representation of the omputational domain (a mesh), whi h in turn is used as input for CFD analysis. Aerodynami design and optimization modies the aerodynami shape by hanging a set of geometri design variables. An obvious hoi e is to use all the surfa e grid points of a wing, but this approa h auses two problems. First, it makes the design spa e very large and this may lead to a highly expensive optimization. Se ond, it may lead to a non-smooth geometry due to the displa ement independen e of a surfa e mesh point; this prob- lem an be solved by the use of a smoothing fun tion, as des ribed by Jameson [41℄. To avoid these problems, most optimization approa hes rely on some form of geo- metri parametrization. Geometri parametrization te hniques an be lassied as analyti al; pie e-wise spline tting; CAD based; and free form deformation (FFD) approa hes. In Se tion 2.1, a review of various geometry parametrization te hniques is presented. Se tion 2.2 des ribes in detail the geometry parametrization te hnique used in this resear h, a novel pie e-wise least-squares tting te hnique [5, 6℄. 2.1 Survey of Geometri Parametrization Te hniques In this review, various te hniques applied to aerodynami optimization are presented; the basi s of these te hniques and their appli ation limitations are dis ussed, in lud- ing referen e to the original papers des ribing them in more detail . 12 2.1.1 Analyti al parametrization The analyti al approa h was rst applied by Hi ks and Henne for airfoil optimiza- tion [35℄. They suggest that weighted sinusoidal displa ement shape fun tions be added to the base geometry to modify the airfoil shape; the weights are the opti- mization design variables. The sinusoidal displa ement shape fun tion is expressed as h(x) = ( sin ( πx ln(0.5) ln(b) ))a 0 < x < 1 (2.1) where a and b are onstants to ontrol the peak lo ation and the width of the sinusoidal displa ement shape fun tion. a = 4 is re ommended for most ases, while b must be between 0 and 1 [45℄. The Class-Shape fun tion-Transformation method (CST) is another analyti al geometry parametrization te hnique presented by Kulfan [49, 48℄. CST parametrizes the airfoil geometry using the following formula y = x0.5(1− x)S (x) + x△ZTE (2.2) where S (x) is the shape fun tion and △ZTE is the nite thi kness of the trailing edge. Kulfan and Bussoletti re ommended using a weighted Bernstein binomial of order n as a shape fun tion S (x) = n∑ i=0 bi n! i! (n− i)!x i (1− x)n−i (2.3) where the weights bi are used as the design variables. Although CST gives non-wavy proles, it is not apable of representing omplex geometries. Mousavi and Nadiragah ompared the impa t of using dierent geometry parametrization te hniques on the optimal wing shape; using CST parametrization gave drag oe ients higher by about 15% than the optimal geometry from a B-spline parametrization [58℄ for a three-dimensional lift- onstrained transoni drag minimization problem. 13 2.1.2 Pie e-wise spline parametrization Bezier urves an be used for airfoil shape parametrization. Obayashi used Bezier urve for aerodynami optimization using geneti algorithm [67℄; he noti ed that the Bezier urve representation fails to represent geometries that gives rooftop pres- sure distributions 2 be ause Bezier urves are always onvex. B-splines, a general- ization of Bezier urves, were found to be more suitable. A ubi B-spline rep- resentation is a very good geometry parametrization te hnique. Better ontrol of the ubi spline representation an be obtained by in reasing the number of splines that represent the airfoil; Li et al. optimized the NACA 0012 at single point and multi-point operating onditions with a lift onstraint using spline representation as geometry parametrization te hnique [50℄. CAD systems typi ally use non-uniform rational B-spline (NURBS) representation for geometry modeling, allowing them to represent any omplex geometry; a detailed dis ussion of NURBS an be found in The NURBS Book [72℄. Mengistu and Ghaly applied su essfully a NURBS parametrization s heme to turboma hine blade aerodynami optimization using a gradient-free method [53℄. While pie e-wise spline parametrization is well-suited for two-dimensional shapes and simple three-dimensional geometries, omplex shapes require a large number of ontrol points, redu ing the ee tiveness of gradient-based optimization te hniques [76℄. 2.1.3 CAD parametrization Computer aided design pa kages have evolved to implement NURBS for geometry representation due the ex ellent properties of NURBS. Linking the CAD and grid generation software an be done using an API that allows a ess to the CAD system's internal interfa e [81, 8℄. However, imposing geometri onstraints is still an obsta- le. Mesh sensitivity al ulation is another obsta le for gradient-based optimization te hniques based of CAD parametrization: analyti al mesh sensitivity with respe t to geometry design variables the NURBS knots an in prin iple be omputed with the use of automati dierentiation of the CAD software but this is not possible without CAD sour e ode and is unlikely to be pra ti al even then, given the size of 2 Distributions for whi h pressure remains almost un hanged over a signi ant hordwise dis- tan e. 14 the CAD ode base. This derivative an also be omputed using nite dieren es, but the risk of poor a ura y still exists, and omputational osts are higher [76, 77℄. 2.1.4 Free form deformation (FFD) Computer graphi s requires large graphi al deformations su h as stret hing, twisting and other surfa e morphing operations; soft obje t animation (SOA) algorithms were developed to help with the geometry morphing required in graphi s animation [86℄. In SOA algorithms, the obje t surfa e is treated as a pie e of rubber and the desired deformation an be obtained by applying loads on it, so geometry morphing an be obtained without a hange in surfa e topology. The surfa e itself an be parametrized using Bezier or B-splines or even NURBS splines. A related approa h, the free form deformation algorithm (FFD) treats the geometry as a void in a box-shaped pie e of rubber. Deformation an be ontrolled by moving ontrol points pla ed on the outer surfa e of the box; the interior of the rubber box with its void is parametrized using a tensor produ t of three spline representations (one in ea h oordinate dire tion). Sederberg and Parry developed an algorithm that uses the FFD on ept with Bezier tri-variate volume representation [78℄. A disadvantage of the FFD method is that it requires large numbers of ontrol points to obtain lo al deformation in the deformed geometry. However, Borrel and Rappoport presented a method to allow lo al shape deformation via FFD by introdu ing a set of ontrol points with onstrained lo al B-splines that an be used to obtain deformation in a radius of inuen e determined by the designer [10℄. 2.1.5 Multidis iplinary aerodynami /stru tural shape optimization using deformations (MASSOUD) The MASSOUDmethod is an analogy to analyti al methods that tries to parametrize the deformations in the geometry rather than the geometry itself. It also uti- lizes SOA algorithms and allows strong lo al deformation ontrol. The MASSOUD parametrization requires few design variables be ause it parametrizes the deforma- tion. Samareh has applied the MASSOUD method to parametrize a simple wing, a wing body blend, and a omplex air raft onguration with su ess [76℄. Nielsen and Anderson su essfully adopted the MASSOUD parametrization s heme for aerody- 15 nami optimization of turbulent ow using unstru tured grids [64℄. 2.2 New Least Squares Parametrization Te hnique A new least-squares based surfa e parametrization method is presented in this se - tion; the airfoil surfa e is parametrized using pie e-wise polynomials whose oe- ients are found by solving a least squares problem [6℄. This se tion also des ribes how to implement a thi kness onstraint with the new parametrization and presents some validation test ases. 2.2.1 Airfoil geometry parametrization In the proposed te hnique, the geometry is parametrized using pie e-wise polyno- mials found by a least-squares t. The parametrization polynomials are ontrolled by a set of ontrol points and satisfy C2 ontinuity at their meeting points. The airfoil upper and lower surfa es are represented using two least square splines for ea h surfa e as shown in Fig. 2.1. The polynomials 3 used are P1(x) = a0 √ x+ a1x+ a2x 2 + a3x 3 0 < x < L1 P2(x) = b0 + b1x+ b2x 2 + b3x 3 L1 < x < L (2.4) where x is the normalized hord-wise position, L1 is hord-wise position that sep- arates the polynomial regions, and L = 1. These polynomials are suitable for an airfoil with a rounded leading edge due to the existen e of the √ x term in P1 (x), whi h gives an innite slope at x = 0. The x and y oordinates of the design ontrol points shown in Fig. 2.1 are used to nd the values of the polynomial oe ients. These polynomials must satisfy ontinuity of value, slope, and urvature at their meeting point x = L1. These onditions an be written as Value ontinuity: P1(L1)− P2(L1) = 0 (2.5) 3 Te hni ally, P1 is not a polynomial be ause of the presen e of the √ x term used to give innite slope and nite radius of urvature at x = 0. However, the label is onvenient and not overly onfusing. 16 Slope ontinuity: P ′ 1(L1)− P ′ 2(L1) = 0 (2.6) Curvature ontinuity: P ′′ 1 (L1)− P ′′ 2 (L1) = 0 (2.7) An additional onstraint should be added on P2(L) to assure zero thi kness at the trailing edge. The above hard onstraints should be stri tly satised by the geom- etry parametrization polynomials; they an be written in matrix form as BP = 0 where B = √ L1 L1 L 2 1 L 3 1 −1 −L1 −L21 −L31 1 2 √ L1 1 2L1 3L 2 1 0 −1 −2L1 −3L21 −1 4 √ L31 0 2 6L1 0 0 −2 −6L1 0 0 0 0 1 L L2 L3 , P = a0 a1 a2 a3 b0 b1 b2 b3 The free parameters are hosen to best approximate the y oordinates of the airfoil shape ontrol points; the x oordinates of the ontrol points are xed. The resulting least square system with onstraints applied an be expressed as min ‖AP − c‖2 subje t to BP = 0 (2.8) where A ontains powers of the x oordinates at the design points so that AP gives the y oordinates of the parametrized shape at the design points and c ontain the a tual y oordinates of the design points. I have used a set of ontrol points that ontrols the airfoil shape fun tions instead of using the oe ients of the shape fun tions to ease onstraints and boundary denitions. To give an example of how this least-squares system is onstru ted, onsider the parametrized airfoil surfa e 17 Figure 2.1: Least square surfa e presentation of RAE 2822 airfoil using two polyno- mials P1 (x) & P2 (x) tted using nine ontrol points. shown in Fig. 2.1. Six ontrol points lies in the region of the polynomial P1 (x), while three points lies in the P2 (x) region. The orresponding least-squares system is √ x1 x1 x 2 1 x 3 1 0 0 0 0√ x2 x2 x 2 2 x 3 2 0 0 0 0√ x3 x3 x 2 3 x 3 3 0 0 0 0√ x4 x4 x 2 4 x 3 4 0 0 0 0√ x5 x5 x 2 5 x 3 5 0 0 0 0√ x6 x6 x 2 6 x 3 6 0 0 0 0 0 0 0 0 1 x7 x 2 7 x 3 7 0 0 0 0 1 x8 x 2 8 x 3 8 0 0 0 0 1 x9 x 2 9 x 3 9 a0 a1 a2 a3 b0 b1 b2 b3 = y1 y2 y3 y4 y5 y6 y7 y8 y9 (2.9) Be ause the onstraint equation, Eq. 2.8, has a zero right hand side, the solution ve tor P must lie in the null spa e of the onstraint equations, i.e. it should be a linear ombination of the null spa e basis of the onstraint equations. The matrix B is full row rank and to nd the null spa e basis of it, QR fa torization an be used: BT = QR 18 Q = [ −→q1 −→q2 −→q3 −→q4 | −→q5 −→q6 −→q7 −→q8 ] = [ Q1 | Q2 ] The ve tors of the [Q2] matrix are unit ve tors forming a basis of the null spa e of matrix B. The solution ve tor P must be a linear ombination of the ve tors of the matrix Q2 P = z1 · −→q5 + z2 · −→q6 + z3 · −→q7 + z4 · −→q8 = Q2z (2.10) Substituting Eq. 2.10 into Eq. 2.8 AP = AQ2z = c (2.11) Solving the least-squares system by singular value de omposition for numeri al sta- bility, z = [AQ2] † c (2.12) where the re tangular matrix [AQ2] † is the pseudo-inverse of [AQ2]. Finally, P = Q2 [AQ2] † c (2.13) Equation 2.13 gives the relationship between the polynomial oe ients P and the y lo ations of the ontrol points c. The sensitivity of the polynomial oe ients with respe t to the y lo ation of the ith ontrol point, whi h is needed to al ulate the mesh sensitivity ∂M/∂D, is the ith olumn of the matrix Q2 [AQ2] † .The dependen y of an airfoil surfa e point −→rb = [ xa ya ] on a design variable Di an be found from ∂−→rb ∂Di = [ 0 ∂a0 ∂Di √ xa + ∂a1 ∂Di xa + ∂a2 ∂Di x2a + ∂a3 ∂Di x3a ] 0 < xa < L1 (2.14) ∂−→rb ∂Di = [ 0 ∂b0 ∂Di + ∂b1 ∂Di xa + ∂b2 ∂Di x2a + ∂b3 ∂Di x3a ] L1 < xa < L 19 where [ ∂a0 ∂Di ∂a1 ∂Di ∂a2 ∂Di ∂a3 ∂Di ∂b0 ∂Di ∂b1 ∂Di ∂b2 ∂Di ∂b3 ∂Di ]T is the ith ve tor of the on- strained pseudo inverse matrix Q2 [AQ2] † . This pro edure an be extended to parametrize airfoil surfa e using any number of pie e-wise polynomials, with an a - ompanying in rease in system size. Leading edge radius and trailing edge thi kness onstraints an be added to the C2 ontinuity requirements forming the onstraint system BP = d, and the ontrol point oordinates are used to onstru t the system AP = {yc}. In pra ti e, all airfoil surfa e points are used in the least squares system to nd the polynomials oe ients. After nding the polynomials oe ients, the designer an sele t a set of ontrol points whi h lies on the parametrized surfa e; the least number of ontrol points is two per polynomial. The less ontrol points used, the less ontrol of airfoil geometry. I have sele ted nine ontrol points at spe i hord- wise x-stations based on my engineering sense; it turns out that my sele ted set of ontrol points were able to produ e sho k free optimal proles as will be shown in the next hapters, however, sele tion of the ontrol points x-stations an be made by formulating a simple minimization problem. In this problem, a set of airfoil geometri data gathered and the obje tive is to minimize the RMS error between the original airfoil surfa es and the parametrized airfoil surfa es, where the design variables are the ontrol points x-stations. 2.2.2 Thi kness onstraint Some air raft fuel tanks are pla ed inside the wing, and the main landing gear of some air raft are stored in the wings after being retra ted. In addition, the wing must have su ient bending rigidity from a stru tural point of view; therefore, the wing thi kness is a onstraint at some hord-wise stations. This subse tion demonstrates how to add a thi kness onstraint to the parametrization. Suppose the airfoil is parametrized using four polynomials, as shown in Fig. 2.2: two for the upper surfa e P1(x) = a0 √ x+ a1x+ a2x 2 + a3x 3 0 < x < L1 P2(x) = b0 + b1x+ b2x 2 + b3x 3 L1 < x < L 20 Figure 2.2: Parametrized airfoil using four polynomials and two for the lower surfa e P3(x) = c0 √ x+ c1x+ c2x 2 + c3x 3 0 < x < L1 P4(x) = d0 + d1x+ d2x 2 + d3x 3 L1 < x < L If the thi kness tc need to be onstrained at a station xc where 0 < xc ≤ L1 , this onstraint will be Ct : P1(xc)− P3(xc) = tc (2.15) Ct : ( a0 √ xc + a1xc + a2x 2 c + a3x 3 c )− (c0√xc + c1xc + c2x2c + c3x3c) = tc Equation 2.15 provides a link between the upper and lower surfa e polynomials; therefore, a oupled least-squares system needs to be onstru ted and solved; the 21 hard onstraint equations are expressed as B 00 B P1(xc) {0}1×4 −P3(xc) {0}1×4 a0 a1 . . . d2 d3 = 0 0 . . . 0 tc (2.16) [Ba] {Pa} = {da} (2.17) The global least-squares system is[ A 0 0 A ] {Pa} = { cu cl } (2.18) [Aa] {Pa} = {ca} The above hard onstraint equations do not have a null right hand side due to the thi kness onstraint; therefore, the solution (polynomials oe ients) does not belong to the null spa e of Ba. The solution pro edure for the onstrained least- squares problem expressed by Eq. (2.16,2.18) is des ribed by Masuda et al [52℄: Apply QR fa torization to get BTa = QR. Let Q = [ Q1 Q2 ] , R = [ R1 0 ] where Q1 ontains the rst 7 olumns of Q and Q2 ontains the rest of the olumns, and R1is the rst 7× 7 sub matrix of R. The hard onstraints an now be written as BPa ≡ RTQTPa = da Let QTPa = [ y z ] . Then Pa = Q [ y z ] = Q1y +Q2z ∴ RTQTPa = [ RT1 0 ] · [ Q1 Q2 ]T Pa = da 7−→ RT1 y + 0 = da ∴ y = R−T1 da (2.19) 22 Be ause R is upper triangular, Equation 2.19 an be solved using forward substitution. The soft onstraints in equation 2.18 an be rewritten as AaPa = Aa [ Q1y Q2z ] = ca (2.20) ∴ AaQ2z = ca −AaQ1y (2.21) ∴ z = [AaQ2] † [ca −AaQ1y] (2.22) Finally, the polynomial oe ients an be found as Pa = Q1y +Q2z: Pa = [ Q1 −Q2 [AaQ2]†AaQ1 ] R−T1 da +Q2 [AaQ2] † ca (2.23) Using the last equation we noti e that the sensitivity of the polynomial oe- ients with respe t to the ith ontrol point y lo ation, yci, is the i th ve tor of the matrix [Q2] [AQ2] † . 2.2.3 Validation ases In this subse tion, testing of the proposed least-squares surfa e parametrization is arried out; parametrization of various types of airfoils is done to show that the proposed geometry parametrization s heme has the required exibility to represent various types of airfoils. These airfoils in lude NACA 4-, 5-, and 6-digit series, laminar ow, supersoni , and super- riti al airfoil se tions. 4 Figures 2.32.10 show the least-squares tting polynomials for various airfoils, while Table 2.1 shows the RMS error in the parametrized geometry for upper and lower surfa es. The last parametrized airfoil, the laminar LV2 5 airfoil of German Aerospa e Center, was parametrized using 10 polynomials, 5 for ea h of its surfa es. Although the RMS error is small (order of 10−5 hord for all parametrization method), as shown in Table 2.1, Fig. 2.11 shows that the dieren e of the surfa e pressure distribution, espe ially at the peak velo ity at the leading edge upper surfa e, is still signi ant. There are many x-stations at whi h LV2 airfoil hanges its urvature, and also a 4 Readers interested in parti ulars for the NACA airfoils are referred to Abbott and von Doen- ho's textbook [36℄. The seminal referen e for the RAE airfoil is Cook et al [18℄. 5 LV2 geometry obtained by personal onta t with German Aerospa e Center (DLR) resear hers. 23 Figure 2.3: Parametrized NACA 0011SC, 'O' are the original airfoil ordinates Figure 2.4: Parametrized NACA 0012, 'O' are the original airfoil ordinates 24 Figure 2.5: Parametrized NACA 6509, 'O' are the original airfoil ordinates Figure 2.6: Parametrized NACA 16006, 'O' are the original airfoil ordinates 25 Figure 2.7: Parametrized NACA 63412, 'O' are the original airfoil ordinates Figure 2.8: Parametrized NACA 644421, 'O' are the original airfoil ordinates 26 Figure 2.9: Parametrized RAE2822, 'O' are the original airfoil ordinates long interval of almost innite urvature value whi h makes it hard for the geometry parametrization s heme to a urately present it and this auses the u tuation in pressure resulted in the parametrized airfoil be ause of la k of a urate presentation of urvature u tuations. Better mat hing an be obtained but this will in rease the number of geometry design variables signi antly. This ase illustrates learly the trade o between a urately representing the geometry (whi h will hange during optimization iteration) and hoosing a reasonable number of design variables; this hoi e must be left for the designer. If the hange in pressure distribution be omes una eptable, or if the pressure distribution requires large number of airfoil ontrol points (design variables) to be a urately represented, parametrization the geometry perturbation is the more attra tive option. The values of RMS error in Table 2.1 are small ompared to the maximum airfoil thi kness value: only the NACA 0011SC has an RMS error of more than 0.2% of maximum thi kness. However, if the RMS error was large, in reasing the number of surfa e parametrization polynomials would redu e the error. 27 Figure 2.10: Parametrized LV2 laminar airfoil, parametrized using 5 polynomials per surfa e, 'O' are the original airfoil ordinates Figure 2.11: LV2 airfoil parametrized using dierent methods with 20 design vari- ables, presented with permission of Brezillon. 28 Airfoil Upper surfa e error Lower surfa e error Maximum Thi kness NACA 0011SC 6.37 · 10−4 6.37 · 10−4 11 · 10−2 NACA 0012 2.35 · 10−4 2.35 · 10−4 12 · 10−2 NACA 6409 9.28 · 10−5 1.24 · 10−4 9 · 10−2 NACA 16-006 1.31 · 10−6 2.03 · 10−5 6 · 10−2 NACA 63-412 1.17 · 10−5 5.63 · 10−5 12 · 10−2 NACA 64-421 1.70 · 10−4 7.50 · 10−5 21 · 10−2 RAE 2822 1.20 · 10−5 3.20 · 10−5 slightly > 12 · 10−2 LV2 1.45 · 10−5 3.22 · 10−5 slightly > 12 · 10−2 Table 2.1: RMS error in dierent parametrized airfoil geometries 29 Chapter 3 MESH MORPHING AND MESH SENSITIVITY 3.1 Mesh Morphing The modi ation of the aerodynami shape during optimization requires a hange of the mesh that presents the shape. This an be done by grid regeneration around the new geometry but this is time onsuming and will hange the dis retization error [51℄. Another strategy is to adapt the old mesh to t the new shape of the airfoil using mesh movement. The tension spring analogy is one of the most widely used mesh deformation strategies for aerodynami optimization. The main idea is to repla e the grid edges by springs with stiness inversely proportional to their length. The boundary points that lie on the airfoil surfa e are moved with displa ement spe ied by the optimizer, far eld points are kept xed, and the interior point displa ements are determined by equilibrium of the spring network [7℄. For large grid displa ements, the linear spring analogy is not robust, and negative area ells an present after mesh morphing. Farhat et. al. improved the tension spring analogy by adding torsional springs at the grid nodes to prevent element ipping [25℄. Ea h edge fa es two angles as shown in Fig. 3.1; the edge stiness is modied to in lude terms with the re ipro al of the sine of these angles. This allows the edge stiness to grow to innity if the angle tends to be zero and therefore prevent element ipping. Another strategy is to modify the mesh by solving a linear elasti ity problem in whi h the boundary displa ements are known [14℄; the element modulus of elasti ity an be the re ipro al of the distan e from the wall, or it an be the re ipro al of the element size [82℄. The later strategy was applied by Stein et al. and the results showed that this method is robust, espe ially for vis ous al ulations. In this ase, 30 Figure 3.1: S hemati drawing of an edge pq and its fa ing angles a and b. the elements of the boundary layer experien ed small geometri al hange while the elements away from the airfoil experien ed larger hanges [80℄. However, the linear elasti ity mesh movement s heme is omputationally expensive ompared to the spring analogy method. The semi-torsional spring analogy method is adopted in this resear h due to its simpli ity and robustness. Consider the edge pq shown in Fig 3.1. The relationship between the for es and the node displa ements when treating this edge as a spring follows Hooke's Law as Fpx Fpy Fqx Fqy = 1lpq ( 1 + 1 sin (θa) + 1 sin (θb) ) −1 0 1 0 0 −1 0 1 1 0 −1 0 0 1 0 −1 upx upy uqx uqy (3.1) where lpq is the length of edge pq, and θa and θb are the angles fa ing the edge. After the assembly of the global stiness matrix, the system of equations that relates grid point displa ement with nodal for es an be written as [ Kii Kib Kbi Kbb ]{ Ui Ub } = { 0 Fb } (3.2) 31 where Ui and Ub are the interior and boundary mesh point displa ements respe tively. We do not need to know the values of boundary nodal for es Fb be ause the boundary points displa ement ve tor Ub is known expli itly: it is the deformation required in the airfoil prole to minimize the obje tive fun tion. Therefore, Eq 3.2 an be written as [ Kii Kib 0 I ]{ Ui Ub } = { 0 Ub } (3.3) Substituting Uki = −−→ rk+1i − −→ rki , Ub = −→rb −−→rboin Eq 3.3 we get[ Kii Kib 0 I ]k{ −−→ rk+1i−→rb } = [ Kii Kib 0 I ]k{ −→ rki−→rbo } + { 0 Ub } (3.4) where −→rbo is the initial position ve tor of the boundary points before mesh morphing. The stiness matri es [Kii] and [Kib] depend on the mesh fa e lengths, whi h hange during the mesh morphing stage; therefore, Eq. 3.4 is a non-linear equation and needs to be solved iteratively. 3.2 Testing Mesh Morphing Farhat et al [25℄ demonstrated the robustness of the semi-torsional mesh movement s heme. In this se tion, testing results are presented to demonstrate that the ur- rent implementation of this s heme shows the same good behavior for several high- deformation ases. The rst test ase is an unstru tured triangular mesh around a NACA 0012. The thi kness of the airfoil is doubled, whi h means that the airfoil surfa e points will translate by several omputational ells in the y-dire tion. Fig- ure 3.2 shows that even with this large displa ement (multiple ells) of the boundary, no ells are inverted and no edges interse ted with another. Displa ement de reases with the distan e from the airfoil surfa e, and there is almost no movement on the symmetry line. The se ond ase tests mesh movement when the outer boundary is hanging. The original mesh is an unstru tured triangular mesh (shown in Fig 3.3a). The outer boundaries are hanged, redu ing the total area by almost 50% and turning the original right-angled orners nearly into usps. Figure 3.3b shows that the semi- 32 Figure 3.2: Mesh movement s heme results of doubling the thi kness of NACA 0012 torsional mesh movement was apable of adapting the mesh in the entire eld without element ipping. 3.3 Mesh Sensitivity For gradient-based optimization, as shown in the next hapter, the gradient ompu- tation requires al ulation of the dependen y of the residual on the design variables ∂R/∂D, whi h in turn depends on the evaluation of mesh sensitivity. The mesh sensitivity tells how mesh points translate in the (x, y) plane with the perturbation of the geometri parametrization ontrol points (design variables). This translation, obviously, depends on the mesh movement s heme. Truong et. al. ompared al- gebrai mesh movement to linear elasti ity mesh movement s hemes to study the impa t of the adopted mesh movement s heme on the nal optimized shape. There was a noti eable dieren e between the nal optimal airfoil shapes for a subsoni ase, although the dieren e in the optimal obje tive fun tion value was of order of dis retization error. For a transoni test ase, the dieren e in the optimal shapes was almost negligible [83℄. The mesh sensitivity with respe t to one of the design variables ∂M/∂Di (that is, 33 (a) Initial mesh (b) Final mesh Figure 3.3: Mesh movement results, large outer boundary deformation of a re tan- gular domain the hange in mesh point lo ations with a hange of a design variable) is al ulated by dierentiating Eq. 3.3:[ Kii Kib 0 I ]{ ∂Ui ∂Di ∂Ub ∂Di } = { 0 ∂−→rb ∂Di } ∴ { ∂M ∂Di } = { ∂Ui ∂Di ∂Ub ∂Di } = [ Kii Kib 0 I ]−1{ 0 ∂−→rb ∂Di } (3.5) where ∂−→rb/∂Di is obtained using Eq. 2.14 and is related to the design variables via the pseudo inverse of the onstrained least-squares system solved in parametrizing the geometry. 34 Chapter 4 GRADIENT CALCULATION USING ADJOINT APPROACH Gradient al ulation plays a key role in gradient-based optimization. The traditional nite dieren e strategy is omputationally expensive as it requires at least as many CFD simulations as the number of design variables to ompute the gradient; ea h of these is the solution to a large non-linear system of equations. The forward strategy to ompute the gradient is less expensive as it requires the al ulation of the ow sensitivity with respe t to the design variables, and uses the omputed ow sensitivity to ompute the gradient; the forward strategy requires solving a number of linear systems equal to the number of design variables to nd the ow property sensitivity with respe t to all the design variables. The adjoint strategy is omputationally heaper; it requires the solution of one linear system whose right hand side is the dependen y of the aerodynami obje tive fun tion on the ow eld properties. Due to its numeri al e ien y and the orresponding redu tion in omputational eort, the adjoint strategy is adopted in this resear h. 4.1 Forward and Adjoint Formulations The obje tive fun tion, F , for aerodynami optimization is a fun tion of the design variables, D, and the ow eld solution at the surfa e points of the boundary ontrol volumes Us F = F (Us,D) (4.1) Us is expressed most onveniently in primitive variables: U = ( ρ u v P )T . Consider, for instan e, the lift and drag for es of a two dimensional airfoil, whi h are perpendi ular and parallel, respe tively, to the in oming ow, whi h is in lined 35 at an angle α to the airfoil hord. These an be evaluated as follows: FL = − {˛ Psnx ds } sinα+ {˛ Psny ds } cosα FD = {˛ Psnx ds } cosα+ {˛ Psny ds } sinα or in dis rete form, FL = − ∑ wsPsnx sinα+ ∑ wsPsny cosα (4.2) FD = ∑ wsPsnx cosα+ ∑ wsPsny sinα (4.3) where Ps is the pressure at surfa e integration point pressure, nx and ny are the unit normal omponents at the surfa e integration point, α is the angle of atta k, and ws is the ar length asso iated with the surfa e integration point. Note that this form uses the dimensional pressure and oordinates gives the dimensional lift and drag for es. The lift and drag oe ients, whi h are the non-dimensional equivalents, are identi al in form but use non-dimensional pressure and oordinates. These dis rete integrals expressed as a fun tion of geometri and ow properties of the ontrol volume su h as the length of ea h fa e and the unit normal at ea h Gauss point. The geometri properties depend on the design variables through the mesh sensitivity, while the ow properties at the Gauss points depend on the ow properties of the ontrol volume itself and its neighbors, whi h in turns depend on the mesh and the boundary shape whi h ultimately depends on the geometri design variables. The gradient of the obje tive fun tion an be obtained by using the hain rule dF dD = ∂F ∂Ubg ∂Ubg ∂U ∂U ∂M ∂M ∂D + ∂F ∂M ∂M ∂D (4.4) where Ubg is the boundary Gauss point ow properties, U is the ontrol volume averaged solution written in primitive variables, and M is the mesh point lo ations. ∂M/∂D is the mesh dependen y on the design variables whi h omputations are presented in details in Chapter 3. The residual of the ow governing equations an be written as a fun tion of the ow eld solution U and mesh geometri design variables D. If we apply the onstraint that the ow solution is onverged regardless 36 of variations in the design variables, we an write dR dD = ∂R ∂M ∂M ∂D + ∂R ∂U ∂U ∂M ∂M ∂D = 0 (4.5) Solving this for the solution sensitivity ∂U ∂D ≡ ∂U ∂M ∂M ∂D , we get ∂U ∂M ∂M ∂D = − [ ∂R ∂U ]−1 · ∂R ∂M ∂M ∂D = − [ ∂R ∂Q ∂Q ∂U ]−1 · { ∂R ∂M ∂M ∂D } (4.6) where the last equality expands the residual Ja obian with respe t to the primitive ow variables into the produ t of the residual Ja obian with respe t to the onserved ow variables (whi h is used in impli it ow solvers) and a hange of variables for onserved to non- onserved. Note that ∂Q ∂U is a blo k diagonal matrix. Substituting Eq.4.6 in Eq. 4.4, we get the forward formulation for gradient omputation dF dD = − ∂F ∂Ubg ∂Ubg ∂U [ ∂R ∂Q ∂Q ∂U ]−1 · { ∂R ∂M ∂M ∂D } + ∂F ∂M ∂M ∂D (4.7) This form of the gradient requires solving as many linear systems as there are design variables in the optimization problem. Taking the transpose of Eq. 4.7, dF dD T = − { ∂R ∂M ∂M ∂D }T · [ ∂R ∂Q ∂Q ∂U ]−T { ∂F ∂Ubg ∂Ubg ∂U }T + { ∂F ∂M ∂M ∂D }T (4.8) where the residual sensitivity to mesh movement an be written using the hain rule as: ∂R ∂M = ∂R ∂AΩi ∂AΩi ∂M + ∑ fa esΩi { ∂R ∂nx dnx dM + ∂R ∂ny dny dM + (4.9) ∂R ∂wi dwi dM + ∂R ∂U fa e ∂U fa e ∂M } We get the adjoint method, presented by A. Jameson [37, 41, 43, 38, 40℄. Now only one linear system solve is required. However, this linear system solve requires 37 expli itly forming the global Ja obian matrix ∂R ∂Q be ause its transpose is required. Be ause the CFD solver used in this thesis an form the global Ja obian matrix expli itly (author?) [57℄, the transpose of the Ja obian an easily be formed as well. The Ja obian from the last GMRES iteration is re-used, so the omputational eort of solving the adjoint problem is redu ed to solving this linear system, a ost on the order of 1% of the CFD simulation omputational eort. To ease the programming eort when hanging the obje tive fun tion, we form three adjoint problems, one ea h for the lift oe ient CL, the drag oe ient CD, and the moment oe ient CM to nd ∂CL/∂D, ∂CD/∂D and ∂CM/∂D, respe - tively. These aerodynami oe ient gradients an be used to evaluate any obje tive fun tion gradient for a fun tion that depends on aerodynami for es: dF dD = f ( dCL dD , dCD dD , dCM dD ) The solution pro edure of the three adjoint problems an be summarized as follows, Using the steady state ow solution, we onstru t the CFD simulation Ja obian matrix ∂R/∂Q expressed in Eq. 1.6. We onstru t ∂R/∂U as: ∂R ∂U = ∂R ∂Q ∂Q ∂U where ∂Q/∂U is the transformation matrix from onservative to primitive ow variables and is blo k diagonal. We onstru t ∂F ∂Ubg ∂Ubg ∂U , where ∂F/∂Ubg is the analyti dependen y of the ob- je tive fun tion on the primitive ow properties at the airfoil surfa e points, and ∂Ubg/∂U is the dependen y of surfa e point primitive ow properties on the ontrol volume average values of the primitive ow properties. The latter is known as a side ee t of solution re onstru tion. We solve the three linear systems [ ∂R ∂U ]T ΨL,D,M = [ ∂CL,D,M ∂Ubg ∂Ubg ∂U ]T to get the adjoint ve tors ΨL,D,M . We onstru t the obje tive fun tion gradient ve tor dCL,D,M dD T = { ∂F ∂M ∂M ∂D }T −{ ∂R ∂M ∂M ∂D }T ΨL,D,M , whi h requires only a ve tor dot produ t for ea h design 38 variable. As an example of how to use dCL,D,M dD to onstru t the gradient of an aerodynami fun tion, onsider the following obje tive fun tion whi h represents a typi al drag minimization fun tion with a lift onstraint applied using a penalty term F = CD + k1 (CL − CLc)2 (4.10) The gradient of the above fun tion an be written as dF dD = dCD dD + 2k1 (CL − CLc) dCL dD (4.11) also ∂F ∂D = ∂CD ∂D + 2k1 (CL − CLc) ∂CL ∂D (4.12) The dis rete (invis id) forms of CL and CD follow from the non-dimensionalization of Eqs. 4.2 and 4.3, where now the pressure and oordinates are non-dimensionalized: CL = − (∑ wsPsnx ) sinα+ (∑ wsPsny ) cosα (4.13) CD = (∑ wsPsnx ) cosα+ (∑ wsPsny ) sinα and their partial derivatives with respe t to the geometry design variables are ∂CL ∂D = − (∑ ∂ws ∂D Psnx ) sinα+ (∑ ∂ws ∂D Psny ) cosα (4.14) − (∑ wsPs ∂nx ∂D ) sinα+ (∑ wsPs ∂ny ∂D ) cosα ∂CD ∂D = (∑ ∂ws ∂D Psnx ) cosα+ (∑ ∂ws ∂D Psny ) sinα+(∑ wsPs ∂nx ∂D ) cosα+ (∑ wsPs ∂ny ∂D ) sinα Using Eq. 4.9 to al ulate ∂R/∂M , we need to ompute the terms ∂AΩi/∂M , ∂nx/∂M , ∂ny/∂M , and ∂wi/∂M , whi h only depend on mesh points' spatial lo a- tions, while the terms ∂R/∂nx and ∂R/∂ny are obtained by dire t dierentiation of the Roe ux, and nally (∂R/∂Uface) · (∂Uface/∂M) is obtained by dierentiating 39 Figure 4.1: S hemati drawing of an element fa e κ, and illustration of its left and right sides the fa e property re onstru tion s heme whi h depends on CFD solver te hnology. The remainder of this se tion will fo us on the derivatives of the residual R, with the derivatives on the geometri terms in the following se tion. Let us onsider the residual ontribution of fa e κ of the ontrol volume Ωi shown in Fig. 4.1. The dis rete form of this edge's residual ontribution to the ontrol volume Ωi is Rκ,Ωi = 1 AΩi J∑ j=1 { 1 2 ( ~FRj + ~FLj − |△F1|j − |△F4|j − |△F5|j ) wj } (4.15) where J = 1 when using one Gauss integration point at the middle of fa eκ for a se ond order a urate s heme. For the fourth order s heme, we have two Gauss integration points, i.e. J = 2, on the fa e κ. wj is the integration weight asso iated with the Gauss point j. 40 To nd ∂Rκ,Ωi/∂AΩi , Eq. 4.15 an be dire tly dierentiated to get ∂Rκ,Ωi ∂AΩi = − 1 (AΩi) 2 ∑ j { 1 2 ( ~FRj + ~FLj − |△F1|j − |△F4|j − |△F5|j ) wj } = − 1 AΩi Rκ,Ωi (4.16) The term ∂Rκ,Ωi/∂nx,y an be written as ∂Rκ,Ωi ∂nx,y = 1 AΩi ∑ j { 1 2 ( ∂ ~FRj ∂nx,y + ∂ ~FLj ∂nx,y − ∂ |△F1|j ∂nx,y − ∂ |△F4|j ∂nx,y − ∂ |△F5|j ∂nx,y ) wj } (4.17) where omponent terms an be expanded using the denition of the Roe ux to give: ∂ ~FR,Lj ∂nx = (ρu)R,L( ρu2 + p ) R,L (ρuv)R,L ({E + p}u)R,L , ∂ ~FR,Lj ∂ny = (ρv)R,L (ρuv)R,L( ρv2 + p ) R,L ({E + p} v)R,L and ∂ |△F1|j ∂nx = uŨn√ Ũ2n ( △ρ− △P ã2 ) 1 ũ ṽ ũ2+ṽ2 2 + ρ̃ 0 △u− nx△U △v − ny△U ũ△u+ ṽ△v − Ũn△U + ∣∣∣Ũn∣∣∣ ρ̃ 0 −△U − nx△u −ny△v −ũ△U − Ũn△u 41 ∂ |△F1|j ∂ny = vŨn√ Ũ2n ( △ρ− △P ã2 ) 1 ũ ṽ ũ2+ṽ2 2 + ρ̃ 0 △u− nx△U △v − ny△U ũ△u+ ṽ△v − Ũn△U + ∣∣∣Ũn∣∣∣ ρ̃ 0 −nx△u −△U − ny△v −ṽ△U − Ũn△v ∂ ∣∣∣△F̃4,5∣∣∣ ∂nx = ( Ũn ± ã ) ũ√( Ũn ± ã )2 (△P ± ρ̃ã△U 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã + ∣∣∣Ũn ± ã∣∣∣ (±ρ̃ã△u 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã + ∣∣∣Ũn ± ã∣∣∣ (△P ± ρ̃ã△U 2ã2 ) 0 ±ã 0 ±ũã 42 ∂ ∣∣∣△F̃4,5∣∣∣ ∂ny = ( Ũn ± ã ) ṽ√( Ũn ± ã )2 (△P ± ρ̃ã△U 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã + ∣∣∣Ũn ± ã∣∣∣ (±ρ̃ã△v 2ã2 ) 1 ũ± nxã ṽ ± nyã h̃o ± Ũnã + ∣∣∣Ũn ± ã∣∣∣ (△P ± ρ̃ã△U 2ã2 ) 0 0 ±ã ±ṽã To nd ∂Rκ,Ωi/∂wj , Eq. 4.15 is dierentiated; we get ∂Rκ,Ωi ∂wj = 1 (AΩi) { 1 2 ( ~FRj + ~FLj − |△F1|j − |△F4|j − |△F5|j )} (4.18) The dependen y of the fa e residual ontribution Rκ,Ωi on the re onstru ted ow properties at the fa e U fa e an found as the sum of the dependen y on the right and left fa e ow properties as ∂Rκ,Ωi ∂U fa e = ∂Rκ,Ωi ∂UR + ∂Rκ,Ωi ∂UL (4.19) = 1 (AΩi) ∑ 1 2 ∂ ~FRj ∂UR + ∂ ~FLj ∂UL − ∂ { |△F1|j + |△F4|j + |△F5|j } ∂UR − ∂ { |△F1|j + |△F4|j + |△F5|j } ∂UL wj The terms in Eq. 4.19 an be obtained using symboli manipulator like Maple® or Matlab®. Another possible approa h is to use an automati dierentiation pa k- age to dierentiate the C, C++, or FORTRAN fun tion that omputes the fa e ux. 43 Figure 4.2: General triangular element with unit normal n̂ on one of its fa es κ. Finite dieren es an also be easily implemented to evaluate ∂RκΩi/∂U fa e ∂Rκ,Ωi ∂U fa e = Rκ,Ωi (UR + ǫ)−Rκ,Ωi (UR − ǫ) 2ǫ + Rκ,Ωi (UL + ǫ)−Rκ,Ωi (UL − ǫ) 2ǫ +o (ǫ)2 (4.20) where ǫ is hosen to be 10−8 in the above entral nite dieren e formula so that the error will be on the order of ma hine zero. The ANSLib CFD ode an use both approa hes; they lead to nearly identi al answers. 4.2 Element and Fa e Geometri Properties Dependen y on Mesh Coordinates Element and fa e geometri properties depend dire tly on the spatial oordinates of the verti es. Fig. 4.2 shows a s hemati drawing of a triangular element used for ell entered nite volume simulations, with its three verti es and one of its three fa es labeled for later referen e. In the next subse tions, the pro edure for evaluating geometri properties like element area, fa e length, and fa e normals will be presented in addition to the mesh oordinate dependen y of these properties. 44 4.2.1 Element area mesh dependen y The element area of the triangular element shown in Fig 4.2 an be obtained as half of the ross produ t of the two ve tors −→r12,−→r13 A = σ (x2 − x1) (y3 − y1)− (y2 − y1) (x3 − x1) 2 (4.21) The verti es' spatial lo ations are related to the airfoil surfa e points via the mesh movement s heme, while the airfoil surfa e points are related to the design vari- ables through the pseudo inverse matrix in the least-squares t during geometry parametrization. It worth mentioning that any general polygon an be split into several triangles and the area of ea h triangle an be al ulated using Eq. 4.21. The element area mesh dependen y an be al ulated by dire t dierentiation of Eq. 4.21 dA dM = ∂A ∂x1 ∂x1 ∂M + ∂A ∂x2 ∂x2 ∂M + ∂A ∂x3 ∂x3 ∂M (4.22) + ∂A ∂y1 ∂y1 ∂M + ∂A ∂y2 ∂y2 ∂M + ∂A ∂y3 ∂y3 ∂M where ∂A ∂x1 = σ (y2 − y3) 2 ∂A ∂x2 = σ (y3 − y1) 2 ∂A ∂x3 = σ (y1 − y2) 2 ∂A ∂y1 = σ (x3 − x2) 2 ∂A ∂y2 = σ (x1 − x3) 2 ∂A ∂y3 = σ (x2 − x1) 2 45 σ = sgn ((x2 − x1) (y3 − y1)− (y2 − y1) (x3 − x1)) 4.2.2 Fa e length mesh dependen y The length of fa e κ in Fig. 4.2, whi h is a general straight fa e in the mesh, an be al ulated as L = √ (x3(M)− x2(M))2 + (y3(M)− y2(M))2 (4.23) The fa e length mesh dependen y an be obtained as dL dM = ∂L ∂x2 ∂x2 ∂M + ∂L ∂x3 ∂x3 ∂M + ∂L ∂y2 ∂y2 ∂M + ∂ L ∂y3 ∂y3 ∂M where ∂L ∂x2 = −x3 + x2√ (x3 − x2)2 + (y3 − y2)2 = − (x3 − x2) L ∂L ∂x3 = x3 − x2√ (x3 − x2)2 + (y3 − y2)2 = (x3 − x2) L ∂L ∂y2 = −y3 + y2√ (x3 − x2)2 + (y3 − y2)2 = − (y3 − y2) L ∂L ∂y3 = y3 − y2√ (x3 − x2)2 + (y3 − y2)2 = (y3 − y2) L The previous relations are for straight fa es whi h exist all over the domain interior, but for higher order s hemes, the boundary fa es are urved to enable higher order ux integration. A boundary element with one urved fa e is shown in Fig. 4.3 The urved fa e length is obtained using numeri al integration, therefore to avoid dierentiating the numeri al integration s heme, ∂L/∂x2,3 and ∂L/∂y2,3 are evalu- 46 Figure 4.3: Boundary element with urved fa e for high order integration s heme ated using nite dieren es as ∂L ∂x2,3 = L(x2,3 + ǫ)− L(x2,3 − ǫ) 2ǫ ∂L ∂y2,3 = L(y2,3 + ǫ)− L(y2,3 − ǫ) 2ǫ (4.24) For a se ond order s heme, only one Gauss integration point exists at the middle of the fa e with wi = L, therefore dwi dM = dL dM For the fourth order s heme, two Gauss integration points are required for ux integration. Figure 4.2.2 shows a s hemati drawing of a urved edge with two Gauss points on it. Table 4.1 shows their integration weights and their parametri lo ations ki on the fa e starting from the point x2, y2. 47 Figure 4.4: S hemati drawing of a urved fa e with two Gauss points used Point 1 Point 2 wi L/2 L/2 ki ( 1 2 − 1√12)L ( 1 2 + 1√ 12 )L Table 4.1: Two point Gauss quadrature rule 48 4.2.3 Fa e normal mesh dependen y The unit normal ve tor of fa e κ that points outward from element Ωi (shown in Fig. 4.2) an be found using the three verti es of the element as n̂ = [ nx ny ] = (~c3 − βr̂23) ‖~c3 − βr̂23‖ (4.25) where ~c3 = [ 2/3x3 − (1/3x1 + 1/3x2) 2/3 y3 − (1/3 y1 + 1/3 y2) ] r̂23 = x3−x2√(|−x3+x2|)2+(|−y3+y2|)2 y3−y2√ (|−x3+x2|)2+(|−y3+y2|)2 β = (x3 − x2) (2/3x3 − (1/3x1 + 1/3x2))√ (|−x3 + x2|)2 + (|−y3 + y2|)2 + (y3 − y2) (2/3 y3 − (1/3 y1 + 1/3 y2))√ (|−x3 + x2|)2 + (|−y3 + y2|)2 The mesh dependen y of the unit normal an be written as dn̂ dM = ∂n̂ ∂x1 ∂x1 ∂M + ∂n̂ ∂x2 ∂x2 ∂M + ∂n̂ ∂x3 ∂x3 ∂M + ∂n̂ ∂y1 ∂y1 ∂M + ∂n̂ ∂y2 ∂y2 ∂M + ∂n̂ ∂y3 ∂y3 ∂M (4.26) where ∂n̂ ∂xi = α1 − α2 ‖~c3 − βr̂23‖2 α1 = ‖~c3 − βr̂23‖ · ∂ ∂xi (~c3 − βr̂23) α2 = (~c3 − βr̂23) · ∂ ∂xi ‖~c3 − βr̂23‖ 49 ∂n̂ ∂yi = α3 − α4 ‖~c3 − βr̂23‖2 α3 = ‖~c3 − βr̂23‖ · ∂ ∂yi (~c3 − βr̂23) α4 = (~c3 − βr̂23) · ∂ ∂yi ‖~c3 − βr̂23‖ i = 1, 2, 3 The terms in the last equations an be found using automati dierentiation of the unit normal expression. 4.3 Fa e Flow Properties Dependen y on The Mesh The evaluation of the ∂U fa e /∂M term in the mesh sensitivity of the residual, using Eq. 4.9, depends on the details of the CFD solver. The CFD solver used in this resear h is the Advan ed Numeri al Simulation Library (ANSLib) whi h is a multi- physi s nite volume solver apable of ondu ting CFD simulations up to fourth or- der a ura y; this solver has been built by Ollivier-Goo h and o-workers [61, 56, 54℄. In ANSLib, the ow properties are assumed to hange within ea h ontrol volume a ording to a linear polynomial for se ond order simulations, or a ubi polynomial for fourth order simulations. These polynomials are found using a least-squares re- onstru tion; the least-squares system depends on the ontrol volume properties like moments of area whi h eventually depend on mesh point lo ations. In the next two subse tions, the re onstru tion and mesh dependen y of ow properties at fa es will be presented. 4.3.1 Fa e ow properties re onstru tion Solution re onstru tion is the key to determine solver a ura y. The solution is assumed to vary linearly within the ontrol volume for se ond order a ura y; for higher order methods, ow properties are assumed to vary a ording to higher order polynomial; variation a ording to ubi polynomial is fourth order a urate. The ow solution re onstru tion method presented in this subse tion follows the method des ribed in Ollivier-Goo h and Van Altena [68℄ and Ollivier-Goo h et al [69℄. 50 For se ond order solution re onstru tion, the primitive variables U = [ ρ u v p ]T are re onstru ted at the fa e Gauss points as UR2 (x, y) = Uref + (x− xref ) ∂U ∂x ∣∣∣∣ ref + (y − yref ) ∂U ∂y ∣∣∣∣ ref (4.27) For fourth order solution re onstru tion, the re onstru tion polynomial takes the form of a third degree Taylor series expansion in two variables: UR4 (x, y) = U R 2 (x, y) + (x− xref )2 2 ∂2U ∂x2 ∣∣∣∣ ref (4.28) + (x− xref) . (y − yref) ∂ 2U ∂x∂y ∣∣∣∣ ref + (y − yref )2 2 ∂2U ∂y2 ∣∣∣∣ ref + (x− xref )3 6 ∂3U ∂x3 ∣∣∣∣ ref + (x− xref)2 (y − yref) 2 ∂3U ∂x2∂y ∣∣∣∣ ref + (x− xref ) (y − yref )2 2 ∂3U ∂x∂y2 ∣∣∣∣ ref + (y − yref)3 6 ∂3U ∂y3 ∣∣∣∣ ref The referen e point ~xref is hosen to be the triangle enter for ell- entered nite volume s heme and the vertex lo ation for vertex- entered s heme. A onstrained least-squares system is onstru ted to nd the re onstru tion poly- nomial oe ient. The hard onstraint is that the re onstru tion polynomial should preserve the omputed ontrol volume average solution. The rest of the equations whi h are satised only in a least-squares sense are obtained by requiring that the average of the re onstru tion polynomial approximates the ontrol volume average solution for neighboring ontrol volumes. Neighbors are hosen based on their topo- logi al distan e from the element, and entire layers of neighbors is added until the number of neighbor elements in the sten il equals or ex eeds 4 for the se ond order s heme and 16 for the fourth order s heme. Figure 4.5 shows a s hemati drawing of an element k and its neighbor layers. The mean onstraint equation is obtained by requiring onservation of the ontrol 51 Figure 4.5: S hemati drawing of rst, se ond and third neighbor layers of triangular element k. volume average by the re onstru tion polynomial and an be written as Uk = 1 Ak ˆ Ak URk dA ≡ Uk ref + xk ∂U ∂x ∣∣∣∣ k ref + yk ∂U ∂y ∣∣∣∣ k ref + (4.29) x2k 2 ∂2U ∂x2 ∣∣∣∣ k ref + xyk ∂2U ∂x∂y ∣∣∣∣ k ref + y2k 2 ∂2U ∂y2 ∣∣∣∣ k ref + ... where xnymk = 1 Ak ¨ Ak (x− xk ref )n (y − yk ref )m dA 52 Mat hing the ontrol volume average in a neighbor j would require that U j = 1 Aj ˆ Aj URk dA ≡ Uk ref + 1 Aj ˆ Aj (x− xk ref ) ∂U ∂x ∣∣∣∣ ref dA+ (4.30) 1 Aj ˆ Aj (y − yk ref ) ∂U ∂y ∣∣∣∣ ref dA+ 1 Aj ˆ Aj (x− xk ref )2 2 ∂2U ∂x2 ∣∣∣∣ ref dA+ 1 Aj ˆ Aj (x− xk ref ) (y − yk ref ) ∂ 2U ∂x∂y ∣∣∣∣ ref dA+ 1 Aj ˆ Aj (y − yk ref )2 2 ∂2U ∂y2 ∣∣∣∣ ref dA+ ... Using (x− xk ref ) = (x− xj ref )− (xk ref − xj ref ) and (y − yk ref ) = (y − yj ref )− (yk ref − yj ref ) , we get U j = 1 Aj ˆ Aj URk dA ≡ Uk ref + x̂k,j ∂U ∂x ∣∣∣∣ k ref + ŷk,j ∂U ∂y ∣∣∣∣ k ref + (4.31) x̂2k 2 ∂2U ∂x2 ∣∣∣∣ k ref + x̂yk,j ∂2U ∂x∂y ∣∣∣∣ k ref + ŷ2k,j 2 ∂2U ∂y2 ∣∣∣∣ k ref + ... where x̂nymk,j = 1 Aj ˆ Aj ((x− xj ref )− (xk ref − xj ref ))n · ((y − yj ref )− (yk ref − yj ref ))m dA = m∑ r=0 n∑ q=0 m! r! (m− r) ! n! q! (n− q) ! (xj ref − xk ref ) q · (yj ref − yk ref )r xn−qj ref ym−rj ref for simpli ity, the term x, yk ref will be written as x, yk, and x, yj ref will be repla ed 53 with x, yj . The resulting onstrained least-squares system an be written as 1 xk yk x 2 k xyk y 2 k · · · 1 x̂k,1 ŷk,1 x̂2k,1 x̂yk,1 ŷ 2 k,1 · · · 1 x̂k,2 ŷk,2 x̂2k,2 x̂yk,2 ŷ 2 k,2 · · · 1 x̂k,3 ŷk,3 x̂2k,3 x̂yk,3 ŷ 2 k,3 · · · . . . . . . . . . . . . . . . . . . . . . 1 x̂k,N ŷk,N x̂2k,N x̂yk,N ŷ 2 k,N · · · · · · U ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . . k = Uk U1 U2 U3 . . . UN (4.32) This onstrained least squares system an be rewritten as min ∥∥∥[ARec](P̃Rec)− (URec − Uk)∥∥∥ satisfying [BRec] (PRec) = ( Uk ) where [ARec] = x̂k,1 − xk, ŷk,1 − yk, x̂2k,1 − x2k, x̂yk,1 − xyk, ŷ2k,1 − y2k, · · · x̂k,2 − xk, ŷk,2 − yk, x̂2k,2 − x2k, x̂yk,2 − xyk, ŷ2k,2 − y2k, · · · x̂k,3 − xk, ŷk,3 − yk, x̂2k,3 − x2k, x̂yk,3 − xyk, ŷ2k,3 − y2k, · · · . . . . . . . . . . . . . . . . . . x̂k,N − xk, ŷk,N − yk, x̂2k,N − x2k, x̂yk,N − xyk, ŷ2k,N − y2k, · · · · · · [BRec] = [ 1 xk yk x 2 k xyk y 2 k · · · ] and PRec = U ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . . k , P̃Rec = ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . . k , ( URec − Uk ) = U1 − Uk U2 − Uk U3 − Uk . . . UN − Uk 54 The least-squares system is solved to obtain the re onstru tion polynomial oef- ients using singular value de omposition method to nd the pseudo inverse of the least-squares oe ient matrix as follows ARec = UΣV T A†Rec = V Σ †UT ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . . k = x̂k,1 − xk, ŷk,1 − yk, x̂2k,1 − x2k, x̂yk,1 − xyk, ŷ2k,1 − y2k, · · · x̂k,2 − xk, ŷk,2 − yk, x̂2k,2 − x2k, x̂yk,2 − xyk, ŷ2k,2 − y2k, · · · x̂k,3 − xk, ŷk,3 − yk, x̂2k,3 − x2k, x̂yk,3 − xyk, ŷ2k,3 − y2k, · · · . . . . . . . . . . . . . . . . . . x̂k,N − xk, ŷk,N − yk, x̂2k,N − x2k, x̂yk,N − xyk, ŷ2k,N − y2k, · · · · · · † · U1 − Uk U2 − Uk U3 − Uk . . . UN − Uk (4.33) Uk = Uk − [ xk yk x 2 k xyk y 2 k · · · ] · ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . . (4.34) The unlimited fa e ow properties an be found at fa e Gauss points ( xfi , yff ) 55 as follows Ufi = Uk + (xfi − xk) ∂U ∂x ∣∣∣∣ k + (yfi − yk) ∂U ∂y ∣∣∣∣ k + (4.35) (xfi − xk)2 2 ∂2U ∂x2 ∣∣∣∣ k + (xfi − xk) (yfi − yk) ∂2U ∂x∂y ∣∣∣∣ k + (yfi − yk)2 2 ∂2U ∂y2 ∣∣∣∣ k + ... To ensure monotoni ity and solution onvergen e, a limiter is applied in order not to reate a new extremum at the ontrol volume fa e; the limited ow properties at the fa e an be written as Ufi = Uk + φk ( (xfi − xk) ∂U ∂x ∣∣∣∣ k + (yfi − yk) ∂U ∂y ∣∣∣∣ k + (4.36) (xfi − xk)2 2 ∂2U ∂x2 ∣∣∣∣ k + (xfi − xk) (yfi − yk) ∂2U ∂x∂y ∣∣∣∣ k + (yfi − yk)2 2 ∂2U ∂y2 ∣∣∣∣ k + ... ) where φk is the limiter value in the ontrol volume k. ANSLib uses two limiter fun tions to al ulate φ: Venkatakrishnan's limiter [84℄ and a higher-order limiter [56℄ designed not to degrade the a ura y of high-order s hemes for smooth ows. We will derive the fa e properties-mesh dependen y for the limited ase; the unlimited ase an be obtained by simpli ation of the limited ase by setting φk = 1. 56 4.3.2 Mesh dependen e of the fa e ow property re onstru tion The term ∂U fa e /∂M for the limited ase an be obtained using dire t dierentiation of Eq. 4.36 ∂Ufi ∂M = ∂Uk ∂M + ∂φk ∂M ( (xfi − xk) ∂U ∂x ∣∣∣∣ k + (yfi − yk) ∂U ∂y ∣∣∣∣ k + (4.37) (xfi − xk)2 2 ∂2U ∂x2 ∣∣∣∣ k + (xfi − xk) (yfi − yk) ∂2U ∂x∂y ∣∣∣∣ k + (yfi − yk)2 2 ∂2U ∂y2 ∣∣∣∣ k + ... ) + φk ( ∂ (xfi − xk) ∂M ∂U ∂x ∣∣∣∣ k + ∂ (yfi − yk) ∂M ∂U ∂y ∣∣∣∣ k + ∂ (xfi−xk) 2 2 ∂M ∂2U ∂x2 ∣∣∣∣ k + ∂ (xfi − xk) (yfi − yk) ∂M ∂2U ∂x∂y ∣∣∣∣ k + ∂ (yfi−yk) 2 2 ∂M ∂2U ∂y2 ∣∣∣∣ k + ... + φk (xfi − xk) ∂ ∂U∂x ∣∣ k ∂M + (yfi − yk) ∂ ∂U ∂y ∣∣∣ k ∂M + (xfi − xk)2 2 ∂ ∂ 2U ∂x2 ∣∣∣ k ∂M + (xfi − xk) (yfi − yk) ∂ ∂ 2U ∂x∂y ∣∣∣ k ∂M + (yfi − yk)2 2 ∂ ∂ 2U ∂y2 ∣∣∣ k ∂M + ... The unlimited fa e ow property-mesh dependen y an be obtained by setting φk = 1 and ∂φk/∂M = 0 in Eq. 4.37. In Eq. 4.37, the term ∂φk/∂M an be obtained by dierentiating the limiter 57 expression, and the terms ∂ ∂ nU ∂xm∂yn−m ∣∣∣ k /∂M an be evaluated as follows: [ARec] { P̃Rec } = { URec } ∂ [ARec] ∂M { P̃Rec } + [ARec] ∂ { P̃Rec } ∂M = ∂ { URec } ∂M = {0} ∴ ∂ { P̃Rec } ∂M = − [ARec]† { ∂ [ARec] ∂M { P̃Rec }} ∂ { P̃Rec } ∂M = − [ARec]† ∂ [ARec] ∂M [ARec] † {URec}(4.38) and ∵ Uk =Uk − [ xk yk x 2 k xyk y 2 k · · · ] · { P̃Rec } ∴ ∂Uk ∂M =− ∂ ∂M [ xk yk x 2 k xyk y 2 k · · · ] · { P̃Rec } − (4.39) [ xk yk x 2 k xyk y 2 k · · · ] · ∂ { P̃Rec } ∂M where ∂ [ARec] ∂M = ∂x̂k,1 ∂M ∂ŷk,1 ∂M ∂x̂2k,1 ∂M ∂x̂yk,1 ∂M ∂ŷ2k,1 ∂M · · · ∂x̂k,2 ∂M ∂ŷk,2 ∂M ∂x̂2k,2 ∂M ∂x̂yk,2 ∂M ∂ŷ2k,2 ∂M · · · ∂x̂k,3 ∂M ∂ŷk,3 ∂M ∂x̂2k,3 ∂M ∂x̂yk,3 ∂M ∂ŷ2k,3 ∂M · · · . . . . . . . . . . . . . . . . . . ∂x̂k,N ∂M ∂ŷk,N ∂M ∂x̂2k,N ∂M ∂x̂yk,N ∂M ∂ŷ2k,N ∂M · · · · · · ∂ { P̃Rec } ∂M = ∂ ∂M ∂U ∂x ∂U ∂y 1 2 ∂2U ∂x2 ∂2U ∂x∂y 1 2 ∂2U ∂y2 . . . 58 and ∂ ∂M {(xfi − xk)m (yfi − yk)n} = m (xfi − xk)m−1 (yfi − yk)n ( ∂xfi ∂M − ∂xk ∂M ) + n (xfi − xk)m (yfi − yk)n−1 ( ∂yfi ∂M − ∂yk ∂M ) ∂x̂nymk,j ∂M = m∑ r=0 n∑ q=0 { m! r! (m− r) ! n! q! (n− q) ! q (xj ref − xk ref ) q−1 · (yj ref − yk ref ) · xn−qj ref y m−r j ref · ( ∂xj ref ∂M − ∂xk ref ∂M ) + m! r! (m− r) ! n! q! (n− q) ! r (xj ref − xk ref ) q (yj ref − yk ref )r−1 · xn−qj ref y m−r j ref · ( ∂yj ref ∂M − ∂yk ref ∂M ) 59 Chapter 5 GRADIENT VALIDATION In this se tion, the omparison of the gradient a ura y evaluated using se ond and fourth order s hemes is arried out. Both se ond and fourth order gradients are evaluated using the adjoint approa h and ompared to their orresponding nite dieren e gradient. We present three test ases: subsoni , non-limited transoni and limited transoni test ases. In all test ases we use 18 design points to parametrize the airfoil geometry. 5.1 Subsoni Test Case In this test ase, the evaluation of the lift oe ient gradient with respe t to the airfoil geometri design ontrol points is presented for both se ond and fourth order s hemes, in luding a omparison with a nite dieren e gradient al ulated using the same order of a ura y in the ow solver. The airfoil used in this test ase is NACA 0012 at subsoni onditions M = 0.5 and α = 2o. Figure 5.1 shows a representative omparison between the eld of pressure sensitivity with respe t to one of the geometry design points omputed using nite dieren es and the solution sensitivity al ulation of Eq. 4.6 for both se ond and fourth order a urate ompu- tations. Agreement is ex ellent in all ases; this is also true for the other design ontrol points. The ex ellent mat hing of the pressure sensitivity when omparing the se ond order and the fourth order results indi ates that the two s hemes will give similar gradient ve tors and similar optimization des ent dire tions for subsoni ow optimization. Table 5.1 shows quantitatively an ex ellent mat hing between the obje tive fun - tion gradient magnitude (less than half a per ent dieren e) and dire tion (less than a half a degree dieren e) when omparing nite dieren es and the adjoint ap- proa h of Eq. 4.8 for both se ond and fourth order s hemes. Figure 5.1 shows good 60 (a) Se ond order sensitivity al ulation (b) Se ond order nite dieren e ( ) Fourth order sensitivity al ulation (d) Fourth order nite dieren e Figure 5.1: The pressure sensitivity with respe t to one of the design ontrol points omputed for subsoni ow over NACA 0012, omparing the sensitivity al ulation of Eq. 4.6 with nite dieren e results. agreement between the obje tive fun tion gradient omponents ∂CL/∂yi, where yi are the y lo ations of the design ontrol points, omputed using adjoint and nite dieren e approa hes for both se ond and fourth order omputations. The maxi- mum error, normalized by the gradient magnitude, is only 0.005 whi h gives a high level of a ura y for the omputed gradient. Taken together, these results imply that the order of dis retization error has little ee t on the omputed gradient ve tor for subsoni ow. 5.2 Transoni Test Case with No Limiter In this test ase, the sensitivity analysis is arried out for NACA 0012 airfoil with M = 0.8 and α = 2◦. The drag oe ient gradient is evaluated without using limiters in the CFD simulation, so overshoot/undershoot at the sho k lo ation is expe ted. Figure 5.3 shows good agreement of the pressure sensitivity omputed using Eq. 4.6 61 Se ond order Fourth order Adjoint FD Adjoint FD Gradient ve tor magnitude 18.605 18.561 18.575 18.655 Angle with se ond order adjoint 0 ◦ 0.247 ◦ 1.117 ◦ 1.330 ◦ Angle with fourth order adjoint 1.117 ◦ 0.917 ◦ 0 ◦ 0.249 ◦ Table 5.1: The magnitude of se ond and fourth order CL gradients and angles be- tween the evaluated gradients for NACA 0012 in subsoni ow. Figure 5.2: CL gradient error in se ond and fourth order s hemes with respe t to the design points normalized by gradient magnitude. 62 (a) Se ond order sensitivity al ulation (b) Se ond order nite dieren e ( ) Fourth order sensitivity al ulation (d) Fourth order nite dieren e Figure 5.3: The pressure sensitivity with respe t to one of the design ontrol points omputed for an unlimited transoni ow around NACA 0012, omparing the sensi- tivity al ulation of Eq. 4.6 with nite dieren e results. and nite dieren e; the gure also shows that for transoni ows, the pressure sensitivity omputed using se ond order and fourth order a urate adjoint s heme are dierent, espe ially near the sho k lo ation. Table 5.2 shows again the ex ellent mat hing between the omputed adjoint and nite dieren e gradient with a dire tion dieren e less than a degree and nearly identi al magnitude. Figure 5.4 shows that for unlimited transoni ow, both se - ond and fourth order adjoint gradients are an ex ellent mat h to the orresponding nite dieren e gradient, with the se ond order s hemes mat hing more losely than the fourth order s heme. Comparison of the se ond and fourth order gradient ve tors shows little dieren e between them (about 4% in magnitude and 2◦ in dire tion). Again, we expe t that the se ond and fourth order s hemes will give similar opti- mization sear h dire tions. 63 Se ond order Fourth order Adjoint FD Adjoint FD Gradient ve tor magnitude 1.895 1.895 1.814 1.813 Angle with se ond order adjoint 0 ◦ 0.252 ◦ 2.236 ◦ 1.997 ◦ Angle with fourth order adjoint 2.236 ◦ 2.354 ◦ 0 ◦ 0.518 ◦ Table 5.2: The magnitude of se ond and fourth order CD gradients and angles between the evaluated gradients for NACA 0012 in an unlimited transoni ow. Figure 5.4: CD gradient error in se ond and fourth order s hemes with respe t to the design points, normalized by gradient magnitude. 64 (a) Se ond order sensitivity al ulation (b) Se ond order nite dieren e ( ) Fourth order sensitivity al ulation (d) Fourth order nite dieren e Figure 5.5: The pressure sensitivity with respe t to one of the design ontrol points omputed for an unlimited transoni ow around NACA 0012. 5.3 Transoni Test Case with Limiter In this test ase, the impa t of using a limiter in the CFD simulation on the a ura y of the omputed gradient for se ond order and fourth order s hemes is studied. Two dierent limiters are used, the Venkatakrishnan limiter [84℄ and the higher order limiter of Mi halak and Ollivier-Goo h [56℄. Figure 5.5 shows very good mat hing of the se ond order pressure sensitivity omputed using adjoint and nite dieren e te hniques with the use of Venkatakrishnan limiter. Mat hing is less good between the pressure sensitivity and nite dieren e results for the fourth order a urate s heme; this lower level of pressure sensitivity mat hing will lead to less a urate gradient values when using the limited fourth order s heme. Table 5.3 shows that with the use of the Venkatakrishnan limiter, the se ond order gradient magnitude is a very good mat h with the orresponding nite dieren e gradient; the larger error in gradient value observed for the fourth order s heme is omparable to the dieren e in magnitude between se ond and fourth order s hemes 65 Se ond order Fourth order Adjoint FD Adjoint FD Gradient ve tor magnitude 1.897 1.895 1.887 1.812 Angle with se ond order adjoint 0 ◦ 3.612 ◦ 2.292 ◦ 3.334 ◦ Angle with fourth order adjoint 2.292 ◦ 5.275 ◦ 0 ◦ 4.873 ◦ Table 5.3: Magnitudes of se ond and fourth order CD gradients and angles between the evaluated gradients for NACA 0012 in Venkatakrishnan limited transoni ow Se ond order Fourth order Adjoint FD Adjoint FD Gradient ve tor magnitude 1.676 1.753 1.576 1.672 Angle with se ond order adjoint 0 ◦ 2.845 ◦ 30.1 ◦ 18.8 ◦ Angle with fourth order adjoint 30.1 ◦ 21.05 ◦ 0 ◦ 18.9 ◦ Table 5.4: The magnitude of se ond and fourth order CD gradients and angles be- tween the evaluated gradients for NACA 0012 using higher order limiter in transoni ow. for the unlimited transoni ase. Also, the dieren e in nite dieren e and adjoint gradient dire tion grows to several degrees with the use of Venkatakrishnan limiter. Table 5.4 shows the same behavior with the higher order limiter of C. Mi halak and C. Ollivier-Goo h. It shows also that the error in the gradient magnitude is larger ompared to the Venkatakrishnan results. The mat hing in adjoint and nite dieren e gradient dire tions is good for the se ond order s heme but it is not for the fourth order s heme with the use of higher order limiter. It is well known that the use of limiters auses onvergen e problem due to the non-dierentiability of the limiting pro edure [55℄. The limiter ae ts both the right and left hand sides of Eq. 1.6. We spe ulate that the poor mat hing between nite dieren e and adjoint gradients for the high-order limiter is related to the non-dierentiability of the limiting pro edure and that this ee t will vary in strength from one limiting pro edure to another. To redu e the inuen e of ∂Ri/∂D in ontrol volume i when al ulating the adjoint gradient, the fourth order Ja obian was modied numeri ally by making the non zero stru ture of the fourth order Ja obian matrix the same as the non zero stru ture of the se ond order Ja obian, and dropping the rest of values in the fourth order Ja obian matrix. The right hand side is still onstru ted with fourth order 66 Mod 4th order Limiter Venkat. HO Gradient s heme Adj. FD Adj. FD Gradient ve tor magnitude 1.735 1.812 1.966 1.672 Angle with mod 4th order Venkat. adjoint 0 ◦ 2.845 ◦ 6.247 ◦ 3.338 ◦ Angle with mod 4th order HO. adjoint 6.247 ◦ 7.256 ◦ 0 ◦ 7.207 ◦ Table 5.5: The magnitudes and angles between the evaluated modied fourth order CD gradients using adjoint, and nite dieren e for NACA 0012 using Venkatakr- ishnan and higher order limiters in transoni ow. a ura y. The above modi ation doesn't ae t the a ura y of the CFD simulation be ause the right hand side remains fourth order a urate. The omputed gradient using this approa h is presented in Table 5.5 and shows a redu tion in the error in the fourth order omputed adjoint gradient, espe ially in the dire tion of the gradient. The limiter used from now on is the Venkatakrishnan limiter as it produ es an adjoint omputed gradient with better mat hing to the nite dieren e gradient ompared with the high order limiter. 5.4 Sensitivity of nite dieren e gradient to perturbation magnitude In the previous se tions, I have used perturbation amplitude ǫ = 10−6 and using the entral dieren e formula to evaluate the nite dieren e gradient. In this se tion, the sensitivity of the evaluated gradient with respe t to ǫ is presented for the hardest ase, transoni ow with limiter. The limiter used in the transoni simulation is Venkatakrishnan's limiter. Table 5.6 shows the norm of the drag gradient for dierent values of perturbation amplitude ǫ; it shows also that de reasing ǫ to be less than 10−6 will not pra ti ally hange the value of the omputed gradient norm. Therefore, ǫ = 10−6 is hosen. 67 Figure 5.6: The normalized CD gradient error in se ond, fourth, and modied fourth order s hemes with respe t to the design points in a limited transoni ow (Venkatakrishnan limiter) Figure 5.7: The normalized CD gradient error in se ond, fourth, and modied fourth order s hemes with respe t to the design points in a limited transoni ow (higher order limiter). 68 ǫ Se ond order ∥∥∥∂CD∂D ∥∥∥ Fourth order ∥∥∥∂CD∂D ∥∥∥ 10−3 1.90237 1.86814 10−6 1.90275 1.87105 10−9 1.09278 1.87114 Table 5.6: Finite dieren e drag gradient sensitivity with respe t to perturbation amplitude ǫ 69 Chapter 6 GRADIENT BASED OPTIMIZATION TEST CASES In this se tion we present four optimization test ases; in all of them we ompare the optimal shape resulting from using se ond and fourth order s hemes. The rst two ases are inverse design problems, one subsoni and the other transoni . In both test ases, a target pressure distribution is obtained using CFD simulation of a parametrized NACA 2412 airfoil and the starting geometry is a NACA 0012 airfoil. The optimizer will try to nd the geometry whose surfa e pressure distribution mat hes the target pressure distribution. Those two test ases are important, as the design spa e has only one solution point at whi h the resulting pressure distribution from the optimal geometry will mat h the target pressure distribution. The third test ase is a transoni drag minimization with no lift onstraint starting from the RAE 2822 airfoil. The obje tive of this test ase is to minimize CD at M = 0.73 and angle of atta k = 2◦. In this test ase a strong sho k wave is formed near mid hord of the initial airfoil geometry and we are seeking a sho k free geometry, or at least a geometry that produ es a mu h weaker sho k wave. Geometri onstraints are applied so the airfoil thi kness will be positive all the way along the airfoil se tion. The fourth test ase repeats test ase three but adds the lift oe ient as a onstraint; in this ase, we ompare the resulting optimal shape with the optimization results of Brezillon and Gauger [11℄ . 6.1 Subsoni Inverse Design In this test ase, the target pressure distribution is obtained for the parametrized NACA 2412 at a subsoni ondition, M = 0.5 and α = 2◦ , using se ond and fourth order CFD simulations. The starting geometry is the parametrized NACA 0012. 70 The obje tive fun tion to be minimized is F = ˛ (PT − Pi)2 dS (6.1) The above obje tive fun tion and its gradient an be expressed in dis rete form as F = ∑ (PT − Pi)2 ws (6.2) dF dxd = ∑ 2 (PT − Pi) (−∂Pi ∂xd ) ws + (PT − Pi)2 ∂ws ∂xd (6.3) where ws is the ar length asso iated with the surfa e Gauss point. Figure 6.1 shows the target pressure distribution of the NACA 2412, the initial pressure distribution of NACA 0012, and the optimized airfoil pressure distribution obtained by the se ond order and fourth order s hemes; both s hemes su essfully rea hed the target pressure distribution. Figures 6.2 and 6.3 show the onvergen e history and gradient norm history for both s hemes. Both s hemes took about the same number of optimization iterations (28 for se ond order and 32 for fourth order) to drop the obje tive fun tion value by eight orders of magnitude. Figures 6.4 and 6.5 show how lose the optimized NACA 0012 is to the NACA 2412 using both se ond order and fourth order s hemes. The error between the target geometry and the optimized prole is larger in the fourth order s heme due to ina ura y of the pressure interpolation s heme used to evaluate the obje tive fun tion; nevertheless the resulting geometry is an ex ellent mat h with the NACA 2412, with an error of less 0.1% of the surfa e movement. 6.2 Transoni Inverse Design In this test ase, the target pressure is obtained using CFD simulation of the ow over a NACA 2412 airfoil at transoni onditions, M = 0.73, α = 2◦. The obje - tive fun tion to be minimized is again the integration of the square of the pressure dieren e between the target pressure and the optimized pressure as expressed in Eq. 6.1. The fourth order optimization is based on the modied fourth order gra- dient evaluation strategy. Figure 6.6 shows the initial, target, and the optimized 71 (a) Se ond order (b) Fourth order Figure 6.1: Subsoni NACA 2412 inverse design pressure distributions for the initial, target, and optimized airfoil proles 0 5 10 15 20 25 30 35 10−12 10−10 10−8 10−6 10−4 10−2 Iteration F Fourth order Second order Figure 6.2: Se ond and fourth order optimization onvergen e history. 72 0 5 10 15 20 25 30 35 10−6 10−5 10−4 10−3 10−2 10−1 Iteration Fi rs t− or de r o pt im al ity Fourth order Second order Figure 6.3: Gradient norm history Figure 6.4: Subsoni inverse design optimal airfoil shapes 73 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 x 10−5 Upper surface Dy x/c D y 4th order 2nd order 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 x 10−5 Lower surface Dy x/c D y 4th order 2nd order (a) Upper surfa e (b) Lower surfa e Figure 6.5: The dieren e between the target prole and the optimized proles, se ond and fourth order pressure distribution. The onvergen e history and gradient norm history are shown in Figs. 6.7 and 6.8; the fourth order s heme is slower to onverge ompared to the se ond order s heme, due to the use of larger number of airfoil surfa e Gauss quadrature points used in the fourth order omputations (double the number used for the se ond order s heme) whi h makes the minimizer slower to onverge. The obje tive fun tion dropped only three order of magnitudes before onvergen e stall. This stall is due to the high non-linearity in the target pressure distribution be ause of the existen e of a strong sho k wave in it. The noise in the pressure sensitivity generated by the sho k wave in the target pressure distribution doesn't allow further onvergen e; however, the gradient magnitude dropped four order of magnitudes from its initial value. Figures 6.9 and 6.10 show that the optimal shapes for the two s hemes dier by less than 10−4 of the hord length on the lower surfa e and less than 10−3 hord length on the upper surfa e (where the strong sho k wave exists); both optimal shapes are in good agreement with the NACA 2412: the maximum deviation is about 5% of maximum surfa e movement. 6.3 Drag Minimization without Lift Constraint In this test ase, minimization of the drag oe ient will be arried out with no lift onstraint applied. The airfoil to be optimized is an RAE 2822 at transoni 74 (a) Se ond order (b) Fourth order Figure 6.6: Subsoni NACA 2412 inverse design pressure distributions for the initial, target, and optimized airfoil proles Figure 6.7: Se ond and fourth order optimization onvergen e history. 75 0 5 10 15 20 25 30 35 40 45 10−4 10−3 10−2 10−1 100 G ra de N or m iteration 2nd order 4th order Figure 6.8: Gradient norm history (a) Upper surfa e (b) Lower surfa e Figure 6.9: The dieren e between the target prole and the optimized proles, se ond and fourth order 76 Figure 6.10: Transoni inverse design optimal airfoil shapes. CL CL optimized airfoil CD CD optimized airfoil Se ond order 0.865 0.765 0.0081 0.0046 6 Fourth order 0.849 0.759 0.0099 0.0047 Table 6.1: Aerodynami oe ients of original and optimized RAE 2822 airfoil at transoni onditions. onditions: M = 0.73 and α = 2◦. Fig. 6.11 shows the initial solution with a strong sho k wave standing at 70% hord. Geometri onstraints are added to insure that there is no interse tion between the airfoil upper and lower surfa es along the airfoil. Figure 6.12 shows the optimal shapes and the optimized pressure elds resulting from using se ond and fourth order s hemes. Figures 6.12 and 6.13 show that the dieren e between the se ond and fourth order optimal proles is notable on both upper and lower surfa es. The surfa e pressure distribution of the original RAE 2822 airfoil and of the optimized airfoil for se ond and fourth order s hemes is shown in Fig. 6.14. Both s hemes su essfully produ ed a similar sho k free pressure distribution but with dierent proles. The original and optimized CD are shown in Table 6.1. A drag redu tion of about 50% is a hieved by both s hemes. The value of the CD of the se ond order optimized airfoil is obtained from fourth order a urate CFD simulation over the optimized se ond order airfoil prole to eliminate the dieren es in dis retization error in the nal results. Figure 6.16 shows the onvergen e history of the optimization. Both s hemes rea hed their optimal solution in about the same number of optimization iterations. 77 Figure 6.11: Pressure ontours of RAE 2822 at Ma h 0.73 and angle of atta k 2 Figure 6.12: Optimized RAE 2822. 78 (a) Upper surfa e (b) Lower surfa e Figure 6.13: Optimization surfa e displa ements of the original RAE2822 surfa es. Figure 6.14: Pressure distribution omparison for the original and the optimized geometries. 79 0 5 10 15 20 25 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 −grade −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−− +grad C D BFGS optimal profile 0 5 10 15 20 25 0 0.002 0.004 0.006 0.008 0.01 0.012 −grade −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−− +grad C D BFGS optimal profile (a) Se ond order response surfa e (b) Fourth order response surfa e Figure 6.15: Obje tive fun tion response surfa e along positive and negative gradient dire tions entered at the optimal prole of RAE 2822 transoni drag minimization without lift onstraint. The number of CFD simulations for the se ond order s heme is 50 ompared to 45 for the fourth order s heme, yet the overall omputational ost of the se ond order s heme is 40% less than the fourth order omputations. The gradient of the obje tive fun tion is not zero at the optimal point; to understand the reasons for that, I have plotted the obje tive fun tion values along positive and negative gradient dire tion at the optimal solution point as shown in Fig 6.15. The se ond order response surfa e shows that the obtained optimal solution is a possible lo al optimum and the shape of the obje tive fun tion shows that the gradient omputed using nite dieren e will not be zero at the optimal point. For the fourth order omputations, the gure shows that the optimizer rea hed a near optimal solution. While a slightly lower value is available,the small dieren e between the obje tive fun tion value at the obtained optimal point and the minimum value along the negative gradient dire tion (on the order of the dis retization error) may ause a line sear h failure due to insu ient de rease in the obje tive fun tion. 6.4 Drag Minimization with Lift Constraint This test ase represents a typi al optimization task required in aerospa e industry, as it is required to minimize the drag while the lift is un hanged. Therefore the lift oe ient CL will be an aerodynami onstraint. The original RAE 2822 geometry 80 Figure 6.16: Se ond and fourth order optimization onvergen e history. will be used as a starting shape; the obje tive fun tion to be minimized is F = CD + 10 · (CL − CLC )2 (6.4) where CLC is the original lift oe ient of RAE 2822 atM = 0.73 and α = 2 ◦ . In this test ase both geometri and aerodynami onstraints are applied to avoid getting a non-feasible airfoil geometry, while not ae ting airfoil lift oe ient. The nonlinear onstraint, the lift oe ient, is added as a penalty term in the obje tive fun tion as shown in Eq. 6.4. Figure 6.17 show the optimized airfoil's pressure distribution for both se ond and fourth order s hemes. Both se ond and fourth order s hemes travel 99% of the way towards their optimal obje tive fun tion after 7 iterations; the se ond order s heme then be omes very slow to rea h its optimal value. The se ond order s heme osts 112 CFD simulations to rea h its optimum ompared to 49 CFD simulations for the fourth order. Figure 6.20 shows a noti eable dieren e between the optimized proles, in luding a notably larger nose radius and less aft hamber for the fourth order s heme. A sho k free geometry is obtained with weak ompression waves on the airfoil upper surfa e. Drag is redu ed by about 50%, as shown in Table 6.2, while the lift oe ient is about 1% higher than its original value. 81 CL CL optimized airfoil CD CD optimized airfoil Se ond order 0.865 0.871 0.0081 0.0047 7 Fourth order 0.849 0.853 0.0099 0.0051 Table 6.2: Aerodynami oe ients of original and optimized RAE 2822 airfoil at transoni onditions. Figure 6.17: Surfa e pressure distribution, of the initial and optimized geometries Figure 6.18: Se ond and fourth order optimization onvergen e history. 82 Figure 6.19: Optimal shapes omparison: se ond order, fourth order, and optimized prole by Brezillon and Gauger ompared with the original RAE2822. (a) Upper surfa e (b) Lower surfa e Figure 6.20: Se ond and fourth order optimized prole surfa e displa ements from the initial shape. 83 0 5 10 15 20 25 0 1 2 3 4 5 6 7 8 −grade −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−− +grad F BFGS Optimal profile 0 5 10 15 20 25 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 −grade −−−−−−−−−−−−−−−−BFGS Optimal −−−−−−−−−−−−−−−−− +grad F BFGS Optimal Profile (a) Se ond order response surfa e (b) Fourth order response surfa e Figure 6.21: Obje tive fun tion response surfa e along positive and negative gradient dire tions entered at the optimal prole of RAE 2822 transoni drag minimization without lift onstraint. At the optimization breakdown point (optimal solution found), the gradient is not zero, possibly due to the existen e of the sho k wave whi h auses high non-linearity in the obje tive fun tion and lift onstraint behavior. Figure 6.21_shows plots of the obje tive fun tion response along negative and positive gradient dire tion. It shows also that the obtained optimal solutions are lo al minima along the gradient dire tion, and shows also that the gradient omputed by nite dieren e is not zero. It is also worth noting that the gradient of the obje tive fun tion at the breakdown point is dominated by the gradient of the lift onstraint. Figures 6.19 and 6.20 show a omparison between the optimized airfoil proles based on se ond and fourth order omputations. Dieren es in shape are espe ially noti eable near the leading edge, and in the reexed part of the lower surfa e near the trailing edge. The se ond order optimized prole Cd is four drag ounts less than the fourth order prole drag oe ient, nearly 10%. The same problem was studied by Brezillon and Gauger [11℄; they used the se ond order solver MEGAFLOW of the German Aerospa e Center and obtained their optimal prole by ontrolling the amber line of the RAE 2822 via 20 ontrol points. Comparing the optimal pressure distribution of our se ond order omputations shown in Fig 6.17 and that of Brezil- lon, Fig 6.22, shows that both solvers tend to a elerate the ow near the leading edge upper surfa e more than the original airfoil shape, followed by a gradual de el- eration till the trailing edge; on the other hand, the surfa e pressure on the lower 84 Figure 6.22: initial and optimal pressure distribution obtained by Brezillon and Gauger [11℄ (presented with permission) airfoil surfa e hanges only slightly from the original surfa e pressure distribution of the RAE2822. This similarity in the optimized pressure distribution suggests that whatever geometry parametrization te hnique is used, the surfa e pressure hanges will be qualitatively similar. The lift penalty fa tor used in this test ase is sele ted su h that violating the lift onstraint by one lift ount (10−2) auses 10 drag ounts (10−3) in rease in the obje tive fun tion value. Table 6.3 shows the ee t of the lift onstraint penalty weighting fa tor on the optimal obje tive fun tion value, as well as the optimal drag and lift oe ients; these results are based on se ond order omputations. 85 Lift penalty weight F opt CD opt CL opt Lift penalty 5 0.0052 0.0047 0.873 3.2 × 10−4 10 0.0051 0.0047 0.871 3.6 × 10−4 20 0.0082 0.0080 0.868 1.8 × 10−4 Table 6.3: Lift penalty weight ee t on Drag minimization of RAE 2822 with lift onstraint 6.5 Mesh Renement Study of a Drag Minimization with Lift Constraint In the last test ase, Table 6.2 shows that the optimal value of CD is dierent when omparing se ond and fourth order s heme results. This subse tion examines the impa t of mesh renement on the optimal value of CD. Three mesh grid sizes are used, with 5000, 11000, and 15000 triangles, respe tively. In all test ases, CLC = 0.84 is used as a lift onstraint value for both se ond and fourth order omputations. Figure 6.23 ompares the dieren e in the optimal CD value for se ond and fourth order s hemes with mesh renement; the plotted value of se ond order CD optimal are obtained using fourth order CFD simulation on the optimal shape obtained using se ond order s heme. 86 Figure 6.23: Optimal CD value with mesh renement The gure shows that the two s hemes tend to rea h the same value of the optimal CD with mesh renement with only one drag ount dieren e (whi h is on order of the dis retization error). This behavior is expe ted as transoni CFD simulations using se ond and fourth order s hemes for the same geometry produ e almost the same surfa e pressure distribution and hen e obje tive fun tion value. The root mean square of dieren e between the se ond and fourth order optimized airfoil proles of drops from 3 · 10−5 on the oarsest mesh to 2 · 10−5 on the nest grid. 87 Chapter 7 PARTICLE SWARM OPTIMIZATION AND A NEW HYBRID OPTIMIZATION METHOD 7.1 Introdu tion to Gradient Free Optimization Gradient based optimization results depend on the optimization starting point. To nd a global optimal solution, gradient independent methods are an obvious an- didate. Perhaps the best-known non-gradient based optimization te hnique is the geneti algorithm [30℄, a biologi ally-inspired approa h whi h an nd the global optimal point of an obje tive fun tion with multiple lo al optima. In geneti al- gorithm optimization, a randomly generated population of andidate solutions are generated; their tness is determined by the value of the obje tive fun tion [29℄. So- lution points with the best tness are sele ted to be parents of the next generation. A new generation of ospring is generated by rossover and mutation of the parents' geneti hara teristi s (that is, the parents' solution point lo ations in the design spa e). This pro edure is repeated until the globally optimal solution is obtained. The number of individuals in ea h generation should not be less than the number of design variables, otherwise poor optimization onvergen e will result. Figure 7.1 shows a ow hart of a simple geneti algorithm optimizer. The dependen y of the population size on the number of design variables and the large number of genera- tions required to onverge to the global optimum ombine to make geneti algorithm impra ti ally expensive for large aerodynami optimization problems. Geneti algo- 88 rithm optimization has been used for aerodynami optimization problems by many resear hers. A geneti algorithm was used for transoni wing drag minimization with a lift onstraint by Gage and Kroo [29℄. Anderson used GA for wing aerody- nami shape optimization with stru tural onstraints [1℄. Jang and Lee used GA to maximize lift to drag ratio of an airfoil using the Euler ow model with the NACA 0012 as a starting geometry [44℄. Oyama et al applied a geneti algorithm with a Navier-Stokes solver for transoni wing optimization [70℄. They also explored the use of fra tal analysis in GA aerodynami optimization [71℄. Parti le swarm optimization is a random sear h te hnique rst introdu ed in 1995 as a te hnique that simulate the behavior of a herd of predators hunting for food [46℄. This te hnique has proved to be faster than geneti algorithm to rea h the global optimal by many resear hers (for example, [34, 59℄); swarm intelligen e has the advantage of not being highly sensitive to problem dimensions [79℄, whi h is not true for evolutionary optimization like geneti algorithm. The parti le swarm optimization method (des ribed in Se tion 7.2) an be made more e ient by hy- bridizing it with sequential quadrati programming, as des ribed in Se tion 7.3. The ee tiveness of this hybrid approa h for un onstrained and onstrained aerodynami optimization is demonstrated in Se tion 7.4. 7.2 Swarm Intelligen e Parti le swarm optimization is a random sear h method that sear hes for the global optimal solution of an obje tive fun tion. It uses a number of swarm parti les that s ans the design spa e looking for the global optimal solution. Ea h swarm parti- le represents a point in the design spa e and it moves through the design spa e a ording to a simple formula. A parti le's velo ity depends on its inertia, its own experien e, and so ial experien e gained by ommuni ating with other swarm par- ti les. The parti le swarm te hnique keeps a re ord of lo ation of the best obje tive fun tion value for ea h parti le ( −→pi ) as well the global best lo ation (−→g ), whi h is the best lo ation among the −→pi 's. Algorithm 7.1 is the general pseudo ode of the parti le swarm te hnique. The velo ity of the parti le in an optimization iteration depends on Its velo ity in the previous iteration (momentum part) 89 Figure 7.1: Geneti Algorithm optimization 90 Algorithm 7.1 Parti le swarm pseudo ode For ea h parti le i in the swarm Initialize parti le position −→xi Set its best position to be its initial position: −→pi = −→xi Initialize its velo ity −→ Vi = ~0 End Set the global best position: −→g = −→pi where F (−→pi ) < F (−→pj ) for all i 6= j Do For ea h parti le i Cal ulate its velo ity −→ Vi Update its position −→xi Cal ulate its obje tive fun tion value (tness) F (−→xi) If the urrent tness value is better than its best tness (F (−→xi) < F (−→pi )) Set the urrent parti le position to be its best position: −→pi = −→xi End if End for Sele t the best parti le obje tive fun tion and position to be the global best if it is better than the stored value While max iteration ount not rea hed or onvergen e riterion not met. 91 Its distan e from its known best value −→pi ( ognitive part) Its distan e from the swarm global best value −→g (so ial part) Figure 7.2 shows s hemati ally how the above three fa tors ae t the parti le ve- lo ity, the velo ity of ea h parti le and its position update are obtained using the following formulas −→ Vi(k + 1) = w −→ Vi(k) + c1 [ −→r1 ] · {−→pi −−→xi(k)} + c2 [−→r2 ] · {−→g −−→xi(k)} (7.1) −→xi(k + 1) = −→xi(k) +−→Vi(k + 1), i = 1..N (7.2) where w = 0.73 is the inertial weight fa tor, c1 = c2 = 1.49 are the ognitive and so ial a eleration fa tors respe tively; these onstants are sele ted a ording to Cler 's onstri tion models [16℄. −→r1 and −→r2 are two random ve tors of size N whose elements are independent and uniformly randomly sele ted from the uniform distribution on the interval [0, 1]. The · operation is a omponent by omponent multipli ation not a dot produ t and this is true from this point through the rest of the hapter. Parti le swarm optimization suers from solution stagnation on e the parti les have onverged to a high quality lo ally optimal solution whi h is not a true globally optimal point [23℄. To avoid this problem, a restart strategy is suggested when premature onvergen e is dete ted [9℄. At ea h restart, swarm parti les start from new positions in the design spa e and onverge to an optimal solution; restarting is arried out m times, and the optimal solution is the best of all global bests. This te hnique is alled the multi-start parti le swarm optimization (MPSO). A drawba k of this me hanism is that restarting may ause repetition of the sear h omputations. Evers and Ben Ghalia suggested a more e ient me hanism to es ape from high quality lo al wells. They re-s atter the parti les around the onvergen e point in a region with radius smaller than the design spa e norm, and large enough some at least some parti les to be outside the basin of attra tion of the lo al optimum [24, 23℄; their s heme is alled the regrouped parti le swarm optimization (RegPSO). The main three elements in their s heme are Prematurity dete tion. 92 Figure 7.2: Velo ity omponents and position update of a parti le Redening upper and lower boundaries for de ision (design) variables. Regrouping the parti les by s attering them on a sear h domain dened by the new upper and lower boundaries. In the following subse tions, the above three elements of the RegPSO strategy will be presented; interested readers are en ouraged to onsult Evers and Ben Ghalia [23℄ for more detail. All onstants presented in the next subse tions were found by Evers and Ben Ghalia using numeri al optimization experiments on standard ben hmark optimization problems. 7.2.1 Premature onvergen e dete tion Convergen e is dete ted when all parti les onverge to an optimal solution and start to lose their momentum; this loss of momentum prevents es aping from the neigh- borhood of a high quality lo al optimal solution. Premature onvergen e an be dete ted when the greatest distan e between a swarm parti le and the global best point −→g falls below a ertain threshold. This was suggested by Van den Bergh and adopted in Ever's work with su ess [9, 23℄. 93 For multidimensional optimization problem of size n andm parti les, the de ision variables −→x = [ x1 x2 x3 · · · xn ]T have initial upper and lower limits −→ x0U =[ x1U x2U x3U · · · xnU ]T and −→ x0L = [ x1L x2L x3L · · · xnL ]T . First, the range is dened as −−−→ range (Ω) = −→xU −−→xL (7.3) Also the normalized swarm radius an be found a ording to δ norm = maxi=1...m ‖−→xi (k)−−→g (k)‖ ‖−−−→range (Ω)‖ (7.4) Convergen e is dete ted when the normalized swarm radius δ norm drops below ε = 1.1 · 10−4. 7.2.2 Redening de ision variables range On e onvergen e is dete ted, the upper and lower limits of the de ision variables for the regrouped swarm are omputed. The initial range of parti le j an be found using its upper limit x0jU and lower limit x 0 jLas rangej ( Ω0 ) = x0jU − x0jL (7.5) Upper and lower limits of the de ision variables are hanged when onvergen e is dete ted in order to s atter swarm parti les in a ball surrounding the onverged solution. Changing the range of ea h de ision variable j in a regrouping stage r is done a ording to the following formula [23℄: rangej (Ω r) = max ( rangej ( Ω0 ) , 1.2 × 104 max i=1...m ∣∣∣∣−−→xr−1i,j −−−→gr−1j ∣∣∣∣) (7.6) When the urrent range of design variables rangej is su iently small, the optimiza- tion problem is onsidered onverged. 7.2.3 Regrouping and position lipping Regrouping parti les in a region surrounding the best position −→g not only moves them away from the lo al minimum but allows them to re over their momentum, 94 whi h had been lost due to onvergen e. Regrouping is done a ording to the fol- lowing formula −→xi = −→g + [−→ri ]T · −−−−−−−→ range (Ωr)− 1 2 −−−−−−−→ range (Ωr) (7.7) where −−−−−−−→ range (Ωr) = [range1 (Ω r) , range2 (Ω r) , . . . , rangen (Ω r)]T . Equation 7.7 requires position lipping to keep the parti les within the on- strained design spa e, i.e xrj,L ≦ xi,j ≦ x r j,U . The upper and lower limits of the omponents of parti le's de ision ve tor an be found as xrj,L = max ( x0j,L, g r−1 j − 1 2 rangej (Ω r) ) (7.8) xrj,U = min ( x0j,U , g r−1 j + 1 2 rangej (Ω r) ) 7.3 Hybrid SQP-RegPSO Te hnique Swarm intelligen e starts with a set of randomly generated position points in the design spa e and denes its initial best lo ation, −→g , as the point of best tness among the swarm parti les; this best point a t as a gravity enter in the design spa e and attra ts other parti les to it; if the initial best point is of high quality, this redu es the total omputational eort to nd the global optimal solution. The te hnique we are proposing here for hybridization of sequential quadrati programing (SQP) and parti le swarm optimization is to use SQP to nd the initial high quality −→g before using RegPSO te hnique. The pseudo- ode of the hybrid s heme is summarized in Algorithm 7.2. 7.4 Optimization Test Cases This se tion will examine the e ien y of the proposed hybrid s heme and ompare it with the original RegPSO te hnique. To do so, a drag optimization that seeks a sho k free prole with CL = 0.4 at M = 0.8 and angle of atta k α = 1.5 ◦ is arried out. Two distin t starting points were sele ted to explore the existen e of a unique 95 Algorithm 7.2 Hybrid SQP-RegPSO pseudo ode For ea h parti le in the swarm Initialize parti le position −→xi Set its best position to be its initial position: −→pi = −→xi Initialize its velo ity −→ Vi = 0 End Set the global best position: −→g = −→pi where F (−→pi ) < F (−→pj ) for all i 6= j, or the SQP optimal solution, whi hever is better Set δ norm = 1 Set t = 1 Do For ea h parti le i Cal ulate its velo ity −→ Vi Update its position −→xi Clip position su h that xL,j ≤ xi,j ≤ xU,j Cal ulate its obje tive fun tion value (tness) F (−→xi) If the urrent tness value is better than its best tness (F (−→xi) < F (−→pi )) Set the urrent parti le position to be its best position: −→pi = −→xi End if End for Sele t the best parti le obje tive fun tion and position to be the global best if it is better than the stored value t = t+ 1 Compute δ norm using Eq 7.4 If δ norm ≤ ε Compute the range of ea h de ision variable range ( Ωrj ) S atter the parti les around −→g a ording to Eq 7.7 Apply position lipping with limits dened a ording to Eq 7.8 End if While (δ norm > ε and t ≤ max iterations). 96 global best 8 ; the rst is to start from NACA 0012, the se ond one is to start from NACA 00083 airfoil. 9 These ases were hosen be ause experiments showed that SQP optimization with these two starting points onverged to distin t optimal solutions with signi antly dierent obje tive fun tion values. All the omputations here are se ond order as we aiming to nd a global optimization method with reasonable omputational ost, not to ompare se ond and fourth order optimal results. 7.4.1 Constrained drag optimization of the NACA 0012 In this test ase, drag optimization of the NACA 0012 at a lift oe ient of CL = 0.4, M = 0.8 and angle of atta k α = 2◦ is presented. The obje tive fun tion an be written as F = CD + 10 · (CL − 0.4)2 (7.9) SQP optimization is used rst to nd a lo al optimum and thus provide an initial high quality −→g as a starting point for the RegPSO s heme. The optimal prole found by SQP optimization with BFGS Hessian approximation was not sho k free. RegPSO is used after SQP and su essfully rea hed a sho k free prole. In this test ase, the swarm size was 10 parti les and the number of de ision variables is 18. The rst phase of the hybrid s heme (SQP optimization) osts 92 CFD simulations while the se ond phase osts 490 CFD simulations whi h makes the total ost 582 CFD simulations to rea h sho k free prole. Using RegPSO alone as an optimization te hnique with a maximum of 1000 CFD simulations does not result in a sho k free prole, as the optimization has not yet onverged to a global optimum. Fig 7.4.1 shows the pressure distribution of the original NACA 0012 at transoni onditions, the SQP optimal prole, the RegPSO optimal shape, and the hybrid sho k free prole. Drag has been redu ed from 150 drag ounts at the initial lo al optimum found by SQP optimization to 8.89 drag ounts at the nal global optimum. Figure 7.4 shows the surfa e pressure distribution of the original NACA 0012 air- foil, and the optimized prole using SQP with BFGS approximate Hessian, RegPSO, 8 In the ase of a nite region in the design spa e for whi h the ow is sho k free, there may be a region of global best solutions diering only in the amount of dis retization error present. 9 That is, an airfoil from the 4-digit NACA series with 8.3% maximum thi kness. 97 NACA 0012 Original BFGS Optimal RegPSO 1000 CFD simulations Hybrid 492 Total CFD simulations Figure 7.3: Pressure Field of the original NACA 0012 and the optimized proles Figure 7.4: Surfa e pressure distribution of original NACA 0012 and the optimized proles 98 Figure 7.5: Airfoil shapes of original NACA 0012 and the optimized proles and hybrid BFGS-RegPSO te hniques. Figure 7.5 shows the dierent optimized air- foil proles resulted from the three optimization te hniques, note that the initial global best point of the hybrid s heme is the SQP optimal shape. To study the fun tion variation between the initial best known shape, whi h results from SQP op- timization, and the nal best shape of the hybrid s heme, the ve tor that onne ts these two points in design spa e was divided into 100 steps and the obje tive fun - tion value at ea h of these points omputed. Figure 7.7 shows that there is a des ent dire tion (the ve tor that onne ts the initial and nal −→g ) at the SQP optimal point but the SQP optimizer an not nd it. Possible reasons for that in lude: 1. An error in the omputed adjoint gradient. 2. The SQP gradient and the ve tor that onne ts the initial and nal −→g are perpendi ular due to the inuen e of numeri al noise and of the behavior of the Euler ow model physi s in the SQP gradient as the des ent dire tion resulting from adjoint omputations tends to minimize the drag by shifting the sho k wave forward so it o urs at lower Ma h number to redu e its strength, while this me hanism to redu e the drag is not stri tly followed by RegPSO s heme . 99 Figure 7.6: Fun tion minimization iterations in the se ond phase of the hybrid s heme (RegPSO) with a NACA 0012 airfoil as the starting point. To determine the real reason, we omputed the adjoint gradient and the nite dif- feren e gradient at the SQP lo al optimum, then omputed the angles between the Adjoint omputed gradient ve tor. Finite dieren e omputed gradient ve tor. Ve tor that onne ts the SQP (lo al) optimum and the hybrid (global) opti- mum. Table 7.1 shows the angles between those three ve tors. The adjoint and nite dieren e gradient dire tions are in ex ellent agreement but are almost perpendi ular to the third ve tor that onne ts lo al and global optima. This observation supports the hypothesis that the inuen e of numeri al noise and Euler ow model physi s prevents nding that des ent dire tion; however, this needs to be veried by another test ase. 100 0 20 40 60 80 100 120 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 BFGS Optimal to Hybrid optimal steps F BFGS optimal Hybrid optimal Figure 7.7: Obje tive fun tion values between SQP and hybrid optimal points (start- ing geometry: NACA 0012). Angle in degrees Adjoint and nite dieren e gradients 1.78 Adjoint gradient and SQP-hybrid ve tor 86.6 Finite dieren e gradient and SQP-hybrid ve tor 86.7 Table 7.1: Angles between adjoint, FD., and SQP-hybrid ve tors 101 7.4.2 Constrained drag minimization of NACA 00083 In this test ase, the same drag optimization problem is studied at the same transoni onditions, but starting from a thinner airfoil. Again, SQP ould not eliminate the sho k wave on the airfoil's upper surfa e, but a sho k free optimal shape is obtained using RegPSO in the se ond phase as shown in Fig 7.8. Figure 7.9 shows a omparison of the surfa e pressure of the original NACA 00083, SQP with BFGS approximate Hessian (lo al) optimum, and hybrid BFGS- RegPSO (global) optimum results. The BFGS optimal prole weakens the sho k strength to redu e the drag, but the hybrid s heme eliminates the sho k wave om- pletely. Figure 7.10 shows a prole omparison of the original NACA 00083, BFGS optimal, and the hybrid optimal prole. The drag oe ient is redu ed from 32 drag ounts at the BFGS optimal solution to 8.89 drag ounts at the hybrid optimal solution. The optimal value of the drag oe ient found by the hybrid s heme in this test ase is the same as the value obtained in the previous ase and in both ases, it is of order of the dis retization error. The total number of CFD simulations is 1080, 80 for the SQP phase and 1000 for the RegPSO phase. Again, to study the obje tive fun tion variation between the BFGS (lo al) opti- mum, and the nal (global) optimum of the hybrid s heme, we took 100 steps on the ve tor that onne ts these points and ompute the obje tive fun tion values. Figure 7.12 shows the variation in the obje tive fun tion value along that ve tor, whi h is not a des ent ve tor at the BFGS optimal point. Also, there is numeri al noise in the obje tive fun tion values at some stations along that ve tor, in addition to an ap- pearan e then vanishing of a sho k wave. As in the previous test ase, we omputed the obje tive fun tion gradient at the initial gBest using adjoint and nite dieren e strategy and al ulated its angle with the ve tor that onne ts the initial and the nal best position ve tors −→g obtained using the hybrid s heme. Table7.2 shows that the ve tor onne t the initial and nal best solutions is almost perpendi ular to the adjoint and nite dieren e gradients; it shows also that adjoint and nite dieren e gradients are a very good mat h. The orthogonality of these ve tors was previously noted in the last test ase as well. This strongly suggests that numeri al noise in the gradient and the inuen e of the ow model physi s prevents nding a des ent 102 NACA 00083 original BFGS Optimal Hybrid Optimal Figure 7.8: Pressure surfa e of NACA 00083 and the optimized proles Angle in degree Adjoint and Finite dieren e 0.27 Adjoint and SQP-hybrid ve tor 86.93 Finite dieren e and SQP-hybrid ve tor 86.87 Table 7.2: Angles between adjoint, FD, and BFGS-hybrid best solutions ve tors dire tion that leads from the SQP lo al optimum to the sho k free global optimum. 7.4.3 Aerodynami and thi kness onstraint drag minimization of NACA 0012 In this test ase thi kness is onstrained to be tc = 0.1 at xtc = 0.35, in addition to the lift onstraint CLc = 0.4 at transoni onditions, M = 0.8 and α = 1.5 ◦ . SQP with BFGS Hessian approximation optimization te hnique is rst used, the obje tive 103 Figure 7.9: Surfa e pressure of NACA 00083 and the optimized proles Figure 7.10: Prole omparison of NACA 00083, BFGS optimal, and hybrid optimal airfoils 104 Figure 7.11: Fun tion minimization iterations in the se ond phase of the hybrid s heme (RegPSO) of NACA 00083 start 0 20 40 60 80 100 120 0 0.005 0.01 0.015 0.02 0.025 BFGS optimal to hybrid optimal steps F BFGS optimal Hybrid optimal Figure 7.12: Obje tive fun tion values between BFGS and hybrid optimal points, starting geometry NACA 00083 105 fun tion to be minimized is F = CD + 10 (CL − 0.4)2 (7.10) subje t tot (xtc) = tc The thi kness onstraint an be ast as an equality onstraint equation using the surfa e parametrization matrixQ2 [AQ1] † ; for this test ase, xc < L1, so the thi kness onstraint equation an be written as t (xtc) = P1(xtc)− P3(xtc) = tc (7.11) ∴ ( a0 √ xc + a1xc + a2x 2 c + a3x 3 c )− (c0√xc + c1xc + c2x2c + c3x3c) = tc Re all that the relationship of polynomial oe ients and the design variables an be found using Q2 [AQ1] † as ai = ∑( Q2 [AQ1] † ) i,j yj (7.12) Equation (7.12) an be used with Eq (7.11) to derive an algebrai linear equation that relates the thi kness onstraint tc to design variables dire tly; this last equation an be used as a onstraint equation added to the optimization problem and this onstraint will be satised using a Lagrange multiplier. For the se ond phase of the hybrid s heme, the thi kness onstraint is satised using a penalty term in the obje tive fun tion; the obje tive fun tion for RegPSO phase an be written as F = CD + 10 (CL − 0.4)2 + (t (xtc)− tc)2 (7.13) The penalty method is used to approximately satisfy the thi kness onstraint term due to the la k of Lagrange multiplier in non-gradient based optimization te hniques. Drag oe ient have been redu ed from 230 to 120 drag ounts in the rst phase, then redu ed to 9 drag ounts using the hybrid s heme, a sho k free prole is obtained as shown in Figs (7.13, 7.14) ; this sho k free prole is obtained with 900 CFD simulations in total. 106 NACA 0012 original Hybrid Optimal Figure 7.13: Pressure surfa e of NACA 0012 and the optimized proles with thi kness onstraint Figure 7.14: Surfa e pressure of NACA 0012 and the optimized prole with thi kness onstraint 107 NACA 00083 original Hybrid Optimal Figure 7.15: Pressure surfa e of NACA 00083 and the optimized proles with thi k- ness onstraint 7.4.4 Aerodynami and thi kness onstraint drag minimization of NACA 00083 This test ase is a repeat of the last one, but the starting geometry is dierent. The starting geometry is NACA 00083 whi h does not satisfy the thi kness onstraint. The SQP optimizer is used rst;it satises the thi kness onstraint after the rst few iterations, then it starts to minimize the drag. Drag is redu ed from 60 drag ounts to 40 drag ounts in the rst phase. The se ond (RegPSO) phase redu ed the drag to 9 drag ounts and produ ed a sho k free optimal shape, as shown in Figs. 7.15 and 7.16. Comparison of the optimized proles of the last two test ases are shown in Fig. 7.17. The two proles are very lose in upper surfa e shape, whi h is tightly onstrained by the need to eliminate the sho k wave; there is more dieren e in the lower surfa e on whi h the ow is subsoni . The last observation suggests that for invis id ow simulations, there is a region of global optimal solution rather than one point global optimum. More pre isely, on e a sho k-free solution is obtained, the obje tive fun tion is too at to distinguish reliably between the drag of sho k-free solutions; this is not surprising, onsidering that the drag for sho k-free invis id ows is due solely to dis retization error. In all test ases, the new hybrid s heme was able to nd a sho k free airfoil prole regardless of the starting points. Only ten parti les were using in the swarm to s an the design spa e, whi h is roughly half of the number of design variables. This is lower than the re ommended population size required by GA to onverge well. The 108 Figure 7.16: Surfa e pressure of NACA 00083 and the optimized prole with thi k- ness onstraint optimization omputational ost is four to nine times the ost of SQP optimization. Using SQP optimization with ten dierent starting geometries might lead to a sho k free prole, but there is no guarantee for that. As the results above show, there is a high han e of getting trapped in lo al minima or suering a breakdown in the optimization pro ess be ause of solution noise. Table 7.3 shows the ee t of the thi kness penalty weight on the optimization results using se ond order a urate omputations. It shows that a thi kness penalty weight of 10 is su ient to satisfy the thi kness geometri onstraint, but it will in rease the value of CD opt and may prevent obtaining a sho k free optimal prole. Figure Shows a omparison of the optimal proles using two dierent thi kness penalty weights; the maximum thi kness of the optimized airfoils were 0.077 and 0.095 respe tively for thi kness penalty of 1 and 10. Therefore it is re ommended to use a thi kness penalty of order of 10 in order to satisfy a maximum thi kness onstraint. Figure 7.18 shows a omparison of the optimal pressure distribution of the two optimal proles; a sho k free pressure distribution is obtained using thi kness penalty weight of 1. The drag oe ient for this prole is 9 drag ounts whi h is 109 Figure 7.17: Prole omparison of the hybrid s heme optimal proles, starting ge- ometries NACA 0012 and NACA 00083 of order of dis retization error. It shows also that with thi kness penalty weight 10, there is a weak ompression wave on the airfoil upper surfa e whi h auses an in rease in the drag oe ient to 17 drag ounts. Thi k. penalty F opt CD opt CL opt t (xtc) Lift penalty Thi k. penalty weight 1 0.00134 0.00088 0.3997 0.077 0.000001 0.00048 10 0.00208 0.0017 0.4005 0.095 0.0000026 0.00025 Table 7.3: Thi kness penalty weight impa t on the optimization results of NACA 00083 using hybrid s heme 110 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1.5 −1 −0.5 0 0.5 1 1.5 x/c − Cp Penalty weight 10 Penalty weight 1 Figure 7.18: Thi kness penalty weight impa t on the Optimal pressure distribution of NACA 00083 using hybrid optimization s heme. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 0.04 x/c y/ c Thickness penalty 1 Thickness renalty 10 Figure 7.19: Thi kness penalty weight impa t on the Optimal prole with NACA 00083 starting geometry using hybrid optimization s heme. 111 Algorithm 7.3 Pseudo ode of RegPSO-SQP hybrid T optimization s heme 1. Use the initial geometry as the initial global best −→g . 2. Randomly generate swarm. 3. Add the initial geometry as a swarm member. 4. Use RegPSO strategy with stopping riterion of max 900 CFD simulations to nd an estimated global best −→g . 5. Use SQP with starting point −→g to nd better global best solution. 7.5 Comparing SQP-RegPSO with RegPSO-SQP optimization strategies I have presented results of the proposed hybrid optimization s heme in whi h SQP optimization is used to obtain a reasonably good upper bound on the global optimum. Based on the results of previous test ases, this s heme seems to be a promising global optimization s heme. Another hybrid s heme in whi h RegPSO pre edes SQP optimization is possible. In this se tion I will ompare the optimization results of the last ve optimization problems using the proposed hybrid s heme with the results of RegPSO-SQP s heme; I will refer to the latter s heme as hybrid T . Algorithm 7.3 shows a pseudo ode of hybrid T optimization. Both s hemes has two phases, SQP phase and RegPSO phase, but they dier in the order of the phases. Tables 7.4 and 7.5 summarize the optimization results for both s hemes. Figures 7.207.24 show a omparison of the obtained pressure distribution of the optimal shape obtained using these s hemes. These omparisons show that the hybrid s heme is better than the hybrid T for all ases; in Case II, the s hemes are nearly identi al in drag oe ient though dierent in nal shape. This was expe ted be ause in hybrid T if the rst (RegPSO) phase ould not rea h the basin of attra tion of the global minimum, the se ond (SQP) phase will not be able to nd it. Tables 7.4, 7.5 show a omparison of Phase I and Phase II optimal obje tive fun tion value for both hybrid and hybrid T for the ve test ases in this hapter. As a hybridization quality measure, the ratio of the optimal obje tive fun tion phase 112 Case Phase I Opt (SQP) Phase II Opt (RegPSO) Phase II/I I 10 0.0150 0.0009 0.0600 II 11 0.0032 0.0009 0.2813 III 12 0.0120 0.0013 0.1083 IV 13 0.0040 0.0013 0.3250 V 14 0.0040 0.0021 0.525 Table 7.4: Optimization results of the hybrid s heme in Phase I, II Case Phase I Opt (RegPSO) Phase II Opt (SQP) Phase II/I I 0.0137 0.0049 0.358 II 0.0019 0.0011 0.579 III 0.0125 0.0125 1.0 IV 0.0017 0.0017 1.0 V 0.0052 0.0052 1.0 Table 7.5: Optimization results of the hybrid T s heme in Phase I, II II to Phase I optimal is used and the smaller this value the more su essful the hybridization sequen e is. The tables show also that the sequen e used in hybrid is better than hybrid T . 113 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 x/c − Cp RegPSO−SQP SQP−RegPSO Figure 7.20: Hybrid and hybrid T optimal pressure distribution for NACA 0012 (Case I). 114 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 x/c − Cp RegPSO−SQP SQP−RegPSO Figure 7.21: Hybrid and hybrid T optimal pressure distribution for NACA 00083 (Case II). 115 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 x/c − Cp RegPSO−SQP SQP−RegPSO Figure 7.22: Hybrid and hybrid T optimal pressure distribution for NACA 0012 with 10% thi kness onstraint and unit thi kness penalty weight (Case III). 116 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 x/c − Cp RegPSO−SQP SQP−RegPSO Figure 7.23: Hybrid and hybrid T optimal pressure distribution for NACA 00083 with 10% thi kness onstraint and unit thi kness penalty weight (Case IV). 117 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 x/c − Cp RegPSO−SQP SQP−RegPSO Figure 7.24: Hybrid and hybrid T optimal pressure distribution for NACA 00083 with 10% thi kness onstraint and ten thi kness penalty weight (Case V). 118 Chapter 8 CONCLUSIONS AND RECOMMENDATIONS 8.1 Contributions and Con lusions This thesis presents the rst study of the use of high order nite volume CFD methods in aerodynami optimization. I have developed the rst nite volume based optimization ode that uses high order Ja obian in gradient omputations using the adjoint te hnique. A new geometry parametrization te hnique has been developed and presented; this new te hnique is based on least squares surfa e tting. To develop an e ient global optimizer, a new hybrid SQP-RegPSO s heme has been developed to nd the global optimal solution in the feasible design spa e. I presented both onstrained and un onstrained aerodynami optimization using sequential quadrati programming with the use of adjoint method to ompute the obje tive fun tion gradient. The omputations were based on higher order nite volume method on unstru tured grids. I took advantage of evaluating the exa t Ja obian matrix to ompute the ow adjoint eld and ow solution sensitivity. The omputed gradient based on the adjoint method is an ex ellent mat h with the orresponding nite dieren e gradient in subsoni ow for both se ond and fourth order s hemes. For transoni ow, the se ond order Ja obian mat hes well with the orresponding order of a ura y nite dieren e gradient; on the other hand, the fourth order Ja obian does not mat h the nite dieren e gradient well when using a limiter. This error in the fourth order gradient is attributed to the high sensitivity of the fourth order Ja obian matrix to the use of the limiter in transoni ow whi h tends to ae t the diagonal dominan e of the Ja obian matrix. To ta kle this problem and to enhan e the diagonal dominan e of the fourth order Ja obian matrix, we used the same non-zero stru ture of the se ond order Ja obian matrix 119 and dropped the rest of the non-zero elements; this tends to in rease the diagonal dominan e of the modied fourth order Ja obian matrix and improves its ondition number. I have developed a new approa h for smooth parametrization of airfoil shapes, based on a least-squares t of the parameters of two polynomials to a set of ontrol points. Compared with using all mesh points as design variables, this approa h redu es the size of the design spa e and eliminates os illations in the shape. We use a semi-torsional spring analogy to deform the mesh grid in the entire ow eld when the surfa e shape evolves during optimization. It is used also to al ulate the mesh sensitivity terms in the adjoint gradient. To test the developed gradient based optimizer; I have used the optimizer to ta kle inverse design problems in whi h only one optimal solution point in the de- sign spa e exists. Additional transoni drag minimization test ases have been pre- sented with and without lift onstraint. As a summary of test ases observations: A subsoni inverse design test ase shows that both se ond and fourth order s hemes rea hed their target geometry after almost the same number of iterations. The same behavior is observed for transoni inverse design test ase. This indi ates that the onvergen e is independent of the order of a ura y of the spatial dis retization. In a drag minimization test ase without lift onstraint, both se ond and fourth order s heme rea hed their optimum shape after almost the same number of optimiza- tion iterations. The dieren e of the resulting optimal airfoil shape is small. Both s hemes redu ed the drag by almost 50% of its original value but the lift also went down. In a drag minimization test ase with lift onstraint, the fourth order s heme was faster to rea h its optimal shape with 8 optimization iterations (whi h ost 49 CFD simulations), while the se ond order s heme took 13 iterations (112 CFD simulations) to rea h its optimum. Both s hemes rea hed their optimum with almost the same wall lo k run time and found nearly identi al airfoil optimal shapes. Based on all the test ases we presented, we on lude that the spatial dis retization error produ es a systemati error in the obje tive fun tion and has a little ee t on the nal optimal shape. A mesh renement study shows that both se ond and fourth order s hemes tend to give the same optimal value for the obje tive fun tion as we rene the mesh. 120 Based on the results of this resear h, I re ommend the use of a high order method for obje tive fun tion value omputations and for a urate predi tion of for es and sho k wave apture in transoni ows, and the use of a se ond order method (but based on the obtained forth order solution) to ompute the obje tive fun tion gradient needed for gradient based optimizer. If the mesh resolution makes se ond and fourth order simulation results similar, using the se ond order s heme for optimization will be e onomi al. I have developed also a new hybrid gradient/non-gradient based optimization te hnique that uses sequential quadrati programming with BFGS Hessian approx- imation te hnique in an initial gradient based optimization phase, followed by the regrouped parti le swarm optimization method in the non-gradient random sear h phase. The SQP optimization phase leads to a high quality initial best-solution point that redu es the total omputational eort of the RegPSO to rea h a sho k free prole. Unlike with the SQP s heme, the transoni drag optimization with lift onstraint of the NACA 0012 and NACA 00083 ta kled with the hybrid s heme lead to sho k free optimal shapes; in both ases, the obje tive fun tion gradient ve tor is almost perpendi ular to the ve tor that onne ts the SQP optimum and the hybrid opti- mum. This explains why the gradient based optimization an not rea h the global optimal, and also shows that its optimal solution depends on the starting geometry. The obje tive fun tion behavior along the ve tor that onne ts the SQP and the hybrid optima shows that numeri al noise may have a harmful impa t on gradient based optimization as well. Transoni drag optimization of NACA 0012 and NACA 00083 have been ta k- led with both aerodynami and geometri onstraints. For RegPSO, the geometri onstraint is satised using penalty method; in both ases, sho k free proles have been obtained. The last test ase (drag redu tion of NACA 00083 with a thi kness onstraint) shows that the hybrid s heme was able to nd a sho k free prole despite in reasing the thi kness during optimization. The omputational expense of the hy- brid s heme is four to nine times more than SQP optimization, but in all ases, a global optimal (sho k free) prole obtained, and in aerospa e industry, this is worth the in rease in omputational eort. 121 8.2 Future Work We intend to apply non gradient-based optimization te hniques like parti le swarm algorithm to study the ee t of CFD s heme order of a ura y on the nal optimal shape in the near future. Another obvious extension of the urrent work is aero- dynami optimization using RANS with turbulen e models. The ow solver should solve the turbulen e equations oupled with the mean ow equations; otherwise, the ow sensitivity values will be less ina urate. Thus, for a three dimensional ow simulation with k − ǫ turbulen e model, the blo k size of the Ja obian matrix is 7× 7. Fully oupled solvers have also been shown to onverge more e iently. For three-dimensional problems, the geometry parametrization should be re- pla ed with a perturbation parametrization (parametrization of the hange in the shape). This will redu e the required number of design variables whi h is very important for three dimensional optimization problems. Importantly, perturbation parametrization also prevents the hange in the initial pressure distribution when parametrizing the geometry itself. Exploring the ee t of the geometry parametrization te hnique used on SQP opti- mization results is needed. We spe ulate that using the right geometry parametriza- tion te hnique may lead to a design spa e in whi h aerodynami onstraints are simply onne ted and an SQP optimizer an nd a globally optimal solution. 122 Bibliography [1℄ Anderson, M.B. The Potential of Geneti Algorithms for Subsoni Wing De- sign. In AIAA Paper 95-3925, presented at the 1st AIAA Air raft Engineering, Te hnology, and Operations Congress, Los Angeles, CA, 1995. [2℄ Anderson, Murray B. Geneti Algorithms In Aerospa e Design: Substantial Progress, Tremendous Potential. In NATO/Von-Karmen-Institute Workshop on Intelligent System, 2002. [3℄ Anderson, W. Kyle, and Bonhaus, Daryl L. Aerodynami Design on Unstru - tured Grids for Turbulent Flows. Te hni al report, NASA Langley Resear h Center, Hampton, Virginia, 1997. [4℄ Avriel, M. Nonlinear Programming Analysis and Methods. Dover Publishing. ISBN 0-486-43227-0., 2003. [5℄ Azab, M. B., and Ollivier-Goo h, C. F. Higher Order Two Dimensional Aerody- nami Optimization Using Unstru tured Grids and Adjoint Sensitivity Compu- tations. In AIAA-2010-1430, 48th AIAA Aerospa e s ien es meeting, Orlando, FL, USA, 4-7Jan 2010, Orlando FL, 2010. [6℄ Azab, M. B., and Ollivier-Goo h, Carl. Constrained and Un onstrained Aerody- nami Quadrati Programming Optimization Using High Order Finite Volume Method and Adjoint Sensitivity Computations. In 49th Aerospa e s ien e meet- ing, Orlando FL, AIAA-2011-0183, 2011. [7℄ Batina, J. T. Unsteady Euler Airfoil Solutions Using Unstru tured Dynami Meshes. AIAA Journal ,, 28 No 8.:13811388, 1990. 123 [8℄ Beall, M. W., Walsh, J., and Shephard, Mark S. A essing CAD Geometry for Mesh Generation. In Pro eedings of the 12th International Meshing Roundtable, Sandia National Laboratories, pages 20033030, 2003. [9℄ Bergh, F., and Engelbre ht, A. A New Lo ally Convergent Parti le Swarm Optimiser. In Conferen e on Systems, Man and Cyberneti s, 2002. [10℄ Borrel, P., and Rappoport, A. Simple Constrained Deformations for Geometri Modeling and Intera tive Design. ACM Transa tions on Graphi s, 13:137155, 1994. [11℄ Brezillon, J. and Gauger, N. R. 2D and 3D Aerodynami Shape Optimisation Using the Adjoint Approa h. Aerospa e S ien e and Te hnology, 8:715727, 2004. [12℄ Broyden, C. G. The Convergen e of a Class of Double-rank Minimization Algo- rithms. Journal of the Institute of Mathemati s and Its Appli ations, 6:7690, 1970. [13℄ Byrd,R. H., Lu, P., and No edal, J. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM Journal on S ienti and Statisti al Comput- ing, 16, 5:11901208., 1995. [14℄ Cavallo, P. A., Hosangadi, A., Lee, R. A., and Dash, S. M. Dynami Unstru - tured Grid Methodology with Appli ation to Aero/Propulsive Flow Field. In AIAA Paper 1997-2310, 2002. [15℄ Chandrashekarappa, P., and Duvigneau, R. Metamodel-Assisted Parti le Swarm Optimization and Appli ation to Aerodynami Shape Optimization. Re- sear h Report RR-6397, INRIA, 2007. [16℄ Cler , M., and Kennedy, J. The Parti le Swarm - Explosion, Stability, and Convergen e in a Multidimensional Complex Spa e. Evolutionary Computation, IEEE Transa tions on, 6(1):5873, 2002. [17℄ Consentino, G., and Holst, T. Numeri al Optimization Design of Advan ed Transoni Wing Congurations. In AIAA 85-0424, 1985. 124 [18℄ Cook, P.H., M Donald M.A., and Firmin, M.C.P. . Aerofoil RAE 2822: Pres- sure Distributions, and Boundary Layer and Wake Measurements. In Experi- mental Data Base for Computer Program Assessment, AGARD Report AR 138. AGARD, 1979. [19℄ Dadone, Andrea, Mohammadi, Bijan, and Petruzzelli, Ni ola. In omplete Sen- sitivities and BFGS Methods for 3D Aerodynami Shape Design. Resear h Report RR-3633, INRIA, 1999. Projet M3N. [20℄ Duvigneau, R. Aerodynami Shape Optimization with Un ertain Operating Conditions using Metamodels. Resear h Report RR-6143, INRIA, 2007. [21℄ Eberhart, R. C., and Kennedy, J. A New Optimizer Using Parti le Swarm The- ory. In Sixth International Symposium on Mi ro Ma hine and Human S ien e, Nagoya, Japan, pp. 39-43., 1995. [22℄ Engelbre ht, A. P. Fundamentals of Computational Swarm Intelligen e. John Wiley and Sons, 2006. [23℄ Evers, G. I., and Ben Ghalia, M. Regrouping Parti le Swarm Optimization a New Global Optimization Algorithm with Improved Performan e Consisten y A ross Ben hmarks. In SMC IEEE International Conferen e, pages 39013908, o t. 2009. [24℄ Evers, George. An Automati Regrouping Me hanism to Deal with Stagna- tion in Parti le Swarm Optimization. Master's thesis, University of Texas-Pan Ameri an, 2009. [25℄ Farhat, C., Degand, C., Koobus, B., and Lesoinne, M. Torsional Springs for Two-Dimensional Dynami Unstru tured Fluid Meshes. Computer methods in applied me hani s and engineering, 163(1-4):231245, 1998. [26℄ Frink, N. T. Three-Dimensional Upwind S heme for Solving The Euler Equa- tions on Unstru tured Tetrahedral Grids. PhD thesis, Virginia Polyte hni In- stitute and State University, 1991. 125 [27℄ Frink, N. T., Parikh, P., and Pirzadeh, S. A Fast Upwind Solver for The Euler Equations on Three-Dimensional Unstru tured Meshes. In AIAA paper 91-0102, 1991. [28℄ Gage, P., and Kroo, I. Development of the Quasi-Pro edural Method for Use in Air raft Conguration Optimization. In AIAA Paper No. 92-4693, Presented at the Fourth AIAA/USAF/NASA/OAI Symposium on Multidis iplinary Analysis and Optimization, 1992. [29℄ Gage, P., and Kroo, I. A Role of Geneti Algorithms in a Preliminary Design Environment. In AIAA Paper No. 93-3933, 1993. [30℄ Goldberg, David E. Geneti Algorithms in Sear h, Optimization, and Ma hine Learning. Addison-Wesley Publishing Company, 1989. [31℄ Goldfarb, D. Fa torized Variable Metri Methods for Un onstrained Optimiza- tion. Mathemati s of Computation, 30(136):796811, 1976. [32℄ Gregg, R.D., and Misegades, K.P. Transoni Wing Optimization Using Evolu- tion Theory. In presented at the AIAA 25th Aerospa e S ien es Meeting, AIAA 87-0520, 1987. [33℄ Hart, Willian E. . Adaptive Global Optimization with Lo al Sear h. PhD thesis, University of California, San Diego, 1994. [34℄ Hassan, R., Cohanim, A. , and de We k,O. A Comparison of Parti le Swarm Optimization and The Geneti Algorithm. In AIAA-2005-1897, page 1897, 2005. [35℄ Hi ks, R. M. , and Henne, P. A. . Wing Design by Numeri al Optimization. Journal of Air raft, 15:407412., 1978. [36℄ Ira H. Abbott and Albert Edward Von Doenho. Theory of Wing Se tions. Dover, 1959. [37℄ Jameson, A. Aerodynami Design via Control Theory. Journal of S ienti Computing, 3:233260, 1988. 10.1007/BF01061285. 126 [38℄ Jameson, A. A Perspe tive on Computational Algorithms for Aerodynami Analysis and Design. Progress in Aerospa e S ien es, 37(2):197243, 2001. doi: DOI: 10.1016/S0376-0421(01)00004-5. [39℄ Jameson, A. Aerodynami Shape Optimization Using The Adjoint Method. VKI Le tur seriese, 2003. [40℄ Jameson, A. Optimum Transoni Wing Design Using Control Theory. In IU- TAM Symposium held in Guttingen, Germany 2-6 September 2002, page 253. Springer Netherlands, 2003. IUTAM Symposium Transsoni um IV: pro eedings of the IUTAM Symposium held in Guttingen, Germany 2-6 September 2002. [41℄ Jameson, A. , and Reuther, J. Control theory based airfoil design using Eu- ler equations. In AIAA/USAF/NASA/ISSMO Symposium on Multidis iplinary Analysis and Optimization, 1994. [42℄ Jameson, A., and Kim, S. Redu tion of The Adjoint Gradient Formula in The Continuous Limit. AIAA paper, 40:2003, 2003. [43℄ Jameson, A., Pier e,N. A., and Martinelli, L. Optimum Aerodynami Design Using The Navier-Stokes Equations. In 35th AIAA Aerospa e S ien es Meeting & Exhibit, Reno, NV, 1997. [44℄ Jang, M., and Lee, J. Geneti Algorithm Based Design of Transoni Air- foils Using Euler Equations. In AIAA Paper 2000-1584, Presented at the 41st AIAA/ASME/ASCE/AHS/ASC Stru tures, Stru tural Dynami s, and Materi- als, 2000. [45℄ Karthik, M. Appli ation Of The Dis rete Adjoint Method To Coupled Multidis i- plinary Unsteady Flow Problems For Error Estimation And Optimization. PhD thesis, Department of Me hani al Engineering, The University of Wyoming, 2009. [46℄ Kennedy, J., and Eberhart, R. Parti le Swarm Optimization. In Pro eedings of IEEE International Conferen e on Neural Networks. IV. pp. 1942-1948, 1995. [47℄ Kroo, I. Drag Due to Lift: Con epts for Predi tion and Redu tion. Annual Review of Fluid Me hani , 33:587617, 2001. 127 [48℄ Kulfan, B. M. Re ent Extensions and Appli ations of the CST Universal Para- metri Geometry Representation Method. In 7th AIAA Aviation Te hnology, Integration and Operations Conferen e, 2007. [49℄ Kulfan, B. M., Bussoletti, J. E. . Fundamental Parametri Geometry Represen- tations for Air raft Component Shapes. In 11th AIAA/ISSMO Multidis iplinary Analysis and Optimization Conferen e, 2006. [50℄ Li, W. , Huyse, L. and Padula, S. Robust airfoil optimization to a hieve drag redu tion over a range of ma h numbers. Stru tural and Multidis iplinary Op- timization, 24:3850, 2002. 10.1007/s00158-002-0212-4. [51℄ Martineau, D. G., Georgala, J. M. A Mesh Movement Algorithm for High Quality Generalised Meshes. In 42nd AIAA Fluid Dynami s Conferen e and Exhibit, AIAA Paper 2004-0614, Reno, NV, 2004. [52℄ Masuda, H., Yoshioka, Y., and Furukawa, Y. Intera tive Mesh Deformation Us- ing Equality-Constrained Least Squares. Computers and Graphi s, 30(6):936 946, 2006. [53℄ Mengistu, Temesgen, and Ghaly, Wahid. Aerodynami Optimization of Tur- boma hinery Blades Using Evolutionary Methods and ANN-Based Surrogate Models. Optimization and Engineering, 9:239255, 2008. 10.1007/s11081-007- 9031-1. [54℄ Mi halak, C. E ient High-Order A urate Unstru tured Finite-Volume Algo- rithms for Vis ous and Invis id Compressible Flows. PhD thesis, The University of British Columbia, 2009. [55℄ Mi halak, C., and Ollivier-Goo h, C. Dierentiability of Slope Limiters on Unstru tured Grids. In Fourteenth Annual Conferen e of the Computational Fluid Dynami s So iety of Canada, 2006. [56℄ Mi halak, C., and Ollivier-Goo h, C. A ura y Preserving Limiter for The High- Order A urate Solution of The Euler Equations. Journal of Computational Physi s, 228(23):86938711, 2009. 128 [57℄ Mi halak, Christopher, and Ollivier-Goo h, Carl. Globalized Matrix-Expli it Newton-GMRES for The High-Order A urate Solution of The Euler Equations. Computers and Fluids, 39:11561167, 2010. [58℄ Mousavi, A., Castonguay, P., and Nadarajah, S. K. Survey of Shape Parame- terization Te hniques and its Ee t on Three-Dimensional Aerodynami Shape Optimization. In Pro eedings of the AIAA 18th Computational Fluid Dynami s Conferen e. AIAA 2007-3837, 25 - 28 June 2007 2007. [59℄ Mouser, C. R., and Dunn, S. A. Comparing Geneti Algorithms and Parti- le Swarm optimisation for an Inverse Problem Exer ise. In Rob May and A. J. Roberts, editors, Pro . of 12th Computational Te hniques and Appli- ations Conferen e CTAC-2004, volume 46, pages C89C101, Mar h 2005. http://anziamj.austms.org.au/V46/CTAC2004/Mous [Mar h 21, 2005℄. [60℄ Nadarajah, S., and Jameson, A. A Comparison of The Continuous and Dis- rete Adjoint Approa h to Automati Aerodynami Optimization. AIAA paper, 667:2000, 2000. [61℄ Nejat, A., and Ollivier-Goo h, C. A High-Order A urate Unstru tured Finite Volume Newton-Krylov Algorithm for Invis id Compressible Flows. Journal of Computational Physi s, 227(4):25822609, 2008. [62℄ Neme , M., and Aftosmis, M. J. Adjoint Sensitivity Computations for an Embedded-Boundary Cartesian Mesh Method. Journal of Computational Physi s, 227(4):27242742, 2008. [63℄ Neme , M., and Zingg, D. W. . Newton-Krylov Algorithm for Aerodynami Design Using the Navier-Stokes Equations. AIAA Journal, 40:11461154, 2002. [64℄ Nielsen, E., and Anderson, K. Aerodynami Design Optimization on Unstru - tured Meshes Using the Navier-Stokes Equations. Te hni al report, NASA Lan- gley Te hni al Report Server, 1998. [65℄ Nielsen, Eri J., and Kleb, Bil. E ient Constru tion of Dis rete Adjoint Op- erators on Unstru tured Grids by Using Complex Variables. In 43rd AIAA Aerospa e S ien es Meeting and Exhibit Reno, Nevada, 2005. 129 [66℄ No edal, J., and Wright Stephen, J. Numeri al Optimization. Springer New York, 2006. [67℄ Obayashi, S. Geneti Algorithm for Aerodynami Inverse Optimization Prob- lems. In Geneti Algorithms in Engineering Systems: Innovations and Appli a- tions, 1995. [68℄ Ollivier-Goo h, C. and Van Altena, M. A High-Order-A urate Unstru tured Mesh Finite-Volume S heme for The Adve tion-Diusion Equation. Journal of Computational Physi s, 181(2):729752, 2002. [69℄ Ollivier-Goo h, C., Nejat, A., and Mi halak, C. On Obtaining and Verifying High-Order Finite-Volume Solutions to the Euler Equations on Unstru tured Meshes. AIAA Journal, 47(9):21052120, 2009. [70℄ Oyama, A., Obayashi, S., and Nakahashi, K. Transoni Wing Optimization Using Geneti Algorithm. In AIAA Paper 97-1854, 13th Computational Fluid Dynami s Conferen e, 1997. [71℄ Oyama, A., Obayashi, S., and Nakahashi, K. Fra tional fa torial design of geneti oding for aerodynami optimization. In AIAA Paper 99-3298, 3rd AIAA Weakly Ionized Gases Workshop, 1999. [72℄ Piegl, L., and Tiller, W. The NURBS Book. Springer, Berlin, 1995. [73℄ Reuther, J., Jameson, A., Farmer, J., Martinelli, L., and Saunders, D. Aerody- nami Shape Optimization of Complex Air raft Congurations via an Adjoint Formulation. In 34rd AIAA Aerospa e S ien es Meeting and Exhibit, no. AIAA 1996-0094, 1996. [74℄ Roe, P. L. Approximate Riemann Solvers, Parameter Ve tors, and Dieren e S hemes. Journal of Computational Physi s, 43(2):357372, 1981. doi: DOI: 10.1016/0021-9991(81)90128-5. [75℄ Saad, You ef, and S hultz, Martin H. GMRES: a Generalized Minimal Resid- ual Algorithm for Solving Nonsymmetri Linear Systems. SIAM J. S i. Stat. Comput., 7:856869, July 1986. 130 [76℄ Samareh, Jamshid A. A Survey of Shape Parameterization Te hniques. In CEAS/AIAA/ICASE/NASA Langley International Forum on Aeroelasti ity and Stru tural Dynami s Williamsburg, VA, Also, NASA/CP-1999-209136, June 1999, pp. 333-343., 1999. [77℄ Samareh, Jamshid A. Geometry and Grid/Mesh Generation Issues for CFD and CSM Shape Optimization. Optimization and Engineering, 6:2132, 2005. 10.1023/B:OPTE.0000048535.08259.a8. [78℄ Sederberg, Thomas W., and Parry, S ott R. Free-Form Deformation of Solid Geometri Models. SIGGRAPH Comput. Graph., 20:151160, August 1986. [79℄ Shi, Y., and Eberhart, R. C. Empiri al Study of Parti le Swarm Optimization. In Pro eedings of Congress onEvolutionary Computation, CEC 99., volume 3, page 1950 Vol. 3, 1999. [80℄ Stein, K., Tezduyar, T., and Benney, R. Mesh Moving Te hniques for Fluid- Stru ture Intera tions with Large Displa ements. J. Appl. Me h., 70:5863, 2003. [81℄ Tautges,Timothy J. The Common Geometry Module (CGM): A Generi , Ex- tensible Geometry Interfa e. In Pro eedings of the 9th International Meshing Roundtable, pages 337348, 2000. [82℄ Tezduyar, T. E., Behr, M., Mittal, S., and Johnson, A. A. Computation of Unsteady In ompressible Flows with The Stabilized Finite Element Methods: Spa e-Time Formulations, Iterative Strategies and Massively Parallel Imple- mentations. New Methods in Transient Analysis, 246:0, 1992. [83℄ Truong, Anh H., Oldeld, Chad A., and Zingg, David W. Mesh Movement for a Dis rete-Adjoint Newton Krylov Algorithm for Aerodynami Optimization. AIAA Journal, 46, issue 7:16951704, 2008. [84℄ Venkatakrishnan, V. Convergen e to Steady State Solutions of The Euler Equa- tions on Unstru tured Grids with Limiters. Journal of Computational Physi s, 118(1):120130, 1995. 131 [85℄ Venter, G., and Sobiesz zanski-Sobieski, J. Multidis iplinary Optimization of a Transport Air raft Wing Using Parti le Swarm Optimization. Stru t Multidis- iplinary Optimization, 26:pp 121131, 2004. [86℄ Watt, A., and Watt, M. Advan ed Animation and Rendering Te hniques. Addison-Wesley Publishing Company, New York, 1992. [87℄ Zhu, C., Byrd, R. H., Lu, P., and No edal, J. Algorithm 778: L-BFGS-B: Fortran Subroutines for Large-S ale Bound-Constrained Optimization. ACM Transa tions on Mathemati al Software (TOMS), 23(4):550560, 1997. [88℄ Zymaris,A.S., Papadimitriou, D.I., Giannakoglou, K.C., and Othmer,C. Feasi- bility Study of Constant Eddy-Vis osity Assumption in Gradient-Based Design Optimization. Journal of Computational Physi s, 229:52285245, 2010. 132
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Aerodynamic optimization using high-order finite-volume...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Aerodynamic optimization using high-order finite-volume CFD simulations Azab, Mohammad Baher 2011
pdf
Page Metadata
Item Metadata
Title | Aerodynamic optimization using high-order finite-volume CFD simulations |
Creator |
Azab, Mohammad Baher |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | The growth of computer power and storage capacity allowed engineers to tackle engineering design as an optimization problem. For transport aircraft, drag minimization is critical to increase range and reduce operating costs. Lift and geometric constraints are added to the optimization problem to meet payload and rigidity constraints of the aircraft. Higher order methods in CFD simulations have proved to be a valuable tool and are expected to replace current second order CFD methods in the near future; therefore, exploring the use of higher order CFD methods in aerodynamic optimization is of great research interest and is one goal of this thesis. Gradient-based optimization techniques are well known for fast convergence, but they are only local minimizers; therefore their results depend on the starting point in the design space. The gradient-independent optimization techniques can find the global minimum of an objective function but require vast computational effort; therefore, for global optimization with reasonable computational cost, a hybrid optimization strategy is needed. A new least-squares based geometry parametrization is used to describe airfoil shapes and a semi-torsional spring analogy mesh morphing tool updates the grid everywhere when the airfoil geometry changes during shape optimization. For the gradient based optimization scheme, both second and fourth order simulations have been used to compute the objective function; the adjoint approach, well known for its low computational cost, has been used for gradient computation and matches well with finite difference gradient. The gradient based optimizer have been tested for subsonic and transonic inverse design problems and for drag minimization without and with lift constraint to validate the developed optimizer. The optimization scheme used is Sequential Quadratic Programming (SQP) with the BFGS approximation of the Hessian matrix. A mesh refinement study is presented for an aerodynamically constrained drag minimization problem to show how second and fourth order optimal results behave with mesh refinement. A hybrid particle swarm / BFGS scheme has been developed for use as a global optimizer. It has been tested on a drag minimization problem with lift constraint; the hybrid scheme obtained a shock free profiles, while gradient-based optimization could not in general. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-10-19 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0072307 |
URI | http://hdl.handle.net/2429/38061 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
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
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
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