Shape Parameterization in
Aircraft Design: A Novel Method,
Based on B-Splines
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
Prof. dr. ir. M.J.L. van Tooren
Samenstelling promotiecommissie:
Rector Magnificus,
voorzitter
Prof. dr. ir. M.J.L. van Tooren,
Technische Universiteit Delft, promotor
Prof. dr. T. B¨
ack,
Universiteit Leiden
Prof. dr. ir. B. Koren,
Universiteit Leiden / CWI Amsterdam
Prof. dr. ir. H.W.M. Hoeijmakers,
Universiteit Twente
Prof. dr. ir. A.W. Heemink,
Technische Universiteit Delft
Dr. ir. G. Carpentieri,
Cardano, Londen
B.M. Kulfan,
Boeing (retired), Seattle, WA
keywords: shape, parameterization, aircraft, design, B-splines,
Class-Shape-Refinement-Transformation, adjoint, euler, optimization
ISBN/EAN: 978-94-6191112-4
Copyright c
⃝ 2012 by Michiel H. Straathof
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-Shape-Transformation (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 Class-Shape-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 opti-mization. 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 to-day. 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 tech-nique. 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
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 optimiza-tion 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 im-provement resulting from the first step, and a smaller imim-provement 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 refine-ment 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 Class-Shape-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 voor-delen 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 vlieg-tuigontwerp 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 opti-malisatie, 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 tech-nieken wordt in meer detail beschreven, waarbij in het bijzonder aandacht wordt besteed aan de Class-Shape-Transformation methode. De Class-Shape-Refinement-Transformation 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 probleem-specifiek is.
Hoofdstuk 4 presenteert een nieuwe manier van omgaan met volumetrische rand-voorwaarden 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 verhou-ding gerealiseerd, maar lokale minima in de ontwerpruimte bleken een potentieel pro-bleem. 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 laat-ste test case liet echter een afwijkend resultaat zien. De eerlaat-ste 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 Bernstein-oppervlak van de shape function niet goed kan omgaan met de ontwerpruimte, wat vervolgens gecompenseerd wordt door het B-spline oppervlak van de refinement
func-tion.
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 test-cases 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 vii 1 Introduction 1 1.1 Background . . . 3 1.2 Current practice . . . 9 1.2.1 Parameterization . . . 9 1.2.2 Aerodynamic analysis . . . 10 1.2.3 Optimization . . . 11 1.3 CSRT method . . . 12 1.4 Thesis overview . . . 15 2 Parameterization 17 2.1 Background . . . 19 2.2 Common methods . . . 21 2.2.1 Discrete approach . . . 21 2.2.2 Analytical approach . . . 21 2.2.3 Geometric approach . . . 222.2.4 Polynomial and spline approaches . . . 24
2.3 CST method . . . 25 2.3.1 CST in two dimensions . . . 25 2.3.2 CST in three dimensions . . . 27 3 CSRT method 31 3.1 B-splines . . . 33 3.1.1 B-spline curves . . . 33 3.1.2 B-spline surfaces . . . 36 3.2 Refinement function . . . 38 3.3 CSRT performance . . . 40 3.3.1 Fitting process . . . 42 3.3.2 Correlation factor . . . 44 3.3.3 Parametric study . . . 44
4 Volume constraints 57
4.1 Typical shape constraints . . . 59
4.2 Volume constraint handling in 2D . . . 60
4.3 Volume constraint handling in 3D . . . 64
5 Design frameworks 69 5.1 Low fidelity design framework . . . 71
5.1.1 Panel methods (VSAERO) . . . 71
5.1.2 Global optimization . . . 78
5.1.3 Response surface methods . . . 79
5.2 High fidelity design framework . . . 83
5.2.1 Euler method . . . 83
5.2.2 Gradient-based optimization . . . 93
6 Test cases 103 6.1 Low fidelity framework: 3D wing . . . 105
6.1.1 Test case definition . . . 105
6.1.2 Test case results . . . 107
6.2 High fidelity framework: NACA0012 . . . 109
6.2.1 Test case definition . . . 111
6.2.2 Test case results . . . 111
6.3 High fidelity framework: ONERA M6 . . . 118
6.3.1 Test case definition . . . 118
6.3.2 Test case results . . . 120
6.4 High fidelity framework: BWB . . . 125
6.4.1 Test case definition . . . 126
6.4.2 Results of the first test case . . . 127
6.4.3 Results of the second test case . . . 131
7 Conclusions 135
A Euler code and meshes 139
B Adjoint validation results 151
Bibliography 157
Acknowledgements 167
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.”
1.1. BACKGROUND
3
1.1
Background
During the early development of aircraft, not much attention was paid to aerody-namic 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 air-craft 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]:
R = H g ηtotln ( W0 Wf ) CL CD (1.1) 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).
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
1.1. BACKGROUND
5
Figure 1.2: The Douglas DC-2 in flight
0.0300 0.0200 0.0100 0 1900 1910 1920 1930 1940 1950 1960 1970 YEAR C DP TURBULENT FLOW SKIN FRICTION BIPLANES MONOPLANE JET ENGINE RETRACTABLE GEAR
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 possi-ble. The first non-linear flow solvers used by Boeing were the well-known FLO27/28
1.1. BACKGROUND
7
Technology Boeing 1965 1970 1975 1980 1985 1990 1995 2000 SST 767 757 KC-135R 737-300 777 737NG 787 Aircraft TA80 TA139/201TA176/217 A230FLEXSTAB A488P582 A502A555 A619 TRANAIR Optimi-zation ZEUS GGNS TRANAIR A588 Boeing Tools Wave Drag Mach Box Woodward Army Fan-In-Wing Contract FLEXSTAB Contract FLO22/ 27/28 A411 PANAIR Cartesian Grid Technology TLNS3D CFL3D TLNS3D-MB OVERFLOW Unstructured Adaptive Grid N-S
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 sim-ulations. 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 poten-tial 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 unstruc-tured 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
Design variables Optimization Aerodynamic analysis results Desired aerodynamic performance achieved? no yes Geometry generation Aerodynamic analysis
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 require-ments 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 air-craft 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
9
1.2
Current practice
1.2.1
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 optimiza-tion. This process is called “shape parameterization” or “geometric representation” and is an essential part of any shape optimization problem. Interestingly, parameter-ization 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 com-plex 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 pos-sess, 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 representa-tions 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 con-straints. A number of different techniques for implementing constraint handling have been investigated in literature. Juhsz and Hoffmann [34] made use of the knot values
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
sepa-rate 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 con-straint 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. Reynolds-Averaged-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 ∂t + v· ∇v ) =−∇p + Fvisc(v) (1.2)
where ρ is the air density, v is the velocity field, p is the pressure field andFviscare 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 ∂t + ¯v· ∇¯v ) + T =−∇¯p + Fvisc(¯v) (1.4) with T = ∂u′2 ∂x ∂u′v′ ∂y ∂u′w′ ∂z ∂u′v′ ∂x ∂v′2 ∂y ∂v′w′ ∂z ∂u′w′ ∂x ∂v′w′ ∂y ∂w′2 ∂z (1.5)
1.2. CURRENT PRACTICE
11
Taking the Navier-Stokes equations and neglecting viscosity leads to the Euler equa-tions, which can be written as:
ρ ( ∂v ∂t + v· ∇v ) =−∇p (1.6)
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:
(1− M2 x)∂ 2Φ ∂x2 + (1− M 2 y)∂ 2Φ ∂y2 + (1− M 2 z)∂ 2Φ ∂z2 − 2MxMy∂ 2Φ ∂x∂y −2MyMz∂ 2Φ ∂y∂z− 2MxMz ∂2Φ ∂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:
(1− M∞2)∂ 2Φ ∂x2 + ∂2Φ ∂y2 + ∂2Φ ∂z2 = 0 (1.9)
Equation (1.9) is aptly called the linearized potential equation. Because of the as-sumptions 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
A
B
C D
Figure 1.7: Example of local maxima
?
Figure 1.8: Example of inability to guar-antee 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
1.3. CSRT METHOD
13
Figure 1.9: Example of a response surface with the black dots representing the sample points 0 0.2 0.4 0.6 0.8 1 −0.05 0 0.05 0.1 0.15 0.2
0 0.2 0.4 0.6 0.8 1 −0.05 0 0.05 0.1 0.15 0.2 Nvar = 4 Nvar = 6 N var = 8 Nvar = 14
Figure 1.11: Capturing a local
defor-mation using Bernstein polynomials for
Nvar = 4, 6, 8, 14 0 0.2 0.4 0.6 0.8 1 −0.05 0 0.05 0.1 0.15 0.2 Nvar = 20
Figure 1.12: Capturing a local
defor-mation using Bernstein polynomials for
Nvar = 20 0 0.2 0.4 0.6 0.8 1 −0.05 0 0.05 0.1 0.15 0.2 Nvar = 4 Nvar = 6 Nvar = 8 Nvar = 14
Figure 1.13: Capturing a local
deforma-tion using a B-spline for Nvar = 4, 6, 8,
14 0 0.2 0.4 0.6 0.8 1 −0.05 0 0.05 0.1 0.15 0.2 Nvar = 20
Figure 1.14: Capturing a local
deforma-tion using a B-spline for Nvar = 20
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 func-tions 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
approxi-mations 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
15
1.4
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.
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.”
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 meth-ods 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 math-ematical descriptions of such products and preferably also easy access to derivatives of such models with respect to predefined design variables. Since many design vari-ables 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 ∂x = ∂J ∂Rf ∂Rf ∂Rs ∂Rs ∂Rg ∂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 vec-tor, domain element, partial differential equation (PDE), discrete, polynomial and spline, CAD-based, analytical, and freeform deformation (FFD).
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 de-formation. The only disadvantage is that the use of the parallelepiped limits the topology of deformation. To alleviate this disadvantage, Sederberg and Parry pro-posed 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 mod-els [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 al-ways 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 com-puter 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 gradient-based 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 three-dimensional geometries can be parameterized in a hierarchical manner allowing free-form 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
Common methods
2.2.1
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 con-nection 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. C2 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 C1 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 pro-cess. 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:
fj(x) = [ sin ( πx log(0.5) log(t1) )]t2 , 0≤ x ≤ 1 (2.2)
where t1and t2can be varied to fit the design problem. Figure 2.2 shows five
Hicks-Henne functions for t2= 5. The weights αjare the design variables. The perturbation
functions fj(x) can now be added to the baseline shape as follows:
y(x) = ybasis+ M ∑ j=1 ¯ αjfj(x) (2.3)
0 0.5 1 0 0.5 1 f 1(x) f2(x) f3(x) f4(x) f5(x)
Figure 2.2: Hicks-Henne functions for t2= 5
Table 2.1: Definitions of typical shape functions
First kind Second kind
Cheby- T0(x) = 1 U0(x) = 1 shev T1(x) = x U1(x) = 2x Tn+1(x) = 2xTn(x)− Tn−1(x) Un+1(x) = 2xUn(x)− Un−1(x) Legen- P0(x) = 1 Q0(x) = 12ln ( 1+x 1−x ) dre P1(x) = x Q1(x) = x2ln ( 1+x 1−x ) − 1 Pn+1(x) = (2n+1)xPn(x)−nPn−1(x) n+1 Qn+1(x) = (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.
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.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 ex-pressed 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 ¯ci will influence the overall
shape of the curve.
y(x) = p ∑ i=0 ¯ cixi (2.4)
A better way to represent a curve would be to use the Bezier representation, which is defined as follows: y(x) = p ∑ i=0 ¯ PiBi,p(x) (2.5)
where ¯Pi 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 ¯Pi, 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 re-quired degree of the Bernstein polynomials increases significantly, making it very inefficient to compute the curve. This follows from the fact that each Bernstein poly-nomial 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=0 ¯ PiNi,p(u) (2.6)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 ¯Pi by a B-spline
control polygon ¯Pi. The parametric variable u is used instead of the Cartesian
co-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:
ζ(ψ) = CN 2N 1(ψ)· S(ψ) + ψ · ζT (2.7)
where ζ = z/c, ψ = x/c and ζT = ∆ζT E/c. ζT is the dimensionless trailing edge
thickness, which is assumed zero from here on for reasons of simplicity. CN 1
N 2(ψ) and
S(ψ) represent the class and shape function respectively. The class function is defined
as follows:
CN 2N 1(ψ) = (ψ)N 1(1− ψ)N 2 (2.8)
For a NACA type round nose and pointed aft end airfoil the C0.5
1.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 ∑ i=0 ¯ PiBi,p(ψ) (2.9)
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2
Figure 2.6: Bernstein polynomial terms for p = 5
Bi,p= ( p i ) ψi(1− ψ)p−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 ¯Pi) minus one. As an example, for six control points
(p = 5) the following terms are found:
B0,5 = (1− ψ)5 B1,5 = 5(1− ψ)4ψ B2,5 = 10(1− ψ)3ψ2 (2.11) B3,5 = 10(1− ψ)2ψ3 B4,5 = 5(1− ψ)ψ4 B5,5 = ψ5
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 ¯Pi 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)
This equation describes the so-called unit C0.5
1.0 airfoil, which is plotted in Fig. 2.7.
Changes in the coefficient vector ¯Pi 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
2.3. CST METHOD
27
0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4Figure 2.7: Unit airfoil for C0.5
1.0 0 0.2 0.4 0.6 0.8 1 −0.4 −0.2 0 0.2 0.4 0.6
Figure 2.8: Airfoil shapes deformed using ¯
Pi = {1 1 1 1 1 1} (dashed line)
and ¯Pi ={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 ¯Pi vector are equal for
both cases. Looking at the figure, this indeed seems to be the case.
S(0) =√2(RLE/C) (2.13)
S(1) = tan β + ∆zT E
c (2.14)
The advantage of using the CST method in combination with Bernstein
polyno-mials is that every possible coefficient vector ¯Pi leads to a realistic, smooth shape.
However, any change in ¯Pi will lead to a global change in geometry, thereby not
allowing for any local modifications.
2.3.2
CST in three dimensions
With a few simple modifications, the CST method can also be used to model three-dimensional shapes. Kulfan [41] gives the three-three-dimensional equivalent of Eq. (2.7) as:
ζ(ψ, η) = ζN(η) + CN 2N 1(ψ)· 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 η = 2yb . The class function
CN 1
N 2 is defined according to Eq. (2.8).
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
0 0.5 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1
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(ψ, η) = px ∑ i=0 py ∑ j=0 ¯ Pi,j· Bi,px x(ψ)· B y j,py(η) (2.16) with Bi,px x(ψ) = ( px i ) ψi(1− ψ)px−i (2.17) and Bj,py y(η) = ( py j ) ηj(1− η)py−j (2.18)
Comparing Eq. (2.16) to (2.9) shows that the coefficient vector ¯Pichanged to a
coeffi-cient matrix ¯Pi,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 ¯Pi,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
2.3. CST METHOD
29
0 0.2 0.4 0.6 0.8 1 0 0.5 1 −0.5 0 0.5Figure 2.10: Unit airfoil in 3D for C0.5
1.0 0 0.5 1 0 0.5 1 −0.4 −0.2 0 0.2 0.4
In conclusion, the CST method combines the advantages of the analytical, geo-metric, 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 Class-Shape-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
B-splines
3.1.1
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=0 ¯ PiNi,p(u) (3.1)
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 if ti≤ u < ti+1
= 0 otherwise (3.2)
and
Ni,p(u) =
(u− ti)Ni,p−1(u)
ti+p−1− ti
+(ti+p− u)Ni+1,p−1(u)
ti+p− ti+1
(3.3)
where ti are the knot values that relate the parametric variable u∈ [0, n − p + 2⟩ to
the control points ¯Pi. The knot values are defined as follows:
ti= 0 if i < p
ti= i− p + 1 if p≤ i ≤ n
ti= n− p + 2 if i > n
(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} (3.5)
N0,1(u) = 1 for u≡ 0 = 0 otherwise N1,1(u) = 1 for u≡ 0 = 0 otherwise N2,1(u) = 1 for 0≤ u < 1 = 0 otherwise N3,1(u) = 1 for 1≤ u < 2 = 0 otherwise N4,1(u) = 1 for 2≤ u < 3 = 0 otherwise N5,1(u) = 1 for 3≤ u < 4 = 0 otherwise (3.6)
Substituting this into Eq. (3.3) gives:
N0,2(u)≡ 0
N1,2(u) = (1− u)N2,1(u)
N2,2(u) = uN2,1(u) + (2− u)N3,1(u)
N3,2(u) = (u− 1)N3,1(u) + (3− u)N4,1(u)
N4,2(u) = (u− 2)N4,1(u) + (4− u)N5,1(u)
N5,2(u) = (u− 3)N5,1(u) (3.7) and again: N0,3(u) = (1− u)2N2,1(u) N1,3(u) = 1
2u(4− 3u)N2,1(u) + 1 2(2− u) 2N 3,1(u) N2,3(u) = 1 2u 2N 2,1(u) + 1 2(−2u 2+ 6u− 3)N 3,1(u) + 1 2(3− u) 2N 4,1(u) N3,3(u) = 1 2(u− 1) 2 N3,1(u) + 1 2(−2u 2 + 10u− 11)N4,1(u) + 1 2(4− u) 2 N5,1(u) N4,3(u) = 1 2(u− 2) 2 N4,1(u) + 1 2(−3u 2 + 20u− 32)N5,1(u) N5,3(u) = (u− 3)3N5,1(u) (3.8)
3.1. B-SPLINES
35
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1Figure 3.1: B-spline basis functions for p = 3 and n = 5
Substituting Eq. (3.8) into (3.1) results in the following expressions for p(u):
p1(u) = (1− u)2P0+ 1 2u(4− 3u)P1+ 1 2u 2P 2, for 0≤ u < 1 p2(u) = 1 2(2− u) 2P 1+ 1 2(−2u 2+ 6u− 3)P 2+ 1 2(u− 1) 2P 3, for 1≤ u < 2 p3(u) = 1 2(3− u) 2 P2+ 1 2(−2u 2 + 10u− 11)P3+ 1 2(u− 2) 2 P4, for 2≤ u < 3 p4(u) = 1 2(4− u) 2P 3+ 1 2(−3u 2+ 20u− 32)P 4+ (u− 3)2P5, for 3≤ u < 4 (3.9)
These functions are shown graphically in Fig. 3.1. Reparameterizing such that 0≤
u≤ 1 on each interval can be done by replacing u by u for p1(u), u by u + 1 for p2(u),
u by u + 2 for p3(u) and u by u + 3 for p4(u):
p1(u) = (u2− 2u + 1)P0+ (− 2 3u 2+ 2u)P 1+ 1 2u 2P 2, p2(u) = 1 2(u 2− 2u + 1)P 1+ 1 2(−2u 2+ 2u + 1)P 2+ 1 2u 2P 3, p3(u) = 1 2(u 2− 2u + 1)P 2+ 1 2(−2u 2+ 2u + 1)P 3+ 1 2u 2P 4, p4(u) = 1 2(u 2− 2u + 1)P 3+ 1 2(−3u 2+ 2u + 1)P 4+ u2P5, all for 0≤ u ≤ 1 (3.10)
From Eq. (3.10) it is clear that the expressions for p2(u) and p3(u) are similar, while
the expressions for p1(u) and p4(u) are unique. It turns out that after
reparameteriz-ing such that 0≤ u ≤ 1 on each interval, there is one expression for the starting curve,
one expression for the interior curves and one expression for the end curve. This is called the uniform B-spline representation and, just like the Bernstein polynomials, it satisfies the partition of unity. It can be expressed in matrix form as follows. For the first curve:
p1(u) = [ u2 u 1] 1 − 3 2 1 2 −2 2 0 1 0 0 PP01 P2 = UMS0P (3.11)
For the interior curves:
pi(u) = 1 2 [ u2 u 1] −21 −2 12 0 1 1 0 PPi−1i Pi+1 = UMSP (3.12)
And for the end curve:
pn−1(u) = 1 2 [ u2 u 1] −21 −3 22 0 1 1 0 PPnn−2−1 Pn = UMSeP (3.13)
These relationships hold for all B-spline curves with p = 3 and n≥ 4. An example
of a B-spline with p = 3 and n = 4 is shown in Fig. 3.2. For B-splines consisting of
higher order Bezier curves (p > 3) the matrices MS0, MS, and MSe can be derived
in a similar fashion.
3.1.2
B-spline surfaces
As was the case for the Bernstein polynomials, B-splines can also be applied in three dimensions. Equation (3.1) is expanded with an extra set of basis functions. This leads to the following parametric equation for the B-spline surface:
p(u, w) = n ∑ i=0 m ∑ j=0 ¯ PijNi,puu(u)N w j,pw(w) (3.14)
The matrix ¯Pij contains the vertices of a control polyhedron. The basis functions
Ni,puu(u) and Nj,pww(w) are defined in the same way as for the two-dimensional case.
Following the same procedure as in Section 3.1.1 leads to the following expression for a B-spline surface:
p(u, w) = UMSiPM
T SjW
T (3.15)
For p = 3, Eq. (3.15) can be expanded as follows for the interior points of the control polyhedron:
3.1. B-SPLINES
37
0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 3 S S 0 S e 0 2 4 p M M M p p 1 p pFigure 3.2: B-spline with p = 3 and n = 4
pi,j(u, w) = 1 2 [ u2 u 1] −21 −2 12 0 1 1 0
PPi−1,j−1i,j−1 PPi−1,ji,j PPi−1,j+1i,j+1
Pi+1,j−1 Pi+1,j Pi+1,j+1
·1 2 −21 −2 12 0 1 1 0 T w 2 w 1 (3.16)
If i = 1 and/or j = 1 then the matrices MSi and/or MSj change to:
MS0 = 1 − 3 2 1 2 −2 2 0 1 1 0 (3.17)
Similarly, if i = n− 1 and/or j = m − 1 then the matrices MSi and/or MSj change
to MSe = 1 2 −21 −3 22 0 1 1 0 (3.18)
