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