Table Of ContentShape Optimization of Low Speed Airfoils
using MATLAB and Automatic Differentiation
Christian Wauquiez
Stockholm 2000
Licentiate’s Thesis
Royal Institute of Technology
Department of Numerical Analysis and Computing Science
Abstract
Thegoaloftheprojectistodevelopan innovative tooltoperform shapeoptimizationoflow
speed airfoils. This tool is written in Matlab, and is constructed by coupling the Matlab
Optimization Toolbox with a parametrised numerical aerodynamic solver.
The airfoil shape is expressed analytically as a function of some design parameters. The
NACA 4 digits library is used with design parameters that control the camber and the
thickness of the airfoil.
The solver has to provide fast and robust computation of the lift, pitching moment and drag
of an airfoil placed in a low-speed viscous flow. A one-way coupled inviscid - boundary
layer model is used.
The inviscid flow is computed with a linear vortex panel method, which provides the lift and
moment coefficients. The boundary layer is computed using an integral formulation : the
laminar part of the flow is computed with a two-equation formulation, and the turbulent part
is solved with Head’s model. An e9-type amplification formulation is used to locate the
transition area. Finally, the drag coefficient is computed using the Squire-Young formula.
In order to be used in optimization, the solver must provide derivatives of the objective
function and limiting constraints with the solution for each set of parameters.
These derivatives are computed by automatic differentiation, a technique for augmenting
computer programs with the computation of derivatives based on the chain rule of differ-
ential calculus. The recent Matlab automatic differentiation toolbox ADMAT is used.
Finallyasanapplication,sampleoptimizationproblemsaresolvedusingtheMatlabOptimi-
zation Toolbox, and the resulting optimal airfoils are analysed.
ISBN 91-7170-520-1l TRITA-NA-0004l ISSN 0348-2952l ISRN KTH/NA/R--00/04--SE
1
Contents
Introduction
1-TheAerodynamicsSolver..................................................7
1.1- Introduction - Overview of the Model
1.2- Airfoil and Flow Parameters
1.3- Inviscid Flow Model
1.4- Boundary Layer Model
2-AutomaticDifferentiation.................................................39
2.1- Introduction
2.2- Method Fundamentals
2.3- Computer Implementation
2.4- ADMAT, Automatic Differentiation Toolbox for Matlab
2.5- Application of ADMAT to the Aerodynamic Solver
3-AirfoilShapeOptimization................................................51
3.1- Definition of the Optimization Problems
3.2- Solving the Optimization Problems
Conclusion................................................................61
References
Acknowledgements
I would like to thank my supervisor, Associate Professor Jesper Oppelstrup, for his interest
and support throughout this work, and Professor Arthur Rizzi for providing information and
for many helpful discussions.
I would also like to thank, for advice and technical support, Associate Professor Ilan Kroo
from Stanford University and Desktop Aeronautics, and Doctor Arun Verma from the
Cornell University Theory Center.
Financial support from NUTEK through the Parallel Scientific Computing Institute, KTH,
and TFR through the Center for Computational Mathematics and Mechanics, KTH, is
gratefully acknowledged.
3
Introduction
The performance of an airfoil can be characterized by three quantities : the lift, moment and
dragcoefficients:Cl,CmandCdrespectively.Theyrepresenttheaerodynamicloadsapplied
totheairfoil.Theactualloadsareproportionaltothecoefficientstimesthesquareoftheflow
velocity.
• Cl corresponds to the force acting on the airfoil in the direction orthogonal to the flow,
which allows an aircraft to fly by compensating its weight.
• Cmcorrespondstothemomentoftheaerodynamicforcewithrespecttothequarterofthe
airfoil chord length. For the equilibrium of an aircraft, the pitching moment of the main
winghastobecompensatedbythemomentofanegative-lifttail.Cmshouldthereforenot
be too large.
• Finally,Cdcorrespondstothecomponentoftheforceintheflowdirection,whichhinders
aerodynamic performances and causes fuel consumption.
Thepresentworkfocusesonshapeoptimizationofairfoilsinlowspeedviscousflows,based
on the analysis ofCl,Cm and Cd.
The airfoils are chosen from the NACA 4 digits library [1], in which the shape is expressed
analytically as a function of three parameters. The library is presented in section 1.2 of this
report.
The formulation used to compute the aerodynamic coefficients is an inviscid - boundary
layer model. The advantage of this kind of approach is that it provides a fast computation of
the flow solution, the disadvantage being that cases with massive flow separation are
impossible to handle. Famous codes based on this formulation include for 2D cases Desktop
Aeronautics’ Panda [9], and Mark Drela’s ISES [10] and Xfoil [12], and for 3D cases Brian
Maskew’s VSAERO [3]. The solver is presented in section 1.
In order to perform optimization, the solver has to provide the derivatives of the defined
objective function and limiting constraints with the solution for a given set of parameters.
The most common method used to compute these derivatives is finite differences. In the
present work, a new approach, called automatic differentiation (AD), is used. Is is a chain-
rule-based technique to compute the derivative of functions defined by computer programs
withrespecttotheirinputvariables,andhasbeeninvestigatedsince1960.In1998,ageneral
AD toolbox for programs written in Matlab, called ADMAT [19], has been developped at
Cornell University by Arun Verma and Thomas F. Coleman. The present work uses this
toolbox. Automatic differentiation is presented in section 2.
Finally,airfoilshapeoptimizationisperformedusingtheMatlabOptimizationToolbox[20],
which uses a Sequential Quadratic Programming algorithm for non linear constrained
problems. The formulation of sample optimization problems and their resolution are
presented in section 3.
5
1- The Aerodynamic Solver
The solver used to provide the objective function as well as the limiting constraints of the
optimization is presented. The solution is computed using an inviscid-boundary layer model,
and consists in the lift, moment, and drag coefficients : Cl, Cm and Cd respectively.
1.1- Introduction - Overview of the Model
The general Navier-Stokes equations are very powerful since they give a complete
description of all possible flow situations. However it is very time-consuming to obtain a
numerical solution using them. In our particular case, the incompressible turbulent flow past
an airfoil, the viscous effects are important only in a small region near the profile. In this
region, the Navier-Stokes equations can be approximated by the so-called boundary layer
equations.Outside,viscouseffectscanbeneglected,andonecanuseaninviscidflowmodel.
Navier-Stokes Euler
V
V
Boundary Layer
The inviscid flow model
Theinviscidpartoftheflowcanbesolvedinseveralways.Afinitedifferencediscretization
ofthesteadyEulerequationonagridaroundtheairfoilcanbeused,asinISES[10],orapanel
method, as in Xfoil [12] and Panda [9]. One advantage of the first option is that by applying
aperiodicboundaryconditiontotheouterboundary,cascadeflowscanbesimulated.Thisis
done in MISES [13], an extension of ISES.
In the present work, a panel method is used. A variety of such methods exists, they differ in
the choice of the singularity used to represent the velocity potential on the airfoil (sources,
doublets or vortices), and by the choice of the Kutta condition, an extra condition that one
must add to the final system of equations in order to obtain a unique solution. Details about
different panel methods can be found in [2]. The present work uses a linear vortex distri-
bution, which gives a good solution accuracy, even with only a few panels. The Kutta
condition chosen is especially well suited for linearly varying singularity distributions.
The inviscid flow solver provides the tangential velocity distribution on the airfoil’s surface
(Ue). The pressure distribution is then computed from the velocity field using the Bernoulli
equation. The lift and moment coefficients, as well as the pressure drag, are calculated by
integrating the pressure over the body surface.
7
1- The Aerodynamic Solver
The boundary layer model
The boundary layer formulation consists of a model for the laminar part of the flow, a
transition criterion, and a model for the turbulent part of the flow. The flow models rely on
one or two differential equations derived from the integration of the Falker-Skan equations
accross the boundary layer thickness, and on additional semi-empirical equations which
close the sytem.
Xfoil and ISES use the same viscous formulation : a set of two differential equations
(integral momentum and kinetic energy shape parameter equations) for both laminar and
turbulent flows, and different closure relations depending on the flow regime. Panda on the
other hand uses simpler formulations : Thwaites’ one equation method [4] for the laminar
part, and Head’s two equations method [5] for the turbulent part. The present work initially
usedsimplemodels,thusThwaites’andHead’s.ButsinceThwaites’modelcannotrepresent
separated flow, Xfoil’s laminar boundary layer model, which can describe thin separated
regions, has been implemented instead.
Concerning the transition criterion, Panda uses Michel’s criterion, and Xfoil uses a more
advanced e9-type formulation. The second method has been chosen after testing.
Details about the reasons of these choices are found in sections 1.4.1 and 1.4.2.
Theresultsprovidedbytheboundarylayersolverincludethedisplacementthickness d * ,the
momentum thickness q , the shape factor H, and the skin friction coefficient Cf. These
quantitiesareusedtocomputethedragcoefficient.Differentmethodscanbeuseddepending
on the coupling between the inviscid flow and the boundary layer.
The coupling
The coupling between the two models is as follows :
• The effect of the boundary layer is that it modifies the shape of the airfoil as seen by the
externalflow.Thisgivestheinviscidflowazeronormalvelocityconditiononaboundary
obtained by adding the boundary layer displacement thickness d * to the airfoil’s surface.
• The boundary layer equations depend on the external tangential velocity distribution.
Ue
Inviscid Flow with
zero normal velocity condition Boundary Layer
on Airfoil +d * d *
Two different approaches are possible for the coupling of inviscid - boundary layer flows.
Two-way coupled computation
In this case, the modification of the boundary where the zero normal velocity condition has
to be met due to the boundary layer thickness is taken into account.
The solution begins with the inviscid flow problem, which produces the velocity field. This
dataisthenfedintotheboundarylayermodelwhichresultsthelocalwallfrictioncoefficient
and the displacement thickness. Then a second iteration is performed, now with modified
8
1.2- Airfoil Shape Parameters
surface geometry. This modification can be obtained by displacing the body panels
according to the local displacement thickness, and the procedure is reiterated until a
converged solution is obtained.
Another way to account for the displacement effects is to modify the boundary condition
instead of the geometry. In this case the normal flow is made non-zero to account for the
effect of d * . This formulation, known as the transpiration velocity concept, states that :
V(cid:215) n = --d----(U d * ) on the airfoil’s surface.
dx e
Using a two-ways coupled model, the total drag is found by adding the friction drag,
obtained by integrating the skin friction coefficient Cf, to the pressure drag, obtained by
integrating the inviscid pressure distribution.
A study of such iterative methods can be found in [14], where it appears that convergence is
not easy to obtain. A more recent approach is the one used in ISES and Xfoil. The
transpiration velocity concept is used, and the entire nonlinear equation set is solved
simultaneously as a fully coupled system by a Newton-Raphson method. This method
provides more robust results than the iterative approach.
One-way coupled computation
In this case the effect of the boundary layer thickness is neglected. One single iteration is
performed : the external tangential velocity is computed by the inviscid model with the
conditionV.n=0 on the airfoil’s surface, and then fed into the boundary layer model.
The drag coefficient is obtained using the Squire-Young formula [8], which in effect
computes the momentum deficit. It is a function of some of the boundary layer results
(momentum thickness and shape factor) at the trailing edge.
Compared to coupled computations, the accuracy of the lift is very good, but not as good for
the drag. However the sensitivity of the results with respect to airfoil shape or flow
parameters is well reproduced. This method is used in Panda, as well as in the present work.
1.2- Airfoil Shape Parameters
The shape of the airfoil can be chosen in the famous Naca 4 digits library [1]. This simple
library is interesting because the shape is expressed analytically as a function of three
parameters, which control the maximum camber, maximum camber location, and maximum
thickness of the airfoil.
y 2412
Maximum Thickness
8412
Maximum Camber
4212
x 4812
4406
Maximum Camber Location 4420
9
Description:The recent Matlab automatic differentiation toolbox ADMAT is used. report. The formulation used to compute the aerodynamic coefficients is an
