Shape Parameterization in Aircraft Design: A Novel Method, Based on B-Splines
Michiel H. Straathof
Shape Parameterization in Aircraft Design: A Novel Method, Based on B-Splines
Proefschrift
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft; op gezag van de Rector Magnificus prof. ir. K.C.A.M. Luyben; voorzitter van het College voor Promoties in het openbaar te verdedigen op vrijdag 3 februari 2012 om 12.30 uur
door
Michiel Hannes STRAATHOF
ingenieur luchtvaart en ruimtevaart geboren te Leiden
Dit proefschrift is goedgekeurd door de promotor: Prof. dr. ir. M.J.L. van Tooren
Samenstelling promotiecommissie: Rector Magnificus, Prof. dr. ir. M.J.L. van Tooren, Prof. dr. T. B¨ack, Prof. dr. ir. B. Koren, Prof. dr. ir. H.W.M. Hoeijmakers, Prof. dr. ir. A.W. Heemink, Dr. ir. G. Carpentieri, B.M. Kulfan,
voorzitter Technische Universiteit Delft, promotor Universiteit Leiden Universiteit Leiden / CWI Amsterdam Universiteit Twente Technische Universiteit Delft Cardano, Londen Boeing (retired), Seattle, WA
keywords: shape, parameterization, aircraft, design, B-splines, Class-ShapeRefinement-Transformation, adjoint, euler, optimization ISBN/EAN: 978-94-6191112-4 c 2012 by Michiel H. Straathof Copyright ⃝ All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording, or by any information storage and retrieval system, without written permission from the author. Cover design: Merijn Straathof Printed in the Netherlands
Summary
Shape Parameterization in Aircraft Design: A Novel Method, Based on B-Splines This thesis introduces a new parameterization technique based on the Class-ShapeTransformation (CST) method. The new technique consists of an extension to the CST method in the form of a refinement function based on B-splines. This ClassShape-Refinement-Transformation (CSRT) method has the same advantages as the original CST method, while also allowing for local deformations in a shape. Most of the work presented in this thesis has been published and presented at conferences [75–80]. Chapter 1 starts with some historical background on aerodynamic shape optimization. It takes the reader from the early WWI fighter planes to Lindbergh’s Spirit of St. Louis, to the iconic Douglas DC-2, and finally to the modern airliners of today. It describes how the emphasis in aircraft design shifted from range and speed to fuel burn and noise. The objective of the work is introduced: to develop a new way of aerodynamically optimizing aircraft shapes, to be able to design novel aircraft geometries. The chapter also gives an overview of current practice in the fields of parameterization, aerodynamic analysis and optimization, which provided the basis of the work presented in this thesis. The chapter ends with a short introduction of the CSRT method. Chapter 2 gives an overview of a number of different parameterization techniques that are currently used in aerodynamic shape optimization. It elaborates on some of these methods, with special emphasis on the Class-Shape-Transformation technique. The Class-Shape-Refinement-Transformation method is explained in detail in Chapter 3. The chapter starts with a theoretical introduction to B-spline curves and surfaces, and their implementation into the refinement function. A parametric study then follows, which investigates the performance of the CSRT method by seeing how well a CSRT curve can be fitted to a given set of airfoil data. A highly non-linear relationship was found between the accuracy of the fit and the number of coefficients used. Furthermore, it was shown that the behavior of the CSRT method is highly problem dependent. Chapter 4 presents a novel way of handling volume constraints in two and three dimensions, making optimal use of the ability of the CSRT method to efficiently model
vi local deformations. With this new constraint handling method, volume constraints can be implemented by applying bounds directly on the design variables. Chapter 5 introduces two different design frameworks (low and high fidelity) that were used to test the performance of the CSRT method within an actual optimization problem. The low fidelity design framework combines the panel method code VSAERO with both local and global optimization algorithms. It makes use of a Kriging response surface to limit computation times. The high fidelity framework employs an in-house Euler/adjoint solver in combination with a trust region reflective algorithm and an active set algorithm. In Chapter 6, the results of a number of test cases are presented. The low fidelity framework was used to optimize a three-dimensional wing. This test case showed some of the capabilities of the low fidelity framework, but also some of its shortcomings. Even though a considerable increase in lift-to-drag ratio could be achieved, local minima in the design space proved to be a potential problem. For all high fidelity cases, a two-step optimization approach was used. In the first “general” step, the Bernstein coefficients of the shape function were varied. In the second “regional” step, the B-spline coefficients of the refinement function were varied to refine the shape. The high fidelity framework was first used to optimize two standard shapes: the NACA0012 airfoil (2D) and the ONERA M6 wing (3D). In both cases, the first optimization step lead to a significant increase in L/D, with a smaller improvement resulting from the refinement step. The last two test cases were performed on a blended-wing-body geometry. The first of these cases followed the trend seen in the previous test cases: a large improvement resulting from the first step, and a smaller improvement resulting from the second step. However, the last test case showed a clear deviation from the trend. The first step only lead to a marginal increase in lift-to-drag ratio, while the refinement step increased L/D by almost 20 %. Apparently, in cases in which the Bernstein surface of the shape function cannot cope with the design space, the B-spline surface of the refinement function can compensate for this. Appendix A shows the structure of the Euler and adjoint codes used in the high fidelity framework. It also gives examples of the unstructured meshes that were used for the two- and three-dimensional cases. Appendix B shows the results of the validation of the implementation of the CSRT method into the adjoint code.
Samenvatting
Parametrisatie van Vliegtuigvormen: Een Nieuwe Aanpak, Gebaseerd op B-Splines Dit proefschrift beschrijft een nieuwe parametrisatietechniek, gebaseerd op de ClassShape-Transformation (CST) methode. Deze nieuwe techniek is een uitbreiding van de CST methode in de vorm van een refinement function, bestaande uit B-splines. Deze Class-Shape-Refinement-Transformation (CSRT) methode heeft dezelfde voordelen als de originele CST methode, met als bijkomend voordeel dat ook plaatselijke vervormingen kunnen worden gemodelleerd. Het grootste deel van dit proefschrift is eerder gepubliceerd en gepresenteerd op conferenties [75–80]. Hoofdstuk 1 begint met de geschiedenis van aerodynamisch vliegtuigontwerp. De lezer wordt meegenomen van de vroege gevechtsvliegtuigen uit WOI naar Lindbergh’s Spirit of St. Louis, naar de beroemde Douglas DC-2 en uiteindelijk naar de moderne lijnvliegtuigen van vandaag de dag. Het hoofdstuk beschrijft hoe de nadruk bij vliegtuigontwerp verschoof van bereik en snelheid naar brandstofverbruik en lawaai. Het doel van dit werk wordt ge¨ıntroduceerd: het ontwikkelen van een nieuwe manier om vliegtuigvormen aerodynamisch te optimaliseren, om het ontwerpen van baanbrekende vliegtuiggeometrie¨en mogelijk te maken. Het hoofdstuk geeft tevens een overzicht van veelgebruikte technieken op het gebied van parametrisatie, stromingsanalyse en optimalisatie, die de basis vormden van het werk in dit proefschrift. Het hoofdstuk eindigt met een korte introductie van de CSRT methode. Hoofdstuk 2 geeft een overzicht van parametrisatietechnieken die veel gebruikt worden bij het aerodynamisch ontwerpen van vliegtuigen. Een aantal van deze technieken wordt in meer detail beschreven, waarbij in het bijzonder aandacht wordt besteed aan de Class-Shape-Transformation methode. De Class-Shape-RefinementTransformation methode wordt in detail besproken in Hoofdstuk 3. Dit hoofdstuk begint met de theorie achter B-splines en de manier waarop deze ge¨ımplementeerd zijn in de refinement function. Dan volgt een parametrische studie, die de prestaties van de CSRT methode onderzoekt op basis van hoe goed een CSRT curve een bepaald vleugelprofiel kan benaderen. Er werd een sterk niet-lineair verband gevonden tussen de nauwkeurigheid van de benadering en het aantal gebruikte vormco¨effici¨enten. Verder werd ook aangetoond dat het gedrag van de CSRT methode zeer probleemspecifiek is.
viii Hoofdstuk 4 presenteert een nieuwe manier van omgaan met volumetrische randvoorwaarden in twee en drie dimensies. Deze nieuwe methode maakt optimaal gebruik van de mogelijkheid van de CSRT methode om effici¨ent plaatselijke vervormingen te modelleren. Met deze nieuwe techniek kunnen volumetrische randvoorwaarden direct worden vertaald naar begrenzingen van de ontwerpvariabelen. Hoofdstuk 5 introduceert twee ontwerpkaders (low en high fidelity), opgezet om de werking van de CSRT methode in een echt optimalisatieprobleem te testen. Het low fidelity kader bestaat uit de paneelcode VSAERO in combinatie met zowel plaatselijke als globale optimalisatie algoritmes. Het maakt gebruikt van een Kriging model om de benodigde rekentijd te beperken. Het high fidelity kader gebruikt een in-house Euler/adjoint code gecombineerd met een trust region reflective algoritme en een active set algoritme. In Hoofdstuk 6 worden de resultaten van een aantal testcases gepresenteerd. Het low fidelity kader werd gebruikt om een drie-dimensionale vleugel te optimaliseren. Deze testcase liet zowel de voordelen als tekortkomingen van het low fidelity kader zien. Er werd weliswaar een behoorlijke verbetering van de lift-weerstand verhouding gerealiseerd, maar lokale minima in de ontwerpruimte bleken een potentieel probleem. Voor alle high fidelity testcases werd een twee-staps strategie gebruikt. In de eerste “algemene” stap werden alleen de Bernstein-co¨effici¨enten van de shape function gevarieerd. In de tweede “regionale” stap werden de B-spline-co¨effici¨enten van de refinement function gevarieerd om de vorm te verfijnen. Het high fidelity kader werd eerst toegepast op de optimalisatie van twee basisvormen: het NACA0012 profiel (2D) en de ONERA M6 vleugel (3D). In beide gevallen leidde de eerste optimalisatiestap tot een significante verbetering van de lift-weerstand verhouding, met nog een kleine verbetering als gevolg van de verfijningsstap. De twee laatste testcases waren gebaseerd op een blended-wing-body geometrie. De eerste van deze cases liet dezelfde trend zien als de vorige testcases: een grote verbetering na de eerste stap en een kleinere verbetering na de tweede stap. De laatste test case liet echter een afwijkend resultaat zien. De eerste stap leidde slechts tot een marginale verbetering van de lift-weerstand verhouding, terwijl de verfijningsstap deze met bijna 20 % liet toenemen. Kennelijk zijn er gevallen waarin het Bernsteinoppervlak van de shape function niet goed kan omgaan met de ontwerpruimte, wat vervolgens gecompenseerd wordt door het B-spline oppervlak van de refinement function. Appendix A geeft de structuur weer van de Euler en adjoint codes die gebruikt zijn in het high fidelity kader. Ook zijn van zowel de twee- als drie-dimensionale testcases voorbeelden opgenomen van de ongestructureerde grids waarop de Euler code is uitgevoerd. Appendix B bevat de resultaten van de validatie van de implementatie van de CSRT methode in de adjoint code.
Contents
Summary
v
Samenvatting 1
2
3
vii
Introduction 1.1 Background . . . . . . . . . . 1.2 Current practice . . . . . . . 1.2.1 Parameterization . . 1.2.2 Aerodynamic analysis 1.2.3 Optimization . . . . 1.3 CSRT method . . . . . . . . . 1.4 Thesis overview . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
1 3 9 9 10 11 12 15
Parameterization 2.1 Background . . . . . . . . . . . . . . . . . 2.2 Common methods . . . . . . . . . . . . . . 2.2.1 Discrete approach . . . . . . . . . 2.2.2 Analytical approach . . . . . . . . 2.2.3 Geometric approach . . . . . . . . 2.2.4 Polynomial and spline approaches 2.3 CST method . . . . . . . . . . . . . . . . . 2.3.1 CST in two dimensions . . . . . . 2.3.2 CST in three dimensions . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
17 19 21 21 21 22 24 25 25 27
. . . . . . . . .
31 33 33 36 38 40 42 44 44 54
CSRT method 3.1 B-splines . . . . . . . . . . 3.1.1 B-spline curves . 3.1.2 B-spline surfaces . 3.2 Refinement function . . . 3.3 CSRT performance . . . . 3.3.1 Fitting process . . 3.3.2 Correlation factor 3.3.3 Parametric study 3.4 Global shape parameters .
. . . . . . . . .
. . . . . . . . .
. . . . . . .
. . . . . . . . .
. . . . . . .
. . . . . . . . .
. . . . . . .
. . . . . . . . .
. . . . . . .
. . . . . . . . .
. . . . . . .
. . . . . . . . .
. . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
x 4
Volume constraints 4.1 Typical shape constraints . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Volume constraint handling in 2D . . . . . . . . . . . . . . . . . . . . 4.3 Volume constraint handling in 3D . . . . . . . . . . . . . . . . . . . .
57 59 60 64
5
Design frameworks 5.1 Low fidelity design framework . . . . 5.1.1 Panel methods (VSAERO) . 5.1.2 Global optimization . . . . . 5.1.3 Response surface methods . 5.2 High fidelity design framework . . . . 5.2.1 Euler method . . . . . . . . 5.2.2 Gradient-based optimization
6
7
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
69 71 71 78 79 83 83 93
Test cases 6.1 Low fidelity framework: 3D wing . . . 6.1.1 Test case definition . . . . . . 6.1.2 Test case results . . . . . . . . 6.2 High fidelity framework: NACA0012 . 6.2.1 Test case definition . . . . . . 6.2.2 Test case results . . . . . . . . 6.3 High fidelity framework: ONERA M6 . 6.3.1 Test case definition . . . . . . 6.3.2 Test case results . . . . . . . . 6.4 High fidelity framework: BWB . . . . 6.4.1 Test case definition . . . . . . 6.4.2 Results of the first test case . 6.4.3 Results of the second test case
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
103 105 105 107 109 111 111 118 118 120 125 126 127 131
Conclusions
135
A Euler code and meshes
139
B Adjoint validation results
151
Bibliography
157
Acknowledgements
167
About the author
169
CHAPTER
1
Introduction
“For a true writer each book should be a new beginning where he tries again for something that is beyond attainment. He should always try for something that has never been done or that others have tried and failed. Then sometimes, with great luck, he will succeed.” - Ernest Hemingway (1899 - 1961)
1.1. BACKGROUND
1.1
3
Background
During the early development of aircraft, not much attention was paid to aerodynamic efficiency. Structural design did not yet allow for cantilever wings, resulting in numerous struts and wires to support the aerodynamic forces. The huge amounts of parasite drag caused by these struts and wires prohibited true optimization of the external shape. This started to change during WWI, when speed and range became important for fighters, bombers and observation aircraft. Since the drag of an aircraft scales with the square of its speed, drag reduction suddenly appeared high on the agenda. Reducing drag also contributes to a larger range, which is illustrated by the Breguet range equation [67]: ( ) H W0 CL (1.1) R = ηtot ln g Wf CD The first term, H, represents the energy contained in the fuel per unit mass, hence it is a measure of fuel quality. Dividing H by the gravitational constant g gives the energy contained in the fuel per unit weight. The next term is the total efficiency of the propulsion system, ηtot . Most improvements in aircraft performance nowadays are actually achieved as a result of better engines. The third term, W0 /Wf is the ratio between the weight at the beginning and end of the cruise leg of the flight profile and is hence a measure of structural efficiency. Using stronger materials, such as carbon fiber composites, this factor can be improved. Finally, the last term is the ratio between the lift coefficient CL and the drag coefficient CD , better known as the lift-to-drag ratio or L/D. The lift-to-drag ratio is a common measure of the aerodynamic performance of an aircraft shape and the main focus of this thesis will be to find means of improving L/D. When speed and range became an issue for aircraft designers, a lot of effort was spent to aerodynamically optimize aircraft shapes. Around the same time, Ludwig Prandtl demonstrated that thick airfoils had superior performance over thin airfoils in the subsonic regime [63]. Earlier, a lack of knowledge about Reynolds effects in wind tunnel testing had led to the wrong conclusion that thin airfoils were better. The advancements in aerodynamic design stemming from Prandtl’s work and that of others eventually led to the development of the Spirit of St. Louis, the aircraft with which Charles Lindbergh performed his famous flight across the Atlantic Ocean in 1927. It could fly non-stop for over 33 hours, covering a distance of almost 6500 km. Looking at Fig. 1.1, a number of aerodynamic design improvements can be seen relative to earlier designs. Firstly, the steel tube fuselage of the Spirit of St. Louis is covered with fabric to allow the air to flow past smoothly. Secondly, the number of struts connecting the wing and landing gear to the fuselage is kept low and they are aerodynamically shaped for low drag. In later versions a cowling was added to the propeller (not present in Fig. 1.1). The end of the 1920s also saw the appearance of wooden monocoque fuselages and wings. Even though this resulted in much cleaner aircraft, a few struts were usually still necessary to support the wing (see for example the Fokker F.VII). The next decade saw the introduction of metal-skin aircraft with the Boeing
4
CHAPTER 1. INTRODUCTION
Figure 1.1: Lindbergh’s 1927 Spirit of St. Louis (bottom right), compared to the 1911 Fokker Spin 3 (top left)
247 (1934) and the Douglas DC-2 (1935). Because the skin of these aircraft could carry part of the loads occurring during flight, external struts were no longer required, making true aerodynamic optimization possible for the first time. Looking at Fig. 1.2, this becomes very apparent. Except for the propellers and the rear landing gear, the entire exterior of the aircraft consists of a smooth aluminum skin. Even the engines are completely covered. Also, the intersection between the wings and the fuselage has been aerodynamically shaped to prevent the air flow from separating. Looking at the wings, it can be seen that they are tapered and swept backwards, which also decreases drag. Still, most of the drag reduction described here can be attributed to a decrease in parasite drag. Figure 1.3 nicely illustrates the history of parasite drag throughout the last century. It can be seen that towards the 1970s, the parasite drag started to approach the value of turbulent flow skin friction drag. The data in Fig. 1.3 already looks fairly impressive, but the following example might help to really put it into perspective. For a medium size subsonic airplane with Boeing 777 technology and a nominal design range of 9000 km: i. a 1 % decrease in drag at a fixed Maximum Takeoff Weight (MTOW) increases the design range by 1.8 % ii. a 1 % decrease in drag for a fixed design range will reduce the MTOW by 800 kg iii. a 1 % decrease in drag for a fixed design range has the same effect as reducing the structural weight by 550 kg
1.1. BACKGROUND
5
Figure 1.2: The Douglas DC-2 in flight
0.0300
BIPLANES
0.0200
CD P RETRACTABLE GEAR
MONOPLANE
0.0100
JET ENGINE
TURBULENT FLOW SKIN FRICTION
0 1900
1910
1920
1930 1940 YEAR
1950
1960
1970
Figure 1.3: History of parasite drag
6
CHAPTER 1. INTRODUCTION
Figure 1.4: The Boeing 787 upon landing
For a supersonic aircraft, the consequences of reducing drag are even greater. A reduction in drag coefficient of 0.0001 (1 drag count) would reduce the airplane gross weight by 4500 kg, save 3500 kg of fuel (per design range trip), and have the same effect as reducing the structural weight by 900 kg. From these examples it can be concluded that seemingly small reductions in cruise drag can have a significant effect on aircraft design. In December 2009, the first flight of the latest airliner to enter production, the Boeing 787, took place. Comparing the B787 of Fig. 1.4 to the DC-2 of Fig. 1.2, it becomes obvious that in an astonishing 75 years, nothing changed in terms of aircraft configuration. However, a number of differences can be distinguished. The wings and tail surfaces of the B787 are very slender and highly tapered, lowering induced and wave drag. Additionally, the nose section of the B787 is more aerodynamically shaped and the landing gear is fully retractable. These differences result in a much higher aerodynamic efficiency compared to that of the DC-2. There are a number of factors that lead to the superiority of the current aircraft compared to the DC-2 from 1935. One is the advancement in aircraft materials. The slender wings of the B787 could simply not have been produced 70 years ago. Other significant factors are the development of swept wings and the introduction of jet engines, which allowed for efficient high speed flight. A final reason is the availability of computer power. The design of the DC-2 was purely driven by the experience of the designers, validated by wind tunnel testing. However, at present, computer algorithms are capable of accurately modeling the viscous airflow around an aircraft, creating the possibility to numerically optimize aircraft shapes. This started in the 1970s and 1980s, with relatively simple computers optimizing simple shapes and has culminated into powerful computers capable of optimizing complete aircraft. Figure 1.5 shows the development of CFD usage at Boeing. The first codes that were used (TA80, TA139-201, TA176/217) were all based on linear flow theory and hence relatively simple. The FLEXSTAB system combined these linear flow solvers with structural and dynamic analysis methods, making aeroelastic analysis possible. The first non-linear flow solvers used by Boeing were the well-known FLO27/28
1.1. BACKGROUND
7
Army Fan-In-Wing Contract Technology
Boeing Tools
Wave Drag Mach Box Woodward
TA80 FLEXSTAB TA139/201 TA176/217 A230
1965
Boeing Aircraft
FLEXSTAB Contract
1970
SST
1975
PANAIR FLO22/ Cartesian 27/28 Grid Technology A411
TLNS3D
Unstructured CFL3D TLNS3D-MB Adaptive Grid OVERFLOW N-S
TRANAIR A555 A488 OptimiP582 A502 TRANAIR A588 A619 zation ZEUS GGNS 1980
767 757 KC-135R
1985
737-300
1990
1995
777
737NG
2000
787
Figure 1.5: Development of CFD usage at Boeing
codes. As the codes became more and more advanced, their complexity increased significantly. As a result, only expert users could set up the models and run the simulations. It wasn’t until the early 1980s, when NASA developed PANAIR (known at Boeing as A502), before more user-friendly codes became available. The A502 solver was less sensitive to specific details of the model and was therefore much easier to use. In fact, the A502 code is still being used today in preliminary design studies. The first user-friendly non-linear solver was called TRANAIR. It was a full potential flow code coupled to a boundary layer solver. The first Euler code to enter the theatre was known as A588, which was also coupled to a boundary layer code. The next logical step was implementing the full Navier-Stokes equations, which was first done in the TLNS3D and CFL3D codes. The first Navier-Stokes code to use unstructured grids was introduced fairly recently, and has been used for designing the Boeing 787. Several interesting papers have been written on the role of Computational Fluid Dynamics (CFD) in aircraft design [29, 31, 83, 84]. Even though CFD gives current aircraft designers a huge advantage, most aircraft still adhere to the standard configuration of a fuselage with wings and a tail. As this configuration has been around for decades, only small, incremental improvements are still achievable. It is often said that the fact that the same configuration has been in use for all this time is proof of its excellence. However, by saying this one ignores the fact that the requirements for aircraft have changed significantly since the 1930s. The emphasis no longer lies on range and speed, but instead on fuel burn and noise. So, what are the odds that with a different set of requirements the optimal solution is the same? It seems that we are just trying to create design variants around a well-established design concept and this cannot be a sustainable way forward if we want to keep improving aircraft performance. In order to improve on this well-established design, incremental changes will not be sufficient. As noted by Green [21], a radical step change will need to be taken to find a design that fulfills current requirements in terms of fuel burn and noise. Such a step change, however, requires more than just a change in geometry. It requires out of the box thinking in all aspects of the design process. The aerodynamic solver needs to be capable of calculating the flow around unconventional aircraft shapes and
8
CHAPTER 1. INTRODUCTION
Design variables
Geometry generation
Aerodynamic analysis
Aerodynamic analysis results
Desired aerodynamic performance achieved?
yes
no Optimization
Figure 1.6: Optimization flow
the optimization algorithm has to be able to handle a more complex design space. In the end, however, the biggest burden will lie on the method of parameterization. If the geometric representation cannot account for the unconventional shapes required, then the whole design process will be extremely limited. Changing the shape of an aircraft is not the only way to improve its aerodynamic efficiency. Techniques such as active flow control (e.g., by boundary layer suction [70] or using plasma actuators [16]) could significantly reduce the drag of an aircraft. However, these methods are usually superimposed on existing designs and therefore will not lead to new aircraft configurations. The objective of the work presented in this thesis was to develop a new way of aerodynamically optimizing aircraft shapes, fulfilling the above mentioned requirements for the design of novel aircraft geometries. The work will include all aspects of the aerodynamic design process (see Fig. 1.6), such as aerodynamic analysis and optimization, but the focus will lie on the parameterization of aircraft shapes. In order to generate a three-dimensional shape, it is possible to simply take a row of two-dimensional shapes and interpolate the points in between. Historically, most aircraft wings have been defined as such an interpolated stack of airfoils. This is quite an easy solution, also considering the production of a wing, which can be done in a similar fashion, with ribs representing the airfoil sections. However, as the wing shape becomes more complex, more airfoil sections have to be defined to describe the wing, making the method less and less efficient. With more sophisticated production techniques and computers available today, it has become feasible to represent the entire three-dimensional shape as a mathematical surface. Various geometry discretization techniques are available to this end, such as polynomials, splines, and NURBS [51]. It is often difficult, however, to use these techniques in the early design phase, because they simply provide too much freedom. Additionally, they often have no physical link to common aerodynamic shapes. As a solution to these problems, a novel parameterization technique was developed, called the Class-Shape-Refinement-Transformation (CSRT) method. This CSRT method, which is described in detail in Chapter 3, consists of a shape function based on Bernstein polynomials and a refinement function based on B-splines, allowing for a two-step optimization approach.
1.2. CURRENT PRACTICE
1.2 1.2.1
9
Current practice Parameterization
In order to optimize any shape, it is necessary to express it with a finite number of variables, preferably as few as possible to minimize computation cost during optimization. This process is called “shape parameterization” or “geometric representation” and is an essential part of any shape optimization problem. Interestingly, parameterization standards differ for different disciplines. In Computer Aided Design (CAD) extensive use is made of non-uniform rational B-splines (NURBS). NURBS curves and surfaces provide a compact and intuitive geometric representation that can be efficiently handled by computer programs. Their main downside is their rather complex definition, which often makes it difficult to generate the same geometry using two different software packages. In Finite Element Analysis (FEA) simple polynomials are often used to prevent this problem. However, these polynomials are not as intuitive and efficient as other means of parameterization. When it comes to Computational Fluid Dynamics (CFD), Kulfan [39] composed a list of desirable features that any geometric representation technique should possess, including features such as smoothness, mathematical efficiency, size of design space, and ability to handle constraints. She then observed that most conventional methods (such as the ones described in Section 2.2) fail to meet all of these criteria and presented a novel technique that did fulfill the complete set of requirements. This Class-Shape-Transformation (CST) technique combines a class function, representing a specific class of shapes, and a shape function, which defines the deviation from the class function. The shape function is defined by Bernstein polynomials and their coefficients form the design variables. Samareh [68] has evaluated a number of shape parameterization techniques and concluded that polynomial and spline representations are the best approach in terms of: i. Ability to handle surfaces ii. Consistency across disciplines iii. Ability to handle large geometry changes iv. Availability of sensitivity derivatives Note the importance of the third item, the ability to handle large geometry changes, for dealing with novel aircraft configurations. Although the CST method fulfills the requirements mentioned by Samareh [68], it has one important shortcoming: it cannot cope efficiently with local shape changes. It appears that Samareh’s list of conditions is not complete and a fifth condition needs to be added: v. Ability to handle local shape changes Local control also provides a larger range of possibilities for handling geometric constraints. A number of different techniques for implementing constraint handling have been investigated in literature. Juhsz and Hoffmann [34] made use of the knot values
10
CHAPTER 1. INTRODUCTION
of the B-spline to create a virtual envelope that the shape cannot cross. A similar technique has been presented by Berglund et al. [6]. Albeit elegant, these techniques only work for external constraints, whereas aircraft constraints are often internal (e.g., volume constraints w.r.t. fuel tanks, passenger cabin, baggage storage). Tuli et al. handled geometric constraints by splitting up the B-spline curve into separate segments. Pourazady and Xu [62] and Celniker and Welch [11] ascribed physical properties, such as stretching stiffness and bending stiffness, to the B-spline curves. Satisfying the constraints was accomplished by minimizing the deformation energy of the curves. This technique was shown to be especially useful in design cases where the physical properties actually play a role. Painchaud-Ouellet et al. [60] used a constraint on the thickness to guarantee a certain internal volume. However, the latter two methods do not allow the geometric constraints to be translated to bounds on the design variables. This translation is, however, desirable, because it limits the required computational time.
1.2.2
Aerodynamic analysis
Apart from the parameterization technique, a flow solver is needed that is appropriate for the design problem at hand. The choice of solver is usually a trade-off between the required complexity of the computed flow and the available resources. ReynoldsAveraged-Navier-Stokes (RANS) solvers use a time-average of the full Navier-Stokes equations and provide the most accurate description of a flow field that can feasibly be used for complete aircraft design applications. The Navier-Stokes equations can be expressed as follows: ( ) ∂v ρ + v · ∇v = −∇p + Fvisc (v) (1.2) ∂t where ρ is the air density, v is the velocity field, p is the pressure field and Fvisc are the viscous forces. In order to derive the RANS equations, the velocity and pressure fields are decomposed into a time-averaged component (¯ v, p¯) and a fluctuating component (v′ , p′ ), as described below. v
=
v ¯ + v′
(1.3)
′
p = p¯ + p
This leads to the following definition of the RANS equations: ( ) ∂¯ v ρ +v ¯ · ∇¯ v + T = −∇¯ p + Fvisc (¯ v) ∂t with ∂u′2
∂x′ ′ ∂u v T= ∂x
∂u′ w′ ∂x
∂u′ v ′ ∂y ∂v ′2 ∂y ∂v ′ w′ ∂y
∂u′ w′ ∂z ∂v ′ w′ ∂z ∂w′2 ∂z
(1.4)
(1.5)
1.2. CURRENT PRACTICE
11
Taking the Navier-Stokes equations and neglecting viscosity leads to the Euler equations, which can be written as: ( ) ∂v ρ + v · ∇v = −∇p (1.6) ∂t Not taking into account viscosity means that the only drag sources that Euler solvers can compute are induced drag and wave drag, unless they are coupled to boundary layer equations. These equations can be derived from the Navier-Stokes equations using an order of magnitude analysis. The assumption of irrotational flow leads to the full potential equation: 2
2
2
2
∂ Φ (1 − Mx2 ) ∂∂xΦ2 + (1 − My2 ) ∂∂yΦ2 + (1 − Mz2 ) ∂∂zΦ2 − 2Mx My ∂x∂y ∂2Φ ∂2Φ −2My Mz ∂y∂z − 2Mx Mz ∂x∂z = 0
(1.7)
where Φ is the velocity potential and M is the local Mach number of the flow. By definition, the velocity components can be derived from the velocity potential using the following relationship: v = ∇Φ
(1.8)
Linearizing the full potential equation finally results in: 2 (1 − M∞ )
∂2Φ ∂2Φ ∂2Φ + + =0 ∂x2 ∂y 2 ∂z 2
(1.9)
Equation (1.9) is aptly called the linearized potential equation. Because of the assumptions made to derive it, its application is limited to the subsonic and supersonic flow regimes. For the low fidelity framework presented in Chapter 6 a panel code was selected that makes use of the linearized potential flow equation. For the high fidelity framework presented in the same chapter, an Euler code was used.
1.2.3
Optimization
Thirdly, an optimization algorithm needs to be selected. When using a gradient-based algorithm, the most important consideration is how to compute the derivatives of the flow variables with respect to the design parameters. Some methods, such as finite differencing and the complex-step method, are easy to implement, but they both require one run of the objective function for each design parameter. A much faster way to obtain the gradients is to make use of the adjoint method. Although this method is fast, a disadvantage of any gradient-based optimization method is that it is only capable of finding local optima. Depending on the design space, this could make the final result highly dependent on the initial geometry, as is illustrated in Fig. 1.7. Additionally, gradient-based algorithms usually cannot handle discrete variables. To overcome these problems, a global optimization method could be employed. A disadvantage of such a method, however, is that it usually requires a high number of function evaluations. Since in aerodynamic optimization function evaluations are
12
CHAPTER 1. INTRODUCTION D
?
A C B
Figure 1.7: Example of local maxima
Figure 1.8: Example of inability to guarantee global optimum
often expensive in terms of computation time, this is highly undesirable. Another problem of global optimization methods is that they cannot guarantee that a global optimum has been found. This is illustrated in Fig. 1.8. A solution to the former drawback is the use of response surfaces [19]. A response surface is essentially a surrogate model of the design space, created using the information from a limited number of sample points, as illustrated in Fig. 1.9 . Evaluating a design point using this surrogate model is much less expensive than performing the actual aerodynamic analysis, hence making global optimization feasible. Choosing the right set of sample points and generating an accurate response surface can be done in many ways. Techniques for achieving this are widely described in literature, e.g., by Forrester et al. [18,19] and Jones et al. [32]. There is also a wide choice in global optimization algorithms. Two well-established techniques are genetic algorithms and particle swarm optimization, see e.g., Green and Cheng [22], Venter and Sobieszczanski-Sobieski [86], and Keane [35].
1.3
CSRT method
The main focus of this thesis will be on a novel parameterization technique called the Class-Shape-Refinement-Transformation (CSRT) method and its implementation into two different design frameworks. The CSRT method is an extension to Kulfan’s Class-Shape-Transformation (CST) method. Figure 1.10 shows an example of a local deformation. In order to model this shape using the CST method, it would be ideal to only increase the order of the Bernstein polynomials that govern the part of the curve around the deformed region. However, each of the basis functions that make up a Bernstein curve are global functions of the same order. Thus, to capture the deformation, the order of the Bernstein polynomials must be increased. Figure 1.11 shows that this is not a very efficient process. The dots represent the local deformation of Fig. 1.10 and the graphs show approximations obtained for different orders of Bernstein polynomials (Nvar -1). The solution for
1.3. CSRT METHOD
13
Figure 1.9: Example of a response surface with the black dots representing the sample points
0.2
0.15
0.1
0.05
0
−0.05 0
0.2
0.4
0.6
0.8
1
Figure 1.10: Example of a local deformation
14
CHAPTER 1. INTRODUCTION 0.2
0.2
0.15
0.15
0.1
0.1 Nvar = 4
0.05
0.05
Nvar = 6 N
0
var
=8
0 N
var
−0.05 0
0.2
0.4
= 14
0.6
Nvar = 20
0.8
1
−0.05 0
0.2
0.4
0.6
0.8
1
Figure 1.11: Capturing a local defor- Figure 1.12: Capturing a local deformation using Bernstein polynomials for mation using Bernstein polynomials for Nvar = 4, 6, 8, 14 Nvar = 20
0.2
0.2
0.15
0.15
0.1
0.1 N
var
0.05
N
var
0
N
var
−0.05 0
0.2
=8
0.4
=4
0.05 =6
0 N
Nvar = 14
0.6
var
0.8
1
−0.05 0
0.2
0.4
0.6
= 20
0.8
1
Figure 1.13: Capturing a local deforma- Figure 1.14: Capturing a local deformation using a B-spline for Nvar = 4, 6, 8, tion using a B-spline for Nvar = 20 14
Nvar = 20 is shown in Fig. 1.12. Even though the deformation itself is modeled relatively well, the rest of the curve is highly oscillatory. In order to allow for local modifications, the CST method was expanded by the author with a refinement function, based on B-splines. Because B-spline basis functions change the shape of the curve only locally, it takes considerably fewer variables to capture a local deformation. Figures 1.13 and 1.14 show B-spline approximations of the deformation from Fig. 1.10 for different values of Nvar . The B-spline approximations are clearly much better than their Bernstein equivalents. For Nvar = 20, the deformation itself is closely followed by the curve and the rest is almost completely flat. By combining Bernstein polynomials with a B-spline curve, the advantages of both types of curves can be utilized. This will be elaborated on in Chapter 3.
1.4. THESIS OVERVIEW
1.4
15
Thesis overview
Chapter 2 will give an overview of existing parameterization methods, with special emphasis on the Class-Shape-Transformation (CST) method in Section 2.3. The Class-Shape-Refinement-Transformation (CSRT) technique will be introduced and thoroughly analyzed in Chapter 3. Chapter 4 will show the rationale behind handling volume constraints with the CSRT method. In Chapter 5 two design frameworks (low and high fidelity) will be introduced that will serve as the basis for the test cases presented in Chapter 6. Finally, the conclusions and recommendations are presented in Chapter 7.
16
CHAPTER 1. INTRODUCTION
CHAPTER
2
Parameterization
“For academic men to be happy, the universe would have to take shape. All of philosophy has no other goal: it is a matter of giving a frock coat to what is, a mathematical frock coat. On the other hand, affirming that the universe resembles nothing and is only formless amounts to saying that the universe is something like a spider or spit.” - Georges Bataille (1897 - 1962)
2.1. BACKGROUND
19
This chapter will focus on parameterization techniques. Section 2.1 will start with an overview of different methods that have been used throughout the past few decades. Subsequently, Section 2.2 describes in more detail some of the parameterization methods that are most common in (preliminary) aircraft design. Finally, Section 2.3 will introduce the Class-Shape-Transformation (CST) technique, which forms the basis of the Class-Shape-Refinement-Transformation (CSRT) method presented in Chapter 3.
2.1
Background
Today, optimization tools are readily available on the market in a wide variety. Their application in the design phase of complex products, however, requires proper mathematical descriptions of such products and preferably also easy access to derivatives of such models with respect to predefined design variables. Since many design variables are commonly geometry related and often the analysis tools invoked during the optimization need geometry and the derivatives thereof, the geometric part of the product model has received much attention in literature. The basic problem is best understood through the following sample equation, applicable to, for example, an objective function related to fluid flows: ∂J ∂Rf ∂Rs ∂Rg ∂J = ∂x ∂Rf ∂Rs ∂Rg ∂x
(2.1)
The optimizer needs the derivative of the objective function, J, with respect to the design variables, x. To calculate this it needs the derivative of the objective function with respect to the field discretization (Rf ), the derivative of the field discretiziation with respect to the product surface model (Rs ), the derivative of the surface model with respect to the geometric model (Rg ), and finally the derivative of the geometric model with respect to the design variables, x. It is very important that the parametric model is both flexible and consistently coupled to the discrete models used by the different analysis tools. Haftka and Grandhi [24] and Ding [15] provided surveys of shape optimization up to 1986. A good survey of parametric models upto the year 2001 for use in combined structural and aerodynamic optimization is given by Samareh [68]. His survey focuses on “the suitability of available techniques for multidisciplinary applications of complex configurations using high-fidelity analysis tools such as computational fluid dynamics and computational structural mechanics”. He concludes that a successful parameterization process must i) be automated, ii) provide consistent geometry changes across all disciplines, iii) provide sensitivity derivatives (preferably analytical), iv) fit into the product development cycle times, v) have a direct connection to the CAD systems used for design, and vi) produce a compact and effective set of design variables for the solution time to be feasible. The parameterization techniques can be divided into eight categories: basis vector, domain element, partial differential equation (PDE), discrete, polynomial and spline, CAD-based, analytical, and freeform deformation (FFD).
20
CHAPTER 2. PARAMETERIZATION
The FFD technique as presented by Barr [4], based on the parallelepiped, is very efficient and easy to implement. This technique is suitable for local and global deformation. The only disadvantage is that the use of the parallelepiped limits the topology of deformation. To alleviate this disadvantage, Sederberg and Parry proposed to use non-parallelepiped objects [72]. Another popular method to define FFD is to use trivariate parametric volumes. Sederberg and Parry used a Bezier volume. Coquillart [12] extended the Bezier parallelepiped to a non-parallelepiped cubic Bezier volume. This idea has been further generalized to NURBS volumes by Lamousin and Waggenspack [42]. Samareh contributed to the improvement of the CAD-based parametric models [68]. The CAD-based approach uses the commercial feature-based solid modeling CAD systems to create the design surfaces. To parameterize an existing model is still a challenging task in today’s CAD systems, and the models created are not always good enough for automatic grid-generation tools. A solution to these challenges was proposed by Samareh and consists of two basic concepts: i) parameterizing the shape perturbations rather than the geometry itself and ii) performing the shape deformation by means of the soft object animation algorithms (SOA) used in computer graphics. The SOA algorithms are used to modify the wing twist and shear distribution. They are modified versions of the algorithms of Watt and Watt [89]. In a later paper Samareh proposed a modified FFD technique [69]. The proposed technique is an alternative shape parameterization technique to a trivariate volume technique. It retains the flexibility and freedom of trivariate volumes for CFD shape optimization, but uses a bivariate surface representation. This reduces the number of design variables by an order of magnitude. The analytical sensitivity derivatives are independent of the design variables and are easily computed for use in a gradientbased optimization. Mousavi presented a survey of shape parameterization techniques in 2007 [52]. This paper applies three shape parameterization techniques: mesh points, B-spline surfaces, and the Class-Shape-Transformation (CST) method [40]. Morris [49] presented a generic domain element shape parameterization method which allows optimization using variable levels of geometric description, ranging from gross planform changes to detailed minor surface changes. The parameterization technique uses radial basis functions (RBFs) to provide an interpolation between a domain element and the surface geometry. Complex, multi-element, two- and threedimensional geometries can be parameterized in a hierarchical manner allowing freeform design of arbitrary geometries with only a few intuitive design variables. The interpolation is also extended to the CFD mesh. Sevant et al. [73] investigated a hierarchical approach to aerodynamic design, using the PDE method for surface generation. With this method he was able to parameterize complex geometries using only a limited number of variables, whilst still maintaining a wide range of possible shapes.
2.2. COMMON METHODS
21
Figure 2.1: The discrete approach
2.2 2.2.1
Common methods Discrete approach
One of the most straight-forward ways of describing a two-dimensional shape is by using the coordinates of a set of points on its boundary and connecting these points with straight lines. This is called the discrete approach and it is illustrated in Fig. 2.1. The discrete approach is very easy to implement, because it allows a direct connection with any type of grid. One major disadvantage of this approach is that a large number of points are required to generate a smooth shape. C 2 continuity is desired for two reasons: to minimize numerical errors in the flow calculation and to prevent separation of the actual flow. However, even reaching C 1 continuity is not possible when describing a curved shape by straight line segments. It can only be approached. Furthermore, once an approximately smooth shape has been reached, it is even more difficult, if not impossible, to maintain this smoothness during the optimization process. This would require a complex set of constraints, thereby nullifying the ease of implementation. Another reason why the discrete approach is hardly ever used to represent wing shapes, is that the round nose of an airfoil is a non-analytic function with infinite derivatives at the nose, and therefore cannot be represented by a discrete function.
2.2.2
Analytical approach
In the analytical approach, a set of perturbation or shape functions is added linearly to a baseline shape. These perturbation functions can be defined either numerically or analytically. An example of the analytical approach with analytically defined shape functions was given by Hicks and Henne [26]. They introduced the Hicks-Henne bump functions, defined as follows: [ ( log(0.5) )]t2 , fj (x) = sin πx log(t1 )
0≤x≤1
(2.2)
where t1 and t2 can be varied to fit the design problem. Figure 2.2 shows five HicksHenne functions for t2 = 5. The weights αj are the design variables. The perturbation functions fj (x) can now be added to the baseline shape as follows: y(x) = ybasis +
M ∑ j=1
α ¯ j fj (x)
(2.3)
22
CHAPTER 2. PARAMETERIZATION f1(x)
1
f2(x)
f3(x)
f (x)
f (x)
4
5
0.5
0 0
0.5
1
Figure 2.2: Hicks-Henne functions for t2 = 5 Table 2.1: Definitions of typical shape functions
Chebyshev
First kind T0 (x) = 1 T1 (x) = x Tn+1 (x) = 2xTn (x) − Tn−1 (x)
Legen-
P0 (x) = 1
dre
P1 (x) = x Pn+1 (x) =
(2n+1)xPn (x)−nPn−1 (x) n+1
Second kind U0 (x) = 1 U1 (x) = 2x Un+1 (x) = 2xUn (x)( − U) n−1 (x) 1+x ln 1−x ( ) 1+x −1 Q1 (x) = x2 ln 1−x
Q0 (x) =
Qn+1 (x) =
1 2
(2n+1)xQn (x)−nQn−1 (x) n+1
The Hicks-Henne functions are just one example of the many analytical shape functions available. Figure 2.3 shows some more examples: Chebyshev polynomials of the first (a) and second (b) kind and Legendre polynomials, also of the first (c) and second (d) kind [3]. Their (recursive) definitions are given in Table 2.1. A disadvantage of the analytical approach is that changing a shape by changing one or more of the perturbation function variables is often non-intuitive, making this method suitable for fully automated optimization processes, but less suitable for “hands on” aircraft design.
2.2.3
Geometric approach
A third approach is to describe a shape using physical parameters, such as leading edge radius, thickness-to-chord ratio or trailing edge angle. Figure 2.4 shows an airfoil generated using a geometric method introduced by Sobieczky, called PARSEC [74]. The airfoil has been described using 11 basic design parameters. This approach is clearly very intuitive. A disadvantage of the method is that
2.2. COMMON METHODS
23
Figure 2.3: Chebyshev polynomials of the first (a) and second (b) kind and Legendre functions of the first (c) and second (d) kind
Figure 2.4: The geometric approach
24
CHAPTER 2. PARAMETERIZATION
Figure 2.5: Airfoil described by a Bezier representation
each shape change has to be translated into a physical property of the airfoil. This is not always straight-forward and might compromise the design freedom. Additionally, the design variables might not be directly related to the properties of the flow, which could complicate the matter further.
2.2.4
Polynomial and spline approaches
The required number of variables needed to generate a smooth shape can be greatly reduced by using a polynomial or spline representation. A polynomial can be expressed in its standard power basis form (Eq. (2.4)), but in that case it will be very difficult to predict how a change in the coefficient vector c¯i will influence the overall shape of the curve.
y(x) =
p ∑
c¯i xi
(2.4)
i=0
A better way to represent a curve would be to use the Bezier representation, which is defined as follows:
y(x) =
p ∑
P¯i Bi,p (x)
(2.5)
i=0
where P¯i is a vector of control points and Bi,p (x) are Bernstein polynomials of degree p, the definition of which is given in Section 2.3.1. Because the resulting curve will closely follow the control polygon defined by the control points of P¯i , this definition is much more intuitive than the power basis form. An example of an airfoil described by a Bezier representation is given in Fig. 2.5. A disadvantage of using Bezier curves is that for more complex shapes, the required degree of the Bernstein polynomials increases significantly, making it very inefficient to compute the curve. This follows from the fact that each Bernstein polynomial has a global influence on the shape of the curve. A solution to this drawback is the use of a string of lower order curves, called a B-spline. A B-spline curve can be described as follows:
2.3. CST METHOD
25
p(u) =
n ∑
¯ i Ni,p (u) P
(2.6)
i=0
Compared to Eq. (2.5), the Bernstein polynomials Bi,p have been replaced by a set of B-spline basis functions Ni,p and the Bernstein coefficient vector P¯i by a B-spline ¯ i . The parametric variable u is used instead of the Cartesian cocontrol polygon P ordinate x. Also note that the number of control points n + 1 is independent of the order p − 1 of the individual Bezier curves. A more thorough definition of B-splines and their basis functions will be given in Section 3.1.
2.3
CST method
The previous section listed a number of parameterization techniques and some of their pros and cons. This section focuses on a technique introduced by Kulfan [41] called the Class-Shape-Transformation (CST) method. This method combines an analytical “class function” with a parametric “shape function”. The class function describes a basic class of shapes and the shape function describes the permutation around this basic shape. This way, characteristic features that are specific to a certain class of shapes, such as the round nose and sharp trailing edge of an airfoil, can be captured in the class function. Consequently, the same shape function can be used regardless of the class of shapes to be optimized.
2.3.1
CST in two dimensions
Symbolically the CST method can be described as follows: N1 ζ(ψ) = CN 2 (ψ) · S(ψ) + ψ · ζT
(2.7)
where ζ = z/c, ψ = x/c and ζT = ∆ζT E /c. ζT is the dimensionless trailing edge N1 thickness, which is assumed zero from here on for reasons of simplicity. CN 2 (ψ) and S(ψ) represent the class and shape function respectively. The class function is defined as follows: N1 N1 CN (1 − ψ)N 2 2 (ψ) = (ψ)
(2.8)
0.5 For a NACA type round nose and pointed √ aft end airfoil the C1.0 (ψ) class function can be used. In other words: C(ψ) = ψ(1 − ψ). For the shape function, Kulfan uses a Bernstein polynomials representation. In accordance with the definition given in Section 2.2.4, the shape function is then given by:
S(ψ) =
p ∑
P¯i Bi,p (ψ)
i=0
The Bernstein polynomial terms are typically defined as:
(2.9)
26
CHAPTER 2. PARAMETERIZATION 1.2 1 0.8 0.6 0.4 0.2 0
0
0.2
0.4
0.6
0.8
1
Figure 2.6: Bernstein polynomial terms for p = 5
Bi,p =
( ) p i ψ (1 − ψ)p−i i
(2.10)
where p is the order of the Bernstein polynomial. Comparing Eq. (2.9) and (2.10) reveals that the order of the polynomial is always equal to the number of control points (or length of the vector P¯i ) minus one. As an example, for six control points (p = 5) the following terms are found: (1 − ψ)5
B0,5
=
B1,5 B2,5 B3,5
= 5(1 − ψ)4 ψ = 10(1 − ψ)3 ψ 2 = 10(1 − ψ)2 ψ 3
B4,5 B5,5
= 5(1 − ψ)ψ 4 = ψ5
(2.11)
Plotting these polynomial terms results in Fig. 2.6. A special feature of Bernstein polynomials is that they provide a partition of unity, i.e., they always add up to one. This is shown by the thick black line in Fig. 2.6. This means that if the coefficients of the vector P¯i are all equal to one and N 1 and N 2 are chosen equal to 0.5 and 1.0 respectively, then Eq. (2.9) will also be equal to one and Eq. (2.7) will simplify to: √ ζ(ψ) = ψ(1 − ψ) (2.12) 0.5 airfoil, which is plotted in Fig. 2.7. This equation describes the so-called unit C1.0 Changes in the coefficient vector P¯i will lead to a variation in shape around the unit airfoil. Figure 2.8 shows the unit airfoil and the airfoil deformed using a coefficient vector of {1 2 0.5 2 0.5 1}. A special feature of the CST method is that the values of the shape function at ψ = 0 and ψ = 1 are directly related to the leading edge nose radius and the boat-tail
2.3. CST METHOD
27
0.4
0.6
0.3 0.4
0.2 0.2
0.1 0
0
−0.1
−0.2
−0.2 −0.4
−0.3 −0.4 0
0.2
0.4
0.6
0.8
1
0.5 Figure 2.7: Unit airfoil for C1.0
0
0.2
0.4
0.6
0.8
1
Figure 2.8: Airfoil shapes deformed using P¯i = {1 1 1 1 1 1} (dashed line) and P¯i = {1 2 0.5 2 0.5 1} (solid line)
angle respectively [41], according to Eq. (2.13) and (2.14). In the example of Fig. 2.8 this means that both the leading edge nose radius and the boat-tail angle should be the same for both airfoils, since the first and last values of the P¯i vector are equal for both cases. Looking at the figure, this indeed seems to be the case. S(0) =
√ 2(RLE /C)
(2.13)
∆zT E (2.14) c The advantage of using the CST method in combination with Bernstein polynomials is that every possible coefficient vector P¯i leads to a realistic, smooth shape. However, any change in P¯i will lead to a global change in geometry, thereby not allowing for any local modifications. S(1) = tan β +
2.3.2
CST in three dimensions
With a few simple modifications, the CST method can also be used to model threedimensional shapes. Kulfan [41] gives the three-dimensional equivalent of Eq. (2.7) as: N1 ζ(ψ, η) = ζN (η) + CN 2 (ψ) · S(ψ, η) + ψ[ζT (η) − ∆αT (η)]
(2.15)
where ζN (η) is the non-dimensional local wing shear and ∆αT (η) is the wing twist, both functions of the non-dimensional semi-span station η = 2y b . The class function N1 CN is defined according to Eq. (2.8). 2 A Bernstein polynomials representation, as explained in the previous section, can also be used in three dimensions. In order to guarantee a smooth surface, Bernstein
28
CHAPTER 2. PARAMETERIZATION
1 0.8 0.6 0.4 0.2 0 0 1 0.5
0.5 1
0
Figure 2.9: 3D Bernstein polynomial terms for px = py = 5
polynomials are necessary in both x- and y-directions. This leads to the following definition of the shape function: S(ψ, η) =
py px ∑ ∑
y x P¯i,j · Bi,p (ψ) · Bj,p (η) x y
(2.16)
i=0 j=0
with ( ) px i ψ (1 − ψ)px −i i
(2.17)
( ) py j = η (1 − η)py −j j
(2.18)
x Bi,p (ψ) = x
and y Bj,p (η) y
Comparing Eq. (2.16) to (2.9) shows that the coefficient vector P¯i changed to a coefficient matrix P¯i,j of size (px + 1) × (py + 1). The Bernstein polynomials now represent a set of surfaces instead of curves, as is illustrated in Fig. 2.9. These surfaces again add up to one everywhere in the domain. Because the class function is only defined in the x-direction, this implies that if all coefficients in the matrix P¯i,j are equal to one, the cross section of the surface will be the unit airfoil as shown in Fig. 2.7 and it will be constant in spanwise direction. This is shown in Fig. 2.10. To show that the surface generated by this method is always smooth, Fig. 2.11 shows the surface resulting from randomly selecting the coefficients of the matrix P¯i,j .
2.3. CST METHOD
29
0.5
0
−0.5 1 0.5
0.2
0
0
0.4
0.6
0.8
1
0.5 Figure 2.10: Unit airfoil in 3D for C1.0
0.4 0.2 0 −0.2 −0.4 1
1 0.5
0.5 0
0
Figure 2.11: Surface resulting from randomly selecting the coefficients of P¯i,j
30
CHAPTER 2. PARAMETERIZATION
In conclusion, the CST method combines the advantages of the analytical, geometric, and spline approaches as mentioned in Section 2.2. A solution to some of the disadvantages of the CST method will be presented in the next chapter.
CHAPTER
3
CSRT method
“The same refinement which brings us new pleasures, exposes us to new pains.” - Edward Bulwer-Lytton (1803 - 1873)
3.1. B-SPLINES
33
This chapter will introduce a B-spline extension to the CST method called the ClassShape-Refinement-Transformation (CSRT) method. Section 3.1 will give an extensive definition of B-splines and B-spline surfaces, as well as a worked out example. The refinement function will be introduced in Section 3.2 and Section 3.3 will present a parametric study investigating the performance of the CSRT method. Finally, Section 3.4 will give the definitions of a number of global shape parameters as used throughout this thesis.
3.1 3.1.1
B-splines B-spline curves
A B-spline is basically a string of lower order curves. Because of its piece-wise nature, B-spline curves are capable of efficiently modeling local deformations of a shape. As was shown in Section 2.2.4, a B-spline is defined as follows [51]: p(u) =
n ∑
¯ i Ni,p (u) P
(3.1)
i=0
Note that the number of control points is independent of the degree of the basis function polynomials. The basis functions are defined iteratively as follows: Ni,1 (u) = 1 =0
ti ≤ u < ti+1
if
otherwise
(3.2)
and Ni,p (u) =
(u − ti )Ni,p−1 (u) (ti+p − u)Ni+1,p−1 (u) + ti+p−1 − ti ti+p − ti+1
(3.3)
where ti are the knot values that relate the parametric variable u ∈ [0, n − p + 2⟩ to ¯ i . The knot values are defined as follows: the control points P ti = 0
if
in
(3.4)
For a curve with six control points (n = 5) and second order basis functions (p = 3), the knot vector according to Eq. (3.4) is: T = {0, 0, 0, 1, 2, 3, 4, 4, 4} This leads to the following values of Ni,1 :
(3.5)
34
CHAPTER 3. CSRT METHOD
N0,1 (u) = 1 =0
for otherwise
u≡0
N1,1 (u) = 1 =0
for otherwise
u≡0
N2,1 (u) = 1 =0
for
0≤u >
ζ 2 ψ1 C(ψ2 ) ζ2 (ψ1 − 1) C(ψ2 )
(4.4)
where C(ψ1 ) and C(ψ2 ) are the values of the class function evaluated at ψ1 and ψ2 respectively. As long as these conditions are satisfied by the box constraint, the fitted curve will be positive and the method described in this section will work. Any box constraint that does not satisfy these conditions will have to be fitted into one that does. As can be seen in Fig. 4.3 there is still some vertical space between the class function and the box constraint. This can be solved by adjusting the control points of the refinement function B-spline such that the gap between the curve and the box is minimized. Any gradient-based optimization algorithm can be used to accomplish this. The result is shown in Fig. 4.5(a). The straight solid line represents the shape function Bernstein curve. It is straight because only the B-spline coefficients have been varied and all Bernstein coefficients are equal to one. The crossed line represents the control polygon of the refinement function B-spline. It is theoretically known which region of the curve is influenced by the fitting process. This is the solid part of the curve. The question is now how to derive the lower bounds on the CSRT coefficients from this “minimum shape”. Changing the coefficients of the shape function has a global effect on the shape of the curve. This means that decreasing the values of the coefficients during the optimization, without changing the refinement function, is very likely to lead to a violation of the box constraint. For this reason, the Bernstein coefficients of the curve in Fig. 4.5(a) form the minimum bounds. The situation for the refinement function is slightly more complicated. If the shape function has not been changed then the B-spline coefficients governing the part of the curve above the box constraint also form the minimum bounds. However, if some of the Bernstein coefficients have been changed, then a situation can occur where some of the B-spline coefficients can be decreased without violating the constraint. Such a situation is shown in Fig. 4.5(b). The question now is what set of constraints can be applied to the B-spline coefficients that allows them to be decreased without violating the box constraint. One way to do this is to divide the values of the B-spline coefficients by the magnitude of the Bernstein polynomial at their specific locations. The number of B-spline control points is usually limited, therefore this does not require a large number of function evaluations. Since the shape of the B-spline curve is generally quite smooth, it follows the control polygon closely enough to assume that in that case the product of shape function and refinement function gives a good approximation of the curve that fits tightly around the box constraint. Figure 4.5(c) shows that after applying this procedure, the curve indeed does not cross the constraint. Finally, it has to be noted that the procedure outlined above does not work exclusively for a quadrilateral box. The
4.2. VOLUME CONSTRAINT HANDLING IN 2D
63
ζ [~]
1.5 Bernstein curve B−spline control polygon Fixed part of CSRT curve Variable part of CSRT curve
1 0.5 0 0
0.2
0.4
0.6
0.8
1
Relative chordwise position, ψ [~]
a)
ζ [~]
1.5 1 0.5 0 0
0.2
0.4
0.6
0.8
1
Relative chordwise position, ψ [~]
b)
ζ [~]
1.5 1 0.5 0 0
c)
0.2
0.4
0.6
0.8
1
Relative chordwise position, ψ [~] Figure 4.5: Fitting the curve around the constraint
64
CHAPTER 4. VOLUME CONSTRAINTS Box constraint dimensions
Does box constraint satisfy Eq. 4.4?
yes
Fit class function through box constraint vertices
Are all Bernstein coefficients equal to 1?
no
no
Fit constraint in box that satisfies Eq. 4.4
Divide B-spline coefficients by Bernstein polynomial
yes
Find B-spline coefficients governing curve above box constraint
Fit curve to box constraint
Bounds on CSRT coefficients
Figure 4.6: Finding the bounds on the CSRT coefficients for a given box constraint
curve could also be fit around more complex shapes. The whole process of finding the constraints on the CSRT coefficients is summarized in Fig. 4.6. A specific characteristic of a B-spline is that its coefficients only influence a specific part of the curve. As a result of this, the positions of the control points in Fig. 4.5 that govern the part of the curve that is above the box constraint, represent their lower bounds. The control points that do not influence the part of the curve above the constraint are free to move up and down during the optimization and hence are unconstrained. It should be noted, however, that the original CST method allows for the leading edge nose radius and the trailing edge boat-tail angle to be expressed analytically in terms of the Bernstein coefficients. In order to retain this useful property, the B-spline coefficients that control the leading and trailing edge regions should remain equal to 1.
4.3
Volume constraint handling in 3D
For dealing with complete aircraft, the constraint handling method has to be extended to three dimensions. As a three-dimensional example, a box constraint is defined as two quadrilateral planes connected by straight lines. Both planes are parallel to the xz-plane and they are located at η = 0 and η = 1. This is shown in Fig. 4.7. For three dimensions, two different class functions need to be determined. One for each side of the box constraint. In between the two sides, the values of a and b (see Eq. (4.1)) are linearly varied between the two extremes. This gives a reasonably good fit around the edges of the box, as is shown in Fig. 4.8. A more accurate, but more laborious way to do it is to determine the values of a and b for each point along the span. Similarly to the two-dimensional case, there is still quite a large space between the surface and the box constraint. This can be eliminated by modifying the positions of the control points above the box constraint. Firstly the x-positions of the control points have to be vertically aligned to the vertices of the box. This can easily be done, like is shown in Fig. 4.9. This figure shows the B-spline grid (red lines), the wing (transparent) and the box constraint (grey).
4.3. VOLUME CONSTRAINT HANDLING IN 3D
65
1 0.8 0.6
ζ [~]
0.4 0.2 0 0
0 0.2
0.2 0.4
0.4
ψ [~]
0.6
0.6 0.8
0.8 1
η [~]
1
Figure 4.7: Three-dimensional box constraint
0.5
0.4
0.3
ζ [~] 0.2
0.1
0 1
η [~]
0.5 0
1
0.8
0.4
0.6
0.2
0
ψ [~]
Figure 4.8: Making the surface go through the edges of the box constraint
66
CHAPTER 4. VOLUME CONSTRAINTS
0 0.1 0.2 0.3 0.4
η [~] 0.5 0.6 0.7 0.8 0.9 1
0
0.2
0.4
0.6
0.8
1
ψ [~]
Figure 4.9: Aligning the B-spline control points with the box constraint
The control points that are allowed to vary in the minimization routine are chosen such that they form a quadrilateral with an equal amount of control points on opposite sides. This is depicted in Fig. 4.10. One vertex of the quadrilateral lies directly above the most forward point of the box constraint, the opposite vertex lies above the most rearward point of the box constraint. After performing a linear optimization routine to minimize the gap between the surface and the box constraint, the wing looks like the one in Fig. 4.11. The way in which the lower bounds on the CSRT coefficients are determined from this “minimum shape” is the same as for the two-dimensional case. If the shape function is unaltered, then the B-spline coefficients that govern the part of the surface above the box constraint form the minimum bounds. If, however, the coefficients of the Bernstein surface of the shape function are modified, then this change should be compensated for by moving the control points of the B-spline surface of the refinement function. This is again similar to the method for two dimensions.
4.3. VOLUME CONSTRAINT HANDLING IN 3D
67
0 0.1 0.2 0.3 0.4
η [~] 0.5 0.6 0.7 0.8 0.9 1 0
0.2
0.4
0.6
0.8
1
ψ [~]
Figure 4.10: The region of the surface fitted to the constraint
0.5
ζ [~]
1 0 1
0.8 0.6 0.8 0.6
0.4 0.4
0.2 0.2
0
0
η [~]
ψ [~]
Figure 4.11: The wing fitted around the box constraint
68
CHAPTER 4. VOLUME CONSTRAINTS
CHAPTER
5
Design frameworks
“I believe that good design is magical and not to be lightly tinkered with. The difference between a great design and a lousy one is in the meshing of the thousand details that either fit or don’t, and the spirit of the passionate intellect that has tied them together, or tried.” - Ted Nelson (1937 - present)
5.1. LOW FIDELITY DESIGN FRAMEWORK
71 Gradients
FD loop
CSRT coefficients VSAERO
Optimization
Shape generation
Mesh generator
Flow solver
Postprocessor
Aerodynamic coefficients
Convergence?
yes
no
Figure 5.1: Low fidelity framework with direct finite differencing
This chapter will introduce two different design frameworks that were used for the test cases presented in Chapter 6. The first framework is a low fidelity optimization tool coupling a panel method program (VSAERO) to a global optimization method (via a surrogate model). The global optimization algorithm employed was particle swarm optimization (PSO). To act as a benchmark, one optimization run was also performed using VSAERO in combination with finite differencing (FD) and an active set algorithm. The second framework is a high-fidelity tool using an in-house Euler solver and a gradient-based optimizer. This framework was first used for a twodimensional case, using FD and an active set algorithm. Subsequently, two threedimensional cases were performed with the same framework, but employing an adjoint solver and a trust region reflective algorithm. The first framework will be described in Section 5.1. The theory behind panel methods in general and VSAERO in particular will be presented, as well as response surface and global optimization methods. Section 5.2 will introduce the second framework, with special emphasis on the theory and validation of the Euler and adjoint solvers and gradient-based optimization algorithms.
5.1
Low fidelity design framework
Figure 5.1 shows a schematic of the low fidelity framework, using finite differencing for sensitivity analysis. The “FD loop” indicated in the figure has to be executed once for each design variable to find the gradients. The flow schematic of the test cases using a response surface is shown in Fig. 5.2. In this scheme, the “sampling loop” (and hence the flow solver) has to be executed once for each sample point used to generate the response surface. Once the response surface has been generated, the “FD/PSO loop” has to be executed once for each design variable in case of finite differencing and once for each particle in case of PSO.
5.1.1
Panel methods (VSAERO)
Panel methods are still widely used in a aerodynamic design problems. In his paper, Hess [25] discusses various theoretical formulations and numerical implementations used to solve the potential flow problem. Most linear panel method programs assume an incompressible, inviscid flow and hence are not capable of capturing all aspects of
72
CHAPTER 5. DESIGN FRAMEWORKS Sample points
Sampling loop
CSRT coefficients VSAERO
Hypercube generation
Shape generation
Mesh generator
Flow solver
Postprocessor
Response surface generation
Aerodynamic coefficients
Gradients / Particles
FD / PSO loop
Response surface evaluation
Optimization
Aerodynamic coefficients
Convergence?
yes
no
Figure 5.2: Low fidelity framework with response surface
a flow field. However, they can be improved by correcting for compressibility and by coupling the potential flow equations to a set of boundary layer equations. For the potential flow cases in this thesis, the Karman-Tsien rule was used for compressibility correction and transition was computed by using the Granville criterion for transition prediction [2]. As a result of these improvements, the lift, induced drag, and friction drag could be determined with reasonable accuracy [75]. Finally, VSAERO is capable of handling complex three-dimensional shapes relatively fast and is thus very suitable for use in any preliminary design application. A number of papers on the theory behind VSAERO have been published by Nathman [53–56].
Flow equations The well known differential form of the continuity equation is: Dρ + ρ∇ · v = 0 Dt Assuming an incompressible fluid reduces Eq. (5.1) to:
(5.1)
∂u ∂v ∂w + + =0 (5.2) ∂x ∂y ∂z Furthermore, for irrotational flow, a velocity potential Φ can be defined such that: ∇·v =
v = ∇Φ
(5.3)
Substituting Eq. (5.3) into (5.2) gives: ∇ · v = ∇ · ∇Φ = ∇2 Φ = 0
(5.4)
Equation (5.4) is Laplace’s equation. Many different control volume types can be envisioned to which Eq. (5.4) can be applied. VSAERO uses the control volume
5.1. LOW FIDELITY DESIGN FRAMEWORK
73
S∞
Φ
P
r
n Φi V∞
SW
SB
Figure 5.3: Control volume
shown in Fig. 5.3. The boundary of this control volume consists of three parts: the body surface SB , the wake surface SW and the outer boundary S∞ . The potential is split into two parts: the potential in the flow field Φ and the potential in the inner region Φi . Applying Green’s theorem to Eq. (5.4) for both the flow field and the inner region results in: ( ) ∫∫ 1 1 ΦP = (Φ − Φi )n · ∇ dS 4π r SB +SW +S∞ ( ) ∫∫ 1 1 − n · (∇Φ − ∇Φi )dS 4π r SB +SW +S∞
(5.5)
This expression describes the velocity potential in any point P expressed in terms of surface integrals. The first integral represents the local jump in potential and the second integral represents the local jump in the normal component of the velocity. Expressed in terms of singularities, the first integral represents the doublet distribution, while the second integral represents the source distribution. Because the surface S∞ is assumed to be infinitely far away from the body, the surface singularities no longer influence the flow field there. Hence the part of the integral evaluated over S∞ reduces to Φ∞P , i.e., the velocity potential at point P due to the freestream. Since the wake is assumed to be thin, no jump in normal velocity is possible across it. This implies that the second integral of Eq. (5.5) evaluated over SW is zero. For
74
CHAPTER 5. DESIGN FRAMEWORKS
the first integral, combining the upper and lower wake elements causes the interior potential term to disappear and what remains is: 1 4π
∫∫ (ΦU − ΦL )n · ∇ SW
( ) 1 dSW r
(5.6)
Equation (5.5) now becomes:
ΦP = +
1 4π
1 4π ∫∫
( ) ∫∫ ( ) 1 1 1 (Φ − Φi )n · ∇ dS − n · (∇Φ − ∇Φi )dS r 4π r SB SB ( ) 1 (ΦU − ΦL )n · ∇ dSW + Φ∞P r SW ∫∫
(5.7)
Boundary conditions To solve the flow, boundary conditions are required that prescribe the velocity on the surface. This can be done either directly by applying the Neumann condition (∇Φ · n specified on the boundary) or indirectly by applying the Dirichlet condition (Φ specified on the boundary). VSAERO makes use of the former. In VSAERO the user has the possibility to specify a certain normal velocity Vn to model for example boundary layer suction. Furthermore, the boundary layer displacement δ ∗ can be modeled by assuming a certain transpiration velocity. Finally, a normal velocity component due to rotation ω can be specified. This implies the following boundary condition: ∇Φ · n = −Vn + (V∞ − ω × rg ) · n +
∂ (V δ ∗ ) ∂s
(5.8)
For the wake a boundary condition is also required. From the fact that the wake surface cannot support a load, it can be derived that the following two conditions hold for the wake panels: V·n=0
and
V · ∇(ΦU − ΦL ) = 0
(5.9)
Singularity model Equation (5.7) does not have a unique solution. There is still an infinite number of combinations of source and doublet distributions that satisfies the flow problem. A number of techniques can be employed to find a unique combination of singularities. It is for example possible to specify the source distribution and solve for the doublet distribution, or vice versa. Another possibility is to apply boundary conditions to the inner flow, either to the velocity or the velocity potential. VSAERO specifies the velocity potential of the inner flow to be equal to the freestream velocity potential, i.e., Φi ≡ Φ∞ . An alternative is to assume zero velocity potential in the inner region, but this method is more sensitive to bad paneling.
5.1. LOW FIDELITY DESIGN FRAMEWORK
75
Numerical procedure The body and wake surfaces are divided into N and NW panels respectively. On each panel there is one point, called the collocation point, where the boundary condition is applied and the velocity potential resulting from the singularity distribution is computed. The influence of the singularity distribution of each panel on each other panel can then be calculated. For each collocation point P the system of equations then becomes: N ∑
Ck µk +
NW ∑
Cl µl +
l=1
k=1
N ∑
Bk σk = 0
(5.10)
k=1
where C and B are the influence coefficients for the doublets and sources respectively. Since the doublet strength in the wake is equal to the difference in doublet strength between the upper and lower trailing edge panels, the influence of the wake element becomes: Ct µt = Ct (µr − µs )
(5.11)
where µr and µs are the upper and lower trailing edge doublet strengths and µt is the wake doublet strength. This coefficient Ct should now be added to the Ck coefficients if the panel is at the trailing edge, such that: Ak = Ck
if panel is not at the trailing edge
Ak = Ck ± Ct
if panel is at the trailing edge
(5.12) (5.13)
Now Eq. (5.10) can be written as: N ∑ k=1
Ak µk = −
N ∑
Bk σk
(5.14)
k=1
This system of N equations and N unknowns can now be solved. Once the singularity distribution is known, the complete velocity potential field can be computed. From this, other quantities, like the velocity and pressure fields, can be derived.
Compressibility In the derivation of Laplace’s equation it was assumed that the flow is incompressible. However, some methods are available to approximate the effects of compressibility on the solution. The best method available in VSAERO for external flows is the KarmanTsien method. This method relates the compressible and incompressible velocities as follows: ( V )inc (1 − λ) V = V∞ V∞ 1 − λ( VV∞ )2inc
(5.15)
76
CHAPTER 5. DESIGN FRAMEWORKS
2 V/V
∞
1 0 2 1 (V / V∞)inc
0
0
0.6
0.4
0.2
1
0.8
M
∞
Figure 5.4: Karman-Tsien compressibility correction
where λ=
M2 √ ∞ 2 )2 (1 + 1 − M∞
(5.16)
These relationships are shown graphically in Fig. 5.4. To summarize, the potential flow calculations result in the incompressible flow solution (V /V∞ )inc . Substituting this solution into Eq. (5.15) and computing λ from Eq. (5.16) gives the flow solution corrected for compressibility, V /V∞ .
Viscosity In order to model all sources of aircraft drag, viscous effects need to be included. Since the potential flow model itself neglects viscosity, it needs to be coupled to a boundary layer method. VSAERO uses a method developed by Thwaites [82], which was later adapted by Curle [13]. A short description of this method is given in this section. First, the integral momentum equation is derived. The momentum equation in x-direction for the boundary layer is given by: ∂u ∂u ∂u −1 ∂p µ ∂ 2 u +u +w = + ∂t ∂x ∂z ρ ∂x ρ ∂z 2
(5.17)
Integrating this equation with respect to z from z = 0 to z = δ gives: ∂ ∂t
∫
∫
δ
δ
udz + 0
u 0
∂u dz + ∂x
∫
δ
w 0
−1 ∂u dz = ∂z ρ
∫
δ 0
µ ∂p dz + ∂x ρ
∫
δ 0
∂2u dz ∂z 2
(5.18)
This equation can be rewritten as (see [2]): d dx
∫
δ 0
dUe (Ue − u)udz + dx
∫
δ
(Ue − u)dz = 0
τw ρ
(5.19)
5.1. LOW FIDELITY DESIGN FRAMEWORK
77
where Ue is the inviscid velocity at the body surface and τw is the shear stress. The displacement thickness δ ∗ and the momentum thickness θ were defined by Lighthill [43] as follows: δ∗ ≡
∫
δ 0
∫
δ
θ≡
u )dz Ue
(5.20)
u u ) dz Ue Ue
(5.21)
(1 − (1 −
0
Using Eq. (5.20) and (5.21), Eq. (5.19) can be rewritten as: dθ 1 dUe τw + (2θ + δ ∗ ) = dx Ue dx ρUe2
(5.22)
Writing Eq. (5.22) in non-dimensional form results in: dθ θ dUe Cf + (H + 2) = dx Ue dx 2
(5.23) ∗
δ w where Cf = 1 τρU 2 is the friction coefficient and H = θ is the shape factor. In the e 2 method of Thwaites an extra term is added to Eq. (5.23), representing streamline convergence: ( ) dθ θ dUe Cf 1 dr + (H + 2) = −θ (5.24) dx Ue dx 2 r dx
This equation is then written in the form: 1 2 (r K/U ′ )′ = L/U r2
(5.25)
2
θ with K = ( c¯/Re )U ′ , L = 2(ℓ − K(H + 2)) and ℓ = ( Uθ )( ∂U ∂y )y=0 . Curle determined the following empirical function for L:
L = F0 (K) − µG0 (K) = 0.45 − 6K + g(K, µ)
(5.26)
F0 and G0 are stored in tabular form. Substituting Eq. (5.26) into (5.25) and integrating gives: ∫ c¯/Re r2 θ2 = (0.45 + g)r2 U 5 dx (5.27) U6 At the start of the process g is chosen to be zero. Equation (5.27) can then be used to calculate θ, after which K and µ can be computed. Once K and µ are known, Eq. (5.26) can be used to determine g. This leads to a new value of θ from Eq. (5.27) and the process can be repeated until it has converged. Once the momentum thickness θ is known, the skin friction coefficient Cf can be calculated from:
78
CHAPTER 5. DESIGN FRAMEWORKS
Cf = 2ℓ/Reθ
(5.28)
The value of ℓ can be found from the following relationship: ℓ2 = F1 (K) − µG1 (K)
(5.29)
where F1 and G1 are stored in tabular form, similar to F0 and G0 . By using the shape factor H the displacement thickness δ ∗ can be calculated. This displacement thickness can be interpreted as the distance by which an external flow streamline is displaced because of the presence of a boundary layer. Multiplying δ ∗ by the inviscid velocity on the body and differentiating this w.r.t. a streamline gives the transpiration velocity required to model the boundary layer. ∂ (V δ ∗ ) (5.30) ∂s The method of modeling viscosity as presented in this section is just one of many possibilities. Several other techniques can be found in literature [23, 33]. Vtransp =
Coupling CSRT method to VSAERO The panel method flow solver used in the low fidelity framework, VSAERO, uses input files to define the flow problem. These input files contain the coordinates of cross-sections of the geometry at specified span wise locations. They also contain instructions on how to build the surface and wake meshes, as well as information about the flow and boundary layer. Because a different input file is required for each flow case, the geometry definition can be performed entirely outside the flow solver. For this work, Matlab was used to build the CSRT geometries and write the VSAERO input files. After performing the flow calculation, VSAERO produces an output file that can again be read by the Matlab algorithm [1]. Note that one of the reasons of introducing the CSRT method, as mentioned in Chapter 1, was to get rid of the idea of defining a wing geometry by just stacking different two-dimensional airfoil sections. However, this is how a wing needs to be defined in VSAERO. Hence after the CSRT geometry has been created, it has to be sliced up into two-dimensional cross-sections again.
5.1.2
Global optimization
Most global optimization methods are not gradient-based, but instead search the design space in a stochastic way. Two examples of search methods are genetic algorithms (GAs) [22] that are based on Darwin’s evolution principle [14] and particle swarm optimization (PSO) [27, 86] which is based on the behavior of a flock of birds or insects. The main advantage of global optimization methods is their capability of finding a global, or near global, solution. Additionally they are able to handle discrete and mixed-integer optimization problems [44], which often occur in preliminary aircraft design. Finally, global optimization methods are considerably easier to program
5.1. LOW FIDELITY DESIGN FRAMEWORK
79
than gradient-based methods, especially when using existing aerodynamic simulation software. The design framework described in this section makes use of particle swarm optimization. PSO is based on swarming theory and models a number of “particles” that migrate towards the optimum solution. These particles decide where to go based on their own fitness (L/D value), as well as the fitness of the swarm as a whole. Each particle at any moment in time represents one design solution. The initial swarm of particles is randomly created within the given bounds, after which both position and velocity of the particles are updated according to: i xik+1 = xik + vk+1 ∆t
i vk+1 = wvki + c1 r1
(pg − xik ) (pi − xik ) + c2 r2 k ∆t ∆t
(5.31)
(5.32)
where w is the inertia of the particle, xik and vki are the position and velocity of particle i at iteration k, r1 and r2 are random numbers between 0 and 1, and c1 and c2 are the “trust parameters”. The first trust parameter is a measure of how much confidence a particle has in itself, while the second parameter indicates its confidence in the swarm as a whole. The parameter pi gives the best position found so far by particle i, while pgk is the best position in the whole swarm at time k. The PSO process and its optimization variables are shown in Fig. 5.5. In order to decide when to stop the optimization process, a stopping criterion is required. Particle swarm optimization offers a number of possibilities, such as: i. Termination tolerance on best solution ii. Termination tolerance on number of iterations iii. Termination tolerance on migration However, no matter what stopping criterion is used, an optimum can never be guaranteed. In gradient-based optimization, the algorithm always moves towards a local optimum, but the Karush-Kuhn-Tucker condition [38] can at least guarantee that the local optimum has been found. A global optimization algorithm has the potential to move towards a global optimum, but it cannot guarantee whether an optimum (local or global) has been found.
5.1.3
Response surface methods
Creating a response surface consists of two main steps. The first step is to choose a sample set that allows the generation of an accurate response surface, while minimizing the number of points required. The second step is to generate the response surface itself. These two steps will be explained in the following sections.
80
CHAPTER 5. DESIGN FRAMEWORKS
Optimization variables Number of particles
Initial swarm
Geometry generation (CSRT)
Number of iterations
Min/max inertia weight
Min/max acceleration coefficients
Update position of each particle
Evaluation
Save best position of each particle
Save best position of the swarm
VSAERO input file generator (MATLAB)
VSAERO
Figure 5.5: PSO scheme
VSAERO output file interpreter (MATLAB)
L/D
81
x3
x2
5.1. LOW FIDELITY DESIGN FRAMEWORK
x1
x2
x3
x3
x1
x1
x2
Figure 5.6: Three-dimensional Latin hypercube
Hypercube sampling process There are various methods available for creating a sample set for response surface generation. The method used in the low fidelity framework is called the hypercube sampling process [19]. The idea behind this method is that when projected onto the axes of each individual variable, the distribution of the sample points should be as uniform as possible. Ideally, no points should overlap in any of the projections. To discretize the problem, the design space is divided up into a number of equal sized boxes (see Fig. 5.6). If the design problem involves only discrete variables, then each box could represent the actual value of the variable. If the variables are continuous, then the value of the variable should merely be contained somewhere within the box. It is now possible to create a set of sample points such that along any of the variable axes only one sample point will be found. Such a set is called a Latin hypercube. An example of a three-dimensional Latin hypercube with 10 sample points is shown in Fig. 5.6.
82
CHAPTER 5. DESIGN FRAMEWORKS
For an n-dimensional design space and m sample points, the number of possible Latin hypercubes is given by (m!)n−1 . This means that the example in Fig. 5.6 is just one out of 1.3 · 1013 possibilities for n = 3 and m = 10. In order to pick the best sample set, some criterion to evaluate uniformity is required. Morris and Mitchell [50] present such a criterion. They define a distance vector d = (d1 , d2 , ..., dm ) (sorted from smallest to largest) in which the distance elements are given by the p-norm of the space: (
dp x(i1 ) , x(i2 )
)
1/p k p ∑ (i1 ) (i ) = xj − xj 2
(5.33)
j=1
Additionally, an index vector J = (J1 , J2 , ..., Jm ) is defined in which Jj is the number of point pairs separated by distance dj . They then derive the following uniformity criterion: Φq (X) =
m ∑
1/q Jj d−q j
(5.34)
j=1
The smaller the value of Φq the more uniform the distribution of the sample points is. In order to find the best sample set X, Eq. (5.34) is minimized using an evolutionary operation algorithm [19].
Generate/evaluate (Kriging) response surface Again, various methods are available to generate a response surface. This section will focus on the Kriging method [32], which was used in this optimization framework. The Kriging method uses the following basis functions to build the response surface: k pj ∑ (i) ψ (i) = exp − θj xj − xj (5.35) j=1
From Eq. (5.35) it follows that each variable has its own value of θj . This value basically determines how far the influence of a specific sample point reaches. In other words, variables that have little influence on the objective function will have a low value of θj , while for variables with a large influence θj will be high. The value of pj typically lies between 1 and 2. The next step is to define the maximum likelihood estimates for the mean value µ and the square of the standard deviation σ 2 . The following definitions can be found in [19]: 1T Ψ−1 y 1T Ψ−1 1
(5.36)
(y − 1µ)T Ψ−1 (y − 1µ) n
(5.37)
µ ˆ=
σ ˆ2 =
5.2. HIGH FIDELITY DESIGN FRAMEWORK
83
where Ψ is the correlation matrix of the sample points and y are the values of the objective function at the sample locations X. Equations (5.36) and (5.37) can be combined into a single concentrated ln-likelihood function as follows: n 1 ln(ˆ σ 2 ) − ln |Ψ| (5.38) 2 2 It is now a matter of finding the values of θ and p for which Eq. (5.38) is minimized. A maximum likelihood estimate for the predicted value yˆ at a specific location x can now be calculated using the following relationship: ln(L) ≈ −
yˆ(x) = µ ˆ + Ψ T Ψ−1 (y − 1ˆ µ)
(5.39)
where Ψ is a vector of correlations between the original and predicted data. For the complete derivation of the above equations, the reader is referred to [19].
Validation of the CSRT response surface In order to test the response surface method, a sample set of 500 analysis runs was performed on a three-dimensional wing using the panel method program VSAERO. Ten design points were randomly generated and evaluated using both VSAERO and the response surface. The result is shown in Fig. 5.7. Considering the fact that only 500 sample runs were performed to recreate a 12-dimensional design space, the results shown in the figure are very good. The difference between the VSAERO and response surface evaluations is never more than 4 %. Additionally, and more importantly, the gradients and minimum and maximum values of the design points coincide. It seems that the response surface generated for this test case is a realistic model of the design space.
5.2
High fidelity design framework
The flow schematic of the high fidelity framework using finite differencing for sensitivity analysis is shown in Fig. 5.8. It is essentially the same as the low fidelity framework of Fig. 5.1, except for the flow solver. Again, the “FD loop” has to be executed once for each design variable. This is no longer the case if the adjoint code, as implemented in the Euler solver, is used to find the gradients of the design variables. This is shown in Fig. 5.9. Contrary to the other design frameworks, this scheme does not contain any loops depending on the number of variables. It only requires one run of the flow solver and one run of the adjoint solver to find the required gradients.
5.2.1
Euler method
The Euler equations are a model of inviscid compressible flow. Solvers based on the Euler equations can predict induced drag and wave drag. Hence, to the extent for which viscosity is not very important, they can be used for aerodynamic analysis and design in transonic flow conditions. When viscosity becomes important suitable
84
CHAPTER 5. DESIGN FRAMEWORKS 18.0 M = 0.10 R e = 2.33 * 10 6 α = 2.00 o
R es pons e s urface V S AE R O
L ift-to-drag ratio, L /D [~]
17.0
16.0
15.0
14.0
13.0
12.0 1
2
3
4
5
6
7
8
9
10
Des ign point
Figure 5.7: Test results of the response surface
Gradients
FD loop
CSRT coefficients Euler code
Optimization
READMESH
PREPROC
DEFORM
FLOW
POSTPROC
Aerodynamic coefficients
Convergence?
yes
no
Mesh
Figure 5.8: High fidelity framework with finite differencing
CSRT coefficients Euler code
Optimization
READMESH
PREPROC
DEFORM
FLOW
POSTPROC
Aerodynamic coefficients
Convergence?
no
ADJOINT
GRADIENT
Gradients
Mesh
Figure 5.9: High fidelity framework with adjoint sensitivity analysis
yes
5.2. HIGH FIDELITY DESIGN FRAMEWORK
85
solvers should be used (e.g., solvers based on the Reynolds-averaged Navier-Stokes equations, or solvers based on the Euler equations coupled with boundary layer methods). The Euler solver used in the high fidelity framework is based on an unstructured finite-volume formulation on the median-dual mesh. The unstructured mesh provides more design flexibility than a Cartesian mesh, such as the one used by Rodriguez [66]. The Euler solver used for this thesis, which has been implemented by Carpentieri, is described in detail elsewhere [8,10]. In the following only a short description is given. Figure 5.8 shows the main modules of the flow solver: the mesh reader, the preprocessor, the deformation module, the flow solver, and the postprocessor. The mesh reader reads the mesh from a mesh-specific format. The preprocessor processes the mesh and creates the median-dual discretization. The mesh deformation algorithm can deform the mesh according to a specified boundary deformation. The postprocessor computes forces along the boundary. In this framework, the CSRT method is applied within the mesh deformation module. The mesh reader reads a mesh of the class function, which is subsequently deformed using the shape and refinement functions. A validation of the implementation of the CSRT method into the Euler code will be presented further on in this section.
Flow equations On a bounded domain Ω with boundary ∂Ω, the Euler equations can be written in integral form as: ∫ I d u dΩ + F(u) · n dΓ = 0 (5.40) dt V ∂V where u is given by: u = [ρ, ρw, ρet ]T
(5.41)
with w the velocity vector. The flux vector is defined as: F(u) = [ρwT , ρwwT + pI, (p + ρet )wT ]T
(5.42)
In order to close the system, a relationship is needed between the total energy et and the pressure p. This is given by the perfect gas equation, which reads: p = (γ − 1)ρ
( ) |w|2 et − 2
(5.43)
Equation (5.40) can be applied to a control volume. In the Euler solver, the control volumes are formed by a median-dual type mesh, which can be created around any unstructured grid. This is shown in Fig. 5.10, where the gray regions indicate the control volumes around the nodes. A more detailed look at one of the control volumes is given in Fig. 5.11. The arrows represent the flux Fij along the control volume interface ∂Vi . The figure also shows the boundary ∂Ω of the complete domain Ω. The
86
CHAPTER 5. DESIGN FRAMEWORKS
Figure 5.10: Three control volumes in the median-dual mesh
∂Ω Ω F ij · n ij
∂ Vi
Vi
Figure 5.11: Detailed look at a control volume
5.2. HIGH FIDELITY DESIGN FRAMEWORK
87
solver makes use of a node-centered approach, meaning that the unknowns are stored at the centers of the control volumes. For the generic control volume i, such as the one in Fig. 5.11, the discretized version of Eq. (5.40) looks as follows: Vi
∑ dui + Fij · nij + Fi · nBi = 0 dt
(5.44)
j=1,Ni
with ∫
∫
nij =
n dΓ
and
nBi =
∂Vij
n dΓ
(5.45)
∂Vi ∩∂Ω
The second integral is only non-zero if the boundary ∂Vi of cell i is part of the domain boundary ∂Ω. Introducing numerical fluxes leads to: Vi
∑ dui + Φij + Φbc i =0 dt
(5.46)
j=1,Ni
The last two terms of Eq. (5.46) are called the residual ri .
Boundary conditions When dealing with a wall or symmetry boundary condition, the velocity component normal to this boundary needs to be zero. As was the case for panel methods (see Section 5.1.1), there is more than one way to accomplish this. The so-called “weak” way is to set the normal velocity to zero in the normal projection of the Euler flux of Eq. (5.46): T Φbc i = [0, pnBi , 0]
(5.47)
Alternatively, one could impose the “strong” boundary conditions directly on the flow variables, i.e.: wi · nBi = 0
(5.48)
Numerical procedure This section explains the numerical procedure used by the Euler solver. Using compact notation, Eq. (5.46) can be written as: D
dU +R=0 dt
(5.49)
where U = [u1 , u2 , ..., uN ]T is the vector of conservative variables and R = [r1 , r2 , ..., rN ]T is the residual vector. The control volumes Vi are stored in the diagonal matrix D. Equation (5.49) is solved by means of an implicit pseudo-time stepping scheme, which can be written as:
88
CHAPTER 5. DESIGN FRAMEWORKS
Dt (Un+1 − Un ) + R(Un+1 ) = 0
(5.50)
where Dt is a diagonal matrix that contains the cell volumes divided by the local time steps, Vi /∆ti . Applying a linear expansion to the second term of Eq. (5.50) gives: R(Un+1 ) ≈ R(Un ) +
∂R n n+1 (U − Un ) ∂U
(5.51)
Substituting Eq. (5.51) into (5.50) gives the equations for the pseudo-time step n: (
∂R Dt + ∂U
)n (Un+1 − Un ) = −Rn
(5.52)
This equation can be solved iteratively by means of a linear solver. It can be written in the form Mz = b, with: ( M=
Dt +
)n
∂R ∂U
, z = Un+1 − Un and b = −Rn
(5.53)
The solver performs linear iteration with Gauss-Seidel preconditioning [9]. The iteration can be written as: k K zk+1 = zk + ∆z, with ∆z = P−1 and RK L = b − Mz M RL
(5.54)
where RL is the residual of the linear system and PM is the preconditioner.
Deformation algorithm The Euler solver also contains a mesh deformation algorithm, which can be employed if a change in geometry is required. This is a very useful feature, since rebuilding the mesh after every optimization step would make the computation times unacceptably large. The deformation algorithm makes use of the spring analogy. With this strategy, an initial solution is created where some of the mesh points on the boundary (e.g., of an airfoil or wing) are displaced. The displacements of the rest of the mesh points are iteratively solved using the following Jacobi iterations: ∑Ni ∆Xm+1 = i
m j=1 kij ∆Xj , ∑Ni j=1 kij
(i = 1, N ; i ∈ / B)
(5.55)
where Ni is the number of nearest neighbors of the node i and kij is the stiffness associated with the edge ij. The stiffness values are chosen to be equal to the inverse of the edge length. When the deformations have been calculated, they are added to the existing mesh as follows: Xnew = X + ∆X
(5.56)
5.2. HIGH FIDELITY DESIGN FRAMEWORK
89
Coupling CSRT method to Euler code As was shown in Section 5.2, the Euler code with adjoint solver consists of seven main executable files: i. READMESH, which reads the mesh file ii. PREPROC3D, which processes the mesh file and creates the median-dual discretization iii. DEFORM, which deforms the mesh according to a boundary deformation specified using the CSRT method iv. FLOW, which runs the flow solver v. POSTPROC, which postprocesses the results vi. ADJOINT, which runs the adjoint solver vii. GRADIENT, which computes the gradients A Matlab batch file was created that executes these files in order. Before it runs DEFORM, it writes a data file that includes the CSRT coefficients of the wing to be modeled. DEFORM then reads this data file and deforms the mesh accordingly. After all executables have been run, the results are presented graphically. Originally, the Euler code made use of (modified) Chebyshev polynomials to model deformations. These polynomials were added to the original wing to form the desired wing shape. This is different from the CSRT technique, in which the deformation functions are multiplied by the original wing to form the new wing shape. To prevent having to extensively modify the code, the following procedure was implemented. First, a mesh is created of the class function (the unit airfoil/wing). Running READMESH loads this mesh into the code. In DEFORM, the CSRT curve is computed by multiplying the (known) class function with the shape and refinement functions. The class function is then subtracted from this curve, leading to the absolute deformation as required by the solver: ∆y(ψ, η) = C(ψ, η) · S(ψ, η) · R(ψ, η) − C(ψ, η)
(5.57)
Validation of CSRT/Euler coupling The framework illustrated in Fig. 5.8 was used to model a NACA0012 airfoil and calculate the flow field for different values of the angle of attack and Mach number. Only the Bernstein polynomials of the shape function were varied to create the NACA0012 airfoil. A high number of shape variables was used (10 for each side), to guarantee an accurate model. Table 5.1 shows the results of this test case, compared to the results calculated by Carpentieri et al. [10] and Viviand [87]. Viviand provided a series of reference cases for inviscid flow field methods.
90
CHAPTER 5. DESIGN FRAMEWORKS
Table 5.1: NACA0012 results compared to Carpentieri et al. [10] and Viviand [87] M∞ [∼]
α [deg]
0.80 0.85 0.95 1.20 1.20
1.25 1.00 0.00 0.00 7.00
cl [∼]
cd [∼]
CSRT
Ref. [10]
Ref. [87]
CSRT
Ref. [10]
Ref. [87]
0.356 0.392 0.000 0.000 0.521
0.349 0.375 0.000 0.000 0.524
0.359 0.353 0.000 0.000 0.521
0.023 0.059 0.109 0.096 0.155
0.022 0.058 0.109 0.095 0.154
0.023 0.055 0.108 0.095 0.154
The table shows that the CSRT model performed similarly (within 5 %) to the code as programmed by Carpentieri. The density and shape of the mesh were approximately the same for all models (∼ 12500 nodes). From this test case it could be concluded that the implementation of the CSRT method into the Euler code was successful. Note that the goal of this verification was not to show that the Euler solver itself is valid, or how well the CSRT method is capable of modeling an airfoil. A thorough validation of the Euler code has been provided by Carpentieri [10] and a parametric study into CSRT performance was presented in Section 3.3.3. The goal of this test case was really to validate the implementation of the CSRT method into the Euler code. The first step in validating the results of the Euler solver for a three-dimensional wing will be a qualitative analysis of the pressure field. Because the goal of this exercise was to validate the correct implementation of the CSRT method into the Euler code (and hence not the performance of the Euler and CSRT codes themselves), a large number of shape variables were used in all directions (10). The geometry selected for this analysis was the ONERA M6 wing [71]. This wing has been specifically designed as a test case for CFD solvers. It has a very simple geometry (constant cross-section) and in transonic conditions a number of typical flow phenomena occur, such as local supersonic flow, shock waves, and turbulent boundary layer separation. Naturally, boundary layer separation cannot be modeled using an Euler solver, but the other phenomena can. What makes this wing even more interesting, is the fact that it is notoriously difficult to model using the CST method. The wing was tested in the ONERA S2MA continuous-flow wind tunnel in France for a number of different flow conditions. Figure 5.12 shows a picture of the ONERA M6 model while it is being tested. For the validation, the wing was modeled and analyzed for two different flow conditions. Figure 5.13(a) shows the pressure distribution on the ONERA M6 wing at a Mach number of 0.84 and an angle of attack of 3.06 degrees, computed using the CSRT method coupled to the in-house Euler solver. Figure 5.13(b) shows the pressure distribution on the same wing in the same conditions, as calculated by Sung and Kwon [81]. The pressure distributions are qualitatively very similar. Both show one shock wave along the leading edge, and another parallel to the trailing edge. The shape of the rest of the isobars also seems to correspond. There appears to be a
5.2. HIGH FIDELITY DESIGN FRAMEWORK
91
Figure 5.12: The ONERA M6 wing is tested in the S2MA wind tunnel
a)
b)
Figure 5.13: Pressure distributions on the ONERA M6 wing from a) the present work and b) Sung and Kwon [81]
difference in the position of the aft shock wave, but that is because both figures have been plotted in different ways. The pressure distribution in Fig. 5.13(b) was plotted on the actual geometry, while the one in Fig. 5.13(a) was plotted on a flat surface. This makes the shock wave in the former figure look more aft. In fact, both shock waves are located around 68 % chord. For a quantitative validation, the results of the ONERA M6 model were compared to the results calculated by Carpentieri [10] and Viviand [87]. Because both references did not document details about the meshes used, the ONERA M6 wing was analyzed using different mesh sizes. Figure 5.14 shows the lift coefficient as a function of the number of mesh elements, as well as the reference values. A trend line is also included in the figure to visualize convergence. The value of the lift coefficient seems to approach the value as mentioned by Carpentieri. Figures 5.15 and 5.16 show the drag coefficient as a function of mesh size, for two different flow conditions. For the first case (M = 0.84, α = 3.06o ), the drag coefficient approaches a value in between those mentioned by Carpentieri and Viviand. For the second case (M = 0.92, α = 0.00o )
92
CHAPTER 5. DESIGN FRAMEWORKS 0.290
Lift coefficient, Cl [~]
0.285 0.280 Carpentieri Viviand CSRT CSRT trend
0.275 0.270 0.265 0.260 0
100000
200000
300000
400000
Number of mesh elements, [~]
Figure 5.14: Lift coefficient as a function of mesh size for M = 0.84, α = 3.06o
0.030
0.030
0.028
0.025 Carpentieri Viviand CSRT CSRT trend
0.023
0.020
Drag coefficient, Cd [~]
Drag coefficient, Cd [~]
0.028
0.026
Carpentieri Viviand CSRT CSRT trend
0.024
0.022
0.018
0.015
0.020 0
100000
200000
300000
0
400000
100000
200000
300000
400000
Number of mesh elements, [~]
Number of mesh elements, [~]
Figure 5.15: Drag coefficient as a function of mesh size for M = 0.84, α = 3.06o
Figure 5.16: Drag coefficient as a function of mesh size for M = 0.92, α = 0.00o
the convergence value of the drag coefficient lies just above that found by Carpentieri. Table 5.2 lists the values of the lift and drag coefficients for both flow cases described above, as computed on the densest mesh, consisting of 64627 nodes and 331345 elements. From these results it was concluded that the implementation of the CSRT method into the Euler code had been successful.
Table 5.2: ONERA M6 results compared to Carpentieri et al. [10] and Viviand [87] M∞ [∼] 0.84 0.92
α [deg] 3.06 0.00
CL [∼]
CD [∼]
CSRT
Ref. [10]
Ref. [87]
CSRT
Ref. [10]
Ref. [87]
0.279 0.000
0.281 0.000
0.284 0.000
0.019 0.025
0.020 0.025
0.017 0.022
5.2. HIGH FIDELITY DESIGN FRAMEWORK
5.2.2
93
Gradient-based optimization
For some of the test cases presented in Chapter 6 the constrained minimization function of Matlab (fmincon) [1] was used. This function will choose an optimization algorithm based on the type of problem and the input provided. If the gradients of the objective function with respect to the design variables are provided, then it will use a trust region reflective algorithm. In most other cases, a medium-scale active set algorithm is employed.
Finite differencing The basis of any gradient-based optimization method is some form of sensitivity analysis. Determining the gradient of the objective function with respect to the design variables can be done in a variety of ways, an overview of which is provided by Newman et al. [58]. The most common way to find the gradient of a function is by finite differencing (FD). In forward finite differencing, the gradient of a function f (x) with respect to a design variable xi is approximated as follows: ∂f (x) f (x + ∆xi ) − f (x) ≈ + O(h) ∂xi ∆x
(5.58)
where ∆xi = ∆x · ˆ ei . Alternative FD methods are backward and central finite differencing. The main advantage of FD is that it is extremely easy to implement. It only requires one base run of the objective function and one additional run for each design variable. This property of FD is, however, also the method’s main disadvantage, since a large number of design variables results in a large number of required runs of the objective function. Performing a large number of function evaluations might be acceptable if the objective function can be computed quickly, but this is usually not the case for aerodynamic flow solvers. Moreover, finite difference derivatives can suffer from truncation and subtractive cancellation errors, which may negatively affect the convergence of the optimization algorithm.
Adjoint method For problems with a large number of design variables, the gradients of the objective function can be determined much faster using the adjoint method, which was pioneered by Jameson [28] and has been used extensively in different optimization frameworks [17, 30, 46, 47, 57]. The workings of the adjoint method are explained briefly in this section as applied to aerodynamic shape optimization. First, let us assume some aerodynamic property J which is a function of the flow variables U and the geometry design variables x: J = J(U, x)
(5.59)
The derivative of J with respect to a specific design variable xi can be written as: dJ ∂J ∂J ∂U = + dxi ∂xi ∂U ∂xi
(5.60)
94
CHAPTER 5. DESIGN FRAMEWORKS
Note that Eq. (5.60) distinguishes between a change in objective function as a result of a variation in the flow solution ∂U and a variation due to the change in geometry ∂xi . In order to solve Eq. (5.59), a relationship between U and xi is needed. Such a relationship is the steady state flow equation, i.e.: R(U, xi ) = 0
(5.61)
Computing the derivative of R with respect to xi gives: dR ∂R ∂R ∂U = + =0 dxi ∂xi ∂U ∂xi
(5.62)
Performing sensitivity analysis using Eq. (5.60) and (5.62) is referred to as solving the “primal problem”. The primal problem, however, has the same disadvantage as FD: it needs to be solved for each design variable xi . Hence for a large number of design variables, it would take a long time to compute all the sensitivities. The adjoint method can be derived by introducing a vector of Lagrange multipliers Λ. Equation (5.62) can be added as a constraint to the sensitivity to obtain: ( ) dJ ∂J ∂J ∂U ∂R ∂R ∂U = + −Λ + dxi ∂xi ∂U ∂xi ∂xi ∂U ∂xi ( ) (5.63) ∂J ∂R ∂J ∂R ∂U = −Λ + −Λ ∂xi ∂xi ∂U ∂U ∂xi The vector of Lagrange multipliers can be chosen to satisfy the following adjoint equation: ∂R ∂J = (5.64) ∂U ∂U Substituting Eq. (5.64) into Eq. (5.63) results in the elimination of the last two terms and hence: Λ
dJ ∂J ∂R = −Λ dxi ∂xi ∂xi
(5.65)
Equations (5.64) and (5.65) represent the “dual problem”. The main advantage of the dual problem over the primal problem is that the former only requires solving as many equations as there are flow variables. For most aerodynamic optimization problems, this number is a lot lower than the number of design variables. Hence, using the dual problem can dramatically decrease the time required to compute the gradients. The numerical procedure used to solve the primal and dual problems is similar to that used for the flow solver, hence that part of the theory will be omitted here. The implicit pseudo-time stepping scheme for the adjoint equations, equivalent to Eq. (5.52), reads: ( )n ( ) ∂R ∂R ∂J n+1 n Dt + (Λ −Λ )=− Λ − (5.66) ∂U ∂U ∂U
5.2. HIGH FIDELITY DESIGN FRAMEWORK
95
Carpentieri describes in much more detail matters such as the derivation of the edgebased assembly, the differentiation of the reconstruction operator, the differentiation of the numerical fluxes, and a two-pass assembly for the adjoint problem [10]. He also presents a verification of the adjoint code.
Adjoint example This section will present an example of how adjoint sensitivity analysis can be applied to an actual design problem [45]. The example deals with steady flow in a duct of cross-section h(x), modeled using the quasi-one-dimensional Euler equations on the interval −1 ≤ x ≤ 1. This is actually an example of the continuous adjoint, but its basic principles are the same as for the discrete adjoint as described in the previous section. The objective function chosen for this problem is the integral of the pressure along the duct: ∫1 J=
P dx
(5.67)
−1
For small perturbations of h and U , the perturbation of the cost function becomes: ∫1 δJ = −1
∂P δU dx ∂U
(5.68)
The quasi-one-dimensional Euler equations for this problem can be expressed as: R(U, h) =
dh d (hF ) − P =0 dx dx
(5.69)
with
ρ U = ρu , ρE
ρu F = ρu2 + p , ρuH
0 P = p 0
(5.70)
where ρ is the density, u the velocity, p the pressure, E the total energy and H the stagnation enthalpy of the flow. One more equation is needed to close the system, which is the equation of state for an ideal gas: [ ] 1 P = (γ − 1)ρ E − u2 , 2
H=E+
p γ p 1 2 = + u ρ γ−1ρ 2
(5.71)
Linearizing Eq. (5.69) with respect to small perturbations in the flow variables (δU ) and in the geometry (δh), gives the following expression for δR: d δR = dx
(
∂F h δU ∂U
) −
dh ∂P dδh d δU − P+ (δhF ) = 0 dx ∂U dx dx
(5.72)
96
CHAPTER 5. DESIGN FRAMEWORKS
Because Eq. (5.72) equals zero, it can be multiplied by a Lagrange multiplier Λ, integrated with respect to x, and subtracted from Eq. (5.68) to yield the following equation: ∫1
∂P δU dx ∂U
δJ = −1
[
∫1 −
ΛT −1
d dx
(
∂F δU h ∂U
)
] dh ∂P dδh d − δU − P+ (δhF ) dx dx ∂U dx dx
(5.73)
Applying integration by parts, the details of which will be omitted here, Eq. (5.73) can be rewritten as follows: [
∫1 ΛT
δJ = −1
∫1 −
δU T −1
] d dδh P− (δhF ) dx dx dx
[
] [ ]1 ∂P T ∂F T dΛ dh ∂P T T ∂F Λ− −h − dx − hΛ δU ∂U dx dx ∂U ∂U ∂U −1
(5.74)
Since the Lagrange multiplier can be chosen arbitrarily, it is chosen such that it eliminates the second and third right-hand side terms of Eq. (5.74), hence: −h
∂P T ∂F T dΛ dh ∂P T Λ− =0 − ∂U dx dx ∂U ∂U
(5.75)
and [ hΛT
∂F δU ∂U
]1 =0
(5.76)
−1
Now the perturbation of the objective function can be expressed as follows: [
∫1 δJ =
Λ −1
T
] dδh d P− (δhF ) dx dx dx
(5.77)
This equation gives a relationship between the perturbation of the geometry δh and the perturbation of the objective function δJ. The perturbation of the flow variables δU does not appear in Eq. (5.77) and as a result the problem can be solved independently of the number of design variables.
5.2. HIGH FIDELITY DESIGN FRAMEWORK
97
Validation of CSRT/adjoint coupling Before the adjoint solver can be used in an optimization test case, it has to be verified that it works in combination with the CSRT method. In order to do this, the gradients of the aerodynamic coefficients (cl , cd , L/D) with respect to the design variables were calculated using both (forward) finite differencing (FD) and the adjoint solver. Two validation cases were run: one two-dimensional case using the NACA0012 airfoil at M = 0.80 and α = 1.00o and one three-dimensional case using the ONERA M6 wing at M = 0.84 and α = 3.06o . The gradients with respect to both the Bernstein and B-spline coefficients were computed. For the two-dimensional case, a step size of 0.01 was used for both the FD and adjoint calculations. For the three-dimensional case, a step size of 0.1 was used. The adjoint solver only outputs the gradients of the lift and drag coefficients with respect to the design variables. However, using this information and knowing the step size used by the solver, the gradient of the lift-to-drag ratio can also be computed. The first step is a finite difference approximation of the L/D gradient: CL ∂( C ) D
∂xi
≈
CL CD (x
+ ∆xi ) −
CL CD (x)
∆x
(5.78)
where again ∆xi = ∆x · ˆ ei . In order to find the gradient of the lift-to-drag ratio, the approximations of the lift and drag coefficient gradients are also needed. They are given below: ∂CL CL (x + ∆xi ) − CL (x) ≈ ∂xi ∆x
(5.79)
∂CD CD (x + ∆xi ) − CD (x) ≈ ∂xi ∆x
(5.80)
From Eq. (5.79) and (5.80) the following relations can be derived: CL (x + ∆xi ) ≈
∂CL ∆x + CL (x) ∂xi
(5.81)
CD (x + ∆xi ) ≈
∂CD ∆x + CD (x) ∂xi
(5.82)
Combining Eq. (5.81) and (5.82) gives: CL (x + ∆xi ) CL (x + ∆xi ) = ≈ CD CD (x + ∆xi )
∂CL ∂xi ∆x + CL (x) ∂CD ∂xi ∆x + CD (x)
(5.83)
This value can be computed using the data available from the solver. Substituting this value into Eq. (5.78) gives the gradient of the lift-to-drag ratio with respect to the design variables. For both cases, the adjoint method proved to be faster than using FD. There was, however, a significant difference between the two- and three-dimensional cases. For the airfoil, FD took about 2 hours to compute the 20 gradients (of either the shape or
CHAPTER 5. DESIGN FRAMEWORKS 100
40
50
20
0 1
2
3
4
5
6
7
8
9
10
FD Adjoint
-50
-100
L/D gradient, [~]
L/D gradient, [~]
98
0 1
2
3
4
5
6
7
8
9
10
FD Adjoint
-20
-40
-150
-60
Coefficient number, [~]
Coefficient number, [~]
Figure 5.17: FD and adjoint gradients of L/D w.r.t. Bernstein coefficients (upper side)
Figure 5.18: FD and adjoint gradients of L/D w.r.t. B-spline coefficients (upper side)
refinement function), while the adjoint solver needed approximately 90 minutes. Even though this difference would save time during an optimization process, it is not very significant. Looking at the three-dimensional case, the strength of the adjoint method becomes apparent. It required approximately 10 hours to compute the 200 gradients using the adjoint, versus over 70 hours for FD. This means that an optimization using 10 design iterations could be performed in approximately 4 days using the adjoint code, while it would take about a month using finite differences. All calculations were performed on a PC with an Intel Core 2 Duo T7500 2.20 GHz processor and 3.50 GB of RAM. The computer was running Windows XP Professional. The results of the validation for the two-dimensional case are shown in Fig. 5.17 and 5.18. They show the L/D gradients of both the Bernstein and the B-spline coefficients for the upper side of the NACA0012 airfoil. The graphs show very good agreement between the FD and adjoint gradients (rBernstein = 0.99964, rB−spline = 0.99892). For the three-dimensional case, the gradients are shown in Fig. 5.19 to 5.22. The FD and adjoint gradients are presented separately as contour plots. There again seems to be good agreement between the FD and adjoint results (rBernstein = 0.95350, rB−spline = 1.00000). More adjoint validation results can be found in Appendix B.
Trust region reflective algorithm This section will explain the basic principles of the trust region reflective algorithm, as used by Matlab’s constrained optimization algorithm fmincon. The main idea of the trust region approach is to approximate the objective function f (x) with a quadratic function q(s), such that the behavior of f (x) is reasonably reflected in a region N around a design point x. A minimization is then performed over the trust region: min{q(s), s ∈ N } s
(5.84)
If f (x + s) < f (x) then the current design point x is updated to x + s. If not, then the trust region N will be shrunk and the minimization of Eq. (5.84) is performed
5.2. HIGH FIDELITY DESIGN FRAMEWORK
3
10
9
2
9
2
8
1
8
1
Coefficient number, [~] (in y−direction)
L/D gradient
0
7
−1
6
−2 5
−3
4
−4
−2 5
2
−6
2
−7
1 1
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
10
Figure 5.19: FD gradients of L/D w.r.t. Bernstein coefficients (upper side)
10
−3
4 3
3
−1
6
−5
2
0
7
3
1 1
3 L/D gradient
Coefficient number, [~] (in y−direction)
10
99
−4 −5 −6 2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
Figure 5.20: Adjoint gradients of L/D w.r.t. Bernstein coefficients (upper side)
10
1
1 L/D gradient
L/D gradient 9
9
0.5
0.5 8
7
0
6 −0.5 5 4
−1
Coefficient number, [~] (in y−direction)
Coefficient number, [~] (in y−direction)
8
7
0
6 −0.5 5 4
−1
3
3
−1.5
−1.5 2
2 1 1
−7
10
2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
10
−2
Figure 5.21: FD gradients of L/D w.r.t. B-spline coefficients (upper side)
1 1
2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
10
−2
Figure 5.22: Adjoint gradients of L/D w.r.t. B-spline coefficients (upper side)
100
CHAPTER 5. DESIGN FRAMEWORKS
again. By defining q as the first two terms of the Taylor approximation of f at x, the trust-region subproblem can be described as follows: { } 1 T min s Hs + sT g such that ∥Ds∥ ≤ ∆ (5.85) 2 where H is the Hessian matrix and D is a diagonal scaling matrix. The Matlab algorithm restricts the trust-region subproblem to a linear, two-dimensional subspace, spanned by the vectors s1 and s2 . The first vector s1 is always in the direction of the provided gradient g. The second vector s2 can be defined either in the direction of steepest descent (Eq. (5.86)) or in the direction of negative curvature (Eq. (5.87)). H · s2 = −g
(5.86)
sT2 · H · s2 < 0
(5.87)
The unconstrained trust-region approach can be summed up as follows: i. Formulate the two-dimensional sub-problem ii. Solve Eq. (5.85) to determine s iii. If f (x + s) < f (x), then set x = x + s iv. Adjust trust region dimension ∆ v. Repeat i-iv until convergence If constraints are present, the situation becomes slightly more complicated. However, the basic outline of the method as described above, does not change. For more background on the algorithm and the way constraints are handled, the reader is referred to literature [1]. In order to decide when to stop the optimization, a number of different stopping criteria can be implemented, such as: i. Termination tolerance on the constraint violation ii. Termination tolerance on the function value iii. Termination tolerance on the design variables Which criterion is most suitable strongly depends on the design problem at hand.
Active set algorithm (SQP) The active set algorithm tries to find a solution to a standard minimization problem subject to the Karush-Kuhn-Tucker (KKT) conditions [38]. This can be described mathematically as follows:
5.2. HIGH FIDELITY DESIGN FRAMEWORK
101
min f (x) subject to:
Gi (x) = 0, i = 1, ..., me Gi (x) ≤ 0, i = me + 1, ..., m ∇f (x∗ ) +
∑m i=1
λi · ∇Gi (x∗ ) = 0
(5.88)
λi · Gi (x∗ ) = 0, i = 1, ..., me λi ≥ 0, i = me + 1, ..., m where f (x) is the objective function, G(x) are the constraints, and λ are the Lagrange multipliers. There are several ways of finding a solution to Eq. (5.88). The Matlab algorithm fmincon makes use of sequential quadratic programming (SQP). It computes an approximation of the Hessian of the Lagrangian function at each major iteration, using a quasi-Newton updating method. Subsequently, a quadratic programming (QP) subproblem is performed, which is then used to find a search direction for a line search procedure. The QP subproblem is given by: min 12 dT Hk d + ∇f (xk )T d ∇gi (xk )T d + gi (xk ) = 0, ∇gi (xk )T d + gi (xk ) ≤ 0,
i = 1, ..., me i = me + 1, ...m
(5.89)
Various ways to solve Eq. (5.89) can be found in literature [59]. The last step is to update the design variable as follows: xk+1 = xk + αk dk
(5.90)
In addition to the stopping criteria mentioned in the previous section for trust region reflective algorithms, the following criteria can be used for the active set algorithm: i. Maximum number of iterations of sequential quadratic programming method allowed ii. Maximum amount of time allowed for the algorithm iii. Constraint violation tolerance on the inner SQP iteration Again, which method is most suitable depends on the specific optimization problem. For more details about the active set algorithm, the reader is again referred to literature [1]. There it is explained for example how the Hessian matrix is updated, how the QP problem is solved and how the problem is initialized.
102
CHAPTER 5. DESIGN FRAMEWORKS
CHAPTER
6
Test cases
“I realized that if I had to choose, I would rather have birds than airplanes.” - Charles Lindbergh (1902 - 1974)
6.1. LOW FIDELITY FRAMEWORK: 3D WING
105
This chapter presents a number of different test cases, based on the design frameworks presented in Chapter 5. All test cases make use of the CSRT parameterization method. In Section 6.1 a three-dimensional optimization case is described, based on a panel method code (VSAERO) in combination with various optimization techniques. An airfoil optimization case using the in-house Euler code with finite differencing (see Fig. 5.8) is presented in Section 6.2. The same code was used with adjoint sensitivity analysis (see Fig. 5.9) to optimize a three-dimensional wing in Section 6.3 and finally a blended-wing-body geometry in Section 6.4. This last geometry was optimized for two different angles of attack: 3 and 2 degrees.
6.1
Low fidelity framework: 3D wing
This test case focused on handling volume constraints in subsonic conditions, using B-splines. The goal of this test case was to show that the three-dimensional CSRT method, including volume constraint handling, can be successfully coupled to a flow solver and to investigate how using the CSRT method influences the optimization problem. The objective function was the lift-to-drag ratio (L/D).
6.1.1
Test case definition
For this test case a symmetric wing was used that consisted of an inner and an outer section, both with equal width (Sin = Sout = 1.0). The outer section had a taper ratio of 0.50 and a variable sweep angle. The wing is shown in Fig. 6.1. Applying sweep is usually done to improve the transonic characteristics of a wing. The reason it was decided to vary the sweep for this subsonic test case was to show that the coupling between parameterization and analysis tool can accommodate such changes in geometry. The inner section had no taper and no sweep. The twist and dihedral of both wing sections were also zero. The wing was described using the threedimensional CSRT method in which only the coefficients of the shape function were design variables. The box constraint consisted of a prismatic beam with constant rectangular cross-section and was only applied to the inner wing section. It could represent for example a fuel tank or other internal volume. An obvious procedure would have been to optimize for constant values of the lift coefficient. However, this would have required a considerably higher number of analysis runs, since the aerodynamic analysis tool could only be run for a specific angle of attack. It was therefore decided to use a constant angle of attack of 2 degrees. Analysis showed that the resulting error was small. As was stated earlier, the wing was described using the three-dimensional CSRT method. However, only one airfoil was used for the entire wing and this airfoil can be described using the two-dimensional CSRT method. This makes explaining the test case and the results more straightforward. Using the two-dimensional analogy, it can be explained how the constraints on the coefficients of the shape function were determined. As explained in Chapter 4, the first step was to find the B-spline coefficients of the refinement function that fitted the airfoil tightly around the box constraint. This resulted in the
106
CHAPTER 6. TEST CASES
cr Λ ct Sin
Sout
Figure 6.1: Definition of the 3D wing
airfoil shown in Fig. 6.2. The values of the Bernstein coefficients describing this airfoil form the lower bounds during the optimization. For this particular test case, these lower bounds were also used as the initial solution. In other words, Fig. 6.2 formed the starting point of the optimization. During the optimization, the refinement function was fixed and the shape function was allowed to vary. For both the upper and the lower side of the airfoil a fifth order Bernstein polynomial was used. Kulfan showed that this is sufficient to approximate a typical airfoil [41]. A fifth order Bernstein polynomial requires six coefficients, so the total number of coefficients governing the airfoil shape would have been twelve. However, to keep the nose of the airfoil round, the first coefficients of both Bernstein polynomials were kept equal. This left as design variables 11 Bernstein coefficients and the outer section sweep angle. It was already mentioned that the lower bounds on the Bernstein coefficients were equal to their initial values. Their upper bounds were set to twice these initial values, which translates to a maximum thickness-to-chord ratio of about 16 %. The reason the thickness had to be constrained in this way is that an inviscid optimizer has no inherent way of doing this. It was thus fully expected that at least some of the coefficients would reach these bounds. This was not a problem, however, since the goal of this test case was merely to show a correct coupling and functioning of parameterization, aerodynamic analysis, and optimizer. Finally, the sweep angle was initialized at 0 degrees and was allowed to vary between 0 and 50 degrees. For higher sweep angles the panel mesh would have become too distorted to give accurate results. The definition of the sweep angle Λ can be found in Fig. 6.1. The taper ratio λ was defined as ct /cr . Three optimization schemes were employed. For the first scheme, the active set algorithm contained in the Matlab function fmincon was used. Every function evaluation required by the optimizer was performed using VSAERO. This scheme will be referred to as “SQP (direct)”. For the second scheme, a response surface was created using Latin hypercube sampling (500 sample points) and the Kriging method (see Section 5.1.3). Matlab’s active set algorithm was used to find the optimum value of the response surface. This second scheme will be referred to as “SQP on response surface”. The last optimization scheme used the response surface from the second scheme, but for this case particle swarm optimization was used instead of SQP. This final case will be referred to as “PSO on response surface”.
6.1. LOW FIDELITY FRAMEWORK: 3D WING
107
Relative position from airfoil centerline, ζ [~]
0.2
0.1
0
-0.1
-0.2 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Relative chordwise position, ψ [~]
Figure 6.2: The airfoil as fitted around the box constraint Table 6.1: Results of the various optimization schemes Test case SQP (direct) SQP on response surface PSO on response surface
6.1.2
CL [∼] 0.726 0.591 0.619
CD [∼] 0.0356 0.0299 0.0313
L/D[∼] 20.4 19.7 19.7
∆L/D[%] 21 17 17
Λ[deg] 10 44 44
Test case results
Table 6.1 presents the results for the three optimization schemes introduced in the previous section. Figure 6.3 presents the shape of the airfoils following from the test cases presented in Table 6.1. The airfoil shapes are represented by the values of their Bernstein coefficients. The graph clearly shows that the direct SQP solution got trapped in a local optimum. Some of the variables stayed at or close to 1.0 while the others went to their extreme values. The other two solutions show a thickness peak on the upper aft part of the airfoil. The actual shape of the optimized profiles can be found in Fig. 6.4(a) to 6.4(c). The initial airfoil is also shown in the plots. The best profile following from the test cases is shown in more detail in Fig. 6.5. Its most striking feature is the aft location of the point of maximum thickness. This can easily be explained by considering the test conditions as implemented in VSAERO. The program was run for subsonic speed (M = 0.1) and a Reynolds number of only 2.33 · 106 . For these conditions the most straightforward way of reducing drag was by postponing the transition point of the boundary layer by moving the region of adverse pressure aft. This is exactly what the optimizer did by moving the point of maximum thickness aft. In fact, comparing the optimized airfoil to the laminar NACA67-110
108
CHAPTER 6. TEST CASES
2.50 Upper side
Lower side
M = 0.10
Re = 2.33 * 106
2.00 α = 2.000 1.50 1.00 0.50 0.00 SQP (direct)
SQP on RS
PSO on RS
Figure 6.3: Optimized airfoils expressed in terms of their Bernstein coefficients
0.1
ζ [~]
0.05
Initial airfoil 0
Optimized airfoil
-0.05
-0.1 0
0.2
0.4
0.6
0.8
1
0.1
0.1
0.05
0.05
ζ [~]
ζ [~]
a)
ψ[~]
0
-0.05
-0.05
-0.1 0
b)
0
0.2
0.4
ψ [~]
0.6
0.8
-0.1 0
1
c)
0.2
0.4
ψ[~]
0.6
0.8
1
Figure 6.4: Optimized airfoils compared to the unit airfoil, (a) SQP (direct), (b) SQP on response surface, (c) PSO on response surface
6.2. HIGH FIDELITY FRAMEWORK: NACA0012
109
Table 6.2: Results of the various optimization schemes for constant lift coefficient Test case SQP (direct) SQP on response surface PSO on response surface
α [deg] 1.5 2.5 2.3
CL [∼] 0.655 0.659 0.661
CD [∼] 0.0327 0.0323 0.0328
L/D[∼] 20.0 20.4 20.1
-1
∆L/D[%] 13 15 14
0.2
-0.5
0.1
0
0
0.5
-0.1
1
Relative position from airfoil centerline, ζ [~]
Pressure coefficient, Cp [~]
M = 0.10 6 Re = 2.33 * 10 α = 2.00o
-0.2 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Relative chordwise position, ψ [~]
Figure 6.5: Solution of the SQP test case, including the box constraint
profile shows that they are very similar. This is shown in Fig. 6.6 and 6.7. Both profiles show laminar flow up to about 80 % of the airfoil. It is not clear why two of the results presented in Table 6.1 show high sweep angles. One would expect an unswept wing to give the best performance in subsonic conditions. However, this does not influence any of the conclusions drawn in this section. As was mentioned in Section 6.1.1, the optimization was performed at a constant angle of attack of 2 degrees. In order to see if this approach was justified, the airfoils were also analyzed for a near constant lift coefficient. The results are shown in Table 6.2. The table also shows the angle of attack required to reach the given lift coefficients. The last column presents the increase in lift-to-drag ratio with respect to the initial solution at a similar value of the lift coefficient (CL = 0.654). As was to be expected, the results differ from the ones presented in Table 6.1. However, the conclusions for this test case remain the same. The three optimization methods perform similarly and all create a good improvement over the initial solution.
6.2
High fidelity framework: NACA0012
This test case modeled an airfoil in transonic conditions using the in-house Euler solver. The goal of this test case was to demonstrate that the solver is capable of
110
CHAPTER 6. TEST CASES
Relative position from airfoil centerline, ζ [~]
0.08 Optimized airfoil NACA67-110
0.06 0.04 0.02 0 0
0.2
0.4
0.6
0.8
1
Relative chordwise position, ψ [~]
Figure 6.6: Optimized airfoil in comparison to the NACA67-110
-1 Optimized airfoil
Pressure coefficient, Cp [~]
NACA67-110
-0.5
0
0.5
1 0
0.2
0.4
0.6
0.8
1
Relative chordwise position, ψ [~]
Figure 6.7: Pressure distribution of the optimized airfoil in comparison to that of the NACA67-110
6.2. HIGH FIDELITY FRAMEWORK: NACA0012
111
modeling shock waves and to demonstrate that using a B-spline refinement leads to lower drag on the airfoil without significantly increasing the computation time.
6.2.1
Test case definition
For this test case the NACA0012 subsonic airfoil was used as an initial guess. Using a subsonic airfoil for a transonic test case guaranteed that the starting point of the optimization was far away from the optimum. This way it could be shown that the optimizer was largely independent on the initial solution. The airfoil was first optimized for a Mach number of 0.80 and an angle of attack of 1.25 degrees. The validation in Section 5.2.1 showed that the Euler solver performs very well for this flow condition. The optimization was performed in two steps. First, a “general” optimization was performed using only the Bernstein coefficients of the shape function as design variables. Second, a “regional” refinement was performed using the B-spline coefficients of the refinement function. The reason for this twostep approach is to exploit the advantageous properties of the CSRT method as much as possible. Because the Bernstein coefficients have a global effect on the airfoil shape, they are very suitable for performing a “general” optimization where the whole airfoil is changed from a subsonic profile to a transonic one. In the second step, a “regional” refinement is made possible by the ability of the B-spline to change its shape locally. This has the potential to result in a further weakening of the shock wave and consequently lower drag of the airfoil. For both steps, 20 shape variables were used, 10 for the upper and 10 for the lower side of the airfoil. The bounds on the Bernstein coefficients were set to [0.1, 0.2], which amounts to about ± 35 % compared to the Bernstein coefficients that describe the NACA0012 airfoil. The B-spline coefficients were allowed to vary between 0.8 and 1.2. Since the B-splines were only meant to locally refine the airfoil shape, ideally these bounds should not be reached. After this initial optimization, another batch of optimization runs was performed for Mach numbers ranging from 0.75 to 0.86 and an angle of attack of 1 degree. This was done to investigate the capability of the algorithm to improve the lift-to-drag ratio for different Mach numbers.
6.2.2
Test case results
Figure 6.8 shows the airfoils resulting from the “general” (Bernstein) and “regional” (B-spline) optimizations, as explained in Section 6.2.1. Both optimizations were performed using Matlab’s active set algorithm. The figure also shows the NACA0012 airfoil that was used as an initial solution. From the figure it is clear that the Bernstein optimization transformed the subsonic NACA0012 profile into a typical transonic airfoil. This is very well illustrated by the pressure distributions over both profiles, which can be seen in Fig. 6.9(a) and 6.9(b). Figure 6.9(a) shows a large shock wave on the upper aft part of the airfoil. In Fig. 6.9(b) it can be seen that this shock wave decreased considerably after the first optimization. Figure 6.9(c) shows the pressure distribution after the second optimization, that employed B-splines. Compared to
112
CHAPTER 6. TEST CASES
0.06
M = 0.80 0 α = 1.25
Relative position from airfoil centerline, ζ [~]
0.04
0.02
0
−0.02
−0.04 Bernstein optimization B−spline refinement NACA0012
−0.06 0
0.2
0.4 0.6 Relative chordwise position, ψ [~]
0.8
1
Figure 6.8: Airfoils resulting from the two-dimensional test case
Table 6.3: Lift and (wave) drag results for M = 0.80, α = 1.25o
NACA0012 Bernstein optimization B-spline refinement B-spline only optimization
cl [∼] 0.356 0.403 0.422 0.292
cdw [∼] 0.0234 0.0025 0.0019 0.0018
L/Dw [∼] 15 161 222 162
Fig. 6.9(a) and 6.9(b), the shock wave has largely disappeared. In Section 6.2.1 it was mentioned that the B-spline coefficients should not reach their bounds of 0.8 and 1.2. With a minimum coefficient of 0.935 and a maximum coefficient of 1.058 they stayed well within these bounds. The conclusions drawn above can be verified by looking at the values of the wave drag, which are presented in Table 6.3. After the B-spline refinement, the wave drag decreased by another 25 %, which is a significant improvement. The computation times for the Bernstein optimization and the B-spline refinement were approximately equal. As a last test case, the NACA0012 airfoil was optimized using B-splines only. Figure 6.9(d) shows the pressure distribution resulting from this optimization. It is clear that the graph is much less smooth than in Fig. 6.9(a), 6.9(b), and 6.9(c). This is likely to have a negative effect on the accuracy of the solver. Looking at Table 6.3 it can be seen that even though the resulting wave drag is similar to the B-spline refinement case, the lift is much lower, resulting in a lower value of the lift-to-drag ratio.
6.2. HIGH FIDELITY FRAMEWORK: NACA0012
a)
b)
c)
d)
113
Figure 6.9: Pressure distributions, (a) for the NACA0012, (b) after the Bernstein optimization, (c) after the B-spline refinement, (d) after the B-spline only optimization
114
CHAPTER 6. TEST CASES
0.1
0
−0.1 NACA0012 (CSRT) 1st derivative * 0.1 2nd derivative * 0.1 NACA0012 (original) 1st derivative * 0.1 2nd derivative * 0.1
−0.2
−0.3
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 Relative chordwise position, ψ [~]
0.8
0.9
1
Figure 6.10: Comparison of the geometrical derivatives of the NACA0012 airfoil
Shape derivatives This section discusses the first and second derivatives of the shapes resulting from each optimization step. Figure 6.10 shows the derivatives of the original NACA0012 airfoil, as well as those of the CSRT model that was used to initialize the optimization. To create this CSRT model, only the Bernstein polynomials of the shape function were varied, which explains the smooth second derivative. Figure 6.11 shows the first and second shape derivatives of the airfoil after the first optimization step. The derivatives of the original NACA0012 airfoil are also shown. The largest change can be seen in the second derivative, which shows an oscillatory behavior. It is still smooth, however, because again only the Bernstein polynomials were varied in this step. Finally, Fig. 6.12 shows the result of the B-spline refinement step. Again, the second derivative shows the most interesting features. At the point where the shock wave had to be reduced, around ψ = 0.7, a small peak was introduced by the refinement function. Interestingly, this discontinuity in the second derivative helped to reduce the strength of the shock wave.
Results at different Mach numbers Figures 6.13 to 6.15 show the lift coefficient, the drag coefficient, and the lift-todrag ratio resulting from the optimization runs performed at different Mach numbers and an angle of attack of 1 degree. Each figure shows the results from the “general” (Bernstein) and “regional” (B-spline) optimizations, hence each point on these graphs represents a different airfoil. The NACA0012 results are also included for reference. In Fig. 6.13 it can be seen that the lift coefficient varies significantly with Mach number. The reason for this is that the objective function of the optimization is the
6.2. HIGH FIDELITY FRAMEWORK: NACA0012
115
0.1
0
−0.1 Bernstein optimization (CSRT) 1st derivative * 0.1 2nd derivative * 0.1 NACA0012 (original) 1st derivative * 0.1 2nd derivative * 0.1
−0.2
−0.3
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 Relative chordwise position, ψ [~]
0.8
0.9
1
Figure 6.11: Comparison of the geometrical derivatives of the airfoil after the Bernstein optimization step
0.1
0
−0.1 B−spline refinement (CSRT) 1st derivative * 0.1 2nd derivative * 0.1 NACA0012 (original) 1st derivative * 0.1 2nd derivative * 0.1
−0.2
−0.3
0
0.1
0.2
0.3 0.4 0.5 0.6 0.7 Relative chordwise position, ψ [~]
0.8
0.9
1
Figure 6.12: Comparison of the geometrical derivatives of the airfoil after the B-spline refinement step
116
CHAPTER 6. TEST CASES Table 6.4: Lift-to-drag ratio results for α = 1.00o
M [∼] 0.75 0.80 0.83 0.84 0.85 0.86
L/Dw [∼]
L/Dw [∼]
initial 66.2 15.6 9.2 7.9 6.6 5.2
“general” 274.4 205.7 144.7 76.6 34.3 19.9
∆[%]
L/Dw [∼]
∆[%]
315 1220 1475 870 418 282
“regional” 275.6 214.3 189.5 140.3 93.3 30.9
0.4 4 31 83 172 55
∑
∆[%]
316 1275 1963 1677 1311 493
0.5
Lift coefficient, c l [~]
0.4
0.3
0.2
0.1 Bernstein optimization B-spline refinement NACA 0012
0.0 0.74
0.76
0.78
0.80
0.82
0.84
0.86
0.88
Mach number, M [~]
Figure 6.13: Lift coefficient versus Mach number
lift-to-drag ratio. Hence for each Mach number the optimizer will find the value of the lift coefficient for which this lift-to-drag ratio is maximum. Figure 6.14 clearly shows the effect of both the “general” and “regional” optimizations. Not only do they considerably decrease the magnitude of the drag, they also postpone the drag divergence Mach number from about 0.70 for the NACA0012 to 0.83 after the Bernstein optimization and to 0.85 after the B-spline refinement. It should be emphasized that for each Mach number a different airfoil was analyzed. Finally, Fig. 6.15 confirms the increase in lift-to-drag ratio after each optimization step. These results have been summarized numerically in Table 6.4. Figures 6.16 and 6.17 show the performance of each optimization step expressed in terms of the relative increase in lift-to-drag ratio compared to the previous solution. Figure 6.18 depicts the combined result of both optimizations. As is clear from these figures, there is a certain Mach range for which the optimizations seem to be most effective. The reason for this lies in the nature of the flow solver and the parameteriza-
6.2. HIGH FIDELITY FRAMEWORK: NACA0012
117
0.08 Bernstein optimization B-spline refinement NACA 0012
Drag coefficient, cd [~]
0.06
0.04
0.02
0.00 0.74
0.76
0.78
0.80
0.82
0.84
0.86
0.88
Mach number, M [~]
Figure 6.14: Drag coefficient versus Mach number
300
Bernstein optimization B-spline refinement NACA 0012
Lift-to-drag ratio, L/D [~]
250
200
150
100
50
0 0.74
0.76
0.78
0.80
0.82
0.84
0.86
0.88
Mach number, M [~]
Figure 6.15: Lift-to-drag ratio versus Mach number
118
CHAPTER 6. TEST CASES
'Global' L/D increase [%]
2000
1500
1000
500
0 0.74
0.76
0.78
0.80
0.82
0.84
0.86
0.88
Mach number, M [~]
Figure 6.16: Relative improvement in lift-to-drag ratio versus Mach number after “general” optimization
tion. Since the only way for the algorithm to reduce the lift-to-drag ratio is to decrease the size and strength of the shock wave, no real gain is expected for Mach numbers lower than about 0.75. As the Mach number increases (0.75 < M < 0.85), both the Bernstein optimization and the B-spline refinement show significant improvements. However, as the Mach number grows even further, a bow shock starts to form on the airfoil and the optimizer is not capable anymore of considerably increasing the performance of the airfoil. This is not a shortcoming of the parameterization method, however. Instead, it demonstrates how the required optimization scheme depends on the flow regime for which the airfoil or wing is to be optimized.
6.3
High fidelity framework: ONERA M6
The goal of this test case was to see whether the optimization strategy used in Section 6.2 would give similar results for a three-dimensional wing.
6.3.1
Test case definition
The ONERA M6 wing, which was used to validate the implementation of the CSRT method into the Euler code in Section 5.2.1, was chosen as the initial solution. The optimization strategy was the same as explained in Section 6.2.1: first a “general” optimization was performed using the shape function, followed by a “regional” refinement using the refinement function. For both steps the angle of attack was kept constant. However, if the value of the lift coefficient resulting from the “general” op-
6.3. HIGH FIDELITY FRAMEWORK: ONERA M6
119
'Local' L/D increase [%]
200
150
100
50
0 0.74
0.76
0.78
0.80
0.82
0.84
0.86
0.88
Mach number, M [~]
Figure 6.17: Relative improvement in lift-to-drag ratio versus Mach number after “regional” refinement
2500
Total L/D increase [%]
2000
1500
1000
500
0 0.74
0.76
0.78
0.80
0.82
0.84
0.86
0.88
Mach number, M [~]
Figure 6.18: Total relative improvement in lift-to-drag ratio versus Mach number
120
CHAPTER 6. TEST CASES
timization differed significantly from the initial value, the angle of attack was changed to give approximately the same value of CL . This new value of the angle of attack was then used for the refinement step. For the two-dimensional test case of Section 6.2 an active set algorithm was employed, using finite differencing to find the gradients of the shape variables. However, for this three-dimensional case, the sharp increase in number of shape variables and time required to solve the flow field renders this approach unfeasible. Instead, adjoint sensitivity analysis was used in combination with a trust region reflective algorithm (see Section 5.2.2).
6.3.2
Test case results
This section will present the results of the ONERA M6 test case. The pressure distributions for each optimization step will be shown, as well as the corresponding wing shapes. Additionally, the shape and refinement functions required to model these wing shapes will be shown. The initial solution will be presented first, followed by the results of the “general” optimization step. Then the results of the “regional” refinement step are discussed and finally the results of an optimization using only the B-spline coefficients will be presented.
Initial solution: the ONERA M6 wing Figures 6.19 and 6.20 show the pressure distribution on the ONERA M6 wing for an angle of attack of 3.06 degrees and a Mach number of 0.84. This is the same pressure distribution that was used for validating the Euler code in Section 5.2.1. It shows supersonic regions and two strong shock waves: one along the leading edge and one parallel to the trailing edge. The shock waves seem to merge near the wing tip. For this condition, the lift and drag coefficients are equal to 0.276 and 0.0237 respectively, resulting in a lift-to-drag ratio of 11.7. The reason why the coefficients differ slightly from those presented in Section 5.2.1 is that for this test case a courser grid was used (20568 nodes, 100069 cells) to speed up the optimization process. The pressure distribution shown in the figures formed the basis for the optimization. Figure 6.21 shows the shape function that was used to generate the mesh of the ONERA M6 wing. The latter is shown in Fig. 6.22. Considering that the undeformed mesh had the shape of the class function (see Section 5.2.1), the shape function represents the ratio between the ONERA M6 wing and the class function.
“General” optimization step For this step the thickness of the ONERA M6 wing was taken as a minimum thickness constraint. In other words, the Bernstein coefficients of the shape function were only allowed to increase. Figures 6.23 and 6.24 show the pressure distribution after the first optimization step. It is clear from both figures that the shock wave along the leading edge has been greatly reduced. Additionally, the pressure field on the rest of the wing has flattened considerably, resulting in a more even distribution of the lift.
6.3. HIGH FIDELITY FRAMEWORK: ONERA M6
121
Pressure coefficient [~] 0 0.2
Chordwise position [m]
−0.2 0.4 −0.4 0.6 −0.6 0.8 −0.8 1 −1 1.2 0
0.2
0.4
Figure 6.19: contours
0.6 0.8 Spanwise position [m]
1
1.2
1.4
ONERA M6 pressure
Figure 6.20: ONERA M6 pressure distribution
0.2
0.18
Shape function
0.25
0.16
0.2 0.14 0.15
1
0.12
0.8 0.1 0
0.6 0.2 0.6 Chordwise direction
0.1
0.4
0.4 0.2 0.8 1
Spanwise direction
0
Figure 6.21: ONERA M6 shape function
Figure 6.22: ONERA M6 mesh
122
CHAPTER 6. TEST CASES
Pressure coefficient [~] 0 0.2
Chordwise position [m]
−0.2 0.4 −0.4 0.6 −0.6 0.8 −0.8 1 −1 1.2 0
0.2
0.4
0.6 0.8 Spanwise position [m]
1
1.2
1.4
Figure 6.23: “General” optimization pressure contours
Figure 6.24: “General” optimization pressure distribution
0.2
0.18
Shape function
0.25
0.16
0.2 0.14 0.15
1
0.12
0.8 0.1 0
0.6 0.2 0.6 Chordwise direction
0.1
0.4
0.4 0.2 0.8 1
Spanwise direction
0
Figure 6.25: “General” optimization shape function
Figure 6.26: “General” optimization mesh
This results in a lift coefficient of 0.381 and a drag coefficient of 0.0276 (L/D = 13.8). In order to make a fair comparison with the initial solution, the lift coefficient should be approximately equal. This was achieved by lowering the angle of attack from 3.06 to 2.00 degrees, leading to a lift coefficient of 0.279 and a drag coefficient of 0.0200. The resulting lift-to-drag ratio of 13.9 means an increase of 19 % compared to the ONERA M6 wing at the same lift coefficient. The shape function belonging to the pressure distributions of Fig. 6.23 and 6.24 is shown in Fig. 6.25. Comparing this shape function to the one in Fig. 6.21, it can be seen that the “leading edge” and “trailing edge” of the shape function have not changed position. This makes sense, since the value of the class function is zero there and as a result deforming the shape function at that location would not have any effect. It is also clear from the figure that the rear of the wing has generally become thicker, especially near the wing tip. It is this thickening that led to the flattening of the pressure field apparent in Fig. 6.24. This also explains why the peak of the shape function lies near the wing tip, since this is where the shock wave was strongest. The resulting wing mesh is shown in Fig. 6.26.
6.3. HIGH FIDELITY FRAMEWORK: ONERA M6
123
2.5 Control point Control polygon B−spline
2 1.5 1 0.5 0 0
0.2
0.4
0.6
0.8
1
Figure 6.27: B-spline example
Pressure coefficient [~] 0 0.2
Chordwise position [m]
−0.2 0.4 −0.4 0.6 −0.6 0.8 −0.8 1 −1 1.2 0
0.2
0.4
0.6 0.8 Spanwise position [m]
1
1.2
1.4
Figure 6.28: “Regional” refinement pressure contours, with new constraints
Figure 6.29: “Regional” refinement pressure distribution, with new constraints
“Regional” refinement step In this step the requirement that the B-spline coefficients were only allowed to increase was found to negatively affect the performance of the method as can be seen in Table 6.5. Only allowing the B-spline coefficients to increase did not leave enough freedom for the refinement function to implement small local changes. Furthermore, Fig. 6.27 illustrates that it is possible for some of the B-spline coefficients to be less than 1, without the B-spline itself becoming smaller than 1 at any point. For this reason, a new constraint strategy was conceived, which gave the refinement function sufficient design freedom while still maintaining the overall thickness of the wing. This time, the bounds on the B-spline coefficients were set to [0.8, 2.0], but with the added constraint that the mean value of all the coefficients had to be greater than or equal to 1. This resulted in the pressure distributions of Fig. 6.28 and 6.29. A further reduction of the shock waves and flattening of the pressure field can be distinguished, leading to a lift coefficient of 0.281 and a drag coefficient of 0.0197. The resulting lift-to-drag ratio of 14.3 means an increase of 3 % from the first optimization step and 22 % compared to the initial solution. Figure 6.30 shows the refinement function resulting from the second optimization
124
CHAPTER 6. TEST CASES
1.15 1.1
Ref nement function
1.3 1.2
1.05
1.1
1
1
1 0.95
0.8
0.9 0.6
0 0.2
0.9
0.4
0.4 0.6
0.2 0.8 1
Spanwise direction
0
Chordwise direction
Figure 6.30: “Regional” refinement function, with new constraints
Figure 6.31: “Regional” refinement mesh, with new constraints
Pressure coefficient [~] 0 0.2
Chordwise position [m]
−0.2 0.4 −0.4 0.6 −0.6 0.8 −0.8 1 −1 1.2 0
0.2
0.4
0.6 0.8 Spanwise position [m]
1
1.2
1.4
Figure 6.32: B-spline only pressure contours
Figure 6.33: B-spline only pressure distribution
step and Fig. 6.31 shows the corresponding mesh. It is clear from these figures that the change in shape is very subtle. In fact, all B-spline coefficients lay between 0.9875 and 1.0625, indicating that it might not even be necessary to implement direct bounds on these coefficients.
B-spline only optimization As a final check, it was investigated what would happen if the B-spline optimization was performed directly on the ONERA M6 wing. The new set of constraints that was used for the refinement step was also used for this case. The result of this optimization is shown in Fig. 6.32 and 6.33. The corresponding lift and drag coefficients (at α = 2.75o ) are 0.283 and 0.0224 respectively. With a lift-to-drag ratio of 12.6, only an 8 % increase in L/D has been achieved, compared to a 22 % increase for the two-step approach. More importantly though, the pressure distribution on this wing is very bumpy and would probably cause the boundary layer to separate early, rendering this solution unusable. Note that this result is very similar to that of the two-dimensional case presented in Section 6.2.
6.4. HIGH FIDELITY FRAMEWORK: BWB
125
1.15 1.1
Ref nement function
1.3 1.2
1.05
1.1
1
1
1 0.8
0.9 0.6
0 0.2
0.95 0.9
0.4
0.4 0.6
0.2 0.8 1
0
Chordwise direction
Spanwise direction
Figure 6.34: B-spline only refinement function
Figure 6.35: B-spline only mesh
Table 6.5: ONERA M6 test case results for similar values of CL (for M = 0.84) Test case ONERA M6 “General” “Regional”, old constraints “Regional”, new constraints B-spline only
α[deg] 3.06 2.00 1.95 2.00 2.75
CL [∼] 0.276 0.279 0.281 0.281 0.283
CD [∼] 0.0237 0.0200 0.0202 0.0197 0.0224
L/D[∼] 11.7 13.9 13.9 14.3 12.6
∆L/D[%] 19 19 22 7
The bumpy nature of the solution is also very apparent from the shape of the refinement function and the wing mesh, as illustrated in Fig. 6.34 and 6.35 respectively. The B-spline coefficients varied significantly more than during the refinement step, their extreme values being 0.8656 and 1.3522.
Summary of results The numerical results of the ONERA M6 test case are summarized in Table 6.5. Looking at Table 6.5, the benefits of the two-step CSRT approach become apparent. Even though the “general” optimization step using Bernstein polynomials produces the biggest improvement in L/D, the “regional” refinement adds an additional improvement of 3 %, which can be very significant. More importantly, this extra 3 % can be reached in a relatively inexpensive way, instead of having to increase the order of all the Bernstein surfaces, which would significantly lengthen the optimization time.
6.4
High fidelity framework: BWB
In addition to the ONERA M6 case presented in Section 6.3, another test case was performed using a geometry closer to that of an actual aircraft. The geometry that was chosen is referred to in literature as the MOB reference aircraft. The MOB
126
CHAPTER 6. TEST CASES
section section section section section section section section section
1 2 3 4 5 6 7 8 9
c[m] 48.0 45.9 41.8 35.7 31.6 27.5 21.4 13.5 9.3
λ[∼] 0.957 0.910 0.853 0.885 0.870 0.780 0.630 0.690 0.420
Λ[deg] 64 64 64 64 64 64 38 38 38
θ0 [deg] 0.0 0.0 0.0 0.2 0.4 0.6 1.1 1.6 2.6
θ[deg] 0.0 0.0 0.2 0.2 0.2 0.5 0.5 1.0 -3.5
Γ[deg] 4.4 4.4 4.4 4.4 4.4 4.4 4.4 5.9 7.4
s(m) 1.0 2.0 3.0 2.0 2.0 3.0 4.5 6.0 14.5
Table 6.6: MOB section data
geometry was developed as part of a European multi-disciplinary optimization project [48] and served as a reference geometry for blended-wing-body (BWB) aircraft in several papers [64,88]. Blended-wing-bodies are widely seen as a very promising future aircraft concept and as such have been the center of many optimization studies [5,61]. The dihedral, taper, twist, and airfoil sections of the BWB are defined for each wing section. The complete definition of the geometry is given in Table 6.6, with each parameter defined as in Section 3.4.
6.4.1
Test case definition
Two optimization cases were carried out. The first case was performed at a Mach number of 0.84 and an angle of attack of 3 degrees. For the second test case, the angle of attack was changed to 2 degrees, while the Mach number remained the same. The baseline geometry used for both test cases is an approximation of the MOB reference aircraft, as presented by Morris [48]. However, the mesh that was read by the mesh reader of the algorithm (see Fig. 5.9), was that of a simple wing with a rectangular planform and a constant cross-section (the class function), as shown in Fig. 6.36. Subsequently, the mesh was deformed in both x- and z-direction using the deformation module to achieve the baseline geometry. The mesh of this geometry is presented in Fig. 6.37. Any shape modification performed during the optimization process is implemented first on the simple wing mesh, after which it is deformed into the blendedwing-body (BWB) geometry. For both optimizations, the same two-step approach was used as in the previous sections. For the first test case, one more optimization was done using B-splines only. After obtaining the results from the different optimization steps, each case was recalculated for similar values of the lift coefficient by varying the angle of attack. In a real design problem, the right approach would have been to keep the lift coefficient constant during the optimization. However, this would require the angle of attack to be an optimization variable, significantly increasing computation time.
6.4. HIGH FIDELITY FRAMEWORK: BWB
Figure 6.36: Shape of the basic mesh as read by the algorithm
6.4.2
127
Figure 6.37: Shape of the baseline geometry
Results of the first test case
This section will present the results of the first test case. The pressure distribution for each optimization step will be shown, as well as the corresponding shape and refinement functions. The initial solution will be presented first, followed by the results of the “general” optimization step. Then, the results of the “regional” refinement step are discussed and finally, the results of an optimization using only the B-spline coefficients will be presented.
Initial solution: the BWB reference geometry Figure 6.38 shows the pressure distribution on the baseline geometry. It shows a strong shock wave on the outboard section of the wing at around two-thirds of the chord. Furthermore, it can be seen that the majority of the wing loading is located on the outboard section. This is common for blended-wing-body designs, since the combination of large relative thickness and high lift of the center section can easily lead to flow separation. Figure 6.39 shows the shape function that was required to create the MOB geometry. As can be seen in Table 6.7, for a Mach number of 0.84 and an angle of attack of 3.00 degrees, the baseline geometry had a lift coefficient of 0.448 and a drag coefficient of 0.0359, leading to a lift-to-drag ratio of 12.5.
“General” optimization step For the first step, the thickness of the baseline BWB geometry was taken as a minimum thickness constraint. In other words, the Bernstein coefficients of the shape function were only allowed to increase. The solution of this “general” optimization step is shown in Fig. 6.40. Compared to Fig. 6.38, the shock wave has decreased in strength and has moved towards the trailing edge, creating a more even lift distribution. For this design, the values of CL and CD were 0.569 and 0.0396 respectively, leading to a value of L/D of 14.5, which is a 16 % improvement over the baseline design.
128
CHAPTER 6. TEST CASES
0.4 0.3
1
0.2 0.25 0.5
−0.2 0.6
−0.4 −0.6
0.4 −0.8 0.2
0.2
0.4 0.3
0.15
0.2 1
0.1 0.6 0.8
0.2
0.3 0.4 0.5 Spanwise position
0.6
0.7
0.4
0.2
0.8
0
Spanwise direction
0
Chordwise direction
Figure 6.38: MOB pressure contours
0.05
0.4
0.6
−1.2
0.2
0.1
0.1
0.8 0 1
−1 Pressure coefficient [~]
0 0
Shape function
Chordwise position
0 0.8
Figure 6.39: MOB shape function
0.4 0.3
1
0.2 0.25 0.5
−0.2 0.6
−0.4 −0.6
0.4 −0.8 0.2
−1 Pressure coefficient [~]
0 0
−1.2
Shape function
Chordwise position
0 0.8
0.2
0.4 0.3
0.15
0.2 1
0.1 0.6 0.8
0.2
0.3 0.4 0.5 Spanwise position
0.6
0.7
0.4
0.8
Figure 6.40: “General” optimization pressure contours
0.05
0.4
0.6 0.2 0.2
0.1
0.1
0.8 0 1
0 Chordwise direction
0
Spanwise direction
Figure 6.41: “General” optimization shape function
The shape function belonging to the pressure distribution of Fig. 6.40 is shown in Fig. 6.41. It shows a large deformation near the trailing edge of the outer section of the wing. This makes sense, since the shock wave was strongest at this position.
“Regional” refinement step For the second step, a slightly more elaborate set of geometrical constraints was used to give the refinement function sufficient design freedom, while still maintaining the overall thickness of the wing. In Section 6.3.2 this was accomplished by constraining the B-spline coefficients of the refinement function such that their mean had to be greater than or equal to 1. However, because of the strong taper of the MOB geometry and the way it is constructed from a rectangular mesh, a deformation of the refinement function near the wing tip leads to a smaller deformation of the actual geometry than a deformation near the wing root. This could lead to the optimizer choosing to make the inner section of the wing thinner, while allowing a smaller thickening of the outer section. In order to prevent this, the constraint was split into two separate constraints: one for the inner section and one for the outer section. For both parts of the wing
6.4. HIGH FIDELITY FRAMEWORK: BWB
129
0.4 1
0.2 1.15
0.6
−0.4 −0.6
0.4 −0.8 0.2
−1 Pressure coefficient [~]
0 0
−1.2
1.1
1.3
−0.2 Ref nement function
Chordwise position
0 0.8
1.2
1.05
1.1 1
1 1
0.9 0.8 0.8 1
0.6 0.8
0.2
0.3 0.4 0.5 Spanwise position
0.6
0.7
0.4
0.8
Figure 6.42: “Regional” refinement pressure contours
0.2 0
Chordwise direction
Figure 6.43: function
0.9
0.4
0.6 0.2
0.1
0.95
0
Spanwise direction
“Regional” refinement
the mean of the B-spline coefficients had to be greater than or equal to 1. Figure 6.42 shows the pressure distribution resulting from the “regional” refinement step. Looking at Fig. 6.43, the refinement function seems rather erratic. However, comparing this solution to the one in Fig. 6.40, a further reduction of the shock wave can be distinguished. Furthermore, the lift has been more equally distributed in spanwise direction, resulting in a lower wing loading of the outer section. The lift and drag coefficients of this refined design were 0.596 and 0.0396 respectively, resulting in a lift-to-drag ratio of 15.1. This means a further 4 % improvement compared to the “general” optimization and a 21 % improvement compared to the baseline geometry.
B-spline only optimization As a final test case, the baseline geometry was optimized using only the B-spline coefficients of the refinement function. The same set of constraints was used as for the refinement step. It was shown in Sections 6.2.2 and 6.3.2 that using only the B-spline coefficients is likely to lead to bumpy and unrealistic solutions. Figures 6.44 and 6.45 show that this was also the case here. Small isolated regions of low pressure can be distinguished on the wing, as well as a double and even triple shock wave on the outer section. So even though there was a 10 % improvement in lift-to-drag ratio, the solution seems to be unfeasible.
Summary of results The results of the first test case have been summarized in Table 6.7. As was mentioned before, all geometries were also analyzed for similar values of the lift coefficient, to see whether the same conclusions hold. Table 6.8 shows the results for each step at CL ≈ 0.53 − 0.54. Looking at the values in the table, it can be seen that not only do the conclusions hold, but the results are actually better than for the cases at constant angle of attack.
130
CHAPTER 6. TEST CASES
0.4 1
0.2 1.15
−0.2 0.6
−0.4 −0.6
0.4 −0.8 0.2
−1
1.1
1.3 Ref nement function
Chordwise position
0 0.8 1.2
1.05
1.1 1
1 1
0.9 0.8 0.8 1
0.6 0.8
Pressure coefficient [~] 0 0
0.1
0.2
0.3 0.4 0.5 Spanwise position
0.6
0.7
0.4
0.2 0.2
0.8
0
0
Chordwise direction
Figure 6.44: B-spline only pressure contours
α[deg] 3.00 3.00 3.00 3.00
CL [∼] 0.448 0.569 0.596 0.505
CD [∼] 0.0359 0.0393 0.0396 0.0368
L/D[∼] 12.5 14.5 15.1 13.7
∆L/D[%] 16 21 10
Table 6.8: MOB test case results for similar values of CL Test case MOB “General” “Regional” B-spline only
α[deg] 3.60 2.80 2.60 3.00
CL [∼] 0.539 0.536 0.531 0.536
Spanwise direction
Figure 6.45: B-spline only refinement function
Table 6.7: MOB test case results Test case MOB “General” “Regional” B-spline only
0.9
0.4
0.6
−1.2
0.95
CD [∼] 0.0456 0.0367 0.0351 0.0393
L/D[∼] 11.8 14.6 15.1 13.6
∆L/D[%] 24 28 15
6.4. HIGH FIDELITY FRAMEWORK: BWB
131
0.4 1
0.2
0.3
0.25
0.8 −0.2 0.6
−0.4 −0.6
0.4 −0.8
0.5 Shape function
Chordwise position
0
0.2
0.4 0.3
0.15
0.2 1 0.1
0.1
0.8
0.2
−1
0 1
0.6 0.8
Pressure coefficient [~] 0 0
0.1
0.2
Figure 6.46: tours
6.4.3
0.3 0.4 0.5 Spanwise position
0.6
0.7
−1.2
0.4
0.2 0.2
0.8
MOB pressure con-
0.05
0.4
0.6 0 Chordwise direction
0
Spanwise direction
Figure 6.47: MOB shape function
Results of the second test case
This section will present the results of the second test case. Again, the pressure distribution and corresponding shape and refinement functions will be shown for each optimization step. The constraints applied to the shape variables were the same as for the first test case. The initial solution will be presented first, followed by the results of the “general” and “regional” optimization steps.
Initial solution: the BWB reference geometry Figure 6.46 shows the pressure distribution on the baseline geometry. Since the only difference with the first test case is a change in angle of attack of 1 degree, this pressure distribution is very similar to that of Fig. 6.38. There is, however, an expected change in the aerodynamic coefficients. For a Mach number of 0.84 and an angle of attack of 2 degrees, the baseline geometry had a lift coefficient of 0.305 and a drag coefficient of 0.0251, resulting in a lift-to-drag ratio of 12.2. Note that this value of L/D is slightly lower than for α = 3o . Naturally, the shape functions of Fig. 6.39 and 6.47, required to generate the initial geometry, are the same.
“General” optimization step The results of the “general” optimization step are presented in Fig. 6.48 and 6.49. Little difference can be distinguished between these graphs and the initial results of Fig. 6.46 and 6.47. The only concernible difference being a slight aftward movement of the shock wave. The values of CL and CD for this solution are 0.368 and 0.0284 respectively, corresponding to a lift-to-drag ratio of 13.0, which is a 7 % improvement over the initial design. This number is much lower than the 16 % improvement that was achieved for the “general” optimization step of the first test case.
132
CHAPTER 6. TEST CASES
0.4 1
0.2
0.3
0.25
0.8 −0.2 0.6
−0.4 −0.6
0.4 −0.8
0.5 Shape function
Chordwise position
0
0.2
0.4 0.3
0.15
0.2 1 0.1
0.1
0.8
0.2
−1
0 1
0.6 0.8
Pressure coefficient [~] 0 0
0.1
0.2
0.3 0.4 0.5 Spanwise position
0.6
0.7
−1.2
0.4
0.2 0.2
0.8
Figure 6.48: “General” optimization pressure contours
0.05
0.4
0.6 0 Chordwise direction
0
Spanwise direction
Figure 6.49: “General” optimization shape function
“Regional” refinement step Figures 6.50 and 6.51 show the results of the “regional” refinement step for the second test case. Looking at Fig. 6.50, a more even distribution of the lift can be distinguished, resulting in a weaker shock wave on the outer section of the wing. The lift and drag coefficients of this design are 0.466 and 0.0299 respectively, leading to a lift-to-drag ratio of 15.6, which is a 20 % improvement compared to the “general” optimization step and a 28 % improvement with respect to the initial geometry. The refinement function of Fig. 6.51 clearly has a much larger amplitude than in the first test case (Fig. 6.43). This seems to be the result of the “general” optimization step not being able to achieve a significant improvement of the objective function. Instead of merely providing a final shape refinement, the “regional” step actually resulted in a larger improvement in L/D than the “general” step. Interestingly, the largest deflection of the refinement function of the second test case occurs in the same region as the largest deflection of the shape function of the first test case (Fig. 6.41). It seems that in the second test case the refinement function compensated for the inability of the shape function to significantly improve the lift-to-drag ratio.
Summary of results Table 6.9 summarizes the results of the second test case for the constant angle of attack (2 degrees) at which the optimizations were performed. Table 6.10 lists the results of the same designs, but for similar values of the lift coefficient (CL ≈ 0.36 − 0.37). Although the relative improvements in lift-to-drag ratio are slightly lower, the conclusions for this second test case remain the same.
6.4. HIGH FIDELITY FRAMEWORK: BWB
133
0.4 1
0.2 1.3
−0.2 0.6
Ref nement function
Chordwise position
0 0.8
−0.4 −0.6
0.4 −0.8 0.2
−1
1.6
1.2
1.4
1.1
1.2 1 1
1 0.8
0.8 0.6
1 0.8
Pressure coefficient [~] 0 0
0.1
0.2
0.3 0.4 0.5 Spanwise position
0.6
0.7
0.4
0.2 0.2
0.8
0 Chordwise direction
Figure 6.50: “Regional” refinement pressure contours
0
Spanwise direction
Figure 6.51: “Regional” refinement function
Table 6.9: Results of second test case Test case MOB “General” “Regional”
α[deg] 2.00 2.00 2.00
CL [∼] 0.305 0.368 0.466
CD [∼] 0.0251 0.0284 0.0299
L/D[∼] 12.2 13.0 15.6
∆L/D[%] 7 28
Table 6.10: Results of second test case for similar values of CL Test case MOB “General” “Regional”
α[deg] 2.40 2.00 1.35
0.8
0.4
0.6
−1.2
0.9
CL [∼] 0.361 0.368 0.368
CD [∼] 0.0286 0.0284 0.0245
L/D[∼] 12.6 13.0 15.0
∆L/D[%] 3 19
134
CHAPTER 6. TEST CASES
CHAPTER
7
Conclusions
“A whole is that which has beginning, middle and end.” - Aristotle (384 BC - 322 BC)
137 This thesis presented a new parameterization technique based on B-splines: the ClassShape-Refinement-Transformation (CSRT) method. A parametric study was presented into the behavior of this CSRT method. It was investigated how well a CSRT curve can be fitted to a set of airfoil data. The results showed a very non-linear relationship between the accuracy of the fit and the number of coefficients used. It was shown that this is a direct result of the behavior of the Bernstein polynomials and B-spline curves, as used in the CSRT method. The Class-Shape-Refinement-Transformation method was successfully coupled to the panel method program VSAERO and an in-house Euler solver. It was shown that including B-splines in the parameterization of the aircraft shape provides the ability to handle volume constraints by applying bounds directly on the design variables. This was not possible with previous parameterization methods based on polynomial and spline representations. The volume constraint handling method as presented in this thesis can be applied equally well in two and three dimensions. Using VSAERO it was shown that the CSRT method can be used to achieve large portions (up to 80 %) of laminar flow on a subsonic wing. The test case also showed that local minima in the design space could pose a problem for the optimization algorithm, possibly requiring the use of non-gradient based methods. In a second test case it was demonstrated that B-splines can also be used as an extra refinement in the shape optimization process. Using an in-house Euler solver, it was shown that this refinement can lead to significant L/D improvements in the order of 25 % to 100 % on the NACA0012 airfoil, by eliminating the shock wave and consequently improving the drag divergence Mach number significantly. A disadvantage of adding the refinement function is that more function variables are required. However, it was shown that the benefits outweigh this disadvantage. A third test case demonstrated a three-dimensional application of the CSRT technique. The lift-to-drag ratio of the ONERA M6 wing was optimized using the same in-house Euler solver as for the two-dimensional case. The sensitivity analysis was performed using an adjoint code integrated in the flow solver. The implementation of the CSRT method in both the flow solver and the adjoint code was validated. A trust region reflective method was chosen as the optimization algorithm. The results of the test case resembled the results of the second test case performed on the NACA0012 airfoil. Using the Bernstein shape function, a significant improvement in lift-to-drag ratio of 19 % was achieved. An extra optimization step using the B-spline refinement function led to an additional improvement of 3 %. With this two-step approach, the lift-to-drag ratio was improved from 11.7 to 14.3. It was also shown that an optimization using only the B-spline coefficients as shape variables can give convergence problems and lead to an unrealistic solution. In a fourth and fifth test case, it was shown that the method previously tested on an airfoil and a simple three-dimensional wing could also be successfully applied to a blended-wing-body design. The two test cases were performed at different angles of attack (two and three degrees) and the same Mach number. For the case at α = 3o , the first optimization step lead to an improvement in lift-to-drag ratio of 24 % (at constant CL ). The refinement step improved L/D by another 3 %, leading to a total increase in lift-to-drag ratio of 28 %. It was again shown that using only the B-spline coefficients
138
CHAPTER 7. CONCLUSIONS
of the refinement function is likely to lead to bumpy and unrealistic solutions. The case for α = 2o showed an interesting deviation from the trend found in the previous test cases. The “general” step resulted in only a marginal improvement in L/D of 3 %, while the “regional” step lead to a significant improvement of 19 % (at constant lift coefficient). Apparently, in cases in which the Bernstein surface of the shape function cannot cope with the design space, the B-spline surface of the refinement function can compensate for this. However, if a large part of the optimization is performed using B-splines, this increases the chance of getting an unrealistic solution such as found in previous test cases. The work presented in this thesis proves that the CSRT method is a very intuitive and effective way of parameterizating aircraft shapes, both in two as well as in three dimensions. The method allows for a two-step approach which has the potential to significantly increase the lift-to-drag ratio of various aircraft shapes. It was also shown that using an adjoint algorithm provides the computational efficiency necessary to perform true three-dimensional shape optimization. A recommendation for future work is to include the lift coefficient in the shape optimization, thereby eliminating the need to optimize for constant angle of attack. This would allow a fairer comparison of the optimization results, albeit at the expense of extra computation time. Another step that needs to be taken is to expand the tool to be able to handle not just airfoils and wings, but also more complex shapes, such as wing-fuselage combinations. Once this is accomplished, a more extensive and more realistic range of test cases can be performed. Additionally, a set of boundary layer equations could be incorporated into the Euler code to be able to capture viscous effects. It would also be interesting to see how applying a constraint on the moment coefficient of the blended-wing-body would influence the results.
APPENDIX A Euler code and meshes
Appendix A contains an overview of the Euler/adjoint code structure. This did not come with the original code and is meant to make implementation easier for future users. The appendix also shows examples of the meshes used for both two-dimensional and three-dimensional cases.
PREPROC
DEFORM
FLOW
POSTPROC
ADJOINT
GRADIENT
Figure A.1: Main program flow of the Euler/adjoint solver, showing the files compiled into each component
READMESH
141
142
APPENDIX A. EULER CODE AND MESHES
USE
read_mesh
CALL
ffa_m
USE
CALL
USE
ffread_m
FFREAD
..
fffsds_m
FFFSDS
..
ffnsds_m
FFNSDS
..
ffrmds_m
FFRMDS
..
..
preprocessor_3d
gridmetrics_3d
DUALGRID
mesh_util math_util
DUALBOUNDARY
mesh_util
ELEMENT_EDGES ELEMENT_FACES VOLUME
math_util
BCNORM
math_util
NODETYPE_EVAL CARLENGTH_EVAL
mesh_util
NDISTONE_EVAL
math_util
REORDERING
math_util
PLOTMATRIX PLOTMESH PLOTDUAL
shape_deform
parameterization
AIRFOIL_PARAMETERS DEFORMATIONS_3
CSRT
CONSTRAINTSEVAL
AIRFOIL_PARAMETERS
MAXDISPEVAL VOLUMEVAL
mesh_util
SPRING
gridmetrics_3d
DUALGRID
AIRF_CONSTRAINT mesh_util
BCNORM
mesh_util
ELEMENT_EDGES
math_util math_util DUALBOUNDARY
postprocess
mesh_util
math_util
ELEMENT_FACES VOLUME
math_util
BCNORM
math_util
math_util
euler_util
Figure A.2: Functional structure of the components READMESH, PREPROC, DEFORM, and POSTPROC
flow_solver
euler_util
math_util
computeresidual
ssor_driver
reconstruction
solver_util
USE
euler_util euler_util
FORCEVAL
euler_util euler_util
LIMITERS
interfaces
math_util
ssor_routines
MUSCL
computefluxes euler_util
FLUX_EDGE FLUX_BOUND
FLUXROE
G_FLUX_EDGE
interfaces
PRECONDITIONING_MRHS
G_FLUX_BOUND
G_FLUX_EDGE
CALL
interfaces
interfaces
USE
MATVECPRODEVAL_MRHS
DIAG_FLUX_EVAL
CALL
G_FLUXBC_FARFIELD
G_FLUXROE_APPROX
euler_util
computefluxes
math_util
euler_util
euler_util
G_FLUXBC_WALL
G_FLUXROE_APPROX
CALL
computefluxes
euler_util
computefluxes
USE
Figure A.3: Functional structure of the component FLOW
RESIDUAL
SSOR
euler_util
LEASTSQUARES_3D
math_util
euler_util
GREENGAUSS
euler_util
LEASTSQUARES_2D
euler_util
WRITECPANDMACH
WRITECPANDMACH_3D
TIMESTEP
euler_util
CHECKPRES
USE
INITIALIZE_FLOW
ALLOCATE_DEALLOCATE_ARRAY
SETTINGS_FLOW
READ_DATA_FLOW
CALL
math_util
euler_util
math_util
euler_util
math_util
euler_util
USE
143
adjoint_solver
math_util
ssor_mrhs_driver
sensitivity
solver_util
reconstruction
adjoint_util
USE
euler_util
USE
euler_util euler_util
LIMITERS
euler_util
euler_util
LEAST_SQUARES_3D_REVERSE
interfaces
interfaces
PRECONDITIONING_MRHS
interfaces
math_util
euler_util
math_util
euler_util
math_util
euler_util
USE
MATVECPRODEVAL_MRHS
DIAG_FLUX_EVAL
G_FLUXBC_FARFIELD
G_FLUXBC_WALL
G_FLUXROE_APPROX
MUSCL
G_FLUXBC_FARFIELD
CALL
G_FLUX_EDGE
G_FLUX_BOUND
G_FLUX_EDGE
CALL
G_FLUXBC_FARFIELD
G_FLUXROE_APPROX
euler_util
computefluxes euler_util
G_FLUXBC_WALL
G_FLUXROE_APPROX
CALL
computefluxes
euler_util
computefluxes
USE
Figure A.4: Functional structure of the component ADJOINT
math_util
ssor_routines
euler_util
LEAST_SQUARES_2D_REVERSE
SSOR_MRHS
euler_util
computefluxes
euler_util
computefluxes
GREENGAUSS_REVERSE
RESIDUAL_REVERSE
RESIDUAL_FREESTREAM
DFORCEVAL_ALPHA
DFORCEVAL_U
TIMESTEP
READ_ARRAY_MESH
math_util
euler_util
LEASTSQUARES_3D
READ_DATA_FLOW
euler_util
GREENGAUSS
LEASTSQUARES_2D
ALLOCATE_DEALLOCATE_ARRAY_SENSITIVITY
SETTINGS_SENSITIVITY
READ_INPUT_SENSITIVITY
INITIALIZE_SENSITIVITY
CALL
math_util
euler_util
math_util
euler_util
math_util
euler_util
USE
144 APPENDIX A. EULER CODE AND MESHES
shape_gradient
euler_util
math_util
DUALGRID
gridmetrics_3d
mesh_util
math_util
mesh_util
BCNORM
VOLUME
ELEMENT_FACES
ELEMENT_EDGES
math_util
math_util
math_util
Figure A.5: Functional structure of the component GRADIENT
DUALBOUNDARY
SPRING
BCNORM
VOLUMEVAL
mesh_util
AIRF_CONSTRAINT
MAXDISPEVAL mesh_util
CSRT
euler_util
AIRFOIL_PARAMETERS
computefluxes
FLUX_EDGE
USE
FLUX_BOUND
DEFORMATIONS_3
AIRFOIL_PARAMETERS
parameterization
interfaces
euler_util
CALL
CONSTRAINTSEVAL
RESIDUAL
euler_util
computeresidual
euler_util
LIMITERS
FORCEVAL
euler_util
LEASTSQUARES_3D
math_util
euler_util
GREENGAUSS
USE
LEASTSQUARES_2D
CALL
solver_util
reconstruction
USE
FLUXROE
MUSCL
CALL
math_util
euler_util
USE
145
146
APPENDIX A. EULER CODE AND MESHES
Figure A.6: Global view of two-dimensional mesh
147
Figure A.7: Detailed view of two-dimensional mesh
148
APPENDIX A. EULER CODE AND MESHES
Figure A.8: Global view of three-dimensional mesh
149
Figure A.9: Detailed view of three-dimensional mesh
150
APPENDIX A. EULER CODE AND MESHES
Figure A.10: Detailed view of three-dimensional mesh
APPENDIX B Adjoint validation results
Appendix B contains the results of the validation of the implementation of the CSRT method into the adjoint code.
100
40
50
20
0 1
2
3
4
5
6
7
8
9
10
FD Adjoint
-50
L/D gradient, [~]
L/D gradient, [~]
153
-100
0 1
2
3
4
5
6
7
8
9
10
FD Adjoint
-20
-40
-150
-60
Coefficient number, [~]
Coefficient number, [~]
Figure B.1: 2D FD and adjoint gradients of L/D w.r.t. Bernstein coefficients (upper side)
Figure B.2: 2D FD and adjoint gradients of L/D w.r.t. B-spline coefficients (upper side)
10
0 1
2
3
4
5
6
7
8
9
10 5
-10 1
-20 FD Adjoint -30
L/D gradient, [~]
L/D gradient, [~]
0 2
3
4
5
6
7
8
9
10
-5 FD Adjoint
-10 -15 -20
-40 -25
-50
-30
Coefficient number, [~]
Coefficient number, [~]
Figure B.4: 2D FD and adjoint gradients of L/D w.r.t. B-spline coefficients (lower side)
1.2
0.5
1.0
0.4
0.8
0.3
0.6 FD Adjoint
0.4
Cl gradient, [~]
Cl gradient, [~]
Figure B.3: 2D FD and adjoint gradients of L/D w.r.t. Bernstein coefficients (lower side)
0.2 FD Adjoint
0.1
0.2
0.0
0.0
-0.1
1 1
2
3
4
5
6
7
8
9
2
3
4
5
6
7
8
9
10
10
-0.2
-0.2 Coefficient number, [~]
Figure B.5: 2D FD and adjoint gradients of cl w.r.t. Bernstein coefficients (upper side)
Coefficient number, [~]
Figure B.6: 2D FD and adjoint gradients of cl w.r.t. B-spline coefficients (upper side)
154
APPENDIX B. ADJOINT VALIDATION RESULTS
0.0
0.04 1
2
3
4
5
6
7
8
9
10
-0.2
0.00 1
2
3
4
5
6
7
8
9
10
-0.6 FD Adjoint
-0.8
Cl gradient, [~]
Cl gradient, [~]
-0.4 -0.04 FD Adjoint
-0.08
-0.12
-1.0 -0.16
-1.2 -1.4
-0.20 Coefficient number, [~]
Coefficient number, [~]
Figure B.8: 2D FD and adjoint gradients of cl w.r.t. B-spline coefficients (lower side)
0.25
0.08
0.20
0.06
0.15
0.04
FD Adjoint
0.10
0.05
Cd gradient, [~]
Cd gradient, [~]
Figure B.7: 2D FD and adjoint gradients of cl w.r.t. Bernstein coefficients (lower side)
FD
0.02
Adjoint
0.00 1
2
3
1
2
3
4
5
6
7
8
9
5
6
7
8
9
10
10
-0.04
-0.05
Coefficient number, [~]
Coefficient number, [~]
Figure B.9: 2D FD and adjoint gradients of cd w.r.t. Bernstein coefficients (upper side)
Figure B.10: 2D FD and adjoint gradients of cd w.r.t. B-spline coefficients (upper side)
0.04
0.04
0.02
0.03
0.00 1
2
3
4
5
6
7
8
9
10
FD Adjoint
-0.02
Cd gradient, [~]
Cd gradient, [~]
4
-0.02
0.00
0.02 FD
0.01
Adjoint
0.00 1
-0.04
2
3
4
5
6
7
8
9
10
-0.01 -0.02
-0.06 Coefficient number, [~]
Figure B.11: 2D FD and adjoint gradients of cd w.r.t. Bernstein coefficients (lower side)
Coefficient number, [~]
Figure B.12: 2D FD and adjoint gradients of cd w.r.t. B-spline coefficients (lower side)
155
3
10
9
2
9
2
8
1
8
1
Coefficient number, [~] (in y−direction)
L/D gradient
0
7
−1
6
−2 5
−3
4
−4
−2 5
2
−6
2
−7
1 1
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
10
Figure B.13: 3D FD gradients of L/D w.r.t. Bernstein coefficients (upper side)
10
−3
4 3
3
−1
6
−5
2
0
7
3
1 1
3 L/D gradient
Coefficient number, [~] (in y−direction)
10
−4 −5 −6 2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
Figure B.14: 3D Adjoint gradients of L/D w.r.t. Bernstein coefficients (upper side)
10
1
1 L/D gradient
L/D gradient 9
9
0.5
0.5 8
7
0
6 −0.5 5 4
−1
Coefficient number, [~] (in y−direction)
Coefficient number, [~] (in y−direction)
8
7
0
6 −0.5 5 4
−1
3
3
−1.5
−1.5 2
2 1 1
2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
10
1 1
−2
Figure B.15: 3D FD gradients of L/D w.r.t. B-spline coefficients (upper side)
10
0.04
6 0.03 5 0.02
4
0.01
3
Figure B.17: 3D FD gradients of CL w.r.t. Bernstein coefficients (upper side)
0.05 0.04
0.02
4
0.01
3
1 1
10
0.06
0.03
−0.01
9
0.07
5
1 1
8
−2
6
2 4 5 6 7 Coefficient number, [~] (in x−direction)
10
7
2
3
9
8
0 2
8
Cl gradient, [~]
Coefficient number, [~] (in y−direction)
0.05
4 5 6 7 Coefficient number, [~] (in x−direction)
9
0.06
7
3
10
0.07
8
2
Figure B.16: 3D Adjoint gradients of L/D w.r.t. B-spline coefficients (upper side)
Cl gradient, [~] 9 Coefficient number, [~] (in y−direction)
−7
10
0 2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
10
−0.01
Figure B.18: 3D Adjoint gradients of CL w.r.t. Bernstein coefficients (upper side)
156
APPENDIX B. ADJOINT VALIDATION RESULTS −3
−3
x 10 10
10
C gradient
L
9
8
8 6
Coefficient number, [~] (in y−direction)
Coefficient number, [~] (in y−direction)
L
9
8
8 7 4
6 5
2
4 0 3
6 7 6
4
5
2
4 0 3
−2
2 1 1
x 10 10
10
C gradient
2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
10
−2
2
−4
1 1
Figure B.19: 3D FD gradients of CL w.r.t. B-spline coefficients (upper side)
2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
−4
10
Figure B.20: 3D Adjoint gradients of CL w.r.t. B-spline coefficients (upper side) −3
−3
x 10 10
10
x 10 10
10
Cd gradient, [~]
Cd gradient, [~] 9
8 7
5
6 5 4
0
Coefficient number, [~] (in y−direction)
Coefficient number, [~] (in y−direction)
9
3
8 7
5
6 5 4
0
3
2
2
1 1
2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
1 1
−5
10
Figure B.21: 3D FD gradients of CD w.r.t. Bernstein coefficients (upper side)
2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
Figure B.22: 3D Adjoint gradients of CD w.r.t. Bernstein coefficients (upper side)
−3
Coefficient number, [~] (in y−direction)
9
2
8
1 6 0.5 5 0 4 −0.5
3
1 0.5 5 0 4
2
−1.5
1 1
8
9
10
Figure B.23: 3D FD gradients of CD w.r.t. B-spline coefficients (upper side)
−0.5
3
1 1
4 5 6 7 Coefficient number, [~] (in x−direction)
1.5
6
2
3
2
7
−1
2
CD gradient
8
1.5
7
x 10 2.5
10
CD gradient
Coefficient number, [~] (in y−direction)
9
−3
x 10 2.5
10
−5
10
−1
2
3
4 5 6 7 Coefficient number, [~] (in x−direction)
8
9
10
−1.5
Figure B.24: 3D Adjoint gradients of CD w.r.t. B-spline coefficients (upper side)
BIBLIOGRAPHY
[1] Anonymous, “Optimization ToolboxT M User’s Guide”, pp. 3.29–3.37, MathWorks, 2010 [2] Anonymous, “VSAERO A Code for Calculating the Nonlinear Aerodynamic Characteristics of Arbitrary Configurations User’s Manual Version 7.1”, Analytical Methods Inc., 2005 [3] Abramowitz, M., Stegun, I.A., “Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables”, Ch. 22, National Bureau of Standards, Applied Mathematics Series, 1965 [4] Barr, A.H., “Global and Local Deformation of Solid Primitives”, Computer Graphics, Vol. 18, No. 3, pp. 21–30, 1984 [5] Berends, J.P.T.J., Tooren, M.J.L. van, and Belo, D.N.V., “A Distributed MultiDisciplinary Optimisation of a Blended Wing Body UAV using a Multi-Agent Task Environment”, Proceedings of the 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Newport, Rhode Island, 1–4 May 2006 [6] Berglund, T., Jonsson, H., and S¨oderkvist, I., “An Obstacle-Avoiding Minimum Variation B-spline Problem”, Proceedings of the 2003 International Conference on Geometric Modeling and Graphics, London, 16–18 July 2003 [7] Buckley, H.P., Zhou, B.Y., and Zingg, D.W., “Airfoil Optimization Using Practical Aerodynamic Design Requirements”, Journal of Aircraft, Vol. 47, No. 5, pp. 1707–1719, 2010 [8] Carpentieri, G., Koren, B., and Tooren, M.J.L. van, “Adjoint-Based Aerodynamic Shape Optimization on Unstructured Meshes”, Journal of Computational Physics, Vol. 224, No. 1, pp. 267–287, 2007 [9] Carpentieri, G., Koren, B., and Tooren, M.J.L. van, “Development of the Discrete Adjoint for a Three-Dimensional Unstructured Euler Solver”, Journal of Aircraft, Vol. 45, No. 1, pp. 237–245, 2008 [10] Carpentieri, G., “An Adjoint-Based Shape-Optimization Method for Aerodynamic Design”, Delft University of Technology, Delft, The Netherlands, 2009 [11] Celniker, G., and Welch, W., “Linear Constraints for Deformable B-spline Surfaces”, Proceedings of the 1992 Symposium on Interactive 3D Graphics, Cambridge, Massachusetts, March 29–April 20 1995 [12] Coquillart, S., “Extended Free-Form Deformation: A Sculpturing Tool for 3D Geometric Modeling”, SIGGRAPH, Vol. 24, No. 4, pp. 187–196, 1990 [13] Curle, N., “Accurate Solutions of the Laminar-Boundary-Layer Equations, for Flows having a Stagnation Point and Separation”, Aeronautical Reseach Council Reports and Memoranda, 1960 [14] Darwin, C., “The Origin of Species by Means of Natural Selection”, John Murray, 1859
160
BIBLIOGRAPHY
[15] Ding, Y., “Shape Optimization of Structures: A Literature Survey”, Computers and Structures, Vol. 24, No. 6, pp. 985–1004, 1986 [16] Enloe, C.L., McLaughlin, T.E., VanDyken, R.D., Kachner, K.D., Jumper, E.J., Corke, T.C., Post, M., and Haddad, O., “Mechanisms and Responses of a Single Dielectric Barrier Plasma Actuator: Geometric Effects”, AIAA Journal, Vol. 42, No. 3, pp. 595–604, 2004 [17] Ferlauto, M., Iollo, A., and Zannetti, L., “Coupling of Optimization and Inverse Problem for Aerodynamic Shape Design”, Proceedings of the 38th Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 10–13 January 2000 [18] Forrester, A.I.J., Bressloff, N.W., and Keane, A.J., “Optimization Using Surrogate Models and Partially Converged Computational Fluid Dynamics Solutions”, Proceedings of the Royal Society A, Vol. 462, No. 2071, pp. 2177–2204, 2006 [19] Forrester, A.I.J., S´obester, A., and Keane, A.J., “Engineering Design via Surrogate Modelling”, John Wiley & Sons, 2008 [20] Fray, J.M.J., Slooff, J.W., Boerstoel, J.W., and Kassies, A., “Inverse Method with Geometric Constraints for Transonic Aerofoil Design”, Numerical Methods in Engineering, Vol. 22, No. 2, pp. 327–339, 1986 [21] Green, J.E., “Civil Aviation and the Environment – the Next Frontier for the Aerodynamicist”, Aeronautical Journal, Vol. 110, No. 1110, pp. 469–486, 2006 [22] Green, M., and Cheng, R., “Genetic Algorithm & Engineering”, John Wiley & Sons, 1997 [23] Hafez, M., Shatalov, A., and Nakajima, M., “Improved Numerical Simulations of Incrompressible Flows Based on Viscous/Inviscid Interaction Procedures”, Computers & Fluids, Vol. 36, No. 10, pp. 1588–1591, 2007 [24] Haftka, R.T., and Grandhi, R.V., “Structural Shape Optimization - A Survey”, Computer Methods in Applied Mechanics and Engineering, Vol. 57, No. 1, pp. 91– 106, 1986 [25] Hess, J.L., “Panel Methods in Computational Fluid Dynamics”, Annual Reviews of Fluid Mechanics, Vol. 22, pp. 255-274, 1990 [26] Hicks, R.M., and Henne, P.A., Wing Design by Numerical Optimization, Journal of Aircraft, Vol. 15, No. 7, pp. 407–412, 1978 [27] Hu, X., and Eberhart, R., “Solving Constrained Nonlinear Optimization Problems with Particle Swarm Optimization”, Proceedings of the 6th World Multiconference on Systemics, Cybernetics and Informatics, Orlando, Florida, 14–18 July 2002 [28] Jameson, A., “Aerodynamic Control via Control Theory”, Journal of Scientific Computing, Vol. 3, No. 3, pp. 233–260, 1988
BIBLIOGRAPHY
161
[29] Jameson, A., and Vassberg, J.C., “Computational Fluid Dynamics for Aerodynamic Design: Its Current and Future Impact”, Proceedings of the 39th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 8–11 January 2001 [30] Jameson, A., Martinelli, L, Vassberg, J., and Cliff, S.E., “Aerodynamic Simulation and Shape Optimization for High Speed Flow”, Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 9–12 January 2006 [31] Johnson, F.T., and Tinoco, E.N., “Thirty Years of Development and Application of CFD at Boeing Commercial Airplanes, Seattle”, Proceedings of the 16th AIAA Computational Fluid Dynamics Conference, Orlando, Florida, 23–26 June 2003 [32] Jones, D.R., Schonlau, M., and Welch, W.J., “Efficient Global Optimization of Expensive Black-Box Functions”, Journal of Global Optimization, Vol. 13, No. 4, pp. 455–492, 1998 [33] Joseph, D.D., Liao, T.Y., and Hu, H.H., “Drag and Moment in Viscous Potential Flow”, European Journal of Mechanics B/Fluids, Vol. 12, No. 1 , pp. 97–106, 1993 [34] Juh´asz, I., and Hoffmann, M., “Constrained Shape Modification of Cubic BSpline Curves by Means of Knots”, Computer-Aided Design, Vol. 36, No. 5, pp. 437–455, 2004 [35] Keane, A.J., “Genetic Algorithm Optimization of Multi-Peak Problems: Studies in Convergence and Robustness”, Artificial Intelligence in Engineering, Vol. 9, No. 2, pp. 75–83, 1995 [36] Khurana, M.S., Winarto, H., and Sinha, A.K., “Application of Swarm Approach and Artificial Neural Networks for Airfoil Shape Optimization”, Proceedings of the 12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Victoria, British Columbia, Canada, 10–12 September 2008 [37] Kolla, M.L., Yokota, J.W., Lassaline, J.V., and Fejtek I., “Curvature Constraints for Airfoil Shape Optimization in Turbulent Flow”, Canadian Aeronautics and Space Journal, Vol. 54, No. 1, pp. 1–7, 2008 [38] Kuhn, H., and Tuckers, A.W., “Nonlinear Programming”, Proceedings of the 2nd Berkeley Symposium, pp. 481–492, University of California Press, 1951 [39] Kulfan, B.M., ““Fundamental” Parametric Geometry Representation for Aircraft Component Shapes”, Proceedings of the 11th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Portsmouth, Virginia, U.S.A., 6–8 September 2006 [40] Kulfan, B.M., “A Universal Parametric Geometry Representation Method CST”, Proceedings of the 45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 8–11 January 2007 [41] Kulfan, B.M., “Universal Parametric Geometry Representation Method”, Journal of Aircraft, Vol. 45, No. 1, pp. 142–158, 2008
162
BIBLIOGRAPHY
[42] Lamousin, H.J., and Waggenspack, W.N., “NURBS-Based Free-Form Deformation”, IEEE Computer Graphics and Applications, Vol. 14, No. 6, pp. 95–108, 1994 [43] Lighthill, M.J., “On Displacement Thickness”, Journal of Fluid Mechanics, Vol. 4, No. 4, pp. 383–392, 1958 [44] Lin, Y-C., Wang, F-S., and Hwang, K-S., “A Hybrid Method of Evolutionary Algorithms for Mixed-Integer Nonlinear Optimization Problems”, Proceedings of the 1999 Congress on Evolutionary Computation, Washington, DC, 6–9 July 1999 [45] Lozano, C., and Ponsin, J., “Remarks on the Numerical Solution of the Adjoint Quasi-One-Dimensional Euler Equations”, International Journal for Numerical Methods in Fluids, 2011 [46] Martins, J.R.R.A., Alonso, J.J., and Reuther, J.J., “A Coupled-Adjoint Sensitivity Analysis Method for High-Fidelity Aero-Structural Design”, Optimization and Engineering, Vol. 6, No. 1, pp. 33–62, 2005 [47] Mavriplis, D.J., “Discrete Adjoint-Based Approach for Optimization Problems on Three-Dimensional Unstructured Meshes”, AIAA Journal, Vol. 45, No. 4, pp. 740–750, 2007 [48] Morris, A.J., “MOB A European Distributed Multi-Disciplinary Design and Optimisation Project”,Proceedings of the 9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Atlanta, Georgia, 4–6 September 2002 [49] Morris, A.M., Allen, C.B., and Rendall, T.C.S., “A Generic Parallel Framework for Shape Parameterisation and CFD-Based Optimisation”, Proceedings of the 26th AIAA Applied Aerodynamics Conference, Honolulu, Hawaii, 18–21 August 2008 [50] Morris, M.D., and Mitchell, T.J., “Exploratory Designs for Computational Experiments”, Oak Ridge National Laboratory, Oak Ridge, Tennessee, 1992 [51] Mortenson, M.E., “Geometric Modeling”, Second Edition, John Wiley & Sons, 1997 [52] Mousvi, A., Castonguay, P., and Nadarajah, S.K., “Survey of Shape Parameterization Techniques and Its Effect of Three Dimensional Aerodynamic Shape Optimization”, Proceedings of the 18th AIAA Computational Fluid Dynamics Conference, Miami, Florida, 25–28 June 2007 [53] Nathman, J.K., “Subsonic Panel Methods – Second (Order) Thoughts”, Proceedings of the AIAA and SAE, 1998 World Aviation Conference, Anaheim, California, 28-30 September, 1998 [54] Nathman, J.K., “Improvement of a Panel Method by Including Panel Warp”, Proceedings of the 42nd AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 5–8 January 2004 [55] Nathman, J.K., and Matarrese, M., “Hybrid Grid (Structured and Unstructured) Calculations with a Potential-Based Panel Method”, Proceedings of the 22nd Applied Aerodynamics Conference and Exhibit, Providence, Rhode Island, 16–19 August 2004
BIBLIOGRAPHY
163
[56] Nathman, J.K., “Induced Drag of High-Aspect Ratio Wings”, Proceedings of the 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 10–13 January 2005 [57] Nemec, M., and Zingg, D.W., “Towards Efficient Aerodynamic Shape Optimization Based on the Navier-Stokes Equations”, Proceedings of the 15th AIAA Computational Fluid Dynamics Conference, Anaheim, California, 11–14 June 2001 [58] Newman, J.C., Taylor, A.C., Barnwell, R.W., Newman, P.A., and Hou G.J-W., “Overview of Sensitivity Analysis and Shape Optimization for Complex Aerodynamic Configurations”, Journal of Aircraft, Vol. 36, No. 1, pp. 87–96, 1999 [59] Nocedal, J., and Wright, S.J., “Numerical Optimization”, Second Edition, Springer, 2006 [60] Painchaud-Ouellet, S., Tribes, C., Tr´epanier, J-Y., and Pelletier, D., “Airfoil Shape Optimization Using a Nonuniform Rational B-Splines Parameterization Under Thickness Constraint”, AIAA Journal, Vol. 44, No. 10, pp. 2170–2178, 2006 [61] Peigin, S., and Epstein, B., “Computational Fluid Dynamics Drive Optimization of Blended Wing Body Aircraft”, AIAA Journal, Vol. 44, No. 1, pp. 2736–2745, 2006 [62] Pourazady, M., and Xu, X., “Direct Manipulation of B-Spline and NURBS Curves”, Advances in Engineering Software, Vol. 31, Nr. 2, pp. 107–118, 2000 ¨ [63] Prandtl, L., “Uber Fl¨ ussigkeitsbewegung bei sehr kleiner Reibung”, Verhandlungen des Dritten Internationalen Mathematiker-Kongresses, Heidelberg, Germany, 8–13 August 1904 [64] Qin, N., Vavalle, A., Le Moigne, A., Laban, M., Hackett, K., and Weinerfelt, P., “Aerodynamic Studies for Blended Wing Body Aircraft”, Proceedings of the 9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Atlanta, Georgia, 4–6 September 2002 [65] Rodgers, J.L., and Nicewander, W.A., “Thirteen Ways to Look at the Correlation Coefficient”, The American Statistician, Vol. 42, No. 1, pp. 59–66, 1988 [66] Rodriguez, D.L., “Response Surface Based Optimization with a Cartesian CFD Method”, Proceedings of the 41st Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 6–9 January 2003 [67] Ruijgrok, G.J.J., “Elements of Airplane Performance”, Second Edition, Delft University Press, 1996 [68] Samareh, J.A., “Survey of Shape Parameterization Techniques for High-Fidelity Multidisciplinary Shape Optimization”, AIAA Journal, Vol. 39, No. 5, pp. 877–884, 2001 [69] Samareh, J.A., “Aerodynamic Shape Optimization Based on Free-Form Deformation”, Proceedings of the 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Albany, New York, 30 August–1 September 2004
164
BIBLIOGRAPHY
[70] Schlichting, H., Gersten, K., “Boundary Layer Theory”, 8th Revised and Enlarged Edition, Springer, 2003 [71] Schmitt, V., and Charpin, F., “Pressure Distribution on the ONERA-M6-Wing at Transonic Mach Numbers”, Experimental Data Base for Computer Program Assessment, AGARD-AR-138, May 1979 [72] Sederberg, T.W., and Parry, S.R., “Free-Form Deformation of Solid Geometric Models”, Computer Graphics, Vol. 20, No. 4, pp. 151–160, 1986 [73] Sevant, N.E., Bloor, M.I.G., and Wilson, M.J., “A Hierarchical Approach to Optimal Aerodynamic Design”, Proceedings of the 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, St. Louis, Missouri, 2–4 September 1998 [74] Sobieczky, H., “Parametric Airfoils and Wings”, Notes on Numerical Fluid Mechanics, Vol. 68, pp. 71–88, 1998 [75] Straathof, M.H., Tooren, M.J.L. van, Voskuijl, M., and Koren, B., “Aerodynamic Shape Parameterisation and Optimisation of Novel Configurations”, Proceedings of the RAeS Aerodynamic Shape Parameterisation and Optimisation of Novel Configurations Conference, London, 27–28 October 2008 [76] Straathof, M.H., Tooren, M.J.L. van, Voskuijl, M., and Vos, R., “Development and Implementation of a Novel Parametrization Technique for Multidisciplinary Design Initialization”, Proceedings of the 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Orlando, Florida, 12– 15 April 2010 [77] Straathof, M.H., and Tooren, M.J.L. van, “Extension to the Class-ShapeTransformation Method Based on B-Splines”, AIAA Journal, Vol. 49, No. 4, pp. 780–790, 2011 [78] Straathof, M.H., and Tooren, M.J.L. van, “Adjoint Optimization of a Wing Using the CSRT Method”, Proceedings of the 29th AIAA Applied Aerodynamics Conference, Honolulu, Hawaii, 27–30 June 2011 (to be published in Journal of Aircraft in 2012) [79] Straathof, M.H., Carpentieri, G., and Tooren, M.J.L. van, “Aerodynamic Shape Optimization Using the Adjoint Euler Equations”, Proceedings of Eurogen 2011: Evolutionary and Deterministic Methods for Design, Optimization and Control, Capua, Italy, 14–16 September, 2011 [80] Straathof, M.H., “Parametric Study of the Class-Shape-Refinement Transformation Method”, to be published in Optimization in 2012 [81] Sung, C., and Kwon, J.H., “Efficient Aerodynamic Design Method Using a Tightly Coupled Algorithm”, AIAA Journal, Vol. 40, No. 9, pp. 1839–1845, 2002 [82] Thwaites, B., “Incompressible Aerodynamics: An Account of the Theory and Observation of the Steady Flow of Incompressible Fluid Past Aerofoils, Wings, and Other Bodies”, Dover Publications, 1987
BIBLIOGRAPHY
165
[83] Tinoco, E.N., “The Changing Role of Computational Fluid Dynamics in Aircraft Development”, Proceedings of the 16th Applied Aerodynamics Conference, Albuquerque, New Mexico, June 15–18 1998 [84] Tinoco, E.N., “An Assessment of CFD Prediction of Drag and Other Longitudinal Characteristics”, Proceedings of the 39th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 8–11 January 2001 [85] Tuli, M., Reddy, N.V., and Saxena, A., “Constrained Shape Modification of B-Spline Curves”, Computer-Aided Design & Applications, Vol. 3, Nos. 1–4, pp. 437–446, 2006 [86] Venter, G., and Sobieszczanski-Sobieski, J., “Particle Swarm Optimization”, Proceedings of the 43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Denver, Colorado, 22–25 April 2002 [87] Viviand, H., “Numerical Solutions of Two-Dimensional Reference Test Cases”, Test Cases for Inviscid Flow Field Methods, AGARD-AR-211, May 1985 [88] Voskuijl, M., La Rocca, G., and Dircken, F., “Controllability of Blended Wing Body Aircraft”, Proceedings of the 26th International Congress of the Aeronautical Sciences, Anchorage, Alaska, 14–19 September 2008 [89] Watt, A., and Watt, M., “Advanced Animation and Rendering Techniques”, Addison-Wesley, Ch. 17, 1992
Acknowledgements
I am very thankful to my promotor Michel van Tooren, without whom my research would not have been possible. A lot of things changed over the past 4 years, but your enthusiasm didn’t. I would also like to thank Mark Voskuijl and Roelof Vos, for always providing me with valuable feedback. This thesis would not have looked the same without it. I am grateful also to Lin Pijpaert and Michiel Haanschoten, who were absolutely indispensable, as well as to all members of my committee for reading this thesis and attending the defense. I performed my research within the framework of CleanEra, which provided an inspirational environment and resulted in a lot of new friends. Thank you Marios, Hui, Fran¸cois, Durk, Gustavo, Ronald, Farid, Chara, Sonell, Marcel and Dipanjay for never making me dread going to the office in the morning. A lot of thanks goes to Albert, Kartik and Roderik, who have always shown interest in what I was doing and have just been great friends for the past 10 years. Thanks also to Alastair, Robin and Joost for coping with me as a flatmate during the first 2 years of my PhD. Thanks, Lize and Merijn, for always being supportive and proud of your brother. That means a lot to me. Thanks also to Noah and Julia, for making home a more lively place. Jona, you have made the last year of my PhD by far the best. Thank you so much for that. Last, but definitely not least, I am eternally grateful to my parents, Hans and Ine. You have always provided me with everything I needed to do the things I wanted to do. And more importantly, you were always there when I needed you. I cannot thank you enough for that.
About the author
Michiel Straathof was born on August 16th 1983 in Leiden, the Netherlands. It didn’t take him long to notice the shiny things moving through the sky above him, and one of his first words was “ieguig”: a rather poor attempt at the Dutch word for airplane. From September 1996 to July 2002 he attended high school at the Stedelijk Gymnasium Leiden, where his interest in science and engineering really became apparent. In 2002, when most of his class mates had decided to become doctors or lawyers, he (and admittedly some of his best friends) went off to pursue a career in engineering. Continuing his childhood fascination he enrolled at the faculty of Aerospace Engineering of Delft University of Technology, where he completed his Bachelor degree in 2005 and his Master degree in 2007. He then got offered a PhD as part of the CleanEra project, where he would work on the aerodynamic design of aircraft. He accepted the offer and performed research at the faculty for 4 more years, culminating in this thesis.
