c
TU Delft, The Netherlands, 2006
AERODYNAMIC SHAPE OPTIMIZATION BY MEANS OF
SEQUENTIAL LINEAR PROGRAMMING TECHNIQUES
Giampietro Carpentieri∗, Michel J.L. van Tooren† and Barry Koren††
Delft University of Technology, Faculty L&R, Kluyverweg 1, 2629 HS Delft, The Netherlands
e-mail: G.Carpentieri@tudelft.nl
Key words: Aerodynamic Shape Optimization, Sequential Linear Programming, Un-structured Meshes, Adjoint equations
Abstract. A Sequential Linear Programming technique, known as the method of centers, is the driver of an aerodynamic shape optimization framework on unstructured meshes. The first order information required by this technique, functional values and their gradi-ents, are computed by a median–dual flow/adjoint solver which is coupled to an analytical shape parameterization. Functional and geometric constraints are easily handled by the algorithm which appears to be very effective in obtaining efficiently near–optimal designs. Shape optimization results are presented for transonic as well as supersonic flows involving appreciable shape deformations.
1 INTRODUCTION
Gradient based shape optimization has received a boost since the introduction of the
adjoint method in CFD.3 Large scale designs involving Finite Volume solvers are
nowa-days possible at a reasonable cost, comparable to that of an additional flow solution for
each design iteration.4 The attractiveness of such benefits in terms of costs have
moti-vated the large effort in the development of Adjoint solvers in the last decade.10,12,14,15
Also, due to their flexibility with respect to the handling of complex geometry, unstruc-tured flow solvers had received particular attention. Among others, median–dual (node centered) formulations were proven to be extremely efficient due to their “edge–based” data structure. Because of this, many median–dual implementations have appeared in the field.6,7,8,9
The present work is part of an effort that is dedicated to the development of an aero-dynamic shape optimization framework for unstructured meshes. The development is horizontal, in the sense that not only CFD aspects are considered but also other aspects such as the shape parameterization and the optimization strategy. The latter are of equal importance for the Shape Optimization practitioner to succeed.
Up to now the capability of the framework is limited to the 2D Euler equations. In
addressed. More in detail, the development and the implementation of the discrete adjoint solver was described. Some results, obtained by the coupling of the flow/adjoint solver with the shape parameterization and a Sequential Quadratic Programming optimization algorithm, were presented. Also, preliminary results from the application of a Sequential Linear Programming technique appeared to be promising.
Therefore, the present paper addresses briefly the flow/adjoint solver and focuses more on the optimization process from an application–oriented point of view. In section 2, after an introduction to the shape optimization problem, the shape parameterization is considered. The Sequential Linear Programming algorithm used to drive the optimization is presented at the end of the section. The algorithm, for the first time applied to a shape optimization problem, is also know as the method of centers. In section 3 the flow/adjoint solver is addressed whereas section 4 deals with the computation of the
gradients. Compared to the last work,1 the gradients are computed with the aid of
Automatic Differentiation. In section 5 the coupling of the whole optimization framework is briefly addressed and, finally, numerical results are presented for shape optimization in transonic as well as supersonic flows.
2 SHAPE OPTIMIZATION
The present work focuses on a shape optimization problem in which the Drag coefficient
(cd) of the airfoil must be reduced and a certain number of constraints must be satisfied.
The constraints are enforced on the Lift coefficient (cl) and on geometric quantities such
as the relative maximum thickness (tcmax) of the section, the nose radius (ru for the upper
side and rl for the lower side, see section 2.1) and the trailing edge angle (θ).
The Lift constraint is considered as an equality, necessary in order not to reduce the Drag coefficient at the expense of the Lift. The relative maximum thickness constraint, considered as inequality, is necessary for avoiding the section to become too thin. The other constraints, on the nose radius and on the trailing edge angle, are included to ensure the feasibility of the shape. For instance, very sharp nose and very thin trailing edge are not desirable from a point of view of manifacturability as well as off–design conditions. Moreover, apart from real life issues, unfeasible shapes can cause the break–down of the
optimization process. In previous work,1,2 without constraints on the nose and on the
trailing edge angle, it was necessary to set bounds (by trial and error) on the design variables in order to avoid unfeasible shapes.
If bounds on the design variables are included, the problem is a bound–constrained optimization which, in general, can be formally stated as:
min f (x) , (1)
gj(x) ≤ 0 (j = 1, ng) ,
hk(x) = 0 (k = 1, nh) ,
xL ≤ x ≤ xU,
indicating a reference value (in the following the subscript R will indicate reference values).
The thickness constraint is defined as g1 = 1 − tcmax/tcmaxR, the upper and lower nose
radius constraints are g2 = 1 − ru/ruR and g3 = 1 − rl/rlR whereas the trailing edge angle
constraint is g4 = 1 − θ/θR. The Lift constraint, which is the only equality constraint
considered here, is defined as h = cl/clR− 1.
The reference values for the constraints have the function of defining the feasible design space. For the objective function, the reference value is only used for the purpose of scaling. In the following, the reference values are considered as initial values multiplied by a constant (the initial values are indicated with the subscript 0).
2.1 Parameterization of the shape
The shape of the airfoil is described using a Chebyshev polynomials parameterization in which the upper and the lower curves are independently approximated. For both curves
the polynomials are made-up of orthogonal basis functions:16
fd(x) = ND X k=0 dkDk(x) , Dk = Tk− Tk+2, (2) Tk(x) = cos(kθ(x)) , θ(x) = cos−1(2√x − 1) , (k ≥ 0 , 0 ≤ x ≤ 1) ,
where ND is the number of basis functions, Dk are the basis functions and dk are the
shape parameters or basis functions coefficients. Those parameters and their number are chosen by means of a minimization procedure which aims at reducing the error between
the curves and the original airfoil coordinates.2 A certain error is achieved in a way that
there are no discrepancies between the flowfield around the original and the approximated airfoil geometry. In the present work the number of basis functions used for each curve
is ND = 11, since this number was found to be effective in approximating many airfoils
accurately.
Therefore, the design variables vector x = [d1, d2, ..., d2ND]
T contains the shape
pa-rameters for the upper and the lower curves (22 design variables in total). Additionally, although not considered in this work, the angle of attack could also be included as de-sign variable. Since the curves are defined analytically, it is possible to calculate the
slope and the radius of curvature at any point (respectively obtained as arctan fd′ and
q
(1 + fd′3)/fd′′). They are used to compute the trailing edge angle θ and the radii of
curvatures ru and rl. Two radii of curvature are computed at the leading edge due to
a lack of continuity of the derivatives at that point. In fact, the derivatives of fd(x) in
x = 0 and x = 1 can only be considered in the limit as x → 0 and x → 1.
2.2 Sequential Linear Programming: the method of centers
A Sequential Linear Programming algorithm has been already applied in a previous
work.1The algorithm5is based upon the linearization of the non–linear optimization
algorithm (e.g. the Simplex method). The solution of the Linear Programming problem is used as a starting point for a new linearization. The process is then iterated until con-vergence to the minimum of the non–linear problem. The algorithm was found effective but slightly inaccurate in terms of the constraints to be satisfied. Such a behavior is not surprising since, as can be seen from Fig. (1.a), the optimum of the Linearized problem is likely to lie outside of the feasible (convex) design space. Therefore, the optimum of the non–linear problem is approached from the unfeasible region.
(a) x0 = linearization point, xc = center of the
circle, xL = optimum of the linearized problem,
xN L = optimum of the non–linear problem
(b) s = move toward the center of the circle, δx = move toward the constraint boundary
Figure 1: Geometric interpretation of the method of centers
In the present work another Sequential Linear Programming algorithm, known as the method of centers, is investigated. The idea behind this algorithm, which is depicted in Fig. (1.a), is to find the largest circle (hypersphere in multi dimension) that fits into the linearized design space. Once the latter has been found, a move in the center of the hypersphere will bring the design in a position which is closer to the optimum. A new linearization may be performed, a new hypersphere may be computed and again the design may be moved to its center. The process may be repeated until convergence to the non–linear optimum.
The interesting feature is that the optimimum is always approached from the feasible region. Moreover, unfeasible designs can be brought into the feasible region simply by computing an hypersphere in which the linearized function is neglected. The largest hypersphere is computed, as can be seen from Fig. (1.b), by means of some geometric
arguments.5 In fact, the hypersphere radius r should be the result of a maximization
under the constraints that df − r ≥ 0 and dj− r ≥ 0 for j = 1, ng with
where dj is obtained considering that, being δx the move to the constraint boundary, one
has gj(x0+ δx) = gj(x0) + ∇gjTδx = 0. If equality constraints have to be included, the
move s towards the center should not violate them, which means hk(x0 + s) = hk(x0) +
∇hkTδx = 0 for k = 1, nh. After some algebra, the maximization problem for the radius
r of the hypersphere may be stated as follows:
max r , (4) ∇fTs + |∇f|r ≤ 0 , ∇gjTs + |∇gj|r ≤ −gj(x0) (j = 1, ng) , ∇hkTs = −hk(x0) (k = 1, nh) , xL− x0 ≤ s ≤ xU− x0, r ≥ 0 .
Equation (4) is a Linear Programming problem for the vector [s, r]T. After the solution
of the problem, the new design located in the center of the hypersphere can be obtained
as xnew = x0+ s. However, as can be seen from Fig. (1.a), one can also take an additional
move in the direction of the steepest descent in order to get closer to the optimum. In
this case xnew = x0+ s − βr∇f/|∇f| where 0 ≤ β ≤ 1. Due to such additional move, also
the line corresponding to the equality constraints in Eq. (4) should be slightly modified
in order to have hk(x0+ s − βr∇f/|∇f|) = 0 for k = 1, nh. Apart from the positive effect
of accelerating the convergence to the optimum, the additional move can increase the risk of violating the constraints.
A negative aspect of this algorithm is the requirement of move limits5for those problems
that are underconstrained. In fact, in such a case the linearization can lead to problems that are unbounded. Unfortunately, this is the case for the problems considered here. Therefore, move limits have to be defined together with a reduction factor. The latter, has to reduce the move limits when the optimum is approached. In practice, when move limits ∆x are introduced the fifth line of Eq. (4) is replaced by −∆x ≤ s ≤ ∆x. One can think the move limits as a box in which the search is limited, an artificial bound for a problem that would be otherwise unbounded.
3 THE FLOW/ADJOINT SOLVER
The Flow solver1 is based upon a median–dual type of discretization7 suitable for
is used to solve the equations: Dt+ g ∂R ∂U !n (Un+1− Un) = −Rn, (5)
where R is the residual vector, U is the vector of conservative variables, Dt is a diagonal
matrix containing the cell volumes divided by their local time steps and ∂R/∂U is ang
approximated Jacobian. The approximations are in terms of the spatial accuracy, which is limited to first order (the reconstruction contribution is neglected), and in terms of the Roe’s fluxes in which the flux Jacobians are frozen when performing the linearization. The solution of the linear system arising at each time step is accomplished using a SSOR scheme implemented using the edge–based data structure of the solver. The scheme is fast and efficient and it has the additional feature of being optionally matrix–free. In spite of the time penalty, due to the recomputation of the Jacobian elements, with the matrix–free options the memory requirements are that of an explicit flow solver.
The discrete adjoint is the preferred approach for computing the adjoint variables, which are the sensitivity of the functionals with respect to the state (the residual vector
R). An exact discrete adjoint has been derived and implemented in a matrix-free fashion.1
Exact means that the differentiation of all the components in the residual vector is made without any approximation. For a MUSCL-like implementation, the development of the discrete adjoint means the differentiation of the fluxes as well as the reconstruction oper-ator. It is useful to consider a vector H containing the numerical fluxes evaluated on each cell interface (edge mid–point). Such a vector is dependent on the reconstructed left and
right states on the interface, indicated respectively with UL and UR. The length of the
three vectors is equal to the numbers of interfaces, E, which in this case correspond to the number of edges in the mesh. The dependency of the residual vector on the conservative variables vector, both of length N equal to the number of nodes, can be expressed as
R(U) = R{H[UL(U), UR(U)]}. Applying the chain rule and transposition one has :
∂R ∂U T = ∂UL ∂U T ∂H ∂UL T + ∂UR ∂U T ∂H ∂UR T! ∂R ∂H T , (6)
where ∂R/∂H is a dummy matrix whose only non–zero elements have values 1 and −1,
∂H /∂UL and ∂H /∂UR are the differentiation of the numerical fluxes and ∂UL/∂U
and ∂UR/∂U are the differentiation of the reconstruction operator. In practice none of
these matrices is formed, but more likely a matrix–vector product involving ∂R/∂UT is
assembled directly on the edges in much the same way as for the residual vector.
The solution process for the adjoint equations is identical to that of the flow equations. The same first order backward Euler scheme is used, with the only difference that the approximated Jacobian is in this case transposed:
It can be argued that this solution method is not necessary for a linear system of equations. However, due to the poor diagonal dominance of the exact Jacobian, it can be an expensive and not always successful approach to use only a linear system solver. On the contrary, the use of the pseudo–time stepping (diagonal dominance is increased due to the approximated
Jacobian and the time step contribution) gives a robust solver for the adjoint equations.11
An additional feature of the pseudo-time stepping and the SSOR scheme adopted for the adjoint equations, is the possibility to solve for multiple right–hand sides. This means that multiple adjoints can be computed in only one-shot. This turns out to be very useful in constrained shape optimization since more than one functional is usually required. For the simultaneous solution of 3 functionals, up to 50% time saving has been noticed.
4 THE COMPUTATION OF THE GRADIENTS
Once the adjoint variables have been computed, the gradient of a functional J with
respect to a shape-parameter xi is obtained as:
dJ dxi = ∂J ∂xi − Λ T∂R ∂xi , (8)
where ∂J /∂xi and ∂R/∂xi can be computed, among other methods,10,12,13 by finite
dif-ference or by means of Automatic Differentiation. The finite difdif-ference approach can be inefficient for a large number of parameters. Moreover, it needs a proper step size in order to reduce the cancellation and the truncation errors. On the contrary, Automatic differentiation relies on basic linearization rules and the exact differentiation of the code
can be derived and executed efficiently.13 More challenging is the derivation of the
ad-joint code with the reverse mode where storage and efficiency issues should be carefully
considered.12,13,14 For this reason, the adjoint code presented before has been hand–coded.
In this work the forward vectorial mode of Tapenade17 has been adopted, which allows
to compute ∂R/∂xi for the complete design vector x in one shot. Time saving is achieved
compared to sequential evaluations since many quantities, having a value which is the same for all parameters, are computed only once.
In order to understand how the differentiated code is derived using Tapenade, the
dependency of the residual vector R on a shape parameter xishould be described first. It is
convenient to start from the mesh coordinates X, which are functions X = X[∆XB(xi)] of
the shape parameter xi. In fact, a change in xi gives a displacement of the boundary ∆XB
through the shape parameterization. In turn, the displacement implies a new repositioning of all the mesh coordinates X. The latter can be achieved using a mesh deforming technique that compute the displacement ∆X to be added to the original mesh points (
Xnew = X + ∆X ). In this work a spring analogy is chosen, where the displacement of
each point ∆Xi is iteratively solved with Jacobi iterations:
where Ni is the number of nearest neighbors of the node i and kij is the stiffness associated
with the edge ij (its numerical value is chosen to be the inverse of the edge length).
Therefore, the boundary displacement ∆XB is iteratively propagated in the whole domain.
The residual vector has a dependency R = R[X, N(X), ∇U(X), Φ(X)] on the mesh points, where N are the grid metrics, integrated normal vectors in the grid and on the boundary, and ∇U and Φ are the least squares gradient and the limiter vector respectively. Applying the chain rule, the complete derivative read
∂R ∂xi = ∂R ∂X + ∂R ∂N ∂N ∂X + ∂R ∂∇U ∂∇U ∂X + ∂R ∂Φ ∂Φ ∂X ! ∂X ∂∆XB ∂∆XB ∂xi (10) where, in practice, each of the derivative corresponds to a differentiated sub–routine.
∂∆XB/∂xi is the differentiation of the shape parameterization routine with respect to the
parameters, ∂X/∂∆XB is the differentiation of the spring analogy routine with respect
to the boundary displacement. ∂∇U/∂X and ∂Φ/∂X are the reconstruction routines differentiated with respect to the mesh points. ∂R/∂X, ∂R/∂N, ∂R/∂∇U and ∂R/∂Φ represent the residual routine which is differentiated with respect to mesh points, metrics, gradients and limiters.
A similar procedure has been implemented for the partial derivative of the functional
∂J /∂xi . In order to check the accuracy of the differentiated code, comparisons with finite
difference has been made for all the differentiated routines. As in the case of the adjoint code, the differentiation of the code in forward mode is exact since no approximations have been made. It is interesting to note that the differentiation of the mesh displacement routine involves an iterative process. The routine has been differentiated properly by Tapenade, as noticed after the checking, without any particular problem.
5 NUMERICAL RESULTS
In order to perform shape optimization, the elements presented above are coupled together in a relatively simpler way. Once the initial geometry, the objective function and the constraints have been defined:
1. the flow solver module is used to compute the flow field U;
2. the adjoint solver module is used to compute the adjoint variables ΛL and ΛD for
the lift and the drag coefficient respectively;
3. a post–processing module receive as input the parameters, the flow and the adjoint
variables. It compute the functionals (cl and cd) and the geometric quantities (
tcmax, ru, rl and θ) together with their gradient ( ∇cl, ∇cd, ∇tcmax, ∇ru and ∇rl )
with respect to the shape parameters x;
5. the mesh is deformed according to the new values of the shape parameters and the process is repeated starting from step 1.
Each element is implemented as a Fortran 90 module at the exception of the Sequential
Linear Programming module. The latter is a script implemented in the Matlab18language,
which use Linear Programming libraries included in the Optimization Toolbox.
The framework has been employed in some shape optimization test cases. A mesh of approximatively 8000 nodes, as shown in Fig. (2.e), has been used for all the tests since different airfoils, starting from an initial airfoil, are obtained trough the application of the shape parameterization and the mesh deforming routines.
For each test case, the initial flow and adjoint field have been computed with a level of accuracy equal to 8 orders of magnitude reductions of the residual norm. During the optimization, both the flow and the adjoint, are converged to reach the same accuracy as the initial computation.
The initial design vector x at the beginning of the optimization is scaled in a way that
x0 = [1, 1, ..., 1]T. The move limits ∆x assume the same value for each component of the
design vector, they are set as ∆x = b x0. Therefore, in the following b will be referred as
the move limit.
5.1 RAE2822
The RAE2822 airfoil is optimized at M∞= 0.73 and α = 2 deg angle of attack. As
can be seen on Fig. (2.c), a shock is present on the upper surface of the airfoil. The reference maximum thickness is equal to the initial maximum thickness of the airfoil
(tcmaxR = tcmax0) whereas the reference values for the nose curvature and the trailing
edge angle are set respectively to 70% and 90% of their initial values (ruR = 0.7 ru0,
rlR = 0.7 rl0and θR = 0.9 θ0). This choice of the reference values means that the thickness
can’t decrease, the trailing edge angle can only decrease 10% of its initial value and the two radii of curvature can only decrease 30% of their initial values. The lift is imposed to remain constant during the optimization, which means that for the equality lift constraint
the reference value is the initial one (clR = cl0). The objective function, which is the drag
coefficient, is scaled by its initial value (cdR= cd0).
At the end of the optimization, as can be seen on Fig. (2.d), a shock–free airfoil is obtained. The constraints, although their trend is not shown, have been satisfied. The objective function, depicted in Fig. (2.a), has a non–smooth trend for a value of the move limit b = 0.1. Probably such value is too large and the design has the tendency to oscillate. In fact, with a move limit b = 0.05 the objective seems to behave smoothly, in spite of the reduced convergence rate.
in proximity of the optimum the oscillations, although mitigated by the reduction factor for the move limit, are unlikely to be suppressed completely.
5.2 NACA0012 transonic
The NACA0012 is optimized at M∞ = 0.73 and α = 2 deg angle of attack. For such
conditions the shape of the NACA0012 airfoil is clearly unsuitable. In fact, Fig. (3.c) shows a strong shock appearing on the upper side of the airfoil, almost at half of the chord. The reference values used for the geometric constraints are different from the
previous case ( ruR = 0.9 ru0, rlR = 0.9 rl0 and θR = 0.9 θ0 ) except for the thickness
which is still not allowed to decrease (tcmaxR = tcmax0).
The move limits have been chosen to be b = 0.3 until the objective had become flat. After, they have been reduced as before. As can be seen on Fig. (3.a), the objective function smoothly reduced 90% with respect to the initial value. The constraints have been satisfied. The final shock–free airfoil is shown in Fig. (3.a). The flow conditions on the bottom side are almost sonic at the quarter of the chord.
5.3 NACA64A410 transonic
The NACA64A410 is optimized at M∞= 0.75 and α = 0 deg angle of attack. A strong
shock, see Fig. (4.c), appears at three quarters of chord on the upper side of the airfoil. The same reference values used for the NACA0012 have been used in this case.
The move limits, with the same values as in the NACA0012 case, are probably large for this case since the objective, see Fig. (4.a), does not seems to behave very smoothly. At the end of the optimization, as in the last cases, a shock–free airfoil has been obtained. However, as can be seen on Fig. (4.d), at this time the optimized shape appears to be rather complex if compared to the initial NACA64A410 shape.
5.4 NACA0012 supersonic
The NACA0012 has also been optimized at supersonic flow conditions, M∞= 1.5 and
α = 2 deg angle of attack. In the supersonic case the shape parameterization and the geometric constraints play a crucial role. It is well known that thin airfoil with sharp noses are likely to produce less wave drag. This is because if the nose is sharp the bow shock tends to be attached and oblique. On the contrary, a rounded nose produce a detached shock which is almost normal just ahead of the nose. The large wave drag which is noticed in this situation is due to the presence of the almost normal shock.
Therefore, an optimization in supersonic regime is expected to produce a sharp airfoil. The situation of a sharp nose cannot be handled properly with the shape parameterization used in this work. Moreover, real life airfoils do not have sharp noses but more likely noses with small radii of curvature. For this reason, a supersonic optimization is a good test case to see if the geometric constraints introduced in this work are effective.
The NACA0012 airfoil has a 12% maximum thickness and a large radius of curvature at the nose. A value of the maximum thickness more suitable for supersonic flow could be 6% and therefore the reference thickness value is imposed to be half of the initial one
(tcmaxR= 0.6 tcmax0). Also, all the other reference geometric quantities has been imposed
to be 10% of their initial values ( ruR= 0.1 ru0, rlR = 0.1 rl0 and θR = 0.1 θ0 ). The latter
means that, for instance, the nose radius curvature can decrease up to 90% of the initial
value. Also in this case, the Lift coefficient must be kept constant (clR = cl0) and the
Drag coefficient,objective function, is scaled by the initial value (cdR= cd0). Move limits
b = 0.3, as in the NACA0012 transonic case, have been used.
The objective function, see Fig. (5.a), smoothly reduced of 70%. Such large reduction in drag is due to the repositioning of the shock, which is almost attached to the nose at the end of the optimization. A magnification of the nose region, see Fig. (5.e), reveals that the shock is at less than 1% of the chord distance from the nose. The latter, appears to be rounded and smooth. The curvature of the nose is now 10% of the initial curvature since both the upper and lower curvature constraints are critical. The same is for the thickness constraint. The three constraints are depicted in Fig. (5.b), which shows how
they approach criticality (g1, g2, g3 → 0). All the other constraints, although not shown,
are satisfied. It is remarkable that also in this case, in spite of the challenging flow conditions and the large deformations involved, the whole framework proved to be effective and robust. Also, it is interesting to see that the final shape resembles that of lenticular airfoils, typically used in the supersonic regime.
6 CONCLUSIONS
Different shape optimization test cases have been successfully addressed by the frame-work presented in this frame-work. The introduction of additional geometric constraints, nose curvature and trailing edge angle, has removed the ambiguity of setting the bounds on the shape parameters in order to avoid unfeasible shapes. In fact, the geometric constraints guarantee to have a smooth and feasible shape according to the prescribed reference val-ues.
ACKNOWLEDGMENTS
This research is supported by the Dutch Technology Foundation STW, applied science division of NWO and the technology program of the Dutch Ministry of Economic Affairs. REFERENCES
[1] G. Carpentieri, M.J.L. van Tooren and B. Koren, “Improving the Efficiency of Aerodynamic Shape Optimization on Unstructured Meshes.” AIAA Paper 2006-298, Reno, NV, Jan. 2006.
[2] G. Carpentieri, M. J. L. van Tooren, M. Kelly and R. Cooper,“Airfoil Optimization Using an Analytical Shape Parameterization” CEIAT Paper 2005–0073, Belfast, UK, Aug. 2005.
[3] A. Jameson, “Aerodynamic Design via Control Theory.” Journal of Scientific Com-puting, 3, 233–260, 1988
[4] A. Jameson, “Optimum Transonic Wing Design Using Control Theory.” Symposium Transsonicum IV, DLR Gottingen, Germany, Sep. 2002.
[5] G.N. Vanderplaats, Numerical Optimization Techniques for Engineering Design, 3th.
Edition, Vanderplaats Research & Development, 2001.
[6] J. A. Desideri and A. Dervieux, “Compressible Flow Solvers using Unstructured Grids.” Lecture series 1988-05, Von Karman Institute for Fluid Dynamics 1988. [7] T.J. Barth, “Aspects of unstructured grids and finite–volume solvers for the Euler
and Navier–Stokes equations.” Lecture series 1991-06, Von Karman Institute for Fluid Dynamics 1991.
[8] V. Selmin and L. Formaggia, “Unified Construction of Finite Element and Finite Volume Discretizations for Compressible Flows.” International Journal for Numerical Methods in Engineering, 39, 1–32, 1996.
[9] V. Venkatakrishnan and D. J. Mavriplis, “Implicit Solvers for Unstructured Meshes.” Journal of Computational Physics, 105, 83–91, 1993.
[10] E. J. Nielsen, “Aerodynamic Design Sensitivities on an Unstructured Mesh Using the Navier–Stokes Equations and a Discrete Adjoint Formulation.” PhD Thesis, Virginia Polytechnic Institute and State University, Blacksburg, Virginia, 1998.
[12] M. B. Giles, M. C. Duta, J.-D. M¨uller and N. A. Pierce, “Algorithm Developments for Discrete Adjoint Methods.” AIAA Journal, 41, 198–205, 2003.
[13] J.-D. M¨uller and P. Cusdin, “On the Performance of Discrete Adjoint CFD Codes
using Automatic Differentiation.” Int. J. Numer. Meth Fluids, 47, 939–945, 2003. [14] B. Mohammadi, “ A New Optimal Shape Design Procedure for Inviscid and Viscous
Turbulent Flows.” International Journal for Numerical Methods in Fluids, 25, 183-203, 1997.
[15] D. J. Mavriplis, “ A Discrete Adjoint–Based Approach for Optimization Problems on Three–Dimensional Unstructured Meshes .” AIAA Paper 2006-50, Reno, NV, Jan. 2005.
[16] A. Verhoff, D. Stooksberry and A. B. Cain, “An Efficient Approach to Optimal Aerodynamic Design Part 1: Analytic Geometry and Aerodynamic Sensitivities.” AIAA Paper 93–0099, Reno, NV, Jan 1993.
[17] http://tapenade.inria.fr:8080/tapenade/index.jsp Software TAPENADE c INRIA
2002, version 2.0
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 OBJ NUMBER OF ITERATIONS b = 0.05 b = 0.1
(a) Objective function (cd/cdR).
0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 RAE2822 OPTIMIZED
(b) Initial and final Mach number distribution.
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
(c) Pressure contours for the RAE2822.
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
(d) Pressure contours of the optimized airfoil.
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
(e) 8K nodes unstructured mesh around the RAE2822.
Figure 2: Optimization of the RAE2822 airfoil at M∞= 0.73 and α = 2 deg angle of attack. At the end
of the optimization: the equality lift constraint is h = −5 × 10−5
, the thickness constraint is g1= −0.011,
the upper radius of curvature and the lower radius of curvature constraints are respectively g2= −0.357
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 OBJ NUMBER OF ITERATIONS
(a) Objective function.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 NACA0012 OPTIMIZED
(b) Initial and final Mach number distribution.
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
(c) Pressure contours of the NACA0012 airfoil.
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
(d) Pressure contours of the optimized airfoil. Figure 3: Optimization of the NACA0012 airfoil at initial condition M∞= 0.75 and α = 2 deg angle of
attack. The final constraints values are h = −6.35 × 10−6
, g1= −1.15 × 10 −3
, g2= −0.011, g3= −3.15
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 OBG NUMBER OF ITERATIONS
(a) Objective function.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 NACA64A410 OPTIMIZED
(b) Initial and final Mach number distribution.
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
(c) Pressure contours of the NACA64A410.
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
(d) Pressure contours of the optimized airfoil. Figure 4: Optimization of the NACA64A410 airfoil at initial condition M∞= 0.75 and α = 0 deg angle
of attack. The final constraints values are h = 6.5 × 10−5
, g1= −2.3 × 10 −3
, g2= −4.415, g3= −0.1684
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 OBJ NUMBER OF ITERATIONS
(a) Objective function.
-1 -0.8 -0.6 -0.4 -0.2 0 0 5 10 15 20 G NUMBER OF ITERATIONS G1 G2 G3
(b) Thickness and nose radius constraints.
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
(c) Pressure contours of the NACA0012.
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1
(d) Pressure contours of the optimized airfoil.
-0.01 -0.005 0 0.005 0.01 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02 0.025
(e) Magnification of the optimized airfoil nose.
Figure 5: Optimization of the NACA0012 airfoil at initial condition M∞= 1.5 and α = 2 deg angle of
attack. The final constraints values are h = 8.8 ×10−5
, g1= −6.5×10−4, g2= −5×10−3, g3= −6×10−3
