Numerical sensitivity analysis for aerodynamic optimization A survey of approaches

更新时间:2023-05-06 03:02:01 阅读量: 实用文档 文档下载

说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。

Review

Numerical sensitivity analysis for aerodynamic optimization:A survey of approaches

Jacques E.V.Peter a ,Richard P.Dwight b,*

a ONERA BP 72–29av.de la Division Leclerc,92322Chatillon Cedex,France b

Aerodynamics Group,TU Delft,P.O.Box 5058,2600GB Delft,The Netherlands

a r t i c l e i n f o Article history:

Received 19November 2007

Received in revised form 21July 2009Accepted 24September 2009Available online 6October 2009Keywords:

Sensitivity analysis

Linearized Navier–Stokes Flow solver linearization Adjoint method

Duality preserving iteration Krylov stabilization

Aerodynamic optimization

a b s t r a c t

The calculation of the derivatives of output quantities of aerodynamic ?ow codes,commonly known as numerical sensitivity analysis,has recently become of increased importance for a variety of applications in ?ow analysis,but the original motivation came from the ?eld of aerodynamic shape optimization.There the large numbers of design variables needed to parameterize surfaces in 3D necessitates the use of gradient-based optimization algorithms,and hence ef?cient and accurate evaluation of gradients.In this context over the last 20years a variety of approaches have been developed to supply these gradi-ents,raising particular challenges that have required novel algorithms.In this paper,we examine the his-torical development of these approaches,describe in some detail the theoretical background of each major method and the associated numerical techniques required to make them practical in an engineer-ing setting.We give examples from our own experience and describe what we consider to be the state-of-the-art in these methods,including their application to optimization of complex 3D aircraft con?gurations.

ó2009Elsevier Ltd.All rights reserved.

Contents 1.Introduction (374)

2.

Sensitivity evaluation methods ..........................................................................................3752.1.Finite differences................................................................................................3752.2.The discrete direct method........................................................................................3762.3.The discrete adjoint method.......................................................................................3762.4.Evaluation of discrete Jacobians....................................................................................3762.5.Derivation of the continuous adjoint equations .......................................................................3772.6.Spatial discretization of the continuous adjoint equations...............................................................3782.7.Discrete versus continuous adjoint .................................................................................3793.

Accuracy of sensitivity computations .....................................................................................3793.1.Fully linearized or frozen turbulence modeling .......................................................................3803.2.Other approximations............................................................................................3803.3.Impact on,and influence of gradient accuracy ........................................................................3814.

Solution strategies for the primal and adjoint equations .....................................................................3814.1.Duality preserving fixed-point iterations.............................................................................3824.2.Krylov stabilization ..............................................................................................3835.

Generalizations of gradient evaluation....................................................................................3845.1.Adjoint mesh deformation ........................................................................................3845.2.Higher-order derivatives using discrete methods ......................................................................

3845.2.1.DD.DD .................................................................................................3855.2.2.DD.AV .................................................................................................3855.2.3.AV.DD .................................................................................................3855.2.4.AV.AV .................................................................................................3855.2.5.Applications ............................................................................................

385

0045-7930/$-see front matter ó2009Elsevier Ltd.All rights reserved.doi:10.1016/0f9759e85ef7ba0d4a733b0ep?uid.2009.09.013

*Corresponding author.

E-mail addresses:jacques.peter@onera.fr (J.E.V.Peter),r.p.dwight@tudelft.nl (R.P.Dwight).Computers &Fluids 39(2010)

373–391

Contents lists available at ScienceDirect

Computers &Fluids

j ou r na l h om e pa ge :w w w.e lse vi e r.c om /lo c at e /c om p?u

id

5.3.Extension to aeroelasticity (386)

6.Examples of sensitivity applications (386)

6.1.Navier–Stokes optimization of a wing–fuselage configuration (387)

6.2.Multi-disciplinary optimization of an engine pylon (388)

7.Conclusion (389)

Acknowledgments (389)

References (389)

1.Introduction

In1987,Sobieszczanski-Sobieski,a specialist in multi-disciplin-ary optimization working at NASA,made a plea to the CFD commu-nity,to extend their codes from simple aerodynamic analysis to sensitivity analysis[1],by which he meant evaluation of the deriv-atives of aerodynamic quantities(typically depending on both geometry and?ow)with respect to some parameterization of the geometry.Such an evaluation is demanding because the?ow itself depends on the geometry over a system of partial differential equa-tions.He pointed out the lead of computational solid mechanics over CFD in this respect,and the signi?cant bene?ts sensitivity analysis can provide to shape design through the use of gradient-based optimization.His arguments were widely heard,and in the following years much effort was invested in ef?cient and accurate evaluation of sensitivities.

However,despite considerable progress over20years,sensitivity analysis has only recently enjoyed widespread use in engineering practice.There are perhaps two principal causes of this:(i)the ques-tionable suitability of gradient-based optimization methods for many engineering problems,and(ii)the lack of availability of cheap, accurate gradients.A close third might be the extra effort involved in setting up gradient-based optimizations,due to the need for smooth mesh deformation,preconditioning of design variables,and the lack of robustness to the failure of any step of the process.This last is not true of,e.g.,genetic algorithms,for which a?ow code failure is sim-ply a member of the population that is not evaluated.

Issue(i)is an unavoidable consequence of the locality of gradient-based methods:in a design space with many local optima they will tend to?nd the nearest(with respect to the starting point),not the best,local optimum.Further there is evidence to suggest that con?g-urations such as multi-element high-lift devices have just such highly oscillatory design spaces[2].This problem may be tackled with hybrid optimization algorithms that combine non-determinis-tic techniques for broad searches,with gradient-based methods for detailed optimization,a signi?cant topic of current research[3,4].

However issue(ii)is perhaps more critical in general,and that with which this article is concerned.Inevitably the availability of sensitivity analysis tools for codes lags behind the codes them-selves,a situation exacerbated by the considerable effort required to linearize complex solvers.A good example is the lack of reliable adjoint solvers for Navier–Stokes problems until recently.Initially the dif?cult and limited adjoint formulation in the continuous con-text slowed development[5],see 0f9759e85ef7ba0d4a733b0eter the stability of iterations for the viscous adjoint in many situations proved prob-lematic[6],examined in Section4.Even today there are few ad-joint solvers that apply linearized turbulence models to large problems–predominantly due to poor stability–instead they rely on constant eddy-viscosity approximations,resulting in an associ-ated gradient error,discussed in Section3.

The situation deteriorates when multi-disciplinary problems are considered,and the sensitivities of coupled multi-code systems are required.Nonetheless these problems are worth overcoming, as when they are available,gradient-based algorithms combined with adjoint gradients are the only methods which can offer rapid optimizations for extremely large numbers of design variables [7,8],as needed for the aerodynamic shape design of air vehicles.

In this work,therefore we consider each of the principle meth-ods of sensitivity evaluation in detail,and discuss their range of applicability and strengths and weaknesses with reference to pub-lished results.We present the derivation of each technique in a uniform notation,and collect references concerning both the his-tory,and recent developments in each method.We also present generalizations to higher-order derivatives and coupled systems, and recent applications to optimization from the authors.

The gradient evaluation problem may be succinctly stated as fol-lows:let JeaTbe a vector of objective functions(in fact functionals of the?ow solution such as lift and drag)and possibly constraints of dimension n J,where a is a vector of design variables of dimension n a which parameterize the problem,and in particular its geometry. Then the aim is to evaluate d J=d a,the gradients of the objective functions and constraints with respect to the design variables,most often for the purposes of optimization using gradient-based meth-ods such as steepest descent and conjugate gradients.

The obvious?rst approach is to apply?nite differences to val-ues of J evaluated at perturbed values of a,but it will be seen in Section2.1that this method has many drawbacks.Other than?-nite differences,techniques may be divided into two broad catego-ries whose relative ef?ciency depends on the values of n J and n a. In the direct approach a linearization of the governing equations is constructed,allowing the sensitivity of the entire?ow with re-spect to each a to be evaluated,given which sensitivities of J are easy and cheap to construct.The dominant cost in this case is the solution of n a linear equations.In contrast in the adjoint ap-proach a dual of the linearized problem is constructed for each J.The solutions of these dual problems allow easy and cheap eval-uation of sensitivities of J with respect to any design variable.The cost of evaluating the full derivative matrix d J=d a is therefore dominated by the cost of n J dual linear problems.In aerodynamics the adjoint approach is almost always preferable,as problems are typically characterized by a low number of cost functions(the inte-grated forces on a body),and a large number of design variables, required to parameterize curves and surfaces in three dimensions; hence n a)n J.

Another distinct classi?cation of sensitivity evaluation tech-niques is into continuous methods,where the governing equations are linearized and adjointed before being discretized,and discrete methods where the linearization and adjointing is performed on the discretization of the non-linear system.The issue of which ap-proach is most suitable under which circumstances has been a sub-ject of much research and contention over the past decade,and we devote Section2.7to the subject.

Perhaps the?rst use of this kind of sensitivity analysis in?uid dynamics was in1973and due to Pironneau,who derived the con-tinuous adjoint of Stokes equations[9]and later of the incompress-ible Euler equations[10],applying the result to pro?le optimization.However it was not until1988that Jameson ex-tended the procedure to inviscid compressible?ows,allowing the optimization of transonic aerofoil pro?les[11],after which point the method became truly widespread.The continuous direct

374J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391

method in contrast has been seldom considered in an aerodynamic context[12].

The discrete direct method was considered as early as1982by Bristow and Hawk[13,14]for panel methods,then in1989by Elbanna and Carlson[15]for the compressible Euler equations,fol-lowed by Taylor et al.[16]and Baysal and Eleshaky[17]in the early 1990s.The later in the same article also considered the adjoint of the discretized equations,which in the same year were considered by Shubin and Frank[18,19]for one-dimensional duct?ow.

Herein,we are concerned purely with the evaluation of gradi-ents;a review of aerodynamic shape optimization techniques may be found in Ref.[20],of multi-disciplinary design in Ref.

[21],and of sensitivity analysis and optimization for complex con-?gurations in Ref.[22].We also note a recent collection of litera-ture in aerodynamic design[23].

The structure of this article is as follows:the various gradient computation methods will be derived in Section2,and their accu-racy under various circumstances discussed in Section3.Iterative solution methods for the resulting equations,which are as ill-con-ditioned as the original?ow equations are described in Section4, with particular reference to duality preserving methods for the ad-joint equations.Section5is devoted to recent generalizations of gradient evaluation methods,namely to higher-order derivatives, adjoint mesh deformation,and multi-disciplinary problems.Final-ly,Section6presents two state-of-the-art gradient-based optimi-zation results using adjoint techniques developed by the present authors.Throughout references are presented to facilitate further reading and highlight points of practical relevance and areas of ongoing research.

2.Sensitivity evaluation methods

This section will derive in some detail the principle methods of gradient evaluation,?rstly?nite differences,followed by discrete and continuous versions of the direct and adjoint methods.We use the convention of uppercase symbols for discrete,and lower-case for continuous quantities,so that W is the solution of the dis-cretized governing equations ReW;XT?0on the mesh X.The discrete adjoint variable will be denoted K.The continuous?ow and adjoint solutions are w and k,respectively.Throughout,the de-sign variables a and objective functions JeaT?JeWeaT;XeaTTare taken to be vector quantities,so that d J=d a generally represents

a matrix of derivatives.

2.1.Finite differences

The application of?nite differences to an entire?ow solver is by far the simplest means of obtaining solution gradients,as it re-quires no modi?cation of the solver itself.As a result it was one of the?rst sensitivity evaluation methods to be used.

To proceed,the numerical?ow solution corresponding not only to a but also to perturbed states atd a and possibly aàd a is cal-culated.For the typical case of d a representing a geometry modi?-cation,this implies a mesh deformation Xeatd aT,and a new?ow solution on the modi?ed mesh satisfying

R Weatd aT;Xeatd aT

eT?0:

An approximation of objective functions derivatives in the direction d a can then be computed using a?nite difference formula,such as the second-order accurate formula

d JeaTd a d a’

1

2

Jeatd aTàJeaàd aT

?

?

1

2

JeWeatd aT;Xeatd aTTàJeWeaàd aT;Xeaàd aTT

? :

e1T

The entire matrix d JeaT=d a may be evaluated at a cost of2?n a

?ow solutions,or if a?rst-order difference is used n at1?ow solu-

tions,making the method impractical for large design spaces.

Another serious disadvantage is that the choice of the step-size

k d a k is critical to the accuracy of the result.If it is too small then

rounding errors become signi?cant;particularly if the convergence

level of the steady state computation is low.For example,if e is the

machine epsilon(the smallest number that may be added to1giv-

ing a number distinct from1),then

k d a k)e k a k;k d a k)e

j JeaTj

d JeaT

d a

;

are weak constraints on k d a k.In practice however J is rarely eval-

uated to machine accuracy,and k d a k must be increased accordingly.

On the other hand too large a value of k d a k invalidates the neglec-

tion of higher-order terms in the Taylor expansion in(1).The choice

is case and parameter dependent,and in practice the only means of

guaranteeing accuracy is to perform a convergence study on k d a k

for each parameter,at considerable cost.

Another dif?culty is that the evaluation of the J s at the two dif-

ferent a s in(1)should be as similar as possible,in order to aid error

cancellation.For example,the mesh topology should not be chan-

ged,implying some means of mesh deformation is required.It may

also be bene?cial to use the same initial solution and number of

iterations in both cases.For more complex con?gurations the

appearance of multiple solutions can also be an issue[24].

All these issues may be alleviated by using a complex?nite dif-

ference formula,such as

d J

d a

d a’I Jeati d aT

? ;

where I represents the imaginary part.As there is no longer any

difference of J s in this expression,it does not suffer from cancella-

tion error,and k d a k may be chosen as machine zero with no loss of

accuracy[25].However,the solver must be modi?ed to accept com-

plex variables throughout,negating the main advantage of?nite

differences.

Finite differences have been used since the1970s in the context of

shape optimization.Early contributions include works of Hicks et al.

[26]and Vanderplaats and Hicks[27]for aerofoil design,and later

wing design[28],using the method of feasible directions in which

the lift was held constant by moving in the design space only normal

to the lift gradient.Here,the?ow was modeled with a small pertur-

bation full potential equation solver.Early contribution also include

articles by Reneaux and Thibert[29,30].As3D optimization was con-

sidered with?nite difference gradients,the number of design

parameters and the cost of computations became severe problems.

In an attempt to reduce the number of parameters needed aerofunc-

tion shapes were introduced for example by Aidala et al.[31],de?ned

to be aerodynamically meaningful parameterizations of an aerofoils.

These might consist of geometry perturbations that control leading

edge pressure,or shock position[32,33].

However,the fundamental limitations of?nite differences have

lead to the investigation of alternative means of gradient

evaluation.

2.2.The discrete direct method

Under the assumption that the discrete residual R is once contin-

uously differentiable with respect to the?ow?eld and mesh in a

neighborhood of WeaTand XeaT,the discrete form of the governing

equations ReW;XT?0may be differentiated with respect to a to give

@R

@W

d W

d a

@R

@X

d X

d a

:e2T

J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391375

This may be regarded as a linear system in unknowns d W=d a, where d X=d a may be evaluated by?nite differences as in the previ-ous section,and the partial derivatives could be evaluated,for example,by hand.The dimension of the system is the number of degrees of freedom in the non-linear equations n W,and it can be re-garded as a linearization of those equations.Of course R is seldom differentiable everywhere,but in practice this rarely causes dif?cul-ties[34].

Given the n a solutions d W=d a the derivatives of J are

d JeaT

a?@J d W

at

@J d X

a;e3T

where again the partial derivatives are in principle easy to evaluate, as J is a known,explicit function of W and X.Hence the2?n a non-linear solutions required for second-order?ow?nite differences have been replaced by one non-linear and n a linear solutions,all of dimension n W,and the dependence on?nite differences has been con?ned to the relatively cheap mesh update procedure.

This method was considered as early as1982by Bristow and Hawk for a subsonic panel method[13,14],and again in1989for the transonic perturbations equations by Elbanna et al.[15].In the early1990s it was applied to the compressible Euler equation by two teams at Old Dominion University;that of Baysal[17,35–37]and that of Taylor and Hou[16].Many of these articles are con-cerned not only with gradient evaluation,but also the linearization of the solver as a tool for investigating small perturbation to the base?ow.On unstructured grids the idea was pursued by New-mann,Taylor et al.from1995onwards[38–40,22].

2.3.The discrete adjoint method

There are many ways to derive the discrete adjoint equations, the one given here is chosen for its similarity to the derivation of the continuous adjoint presented in the following sections.Let the direct linearization(2)be premultiplied by an arbitrary line vector K T of dimension n W,so that

K T @R d W

atK

T

@R d X

a

?0;8K2R n W:

Adding this expression to(3)

d JeaTd a ?

@J

@X

d X

d a

t

@J

@W

d W

d a

tK T

@R

@W

d W

d a

tK T

@R

@X

d X

d a

;8K2R n W;

e4T

and factorizing

d JeaTd a ?

@J

@W

tK T

@R

@W

d W

d a

t

@J

@X

d X

d a

tK T

@R

@X

d X

d a

;8K2R n W;

isolates the term d W=d a,which may be eliminated by choosing the arbitrary vector K to satisfy

@J

tK T @R

?0or equivalently

@R

T

K?à

@J

T

;

e5T

the adjoint equation,a linear system in unknowns K the adjoint solu-tion,with respect to the objective function J.Given K the sensitivi-ties may be written

d JeaT

a?@J d X

atK

T

@R d X

a

:

The critical point is that,because a does not appear in(5),that lin-ear system must only be solved once for each J.Hence the full ma-trix d J=d a may be evaluated at a cost of n J linear system solutions, substantially independent of n a.Note that the issues described in Section 2.2with regard to evaluation of the partial derivatives above also apply here.

Perhaps the?rst application of this method was given by Shu-bin and Frank in1991for a quasi one-dimensional nozzle?ow using the compressible Euler equations[18,41,19],and was de-noted there the‘‘implicit gradient approach”as contrast to the di-rect approach.Baysal et al.also recognized its potential,and offered it as an alternative to the direct approach when n J)n a [17,35–37].Later examples of the discrete adjoint are given in more speci?c contexts elsewhere in this article.

2.4.Evaluation of discrete Jacobians

It was simply stated in the previous sections that the partial derivatives in the above expressions,in particular the Jacobian @R=@W,are straightforward to construct.This,while strictly true, misrepresents the substantial effort involved.Although R is typi-cally a known,explicit function of W,it is also typically extremely complex in any non-trivial code,containing?ux-functions,recon-struction operators,characteristic boundary conditions,turbulence models,etc.,which must all be differentiated and reimplemented. Furthermore,the accuracy of the resulting gradients depends on the completeness and exactness of the linearization;unlike the Jacobian required for an implicit iterative solution scheme for example,which may be simpli?ed arbitrarily without affecting the?nal result of the computation,only its ef?ciency.By hand,per-forming the full linearization is a very considerable,but not com-pletely impractical undertaking[42].Nonetheless it is sometimes sensible to make simplifying assumptions by neglecting the differ-entiation of certain terms[43],thereby introducing an error into the gradients which must be examined numerically,see Section3.

A further complication is that the Jacobian,despite being a sparse matrix,has typically(depending on the discretization sten-cil size and the number of?ow equations)two orders of magnitude more non-zero entries than a solution vector,and if stored explic-itly would represent a memory bottleneck in most codes[7].This storage may be avoided by evaluation of Jacobian-vector products either entirely or partially on-the-?y,and applying solution meth-ods that require only these matrix-vector products.Of course the time required for the evaluation is also important,and partially on-the-?y evaluations offer an effective compromise.

An example of the memory costs involved,from the authors’own experience is given in Table1taken from[6],for two lineari-zations of the unstructured?nite volume solver TAU,with explicit time-integration.The?rst implementation stores the full Jacobian explicitly,and the second uses a partially on-the-?y evaluation in which terms of the spatial discretization with a stencil of immedi-ate neighbors only are pre-calculated and stored,whereas those involving next-neighbors are evaluated as needed.The initial non-linear code is provided for reference,and is extremely mem-ory ef?cient,able to perform a computation on a grid of1.2million points within1GB of memory.With the full Jacobian the code re-quires13.5times as much memory,which can be reduced to4 times with the on-the-?y implementation.The CPU cost of a single Table1

The memory requirements and run-time of the standard code and two linearized versions of the code for a representative3D(hybrid120k points)grid.

Standard TAU Full Jac.On-the-?y Jac.

Memory(bytes)100M1350M400M

Factor increase?1.0?13.5?4.0

Points?tting in1GB1:2?10689?103300?103

Time for a residual eval.?1.0?1.7?0.6

Time for the Jacobian eval.–?28.0?3.0

376J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391

linear residual evaluation (Jacobian-vector product),and the cost of building the Jacobian (or partial Jacobian)is given in the ?nal two rows of Table 1.In all cases the costs are of similar order to those of the non-linear code.

A technique which avoids the implementation and memory problems associated with the per-hand approach is using ?nite dif-ferences to evaluate Jacobian-vector products directly.As is the case for ?nite differences applied to the cost function the choice of step-size is vital to the accuracy.Evaluation of the transpose Jacobian-vector products needed by the adjoint is more dif?cult and was investigated by Anderson et al.[25],who used complex differences to maintain accuracy.They found however the cost of a single transpose Jacobian-vector product to be equivalent to about 12non-linear residual evaluations,an order of magnitude more expensive than the per-hand approach.

Outside the context of this article is the rapidly developing ?eld of automatic differentiation (AD)[44,45],which aims to provide tools (compilers or code generators)for automatically linearizing and adjointing computer codes,thereby reliving the development burden involved.We comment only that at the present time we do not know of any applications of such methods to large problems in aerodynamics,possibly due to the fact that the resulting codes tend to require very signi?cantly greater computational resources than hand-differentiated counterparts,despite the existence of good theoretical upper bounds on the costs [44].A more realistic immediate goal for AD could be in assisting hand-differentiation,for example by automatically differentiating and transposing indi-vidual routines [46].

2.5.Derivation of the continuous adjoint equations

In this approach,the adjoint of the continuous governing equations with respect to a given objective function is derived,before being discretized.As already mentioned,the ?rst appear-ance was due to Pironneau [9],with Jameson providing the ?rst treatment for compressible ?ows [11],which is similar to that given here.For a more detailed introductory treatment see Aber-gel and Temam [47],Giles and Pierce [48],and Jameson [49].The approach has since been reproduced by a variety of authors [50].

It is no longer easy to present the theory independently of the particular equations considered,therefore we consider the 2D Eu-ler equations in body-?tted coordinates in two dimensions.At the end of this section,references for more general cases are presented.It is assumed that the problem in physical space with a body-?t-ted structured grid,can be transformed into computational coordi-nates en ;g T,see Fig.1,in such a way that the transformation of D xy to D n g is direct,that D n g is a rectangular domain ?n min ;n max ??g min ;g max and that n min corresponds to the pro?le surface.Note that whereas the coordinate change operator depends on a ,D n g itself does not.

Let K be the determinant of the Jacobian of the coordinate transformation

K en ;g T?det

@x

@x g @y @n

@y @g

!

representing the change in size of a volume element under the transformation.Then the Euler equations in the computational coordinates are

@F eW Tt

@G eW T

g

?0;e6T

where

W ?K q

q u q v q E 0

B B B @1

C C C A ;F eW T?K q U q Uu tp @n @x q U v tp @n @y q UH 0B B B B @1C C C C A ;G eW T?K q V q Vu tp @g @x q V v tp @g @y q VH

0B B B B @1

C C C C A

:The slip-wall boundary condition is U ?0on n ?n min ,and a far?eld

condition is applied to the n max boundary.The objective function formulated in the new coordinate system is

J ea T?Z

n min

J 1ew T

@y

@g d g t

Z

D n g

J 2ew TK en ;g Td n d g ;

e7T

where the domain of integration is now independent of a .

The continuous adjoint equations can now be derived as fol-lows:?rst write the ?rst variation of the ?ow equations with re-spect to the design parameter a .Referring to (6)there are two distinct types of variation:the ?ux terms f ew Tand g ew Tvary with a ,because the ?ow changes in the transformed coordinate space when the shape changes,and independently all metric terms also depend on a :

f ew T!f ew Tt@f d w

a

d a ;

@x g !@x g t@2x g a

d a :Introduc

e a ew T?d

f ew T=d w and b ew T?d

g ew T=d w the ?ux Jacobi-ans,then the linearized equation corresponding to (6)is

@

a ew T@y g à

b ew T@x g d w a

&'

t@@g àa ew T@y @n tb ew T@x @n d w

d a &'

t@@n f ew T@2y @g @a àg ew T@2x @g @a ()t@@g àf ew T@2y @n @a tg ew T@2x @n @a ()?0;

e8T

and similarly for (7)

:

J.E.V.Peter,R.P.Dwight /Computers &Fluids 39(2010)373–391377

d JeaTd a ?

Z

n min

dJ

1

ewT

dw

d w

d a

@y

@g

tJ1ewT

@2y

@g@a

!

d g

t

Z

D n g

d J

2

ewTd w

a Ken;

gTtJ

2

ewT

@Ken;gT

a

d n d g:e9T

Before continuing it is convenient to introduce notation for the?ux Jacobian in the computational mesh directions n and g,

a1ew;n;gT?aewT@y

@g

àbewT

@x

@g

;a2ew;n;gT

?àaewT@y

tbewT

@x

:e10T

The idea behind the following procedure is to add to(9)the in-ner product of the linearized governing equations with an arbitrary four-component Lagrange multiplier k,analogously to the discrete case in Section2.3.Then we search for a condition on k for the gra-dient to be independent of the d w=d a terms.In this case we as-sume that the?ow and adjoint solutions,w and k,are once continuously differentiable with respect to the computational coordinates,i.e.,w;k2C1eD n gT4.Note that this is a very different restriction to that necessary in the discrete case.Proceeding from (8)we have that

Z D n g k T

@

@n

a1ew;n;gT

d w

d a

t

@

@g

a2ew;n;gT

d w

d a

&'

d n d g

t

Z

D n g k T

@

@n

fewT

@2y

@g@a

àgewT

@2x

@g@a

! (

t@

gàfewT

@2y

atgewT

@2x

a

!)

d n d g?0:

Using integration by parts,and the fact that the?ow sensitivity and coordinate derivatives such as@2y=@n@a are taken to be zero at the far?eld we have

à

Z

D n g @k T

@n

a1ew;n;gT

d w

d a

d n d gà

Z

D n g

@k T

@g

a2ew;n;gT

d w

d a

d n d g

à

Z

D n g @k T

@n

fewT

@2y

@g@a

àgewT

@2x

@g@a

!

à@k T

gàfewT

@2y

atgewT

@2x

a

!

d n d gt

Z

n min

k T a1ew;n;gT

?d w

d a

d gt

Z

n min

k T fewT

@2y

@g@a

àgewT

@2x

@g@a

!

d g

?0:

Adding this expression to(9)and extracting terms multiplying d w=d a we obtain,from the volume and surface integrals, respectively:

d J2ewTKen;gTà@k T a1ew;n;gTà@k T

g

a2ew;n;gT?0;on D n;g;

k T a1ew;n;gTtd J1ewT

d w @y

@g

?0;on n min;

8

<

:

e11Tthe continuous adjoint equations and boundary conditions.

One very signi?cant feature of this problem is that not all objec-

tive functions J

1

lead to a well posed adjoint boundary condition, because the?ux Jacobian a1is singular at a slip wall[11].For the compressible Euler equations described here functions of pressure are admissible,which is fortunate given that integral forces on a pro?le are thereby allowed.On the other hand for the Navier–Stokes equations there is no clear way of accounting for wall shear-stresses,and hence viscous drag[51].Given a solution of

(11)the gradients of J may be written

d JeaT

a?

Z

n min

J

1

ewT

@2y

g a d

gt

Z

n min

k T fewT

@2y

g aàgewT

@2x

g a

!

d g

à

Z

D n g

@k T

fewT

@2y

g aàgewT

@2x

g a

!

d n d g

à

Z

D n g

@k T

@g

àfewT

@2y

@n@a

tgewT

@2x

@n@a

!

d n d g

t

Z

D n g

J

2

ewT

@Ken;gT

@d a

d n d g:e12T

For the extension to the Navier–Stokes equations,also derived in curvilinear coordinates,see the seminary articles of Jameson et al.of1997[5],while a more discursive treatment is given in Ref.[52],which also includes considerations related to the use of the thus obtained gradients in shape optimization.

For?nite volumes on unstructured meshes such coordinate transforms as described above are not used,and a formulation in physical coordinates is necessary.One approach was?rst pub-lished by Anderson and Venkatakrishnan in1998[51],and1year later by Hiernaux and Essers[53,54].Also of interest is an early ad-joint formulation of the thin shear-layer equations in physical coordinates[55].

All theory presented in this section has been under the assump-tion of continuously differentiable?ow and adjoint solutions,and therefore is only valid for?ows without shocks.Given a disconti-nuity which moves due to some parameter perturbation,the change in the solution between the old shock location,and the new one is not small at all.To alleviate this restriction the integrals involved the derivation above may be split into parts downstream and upstream of a shock,before applying the integration by parts formula and imposing the Rankine–Hugoniot condition on the boundary.The extension of the adjoint to discontinuous?ows is a dif?cult topic,and the reader is referred to works by Iollo et al. [56,57],Cliff et al.[58],Giles and Pierce[59]and Gunburger[60], and most especially Pironneau et al.[61–63].

However despite the theory,it should be noted that dif?culties with discontinuous solutions have only been seen to occur numer-ically for special cases,and in transonic design problems they are not an issue.

Finally,we note that the continuous direct formulation,embod-ied by(8),rarely occurs in the gradient evaluation literature.Pelle-tier et al.[64,65]applied it to incompressible?ows,and the one of few compressible references is of Sharp and Sirovitch for hyper-sonic pro?le optimization[12].

2.6.Spatial discretization of the continuous adjoint equations

In the discrete case after derivation of the adjoint one is?nished in the sense that the method is completely prescribed,the spatial discretization of the linear problem having been directly derived from that of the non-linear problem.In the continuous case it re-mains to discretize the adjoint equation(11),and the expression for the gradient(12).

Given the close relationship of the non-linear and adjoint prob-lems,it is reasonable for discretization techniques of the original problem to be applied to the adjoint.In the case of compressible ?ows for example,upwind or arti?cial viscosity?ux formulations are appropriate,and most authors use similar schemes for?ow analysis and gradient computation,and even some of the same subroutines[49].As a result the discretization expression and implementation tends to be simpler than in the discrete case, where derivatives of the discretization are required.

378J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391

For completeness we present a single example discretization of (11)(and with zero volume component to the objective function, J

2

?0),taken from Ref.[11],and used subsequently without mod-i?cation by Brezillon and Gauger[66].A central?ux with arti?cial viscosity is used,on cells in a structured grid denoted by their indi-ces i;j:

K i;j @k i;j

@t

?

1

2

^a

1

ek it1;jàk ià1;jTt^a2ek i;jt1àk i;jà1T

àá

td

it1

2

;j

àd

ià1

2

;j

td

i;jt1

2

àd

i;jà1

2

;

where d are diffusion?uxes on the control volume faces it1;j,etc., and^a1and^a2are discretized versions of the?ux Jacobians(10), with aewTand bewTevaluated at cell centers and

@x @n ’

1

2

@x

@n

ià1;j

t

@x

@n

it1;j

!

;

@y

@g

1

2

@y

@g

i;jà1

t

@y

@g

i;jt1

!

;

etc.,an average of metrics on each side of the cell.Note that because of this last approximation the?ux through a face in one direction is in general different to that in the other,and as a consequence the discretization is non-conservative,corresponding to the fact that the adjoint equations are also not conservative.This latter point is easy to overlook given their origin in a conservation equation,and initially it is tempting to try to impose a?ux balance using the mod-i?ed?uxe@f=@wTT k averaged on cell faces,which naturally fails.

The dominant contributors to the continuous adjoint in aerody-namics must be considered to be Jameson and coworkers,and two detailed descriptions of their approach to discretization may be found in Refs.[5,67].For the convective terms they use a Jameson–Schmidt-Turkel,Roe,or H-CUSP?ux.Viscous?uxes use a cell-ver-tex approximation of velocity and temperature gradients on a dual mesh.The?ow and adjoint equations are separately advanced to a steady state using a Runge–Kutta scheme speci?cally tuned for each system,respectively,and a multigrid procedure.Other signif-icant contributions have been made by Giles and Pierce[55,48], Weinerfelt and Enoksson[68]and Brezillon et al.[66,50].The continuous adjoint on unstructured grids has been pursued by Anderson and Venkatakrishnan[51],Soto et al.[69],Jameson et al.[70],Castro et al.[71],and Widhalm et al.[72]among others.

2.7.Discrete versus continuous adjoint

The discrete and continuous formulations of the adjoint prob-lem clearly represent two very different ways of doing essentially the same thing,and the question naturally arises:which method should be preferred in a given situation?This issue proved quite contentious in the early days of adjoint development,where sev-eral groups promoted their own approaches,and the accuracy and range of validity of each approach was still unclear.Since that time it has been shown by a variety of authors that both the dis-crete and continuous formulations give more than suf?cient accu-racy in practice,and in fact the dominant error source in a typical adjoint gradient evaluation is likely to be the?nite difference approximation of the metric sensitivities d X=d a[6].1As a conse-quence the question of which method is most appropriate has be-come largely an implementational one.

Several authors have dealt with the issue directly by comparing codes based on continuous and discrete formulations,notably Giles and Pierce[48]and Nadarajah and Jameson[73].Among the points raised are the following:–The discrete adjoint gives the exact gradient of the discrete objective function(agreeing with the results of?nite differ-ences),whereas the continuous adjoint gives an approximation to the continuous gradient based on some alternative discretiza-tion.Discrete gradients are therefore bene?cial for an optimizer, which receives gradient information consistent with objective function evaluations.Also the agreement with?nite differences eases debugging and veri?cation of the discrete code.

–A discrete adjoint code may be derived in a purely algorithmic fashion from an existing code,as demonstrated by the existence of automatic differentiation tools[44],whereas a continuous for-mulation often requires domain-speci?c insight to perform the discretization.

–Because of the close relationship of the discrete adjoint and the ?ow discretization,it is often possible to modify the iterative sol-ver of the?ow problem to give an iterative solver for the adjoint problem that converges at the same asymptotic rate,see Section4.–The continuous code is often considerably simpler than the dis-crete in terms of operation count and memory requirements,as well as implementationally.However it is quite possible to build an ef?cient discrete code without explicit storage of@R=@W(at a cost of some small inaccuracy in the gradients),mitigating the ?rst two points[6].

–Mathematical and physical understanding of the adjoint prob-lem is aided by consideration the continuous approach,as exem-pli?ed by analytic adjoint solutions for1D Euler problems[74].–As already noted in Section2.5,the continuous equations are only well posed for certain objective functions,most seriously impeding the evaluation of shear–stress drag gradients in vis-cous?ow.In order words boundary conditions do not exist for every objective function in the continuous case,whereas they are automatic in the discrete case.

–Additional source terms in turbulence models complicate the derivation of the continuous adjoint.To the best of our knowl-edge no such treatment has yet been published.

–Generalization to higher-order derivative evaluation is much easier in the discrete case,as addressed in Section5.2.The cor-responding continuous formulation is too involved to be addressed here.

There have also been efforts towards hybrid schemes[73],where the discrete informs the discretization of the continuous adjoint,or the discrete is used to extend the set of allowed objective functions of the continuous,or the discrete is used for the turbulence equations and the continuous for the mean-?ow equations.However in much such situations it is easier to implement a fully discrete code.

One pertinent question is:assuming smooth solutions,do the two approaches converge to the same result as mesh spacing tends to zero?It turns out that a consistent discretization of the?ow equations is not suf?cient for its adjoint be a consistent discretiza-tion of the continuous adjoint equations.A discretization that does have this property is known as adjoint consistent,and for?nite ele-ment methods corresponding conditions on the discretization have been derived[75].

3.Accuracy of sensitivity computations

In many circumstances it is much easier and cheaper to obtain some approximation of the solution of the adjoint problem,and hence the gradient,than it is to obtain the exact(discrete)gradient. The question then arises:how accurate must gradients be for the application under consideration?In optimization,steepest descent and(restarted)conjugate-gradient algorithms have proven very robust to poor quality gradient data[43].This is not surprising: provided the gradient direction is not more than90°incorrect

1While these may be evaluated to high accuracy if great care is taken and a

parameter study performed on the step-size for each design variable,in practice such

attention to detail is computationally expensive and dif?cult to automate.

J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391379

these methods will produce an improvement in the objective func-tion at each iteration.BFGS on the other hand attempts to con-struct a local quadratic model of the design space using the gradient,and this process compounds gradient error.Also if gradi-ents are to be used in a response surface method,for example,gra-dient-enhanced Kriging,then poor gradients will lead to a poor response surface[76].However that still leaves many applications where only good approximations to gradients are necessary.

This section principally concerns the discrete adjoint,using sim-pli?cations of the Jacobian,and their effect on the accuracy of resulting gradients.Given an exact and complete linearization of the discrete?ow solver excellent agreement with?nite differences is regularly seen in the literature[77,43,48].However,because of the effort required to construct this linearization for all but the most modest?ow solvers,it is perhaps not surprising that in a large proportion of references some simplifying approximation is made,thereby affecting the accuracy of the gradient.In the contin-uous adjoint,approximations are also made due to the dif?culty of the treatment of turbulence models,which are not always formulated as continuous equations.In fact the neglection of the linearization of turbulence models,commonly denoted the frozen eddy-viscosity approximation,is so widely used that it may be con-sidered standard,and is treated in detail in Section3.1.Several other common approximations are discussed brie?y in Section3.2.

Apart from these explicitly added sources of inaccuracy there are more fundamental and unresolved issues regarding the mesh con-vergence of adjoint gradients.Giles et al.[78]considered a case with a strong shock which wandered between mesh points as a design parameter was varied.The lift of the aerofoil was seen to be some-what dependent on the position of the shock within the local mesh cell,as a result of which the design parameter-lift curve took on a slightly scalloped shape.As the mesh was re?ned the amplitude of the scalloping was reduced as expected–but as it always took the same shape the amplitude of its gradient was not reduced.Fortu-nately the effect was not seen for cases with weaker shocks,and has not appeared elsewhere,although we are not aware of any other attempts to demonstrate grid convergence of gradients,adjoint or otherwise.

3.1.Fully linearized or frozen turbulence modeling

Though frozen eddy-viscosity is increasingly the norm,there are indeed many authors who have linearized turbulence models by hand or using AD.To our knowledge three classes of models have been linearized in the literature:

(i)The algebraic model of Baldwin–Lomax was differentiated by,

for example,Le Moigne and Qin[79]and Kim et al.[80],both of whom used their method for pro?le optimization.Pham[81] adjointed the algebraic Michel model for evaluating sensitiv-ities in three dimensions in turbomachinery.

(ii)The one-equation model of Spalart–Allmaras was linearized by under others Giles et al.[82],and the team of Anderson, Nielsen,and Bonhaus at NASA Langley[83–85]who com-pared resulting gradients with frozen eddy-viscosity gradi-ents[77].The complete linearized code was applied to the 3D optimization of an isolated wing.Further examples include Nemec and Zingg[86]and Brezillon and Dwight

[87]who both presented pro?le optimizations,in the latter

case with multiple design points and constraints.

(iii)The two-equation transport models k— ,k—x SST.and Wil-cox k—x have been differentiated and applied to sub-and transonic pro?le design,as well as high-lift pro?le and set-ting optimization by Kim et al.[88,89].The former model has also been adjointed in the context of turbomachinery by Renac et al.[90,91].

There is however a notable lack of linearized transport equation

turbulence models applied to con?gurations signi?cantly more

complex than isolated wings with fully attached?ow[77]or2D

high-lift pro?les[89].We suspect this is consequence,not of the

dif?culty of performing the linearization itself or accounting for

the coupling with the mean-?ow,but of the problems associated

with the solution of the resulting linear system,which may be

exceptionally poorly conditioned[42].

In any case in a very large number of articles the turbulence is

?xed.In both the continuous and discrete cases this is achieved

simply by adjointing only the mean-?ow equations and treating

turbulent quantities appearing therein as constants.For one-equa-

tion models these are just the eddy-viscosity,for two-equation

models the turbulent kinetic energy k also appears in the internal

energy and hence in the pressure.Hence in addition simplifying

the solution of the system,the number of equations to be solved

is reduced.Notable examples of this approach are due to Jameson

et al.[92],Soemarwoto[93],and Valentin[94].

There are additionally several authors that are principally inter-

ested in AD,and the linearization of turbulence models as a second-

ary issue:Hou et al.[95]for Baldwin–Lomax,and Mohammadi[96]

for a k— model.As a related point of interest,Bischof et al.[97]used an AD tool in forwards mode to determine the sensitivity of?ow over

a backward-facing step(particularly the re-attachment point)to

empirical parameters of several turbulence models,the idea being

that a large sensitivity to an experimentally determined parameter

is a weakness of the model.

3.2.Other approximations

After frozen eddy-viscosity,perhaps the second most common

approximation is that of neglecting the derivatives of arti?cial dif-

fusion coef?cients in convective terms when using a centered con-

vective?ux approach[98,43,85],which has almost no effect on the

gradient while allowing a considerable simpli?cation of the corre-

sponding linearized terms.

A much stronger approximation is the use of?rst-order convec-

tive?uxes in the linear calculation[99–103],which has the advan-

tage of reducing the?ll-in of the system matrix to immediate mesh

neighbors only,and reducing the stiffness of the problem

dramatically.

In some articles authors have linearized a simpli?ed viscous

?ux[104,105],and in a very limited number of articles neglected

viscous terms completely[67,106],although this latter has several

serious complications(no-slip boundary conditions must be re-

placed,and a separate Euler mesh is necessary for the adjoint).

Perhaps the most extreme‘‘approximation”–which can barely

be called adjoint at all–was suggested by Mohammadi[107–109],

who proposed ignoring all aerodynamic contributions to the gradi-

ent on the basis that in some very particular cases metric terms

dominate,i.e.

d JeaT

a?

@J d X

at

@J d W

a’

@J d X

a;e13T

thereby completely obviating calculation of the aerodynamic ad-joint.Mohammadi suggested this might be useful in the case where the objection function is an integral of a boundary quantity,and the design parameters deform the shape of that boundary.However at best a gradient vector qualitatively resembling the correct vector will be obtained.Subsequently Kim et al.[110]compared the idea with a full gradient for the design of supersonic transport high-lift devices,and discovered that for their case,while giving reasonable results at the leading edge it performed poorly at the trailing edge, where aerodynamic variations were stronger.

380J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391

3.3.Impact on,and in?uence of gradient accuracy

Given the widespread use of adjoint approximations there have been relatively few published studies examining the effect on the gradient of various approximations.Nielsen and Anderson [103,77]examined the effect of?rst-order?ux and frozen eddy-viscosity approximations for a high-lift con?guration,and ob-served errors in gradients of several design variables of typically 90%with the former and30%with the later method.On this basis they rejected these approximations.

However a pertinent question in this context is:how important is the accuracy of the gradient to the optimization speed and accu-racy?Kim et al.[89,111]tackled this question for the frozen eddy-viscosity assumption by performing the same optimization using gradients from an adjoint with and without the approximation. They observed signi?cant discrepancies in the gradient that for some cases resulted in signi?cantly different optimal geometries when using a BFGS optimizer,which builds an approximation to the Hessian from the gradients.

A similar study[43]that considered?ve separate adjoint approx-imations,showed on the other hand that,provided a robust opti-mizer is used(in this case conjugate gradients(CG)),then even very poor gradients result in almost identical optima to highly accu-rate gradients.This effect is shown in Figs.2and3reproduced from Ref.[43].The?rst?gure shows the variation in gradient for a tran-sonic RAE2822with20shape design variables,for the adjoint approximations considered,where it is evident that,e.g.,frozen eddy-viscosity introduces a large error.The second shows the con-vergence of drag minimization optimizations using these gradients: in every case the same optimum is reached to within an accuracy of 4%.Similar results hold for high-lift optimization,and–consistent with the results of Kim–the use of the more gradient sensitive BFGS damages this independence.This result perhaps explains in part the use of and reliance on these approximations,though inaccurate gra-dients do complicate and decrease the robustness of the design iter-ation,and there is some problem dependence[8].

4.Solution strategies for the primal and adjoint equations

The discrete adjoint and direct equations are both linear prob-lems,so it is somewhat counter-intuitive that–for large test cases –they are as hard to solve as the original?ow equations and often signi?cantly harder.This is due to the demands of stability:itera-tions applied to linear problems never stall as is common for?ow problem,if they do not converge asymptotically they diverge,com-bined with the dif?culty of CPU and memory ef?cient evaluation of the linearized and adjoint residual.

From the simplest perspective the(discretized)continuous or discrete adjoint equations are nothing other than sparse linear

J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391381

systems,and so in principle amenable to the standard solution tech-niques for such problems,for example,Jacobi,Gauss–Seidel or incomplete lower–upper (ILU)factorization preconditioned Krylov methods.Indeed among authors who ?rst considered discrete sensitivity analysis,direct inversion was possible –because of the small size and banded structure of the Jacobians considered –and popular [17,35,36].The straightforward application of ILU precondi-tioned GMRES to the adjoint problem has been advocated by Nemec and Zingg [86],who apply it in a similar manner within a Newton method for the ?ow problem.This technique can be extremely ef?-cient,giving a method that for 2D cases provides a solution in about 5%of the CPU time required for the ?ow calculation if ILU(4)is used [43].However the memory requirements can be 10–30times those of the original solver due to the storage of the factorized matrix and the Krylov basis,and the effectiveness of ILU as a preconditioner has been seen to decrease rapidly with problem size and stiffness –most particularly in three dimensions –requiring an increase of the ?ll-in level to impractical values [43,86].

These dif?culties have been alleviated somewhat by several authors including Newman et al.[112],Nielsen et al.[103]and Ne-mec and Zingg [113],who all applied ILU (usually ILU(0))to a linear-ization of a ?rst-order accurate ?ow discretization.The idea is that if

@R @W

1st W n t1

àW n ?àR eW n T;

is an effective iterative method for the non-linear problem,where @R =@W 1st is some approximate simpli?ed Jacobian,then

@R @W

1st "

#T

K n t1àK n ?@J @W T à@R @W

T K n ;

e14T

should be an effective method for the adjoint problem.The use of the ?rst-order Jacobian on the LHS reduces the stiffness of the linear system to be solved at each iteration so that ILU(0)rather than ILU(4)or higher is suf?cient.(As an aside an interesting paper on the conditioning of @R =@W for upwind and central convective ?ux discretizations with and without low Mach number preconditioning is due to Baysal and Eleshaky [17].)

The general principle is that because the ?ow and adjoint prob-lems are closely related in a very particular way,if an iterative method is effective for the non-linear problem,then some suitable ‘‘adjointing”of that method should be effective for the adjoint problem.In the example above this adjointing is performed simply by the transpose operation on the LHS of (14).The principle was formalized by Giles [114],and is discussed in detail in the follow-ing section.It is of particular importance because the ef?cient solu-tion of ?ow problems has been a subject of intense research since the beginnings of CFD,and a wide variety of effective algorithms have been developed such as multigrid,Runge–Kutta with residual smoothing,and approximately factored implicit methods,all of which may be adjointed and applied to the adjoint problem with the expectation of identical convergence behavior.Increasingly therefore variants of iterations used in the non-linear case are ap-plied to the adjoint case in practice.

This is not quite the end of the story however,as there are regu-larly situations in practice where the original iterative method for the ?ow equations stalls,implying divergence for the corresponding adjoint iteration.This may be remedied by reintroducing a Krylov method as a form of stabilization,see Section 4.2.With these tech-niques the solution of the adjoint problem can be made as robust and ef?cient as that of the original ?ow problem,and two examples are given.

4.1.Duality preserving ?xed-point iterations

An insight due to Giles [114,82,78]based on work by Christianson [115]instructs the construction of dual iterations that give rates of

convergence of the adjoint problem identical to those of the non-lin-ear problem.The method has since been widely adopted [116,98,6].In essence,if a particular ?xed-point iteration has a given asymptotic rate of convergence for a given non-linear problem,the same iteration applied to the corresponding linearized problem will have the same rate of convergence,assuming that the lineari-zation is exact,complete,and based on the fully converged non-linear solution.In practical applications this is a very desirable property:given a successful non-linear computation a linear com-putation can be immediately performed with the same iterative method (and CFL number),and convergence within a ?xed number of iterations is guaranteed.It will now be seen how this property may be extended to the adjoint equations.

A general ?xed-point iteration,including, e.g.,Runge–Kutta,multigrid and implicit schemes,may be written

M eW n t1àW n T?àR eW n T;

e15T

for some preconditioning matrix M .By rearranging (15),and linear-izing about the exact solution f W

,we have W n t1?eI àM à1A TW n tg ef W

TtO k W n àf W k 2;e16T

where A ?@R =@W is the Jacobian evaluated at the exact solution and g eáTis some function independent of n .A close relation to the linearized problem is apparent,and informed by this we consider the ?xed-point iteration (FPI)

M eH

n t1

àH n

T?@R @X d X

d a

àA H n ;

e17T

or rearranged:

H

n t1

?eI àM à1A TH n tM

à1

@R @X d X

d a

:

e18T

The coef?cients of W n and H n in these two iterations are identical,so we can conclude that –for a suf?ciently converged non-linear iteration –the errors reduce at the same rate;i.e.,the asymptotic convergence behavior is identical.The rate of convergence itself is given by the dominant eigenvalue of eI àM à1A T.

To achieve the same for the adjoint equation consider the FPI

M T eK

n t1

àK n

T?

@J

!T

àA T K n ;

e19T

so that

K

n t1

?eI àM àT A T TK n tM

àT

@J

!T

:e20T

where M àT denotes the inverse transpose of M .Since any real ma-trix and its transpose have identical eigenspectra,the asymptotic convergence rate of (20)will be the same as that of (18).For implicit methods constructing the adjoint FPI is easy when M is explicitly available [117];Runge–Kutta may be treated by ?rst constructing M as a matrix,transposing,and then deconstructing [118].

A further convergence statement may be made,by requiring that H 0?K 0?0and expanding (18)and (20)to the N th iteration.Then we have,respectively

H

N t1

?àX

N n ?0

eI àM à1

A Tn

M à1

@R d X

a

;

e21TK N t1

?àX N n ?0

eI àM àT A T Tn M àT

@J !T ;

e22T

given which it quickly follows that

K N t1

h i T @R @X d X d a ?@J @W

!T

H N t1;

e23T

i.e.,that the primal and adjoint results for d J =d a will be identical not

only after the equations have fully converged,but also at every

382J.E.V.Peter,R.P.Dwight /Computers &Fluids 39(2010)373–391

intermediate step.This relation is useful in verifying the correctness of the adjoint FPI and discretization[118,98].As a numerical exam-ple of these results see Fig.4reproduced from[6].The gradients produced by the primal and adjoint codes are identical to machine accuracy at every iteration,as well as the convergence rates of both methods being identical.

4.2.Krylov stabilization

Despite the guarantees regarding convergence provided by the theory of the previous section for both the primal and adjoint equa-tions,for problems involving turbulence modeling there are regu-larly situations in which it is possible to obtain a reasonably converged solution of the non-linear equations,but not of the corre-sponding linear equations.This can occur for two reasons,both of which violate the parity of(15)and(17).Firstly there could a dis-crepancy between the linear and non-linear problem,for example, due to a Jacobian approximation.This effect leads to some promi-nent authors to advocate the use of only the exact linearization when constructing the adjoint[116,119].However,there is a second dif?-culty:the FPI applied to the non-linear problem might not converge asymptotically itself,but rather stall in some limit cycle,which in fact is almost the standard situation in complex3D applications.This behavior does not necessarily indicate a de?cient iterative method, but may often be traced to physical instationary behavior that is not correctly captured by the stationary?ow approximation.

In such a position one simple solution is to add extra arti?cial dissipation to the adjoint problem which is not present in the?ow problem.This may be done explicitly[104,91],or indirectly by lin-earizing only?rst-order convective?uxes as already discussed [99,43].The remarks of Section 3.3regarding the accuracy of resulting gradient apply,and the technique tends only to delay the point at which the iteration starts to diverge rather than removing the instability altogether.

A preferable alternative is to apply a Krylov method.Such meth-ods have been in use for some time in the context of linearized?ow solvers[120,113],but only recently have they come to be regarded as a stabilization technique for an existing iteration,rather than simply as a generic linear solver.This perspective was?rst advanced in the context of sensitivity analysis in CFD by Campobasso and Giles [121,122].Campobasso showed on the basis of eigensystem analysis that the number of diverging modes in the linearized problem is gen-erally quite small,of the order of10s,even for problems with mil-lions of degrees of freedom.The standard iterative method(which may be a duality preserving iteration in the adjoint case)damps all other modes more-or-less effectively,but the presence of even one diverging mode is fatal.Krylov methods identify the most strongly diverging modes of a system and treat them specially,in the case of GMRES[123]by directly minimizing the residual on the subspace spanned by those modes,in the case of the recursive projection method(RPM)[124,125]by using a Newton iteration.Provided all diverging modes are so treated,the resulting iteration will be stable.

Campobasso chose RPM over the more commonly used GMRES, as it has the property of only adding a vector to the Krylov basis(to be specially treated)when an instability is detected,thus minimiz-ing the size of the basis required for stability and thereby the mem-ory overhead.Note that for many Krylov methods(including GMRES and RPM)convergence behavior is a function only of the distribution of the eigenvalues of the system matrix,and therefore transposing the matrix has no effect on convergence.Hence such Krylov methods are in some sense automatically duality preserv-ing.This approach has been applied since by several groups[7] including the current second author[6].

The effect on convergence of applying a Krylov solver can be quite dramatic,as shown in Fig.5for adjoint high Reynolds num-ber viscous?ow about a DLR-F6transport aircraft wing-body con-?guration.The short dotted line shows the rapidity with which the duality preserving iteration diverges without stabilization.With RPM activated,initially the calculation proceeds without modi?ca-tion until a diverging mode is identi?ed.When found this mode is added to the solution subspace P,which will be treated with a Newton method,while the complementary subspace Q will be treated with the standard iteration as before.Shown in Fig.5are the?ow residuals on each of the subspaces P and Q,as well as the dimension of P and the convergence of an integral quantity de-rived from the adjoint solution,d C D=d a.The jumps in the residuals correspond to the identi?cation of a diverging mode,and subse-quent modi?cation of P and Q.As can be seen only12modes are needed in P to bring the calculation to machine accuracy.

Furthermore,if after a while RPM does not identify any more diverging modes,it will add modes which converge most slowly

J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391383

to P,thereby improving the asymptotic convergence rate,an effect which can be observed in the difference in convergence rate before and after the2000th iteration in Fig.5.

5.Generalizations of gradient evaluation

Up until now we have considered the evaluation of?rst deriva-tives of quantities only dependent on the?ow.We discuss now three generalizations:extension of the adjoint procedure to include the mesh deformation,to higher-order derivatives,and to coupled structure–?uid systems.The extension to unsteady?ows,including optimization of time-dependent problems,has recently been per-formed by Mani and Mavriplis[126],and will not be discussed here.

5.1.Adjoint mesh deformation

Through the use of adjoint?ow equations the cost of a gradient evaluation for(typically)50design variables or more becomes dominated by the cost of deforming the mesh for each design var-iable,necessary for the evaluation of the metric sensitivities by?-nite differences.In the continuous case a partial solution to this is given by approximate surface formulations of the volume integral for the gradient,?rst proposed by Weinerfelt and Enoksson[68], which require only deformation of the surface,but which come with a loss of accuracy.In the discrete case there is a more elegant and exact solution due to Nielsen and Park[119]who suggested constructing additionally the adjoint of the mesh deformation equations.

Let the mesh X satisfy DeX;aT?0.Differentiating this equation with respect to a,multiplying by a mesh deformation adjoint var-iable C,and adding to(4)from Section2.3,which already contains the linearization of the?ow equation gives:

d JeaTd a ?

@J

@X

d X

d a

t

@J

@W

d W

d a

tK T

@R

@W

d W

d a

tK T

@R

@X

d X

d a

tC T

@D

@X

d X

d a

tC T

@D

@a

8K;C2R n W

As before,d W=d a can be eliminated by factorizing it out and setting its coef?cient to zero using K,giving the familiar?ow adjoint equa-tion.However now the same may be done for d X=d a using C,giving

@D T

C?à@J T

à

@R T

K;

where K is already known,an adjoint mesh deformation equation. Note that the additional partial derivatives@R=@X and@J=@X must be available.Similar dif?culties apply as for the evaluation of @R=@W,with the signi?cant difference that the metric sensitivities must only be evaluated once,and hence ef?ciency is of reduced importance.

Given C the objective function derivative becomes

d JeaTd a ?C

@D

@a

T

:

Hence n J adjoint mesh deformations must be performed rather than n a direct mesh deformations.The gradient evaluation is almost completely independent of the number of design variables,allowing optimization of3D con?gurations using all surface mesh points as design variables[119].Mavriplis[7]did the same for a3D transport aircraft wing-body design problem,using this technique among others to perform a complete Navier–Stokes drag minimization optimization in about6h on a small cluster of eight processors. The idea has also recently been used by Jakobsson and Amoignon [127],who applied a single volume spline for both mesh deforma-tion and geometry parameterization.Differentiating the spline al-lowed construction of the gradients with respect to the design variables directly.

Alternatively the gradient computation can be split in two parts, making the aerodynamic gradient post-processing and adjoint equation resolution fully independent of the geometric design pro-cess.This technique is a consequence of the relation

d JeaT

d a

?

@J

@X

d X

d a

tK T

@R

@X

d X

d a

?

@J

@X

tK T

@R

@X

d X

d a

;

whereby the two factors on the right-hand side may be calculated independently.The drawback of the method compared to more standard discrete gradient computation method is that the matrix @R

@X

is needed explicitly,rather than just the product@R

@X

d X

d a

,which can be easily estimated by?nite differences.

5.2.Higher-order derivatives using discrete methods

Higher-order derivatives have applications in optimization algorithms based on Newton methods,where a root of the gradient of J is sought,and hence second-derivatives of J are required. There are of course many applications of higher-order Taylor expansions of?ow solvers,for example,response surface construc-tion.Another use of great recent interest is in uncertainty analysis, where it is often desirable to compute the probability distribution of an output variable given those of several input 0f9759e85ef7ba0d4a733b0eing a ?rst-order expansion the mean and variance of the output distribu-tion may be modeled,but second-order is necessary to account for skewness,third-order for kurtosis,and so on.

Here a generalization of the adjoint treatment to second-order derivatives only is considered,the associated development,accu-racy,and ef?ciency problems already being very signi?cant,and becoming much worse for higher derivatives.

As in Section2,we consider the effort(in terms of linear system solutions required)to evaluate the full tensor of second-derivatives for several approaches.Whereas there the choice was between di-rect or adjoint,here there are four possibilities,corresponding to the two methods of solving for the?rst and second derivatives. These were?rst investigated by Sherman et al.[128],and we fol-low his convention of denoting them DD.DD,DD.AV,AV.DD,and AV.AV,respectively,where DD stands for‘‘direct differentiation”and AV‘‘adjoint vector”.The objective functions,?ow solutions and discrete?ow equations are now assumed to be locally twice differentiable with respect to their arguments.

Because of the complexity of the vector algebra we introduce i,j, to index degrees of freedom of the?ow solution,m,n those of the mesh,and a,b,c to index entries in the objective function and design variable vectors,and use summation convention.Then(3)becomes

d J

c

a a?

@J

c

m

d X m

a at

@J

c

i

d W i

a a;

and second derivatives are

d2J c

a a a b?

@2J

c

m n

d X m

a a

d X n

a bt

@J

c

m

d2X m

a a a b

t

@2J

c

@W i@X m

d X m

d a a

d W i

d a b

t

d X m

d a b

d W i

d a a

t

@2J

c

@W i@W j

d W i

d a a

d W j

d a b

t

@J

c

@W i

d2W i

d a a d a b

;e24T

so that there are then n aen at1T=2second-order?ow derivatives d W=d a d a required to evaluate the second-derivative tensor.

384J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391

5.2.1.DD.DD

The most direct approach is to solve (2)for each d W =d a term,and derive the corresponding second-order equation by differenti-ating (2)with respect to a :

@R @W i

d W i d a a ?à@R @X m d X m

d a a ;

@R i d 2

W i a a a b ?à

@2R m n d X m a a d X n

a b

à@R @X m d 2

X m d a a d a b à@2

R @X m @W i d X m d a a d W i d a b t

d X m d a b d W i d a a

à

@2R @W i @W j d W i d a a d W j

d a b

;

e25T

giving a total of n a en a t1T

2

tn a linear systems to be solved.Once found the ?ow derivatives may be substituted into (24)to obtain the de-sired gradient.

At this point it is worth examining the effort involved in applying this method.First Eq.(25)must be constructed requiring the evalua-tions of terms of the form @R =@W @W ,etc.and d X =d a d a .As before,the former are –in principle –easy to evaluate,R and J being known,ex-plicit functions of their arguments.However,it is now hard to see how the derivation could be performed by hand,and automatic dif-ferentiation tools seem to offer the only practical approach.Note also that explicit storage of the tensor @R =@W @W in memory is likely to be impossible for all but the smallest problems.As for the metric deriv-atives,evaluation using ?nite differences is extremely step-size sen-sitive,and will likely lead to large errors.In addition the mesh must be deformed O en 2a Ttimes,though this may be avoided by adjointing the mesh deformation equation as in Section 5.1.

Given these dif?culties the numerical solution of the resulting linear systems is a very secondary problem.Though in DD.DD the number of systems is quadratic in n a we can exploit the fact that they all have the same system matrix,and possibly factorize @R =@W if the problem is small.In the following it will be shown how the AV formulations can serve to reduce the number of sys-tems to be solved,but they do nothing to help the signi?cant prob-lems in formulation;nonetheless these methods have been implemented and applied to problems in aerodynamics [128].5.2.2.DD.AV

The second method consists in adding to (24)the product of Eq.(25)by an arbitrary adjoint vector C in a similar procedure to that in Section 2.3.C is then chosen to cancel the second-order ?ow sensitivity d W =d a a d a b :

d 2J d a a d a b ?@2J @X m @X n d X m d a a d X n d a b t

@J @X m d 2

X m

d a a d a b

t@2J @W i @X m d X m d a a d W i d a b td X m d a b d W i

d a a

t

@2J i j d W i a a d W j a b t@J i d 2

W i

a a a b

tK T

@R @W i d 2

W i d a a d a b tK

T

@2R @X m @X n d X m d a a d X n d a b

tK T @R @X m d 2

X m d a a d a b tK T

@2R @X m @W i d X m d a a d W i d a b t

d X m d a b d W i d a a

tK T

@2R i @W j d W i d a a d W j

d a b

:

e26T

The adjoint equation to solve turns out to be the same as for the ?rst-order gradient computation,

@R @W T K ?à@J

@W

T

;where K is independent of a again.Having solved the n a direct equations for the d W =d a terms,and n J adjoint equations for the K ,i.e.,n a tn J linear systems,the derivative tensor is

d 2J d a a d a b ?@2J @X m @X n d X m d a a d X n d a b t

@J @X m d 2

X m

d a a d a b

t@2J @W i @X m d X m d a a d W i d a b t

d X m d a b d W i d a a

t@2J

@W i @W j ?d W i d a a d W j d a b

tK T @2R @X m @X n d X m d a a d X n d a b tK T

@R @X m ?d 2

X m a a a b tK T @2R m i d X m a a d W i a b t

d X m a b d W i a a tK T

@2R @W i @W j d W i d a a d W j

d a b

;

e27T

where no d W =d a d a terms appear.

5.2.3.AV.DD

We consider AV.DD for completeness,despite the fact that it is always more costly than DD.AV.The starting point for AV.DD and AV.AV is the familiar adjoint expression

d J a a ?@J d X a a tK T @R d X

a a

;@R T K ?à@J T

;which de?nes K as a function of W ea Tand X ea T,and hence a .Differ-entiating both expressions with respect to a design parameter a b

therefore gives

d 2J d a a d a b ?@2J @X i @X j d X i d a a d X j d a b t@2J @X i @W j d X i d a a d W j d a b t

@J @X i d 2

X i

d a a d a b

t

d K T

a b !

@R i d X i a a tK T @2R @X i @X j d X i d a a d X j d a b t@2

R @X i @W j d X i d a a d W j d a b t

@R @X i d 2X i

d a a d a b

!

;e28T

and

@R i @W l d K i d a b ?à@2R i @W l @W j d W j d a b K i à@2R i @W l @X j d X j

d a b

K i

à@2J l j d W j a b à@2J l j d X j

a b

:e29T

Required are therefore n a ?ow derivatives,n J adjoint solutions,and

n J ?n a adjoint derivatives,giving n a tn J tn a n J linear problems.5.2.4.AV.AV

The ?nal method of [128]extends AV.DD by introducing a sec-ond adjoint variable C .The product of C with (29)is added to (28),allowing the elimination of the derivatives of the ?rst adjoint var-iable K ,if C satis?es

@R C ?à@R @X

a

:But this is exactly the direct equation,so C ?d W =d a ,and the expres-sion for the derivative reduces to that of DD.AV (27).Again only the original n a tn J systems must be solved,as for DD.AV.The four for-mulations differ predominantly in the number of linear systems to be solved,and may be chosen on the basis of n a and n J .For typical problems we therefore expect DD.AV/AV.AV to be most appropriate.5.2.5.Applications

Given the considerable obstacles to the evaluation of second derivatives,it is perhaps surprising the techniques described here

J.E.V.Peter,R.P.Dwight /Computers &Fluids 39(2010)373–391385

have been implemented for aerodynamic?ow codes at all.In fact several groups have independently performed the necessary work, though it is telling that all used automatic differentiation in some form.Sherman,Taylor et 0f9759e85ef7ba0d4a733b0eing ADIFOR[129],successfully tested DD.DD for two2D viscous aerofoils,one with turbulence [128].Much more recently they applied both DD.AV and AV.AV to inviscid?ow about a NACA0012[130].The small test cases and the relatively slow pace of progress are perhaps a testament to the current limitations of AD codes(adjointing or backward mode is signi?cantly harder than simple differentiation or forward mode for AD tools[44]).

Based on a purely theoretical derivation corresponding to(DD)n by Guillaume and Masmoudi[131],Aubert et al.applied AD to eval-uate the partial derivatives of R and J,and applied the resulting code to the design of2D and3D low speed gas mixing devices[132].By restricting design variables to aerodynamic values at the injection boundary,derivatives with respect to the mesh could be avoided. In further work the code was coupled to a genetic algorithm,the idea being to provide cheap objective function evaluation based on a sec-ond-order Taylor expansion near to known values[4].The problem considered was optimization of a turbomachinery blade.

Most recently Ghate and Giles[46]constructed the Hessian of a simple code via DD.AV,including mesh sensitivities,by applying automatic differentiation extremely selectively.Rather than apply-ing TAPENADE[133]to the entire solver in one shot,it was used to differentiate individual?ux evaluation routines and boundary con-ditions.The remaining summation loops required to build the full derivatives of the residual were then constructed by hand,being substantially identical to the loops of the?ow solver itself.This mode of use of AD has the potential to scale well to the much larger codes in use in industry,as well as to third-and higher-order derivatives.

5.3.Extension to aeroelasticity

The need for multi-disciplinary modeling of aircraft is well known,and is follows that multi-disciplinary optimization(MDO) is also necessary.Even if a purely aerodynamic objective function is considered,a wing optimized to be shock free in an unloaded(or jig)shape is unlikely to remain so when bent and twisted under a load.Mono-disciplinary optimization of the wing under a?xed structural deformation has also been shown to be insuf?cient [134,135].The issues become even more complex when objective functions are also multi-disciplinary.Consider the trade-off be-tween low aerodynamic drag–which requires a thin wing–and low structural weight which requires a thick wing(at least near the root).

Extension of the direct discrete gradient formulas for multi-disci-plinary problems were presented as early as1990by Sobieszczan-ski-Sobieski[136].Aeroelastic optimization was?rst carried out using?nite differences gradients,Haftka[137],Bowman et al. [138],Friedmann[139],Barthelemy et al.[140]and Giunta and Sob-ieszczanski-Sobieski[141].Later on,in order to reduce the cost of gradient evaluation,direct and adjoint methods for the multi-disci-plinary equations were considered,?rst for2D con?gurations by Ghattas and Li[142]and Moller and Lund[143],then for3D con?g-urations by Maute et al.[144],Hou and Satyanarayana[145],Maute et al.[146],Martins et al.[147],and Fazzolari[135].

These multi-disciplinary discrete sensitivity equations are brie?y summarized in the remaining part of this subsection for an aeroelastic problem.In addition to existing notation introduce the surface displacement?eld U,and the structure mesh X S. Meshes X and X S are taken to match at solid walls and represent the jig-shape.This framework is the most widely used;see Farhat et al.for an alternative approach[148,144].

The objective function is taken to be JeaT?Jea;XeaT;X SeaT; W;UT,and the state equation reads Rea;XeaT;X SeaT;W;UT?eR AeXeaT;W;UT;R SeX SeaT;W;UTTT?0; where R A and R S are aerodynamic and structural parts,respectively. In practical aeronautic applications the dependency of R S on W cor-responds to the aerodynamic stresses applied at the solid wall, whereas the dependency of R A on U corresponds to the displace-ment of the aeroelastic structure.De?ne Y?eW;UT,so that differ-entiating the state equation gives

d R

d a

?

@R

@Y

d Y

d a

t

@R

@X

d X

d a

t

@R

@X S

d X S

d a

?0;

which may be separated into structural and aerodynamic parts

@R A

@W

@R A

@U

@R S@R S

!

d W

d a

d U

a

!

t

@R A

@X

@R S

!

d X

d a

t

@R A

S

@R S

S

!

d X S

d a

?0;

and resolution for d Y=d a allows the computation of objective func-tion gradients as

d JeaT

a?

@J d Y

at

@J d X

at

@J

S

d X S

a:

Yet again,for each state equation an adjoint vector is introduced which multiplies the linearization of that equation,and contributes to the expression for the objective function derivative,in order to factorize out the dif?cult terms:

d J

a?

@J d Y

at

@J d X

at

@J

S

d X S

atK

T

@R d Y

at

@R d X

at

@R

S

d X S

a

!

?

@J

@X

tK T

@R

@X

!d X

d a

t

@J

@X S

tK T

@R

@X S

!d X

S

d a

t

@J

@Y

tK T

@R

@Y

!d Y

d a

;

where K??K A;K S T corresponding to the decomposition of R into R A and R S,satis?es

@R A

@W

T@R

S

@W

T

@R A

@U

T@R

S

@U

T

@

1

A K A

K S

?

à@J

@W

T

à@J

@U

T

@

1

A:

In practice the resolution of this multi-disciplinary system uses an iteration on lagged mono-disciplinary solutions

@R A T KepT

A

?à@J Tà@R S T Kepà1T

S

@R S T KepT

S

?à@J Tà@R A T Kepà1T

A

8

<

:;

converging to the solution of the coupled problem.As for the mono-disciplinary case n J such systems must be resolved to obtain the complete derivative matrix.

6.Examples of sensitivity applications

In closing we brie?y present some applications of sensitivity analysis to optimization.For the period up to1999the review of Newmann et al.[22]already mentioned,provides a comprehensive overview of applications in gradient-based design for aerodynam-ics,so we do not repeat that content here.For the time since we give a limited number of milestones results,to provide an impres-sion of the progress made.

A good test of the progress in gradient-based methods over the last ten years is the high-?delity coupled structure–?uid optimiza-tion of a transport aircraft wing,where the objective is to optimize a multi-disciplinary function–such as range–subject to a large number of both aerodynamic and structural design variables and constraints(including perhaps wing planform variation),necessi-tating an adjoint formulation.This is clearly a problem of great engineering relevance,which has to date,and to the authors’knowledge,remained unsolved.

386J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391

By1999several groups had performed purely aerodynamic shape optimization of transport aircraft wings[149–151],whereby typically only pro?le design parameters were used,as without struc-tural penalties wing planform optimization tends to produce wings with unrealistically large span.Starting in2000aeroelastic designs began to be performed,but always with a predetermined structure and planform[144,147,135].In2004,Leoviriyakit and Jameson [152,153]performed wing planform optimization using a semi-empirical analytical model for wing weight,taking into account aerodynamic loading and planform shape,but without elastic defor-

mation.Finally in2007Jameson et al.[154]brought previous work together with an aeroelastic,multi-objective,planform optimiza-tion,but without modi?cations to the wing structure.

In other areas,in2003Campobasso et al.[155,156]performed an unsteady aerodynamic optimization of turbine blades using an adjointed frequency-domain code to avoid adjointing full time-accurate 0f9759e85ef7ba0d4a733b0etely Mani and Mavriplis[126]have adjointed a full unsteady?ow solver,and used it to minimize the time-dependent drag pro?le of a2D aerofoil without any loss of lift,for example.For this the adjoint of the volume mesh deforma-tion algorithm was also required.We have also already mentioned the?rst use of mesh sensitivities in adjoint based optimization in 2005by Nielsen et al.[119,7].

Finally,we present three recent applications of gradient-based optimization from the authors’personal experience.The?rst uses a continuous adjoint on a structured grid in inviscid?ow,the remaining two consider the Navier–Stokes equations and adjoints on unstructured and structured grids,respectively.

6.1.Navier–Stokes optimization of a wing–fuselage con?guration

The only reason to employ an adjoint code is to evaluate sensitiv-ities with respect to a large number of design variables rapidly, which raises the question:are a large number of variables really nec-essary in design?This issue was partially addressed by Brezillon et al.[8]for the optimization of the DLR-F6transonic wing–fuselage con?guration.The wing was parameterized using a free-form defor-mation(FFD)control box whose node coordinates formed the design variables.Fig.6shows such a box for a very low number of parame-ters;more detailed control is possible simply by adding nodes.The wing thickness is held constant by coupling nodes on the upper and lower sides of the box.A structured mesh is used in order to re-duce the number of mesh points required for the required accuracy.

Gradients were computed using a discrete adjoint solver imple-mented in the unstructured DLR Tau-Code[157]using a duality preserving LU-SGS smoothed multigrid method.Frozen eddy-vis-cosity was used,and for stabilization the recursive projection method was applied[6].

The optimization is drag minimization at constant lift.Satisfy-ing the lift constraint requires sensitivities of the lift as well as the drag,hence two adjoint calculations are performed per gradi-ent evaluation,and the angle-of-attack a is chosen to achieve the desired lift.The optimization is performed using conjugate gradi-ents for two different parameterizations,one coarse with nine vari-ables including twist,and one?ne with42variables.In Table2, reproduced from Ref.[8],the results are given including drag reduction for the wing and fuselage independently.

What is apparent is that the coarse parameterization achieved a drag reduction purely on the fuselage by means of reducing the an-gle-of-attack.The drag on the wing was actually increased as the optimizer attempted to preserve the lift.In contrast the?ne parameterization allowed a signi?cant improvement on the

wing, Fig.6.DLR F6wing–fuselage showing parameterization of the wing with a free-form deformation box.Black circles denote nodes used as design variables.

Table2

Drag improvement for coarse and?ne parameterized optimizations.

Parameters942

D ae Tà0.30à0.29

C L0.5000.500

D C D?10à4Totalà1.4à17.7

D C D?10à4Wing+1.2à16.1

D C D?10à4Fuselageà2.6à1.6

J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391387

and a total drag reduction of18counts.Because there is no wing-body fairing this con?guration always has a region of recirculation in the corner where the wing and fuselage meet,and this separa-tion is the dominant source of removable drag.In Fig.7,the region of separation is plotted for the two optimized con?gurations with the original separation boundary marked.

Although the result of any such study must be very dependent on the method of parameterization,and an informed choice may reduce the number of variables required signi?cantly,it seems that –at least in general–there are substantial advantages to large numbers of variables,and hence adjoint is essential.

6.2.Multi-disciplinary optimization of an engine pylon

An optimization much closer to a concrete engineering applica-tion was performed by Mouton et al.[158,159]who coupled aero-dynamic and structural optimization of the outboard pylon of a generic civil transport aircraft within the European project VIVACE, see Fig.8from Ref.[158]for the computational grid.In order to incorporate advanced CFD in the form of the adjoint RANS equa-tions into this multi-disciplinary optimization without building a structure–?uid coupled adjoint as in Section5.3,they followed a two-level parameterization strategy where three design variables were identi?ed as having a strong coupling role in the problem–namely the vertical and horizontal position of the pylon and its width–the remaining variables of pylon fairing and structural frame shape being almost discipline speci?c(the latter due to the free space available under the fairing).A coarse sampling of the strongly discipline coupling variables was made,and a disci-pline speci?c optimization performed at each sample point.These individual results were fed back into the global optimization prob-lem using a Kriging method.

The individual aerodynamic optimizations were performed using a discrete adjoint solver implemented in the block-structured elsA code using an adjointed implicit method with added dissipation for stabilization,and frozen eddy-viscosity[91].The pylon was parameterized using19design variables and the optimizations were performed with a BFGS method.The relatively large number of de-sign variables and repeated optimizations necessitated a rapid opti-mization strategy,hence the use of adjoint gradients in this

case.

0f9759e85ef7ba0d4a733b0eparison of region of corner separation on the F6before and after

optimization.

Fig.8.Surface mesh for the generic transport aircraft with close up of the pylon.

388J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391

The convergence of a typical adjoint computation involved in this optimization is shown in Fig.9.At regular intervals a selection of gradients obtained from the non-fully converged adjoint solu-tion are plotted in order to evaluate the accuracy of the transient adjoint solution.In this case after500multigrid cycles no signi?-cant further variation in gradients is seen and the computation is halted.For comparison gradients obtained from?nite differences are also plotted and relatively good agreement is seen.

Finally the optimization over all structure and shape variables parameterizing the pylon resulted in a total drag reduction of1:3 counts.

7.Conclusion

A wide variety of gradient computation approaches have been investigated in the past20years.Initially continuous adjoint ap-proaches dominated,being readily constructed for the Euler equa-tions,not requiring linearization of large complex discretizations, and solved with more-or-less the same techniques as the?ow equa-tions.More recently,with wider use of RANS codes,the most practi-cal algorithms for optimization,have been discrete–most often adjoint–approaches,where the Jacobian is obtained by hand,and where the resulting system is solved using duality preserving itera-tions stabilized with Krylov methods.For smaller problems today the Jacobian may be obtained with AD,and solved with generic lin-ear system solution methods,such as ILU preconditioned Krylov methods.

The situation is unlikely to remain so for long however.The trend in simulation today is towards increased coupling and com-plexity of physical modeling,and there is a demand for not just ?ow sensitivities,but mesh and parameter sensitivities.Lineariz-ing a code by hand in all these variations is likely to become increasingly impractical,whereas constructing an accurate contin-uous formulation is likely to be impossible.Hence,if corresponding sensitivity codes are needed they must be obtained with some greater-or-lesser reliance on AD.

However a common mode of application of AD,of applying indiscriminately to an entire?ow code,results today in large run-times and memory requirements,which nonetheless might be solved by future developments in AD compilers.However,more fundamentally,such usage does not take into account the need for extra stabilization of the linear iteration,done with Krylov meth-ods and Jacobian reduction today.

It is our belief that,at least in the medium term,industrially rel-evant linearized codes will be developed by using AD to differenti-ate individual non-linear routines,which are then assembled by an expert,in the model of Ghate and Giles[46],using the many tricks learnt in the past20years.

Acknowledgments

The authors express their gratitude to their colleagues:Jo?l Brezillon for permitting the use of results and?gures from Ref.

[8],and for his many useful comments,as well as Meryem Marc-elet,Chi-Tuan Pham,and Florent Renac for their useful comments and contributions to Sections3.1and6.Daniel Destarac also pro-vided invaluable corrections in the?nal stages of preparation. References

[1]Sobieszczanski-Sobieski J.The case for aerodynamic sensitivity analysis.

Technical Report CP2457,NASA;1987.

[2]Brezillon J,Wild J.Evaluation of different optimization strategies for the

design of a high-lift?ap device.In:Evolutionary and deterministic methods for design.Munich:EUROGEN;2005.

[3]Lombardi G,Mengali G,Beux F.A hybrid genetic based optimization

procedure for aircraft conceptual analysis.Optim Eng2006;7(2):151–71. [4]Kelner V,Grondin G,Léonard O,Moreau S.Multi-objective optimization of a

fan blade by coupling a genetic algorithm and a parametric solver.In: Proceedings of EUROGEN,Munich;2005.

[5]Jameson A,Martinelli L,Pierce N.Optimum aerodynamic design using the

Navier–Stokes equations.Theor Comput Fluid Dyn1998;10(1):213–37. [6]Dwight R,Brezillon J,Vollmer D.Ef?cient algorithms for solution of the

adjoint compressible Navier–Stokes equations with applications.In: Proceedings of the ONERA-DLR aerospace symposium(ODAS),Toulouse;

2006.

[7]Mavriplis D.Discrete adjoint-based approach for optimization problems on

three-dimensional unstructured meshes.AIAA J2007;45(4):740–52.

[8]Brezillon J,Brodersen O,Dwight R,Ronzheimer A,Wild J.Development and

application of a?exible and ef?cient environment for aerodynamic shape optimisation.In:Proceedings of the ONERA-DLR aerospace symposium(ODAS), Toulouse;2006.

[9]Pironneau O.On optimum pro?les in Stokes?ow.J Fluid Mech1973;59(1):

117–28.

[10]Pironneau O.On optimum design in?uid mechanics.J Fluid Mech1974;

64(2):97–110.

[11]Jameson A.Aerodynamic design via control theory.J Sci Comput1988;

3(3):233–60.

[12]Sharp H,Sirovitch L.Constructing a continuous parameter range of

computational?ows.AIAA J1989;27(10):1326–31.

[13]Bristow D,Hawk J.Subsonic panel method for the ef?cient analysis of

multiple geometry perturbations.Technical Report CR3528,NASA;1982. [14]Bristow D,Hawk J,Subsonic panel method for designing wing surface from

pressure distribution.Technical Report CR3713,NASA;1983.

[15]Elbanna H,Carlson L.Determination of aerodynamic sensitivity coef?cients in

the transonic and supersonic regimes.In:AIAA paper series,Paper89-0532;

1989.

[16]Taylor III A,Korivi V,Hou G.Sensitivity analysis applied to the Euler

equations:a feasibility study with emphasis on variation of geometric shape.

In:AIAA paper series,Paper91-0173;1991.

[17]Baysal O,Eleshaky M.Aerodynamic design sensitivity analysis methods for

the compressible Euler equations.J Fluids Eng1991;113(4):681–8.

[18]Shubin G,Frank P.A comparison of implicit gradient approach and the

variational approach to aerodynamic design optimization.Technical Report AMS-TR-163,Boeing Computer Service,Applied Mathematics and Statistics, June;1991.

[19]Frank P,Shubin G.A comparison of optimisation-based approaches for a

model computational aerodynamics design problem.J Comput Phys 1992;98:74–89.

[20]Dulikravitch G.Aerodynamic shape design optimization:status and trends.J

Aircraft1992;29(6):1020–6.

[21]Sobieszczanski-Sobieski J,Haftka R.Multidisciplinary aerospace design

optimization:survey of recent developments.Struct Optim1997;14(1): 153–60.

[22]Newman III J,Taylor III A,Barnwell A,Newman P,Hou G.Overview of

sensitivity analysis and shape optimization for complex aerodynamic con?gurations.J Aircraft1999;36(1):97–116.

[23]Uthup B,Koruthu S-P,Sharma R-K,Priyadarshi P,editors.Recent trends in

aerospace design and optimization.New Delhi:Tata-McGraw Hill;2005. [24]Haftka R.Sensitivity calculations for iteratively solved problems.Int J Numer

Methods Eng1985;21:1535–46.

[25]Anderson W,Newman J,Whit?eld D,Nielsen E.Sensitivity analysis for the

Navier–Stokes equations on unstructured meshes using complex variables.

AIAA J2001;39(1):56–63.

J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391389

[26]Hicks R,Murman E,VanderPlaats G.An assessment of airfoil design by

numerical optimization.Technical Report TMX3092,NASA;1974.

[27]VanderPlaats G,Hicks R.Numerical airfoil optimization using a reduced

number of design coordinates.Technical Report TMX73151,NASA;1976. [28]Hicks R,Henne P.Wing design by numerical optimization.In:AIAA paper

series,Paper77-1247;1977.

[29]Reneaux J.Méthode de dé?nition de pro?ls par optimisation numé0f9759e85ef7ba0d4a733b0e

Rech Aerospatiale1984;5:303–21.

[30]Reneaux J,Thibert J.The use of numerical optimization for airfoil design.In:

AIAA paper series,Paper85-5026;1985.

[31]Aidala P,Davis W,Mason W.Smart aerodynamic optimization.In:AIAA paper

series,Paper83-1863;1983.

[32]Reuther J,Cliff S,Hicks R,van Dam C.Practical design optimization of wing/

body con?gurations using the Euler equations.In:AIAA paper series,Paper 92-2633;1992.

[33]Hager J,Eyi S,Lee K.Design ef?ciency evaluation for transonic airfoil

optimization:a case for Navier–Stokes design.In:AIAA paper series,Paper 93-3112;1993.

[34]Matsuzawa T,Hafez M.Optimum shape design using adjoint equations for

compressible?ows with shock 0f9759e85ef7ba0d4a733b0eput Fluid Dyn J1998;7(3):74–89.

[35]Baysal O,Eleshaky M,Burgreen G.Aerodynamic shape optimization using

sensitivity analysis on third-order Euler equations.In:AIAA paper series, Paper91-1577;1991.

[36]Baysal O,Eleshaky M.Aerodynamic design optimization using sensitivity

analysis and computational?uid dynamics.AIAA J1992;30:718–25.

[37]Baysal O.Flow analysis and design optimization methods for nozzle

afterbody of a hypersonic vehicle.Technical Report CR4431,NASA;1992. [38]Newman III J,Taylor III A,Burgreen G.An unstructured grid approach to

sensitivity analysis and shape optimization using the Euler equations.In: AIAA paper series,Paper95-1646;1995.

[39]Newman III J,Taylor III A.Three dimensional aerodynamic shape sensitivity

analysis and design optimization using the Euler equations on unstructured grids.In:AIAA paper series,Paper96-2464;1996.

[40]Taylor III A,Newman III J,Barnwell R.Aerodynamic shape sensitivity analysis

and design optimization of complex con?gurations on unstructured grids.In: AIAA paper series,Paper97-2275;1997.

[41]Shubin G.Obtaining cheap optimization gradients from computational

aerodynamics codes.Technical Report AMS-TR-164,Boeing Computer Service,Applied Mathematics and Statistics,June;1991.

[42]Dwight R,Brezillon J.Adjoint algorithms for the optimization of3d turbulent

con?gurations.In:Proceedings of the15th STAB Symposium;2006.

[43]Dwight R,Brezillon J.Effect of approximations of the discrete adjoint on

gradient-based optimization.AIAA J2006;44(12):3022–71.

[44]Griewank A.Evaluating derivatives,principles and techniques of algorithmic

differentiation,No.19in Frontiers in Applied Mathematics,SIAM, Philadelphia;2000,ISBN08-987-1451-6.

[45]Corliss G,Le Faure C,Griewank A,Hascoét L,Naumann U,editors.Automatic

differentiation:from simulation to optimization,computer and information science.Springer;2001.

[46]Ghate D,Giles M.Hessian calculation using automatic differentiation.In:

AIAA paper series,Paper2007-4059;2007.

[47]Abergel F,Temam R.On some control problems in?uid mechanics.Theor

Comput Fluid Dyn1990;1(1):303–25.

[48]Giles M,Pierce N.An introduction to the adjoint approach to design.In:

Proceedings of ERCOFTAC workshop on adjoint methods;1999.

[49]Jameson A.Optimum aerodynamic design via boundary control.In:AGARD-

FDP-VKI special course;1994.

[50]Gauger N,Brezillon J.Aerodynamic shape optimization using the adjoint

method.J Aero Soc India2002;54(3):247–54.

[51]Anderson W,Venkatakrishnan V.Aerodynamic design optimization on

unstructured grids with a continuous adjoint 0f9759e85ef7ba0d4a733b0eput Fluids 1998;29:443–80.

[52]Jameson A.Aerodynamic shape optimization using the adjoint method.In:

VKI course in optimization;2003.

[53]Hiernaux S,Essers J-A.An optimal control theory based algorithm to solve2d

aerodynamic shape optimisation problems for inviscid and viscous?ows.In: Proceedings of the RTO-AVT symposium on aerodynamic design and optimisation of?ight vehicles;1999.

[54]Hiernaux S,Hessers J.Aerodynamic optimization using Navier–Stokes

equations and optimal control theory.In:AIAA paper series,Paper99-3297;1999.

[55]Giles M,Pierce N.Adjoint equations in CFD:duality,boundary conditions and

solution behaviour.In:AIAA paper series,Paper97-1850;1997.

[56]Iollo A,Salas M,Ta’asan S.Shape optimization governed by the Euler

equations using an adjoint method.Technical Report93-78,ICASE;1993. [57]Iollo A,Salas M.Contribution to the optimal shape design of two-dimensional

internal?ows with embedded shocks.J Comput Phys1996;125:124–34. [58]Cliff E,Heikenschloss M,Shenoy A.On the optimality system for a1d Euler

?ow problem.In:AIAA paper series,Paper96-3993;1996.

[59]Giles M,Pierce N.On the properties of solutions of the adjoint Euler

equations.In:Proceedings of the6th ICFD conference on numerical method for?uid dynamics,Oxford;1998.

[60]Gunzburger M.Sensitivities,adjoint and?ow optimization.Int J Numer

Methods Fluids1999;35:53–78.

[61]Bardos C,Pironneau O.Derivatives and control in the presence of shocks.

Comput Fluid Dyn J2003;11(4):383–92.[62]Mohammadi B,Pironneau O.Optimal shape design for?uids.Ann Rev Fluids

Mech2004;36:255–79.

[63]Pironneau O.Shape sensitivities and design for?uids with shocks.Int J

Comput Fluid Dyn2003;17(14):235–42.

[64]Pelletier D,Turgeon E,Lacasse D,Boorggaard J.Adaptivity,sensitivity and

uncertainty:towards standards in CFD.In:AIAA paper series,Paper2001-0192;2001.

[65]Pelletier D,Turgeon E,Etienne S,Boorggaard J.Reliable sensitivity analysis via

an adaptive sensitivity equation method.In:AIAA Paper Series,Paper2002-2578;2002.

[66]Brezillon J,Gauger N.2D and3D aerodynamic shape optimization using the

adjoint approach.Aerospace Sci Technol J2004;8(8):715–27.

[67]Jameson A,Alonso J,Reuther J,Martinelli L,Vassberg J.Aerodynamic shape

optimization techniques based on control theory.In:AIAA paper series,Paper 98-2538;1998.

[68]Weinerfelt P,Enoksson O.Numerical methods for aerodynamic optimization.

In:Proceedings of the8th international symposium on CFD,Swansea,Wales;

1999.

[69]Soto O,L?hner R,Yang C.An adjoint-based design methodology for CFD

problems.Int J Numer Methods Heat Fluid Flow2004;14(6):734–59.

[70]Jameson A,Sriram,Martinelli L,Haimes B.Aerodynamic shape optimization

of complete aircraft con?gurations using unstructured grids.In:AIAA paper series,Paper2007-4060;2007.

[71]Castro C,Lozano C,Palacios F,Zuazua E.A systematic continuous adjoint

approach to viscous aerodynamic design on unstructured grids.In:AIAA paper series,Paper2006-0051;2006.

[72]Widhalm M,Ronzheimer A,Hepperle 0f9759e85ef7ba0d4a733b0eparison between gradient-free

and adjoint based aerodynamic optimization of a?ying wing transport aircraft in the preliminary design.In:AIAA paper series,Paper2007-4060;

2007.

[73]Nadarajah S,Jameson A.Studies of continuous and discrete adjoint

approaches to viscous automatic aerodynamic shape optimization.In:AIAA paper series,Paper2001-2530;2001.

[74]Giles M,Pierce N.Analytic adjoint solutions for the quasi-one-dimensional

Euler equations.J Fluid Mech2001;426:327–45.

[75]Hartmann R,Derivation of an adjoint consistent discontinuous Galerkin

discretization of the compressible Euler equations.In:Lube G,Rapin G, editors.Proceedings of the BAIL conference,Goettingen;2006.

[76]Dwight R,Han Z-H.Ef?cient uncertainty quanti?cation using gradient-

enhanced Kriging.In:AIAA paper series,Paper2009-2276;2009.

[77]Nielsen E,Anderson W.Aerodynamic design optimization on unstructured

meshes using the Navier–Stokes equations.AIAA J1999;37(11):185–91. [78]Giles M,Duta M,Muller J-D,Pierce N.Algorithm developments for discrete

adjoint methods.AIAA J2003;41(2):198–205.

[79]Le Moigne A,Qin N.A discrete adjoint method for aerodynamic sensitivities

for Navier–Stokes.In:Proceedings of CEAS Cambridge;2002.

[80]Kim C,Kim C,Rho O,Lee S.Aerodynamic sensitivity analysis for the Navier–

Stokes equations.In:AIAA paper series,Paper99-0402;1999.

[81]Pham C.Linéarisation du?ux visqueux des equations de Navier–Stokes et de

modéles de turbulence pour l’optimisation aérodynamique en turbomachines.

Ph.D.thesis,L’Ecole Nationale Supérieure d’Arts et Métiers,September;2006.

[82]Giles M,Duta M,Muller J.Adjoint code developments using the exact discrete

approach.In:AIAA paper series,Paper2001-2596;2001.

[83]Anderson W,Bonhaus D.Aerodynamic design on unstructured grids for

turbulent?ows.Technical Report TM112867,NASA;1997.

[84]Anderson W,Bonhaus D.Airfoil design optimization on unstructured grids for

turbulent?ows.AIAA J1999;37(2):185–91.

[85]Nielsen E,Anderson W.Recent improvements in aerodynamic design

optimization on unstructured meshes.AIAA J2002;40(6):1155–63.

[86]Nemec N,Zingg D.Towards ef?cient aerodynamic shape optimization based on

the Navier–Stokes equations.In:AIAA paper series,Paper2001-2532;2001.

[87]Brezillon J,Dwight R.Discrete adjoint of the Navier–Stokes equations for

aerodynamic shape optimization.In:Evolutionary and deterministic methods for design.Munich:EUROGEN;2005.

[88]Kim C,Kim C,Rho O.Sensitivity analysis for the Navier–Stokes equations with

two equations turbulence models.AIAA J2001;39(5):838–45.

[89]Kim C,Kim C,Rho O.Effects of constant eddy viscosity assumption on

gradient-based design optimization.In:AIAA paper series,Paper2002-0262;

2002.

[90]Renac F,Pham C-T,Peter J.Sensitivity analysis for the RANS equations

coupled with a linearized turbulence model.In:AIAA paper series,Paper 2007-3948;2007.

[91]Peter J,Mayeur J.Improving accuracy and robustness of a discrete direct

differentiation method and discrete adjoint method for aerodynamic shape optimization.In:Proceedings of ECCOMAS,Egmond ann Zee;2006.

[92]Jameson A,Pierce N,Martinelli L.Optimum aerodynamic design using

the Navier–Stokes equations.In:AIAA paper series,Paper97-101;

1997.

[93]Soemarwoto B.The variational method for aerodynamic optimization using

the Navier–Stokes equations.Technical Report97-71,ICASE,December;

1997.

[94]Valentin V.Optimisation aérodynamique3d des aubages dans les

turbomachines axiales multi-etages.Ph.D.thesis,UniversitéParis6,June;2002.

[95]Hou G,Maroju V,Taylor A,Korivi V,Newman P.Transonic turbulent airfoil

design optimization with automatic differentiation in incremental iterative form.In:AIAA paper series,Paper97-0101;1995.

390J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391

[96]Mohammadi B.A new optimal shape design procedure for inviscid and

viscous?ows.Int J Numer Methods Fluids1997;25(2):183–203.

[97]Bischof C,Buecker H,Rasch A.Sensitivity analysis of turbulence models using

automatic differentiation.SIAM J Sci Comput2005;26(2):510–22.

[98]Mavriplis D.Multigrid solution of the discrete adjoint for optimization

problems on unstructured meshes.AIAA J2006;44(1):42–50,Paper0001-1452.

[99]Beux F,Dervieux A.Exact-gradient shape optimization of a2d Euler?ow.

Finite Elem Anal Design1992;12:281–302.

[100]Beux F,Dervieux A.A hierarchical approach for shape optimization.Technical Report1868,INRIA;1993.

[101]Marco N,Beux F.Multilevel optimization:Application to one-shot shape optimum design.Technical Report2068,INRIA;1993.

[102]Marco N,Dervieux A.Multilevel parameterization for aerodynamical optimization of3d shapes.Technical Report2949,INRIA;1996.

[103]Nielsen E,Anderson W.Aerodynamic design optimization on structured meshes using the Navier–Stokes equations.In:AIAA paper series,Paper98-4809;1998.

[104]Peter J,Drullion F,Pham C-T.Contribution to discrete implicit gradient and discrete adjoint method for aerodynamic shape optimization.In:Proceedings of ECCOMAS,Jyvaskyla;2004.

[105]Meaux M,Cormery M,Voizard G.Viscous aerodynamic shape optimization based on the discrete adjoint state for3d industrial con?gurations.In: Proceedings of ECCOMAS,Jyvaskyla;2004.

[106]Kuruvila G,Narducci R,Agrawal S.Development and application of TLNS3d-Adjoint:a practical tool for aerodynamic shape optimization.In:AIAA paper series,Paper2001-2400;2001.

[107]Mohammadi B.Dynamical approaches and incomplete gradients for shape optimization and?ow control.In:AIAA paper series,Paper99-3374;

1999.

[108]Medic G,Mohammadi B,Petruzzelli N,Stanciu M,Hecht F.3d optimal shape design for complex?ows:applications to turbomachinery.In:AIAA paper series,Paper99-0182;1999.

[109]Mohammadi B.Flow control and shape optimization in aeroelastic con?gurations.In:AIAA paper series,Paper99-0182;1999.

[110]Kim H-J,Obayashi S,Nakahashi K.Flap-de?ection optimization for transonic cruise performance improvement of supersonic transport wing.J Aircraft 2001;38(4):709–17.

[111]Kim C,Kim C,Rho O.Feasibility study of the constant eddy-viscosity assumption in gradient-based design optimization.J Aircraft2003;40: 1168–76.

[112]Newman III J,Taylor III A,Hou G-W,Jones H.An approximately factored incremental strategy for calculating consistent discrete aerodynamic sensitivity derivatives.J Comput Phys1994;113:336–46.

[113]Nemec M,Zingg D.A Newton–Krylov algorithm for aerodynamic design using the Navier–Stokes equations.AIAA J2002;40(6):1146–54.

[114]Giles M.On the iterative solution of the adjoint equations.In:Automatic differentiation:from simulation to optimization;2001.p.145–52.

[115]Christianson B.Reverse accumulation and attractive?xed points.Optim Methods Software1994;3(4):311–26.

[116]Nielsen E,Lu J,Park M,Darmofal D.An exact dual adjoint solution method for turbulent?ows on unstructured grids.In:AIAA paper series,Paper2003-0272;2003.

[117]Dwight R.Ef?ciency improvements of RANS-based analysis and optimization using implicit and adjoint methods on unstructured grids.Ph.D.thesis, School of Mathematics,University of Manchester;2006.

[118]Giles M.On the use of Runge–Kutta time-marching and multigrid for the solution of the steady adjoint equations.In:AD2000conference in Nice;

2000.

[119]Nielsen E,Park 0f9759e85ef7ba0d4a733b0eing an adjoint approach to eliminate mesh sensitivities in aerodynamic design.AIAA J2005;44(5):948–53.

[120]Nielsen E.Aerodynamic design sensitivities on an unstructured mesh using the Navier–Stokes equations and a discrete adjoint formulation.Ph.D.thesis, Virginia State University;1998.

[121]Campobasso M,Giles M.Stabilization of a linear?ow solver for turbomachinery aeroelasticity by means of the recursive projection method.AIAA J2004;42(9):1765–74.

[122]Campobasso M,Giles M.Stabilizing linear harmonic?ow solvers for turbomachinery aeroelasticity with complex iterative algorithms.AIAA J 2006;44(5):1048–59.

[123]Saad Y,Schultz MH.GMRES:a generalized minimum residual algorithm for solving non-symmetric linear systems.SIAM J Sci Stat Comput 1988;7(3):856–9.

[124]Schroff G,Keller H.Stabilization of unstable procedures:the recursive projection method.SIAM J Numer Anal1993;30(4):1099–120.

[125]Keller H.RPM:a remedy for instability.In:Estep D,Tavener S,editors.

Collected lectures on the preservation of stability under discretization,SIAM Proceedings in applied mathematics,vol.109;2002.p.185–96.

[126]Mani K,Mavriplis D.An unsteady discrete adjoint formulation for two-dimensional?ow problems with deforming meshes.In:AIAA paper series, Paper2007-0060;2007.

[127]Jakobsson S,Amoignon O.Mesh deformation using radial basis functions for gradient-based aerodynamic shape 0f9759e85ef7ba0d4a733b0eput Fluids2007;36: 1119–36.[128]Sherman L,Taylor III A,Green L,Newman A,Hou G,Korivi V.First-and second-order aerodynamic sensitivity derivatives via automatic differentiation with incremental iterative methods.J Comput Phys1996;129:307–31.

[129]Bischof C,Carle A,Corliss G,Griewank A,Hovland P.The ADIFOR2.0system for the automatic differentiation of Fortran77programs.Sci Program 1992;187:1–29.

[130]Taylor III A,Green L,Newman P,Putko M.Some advanced concepts in discrete sensitivity analysis.AIAA J2003;48(7):1224–9.

[131]Guillaume P,Masmoudi 0f9759e85ef7ba0d4a733b0eputation of high order derivatives in optimal shape design.J Numer Math1994;67:231–50.

[132]Aubert S,Tournier J,Rochette M,Blanche J,N’Diaye M,Melem M,et al.

Optimisation of a gas mixer using a new parametric solver.In:Proceedings of ECCOMAS CFD,Swansea,Wales;2001.

[133]Hasco?t L.TAPENADE:a tool for automatic differentiation of programs.In: Proceedings of4th European congress on computational methods,ECCOMAS 2004,Jyvaskyla,Finland;2004.

[134]Wakayama S.Lifting surfaces design using multidisciplinary optimization.

Ph.D.thesis,Stanford University,June;1994.

[135]Fazzolari A.An aero-structure adjoint formulation for ef?cient multidisciplinary wing optimization.Ph.D.thesis,Technical University of Braunschweig,September;2005.

[136]Sobieszczanski-Sobieski J.Sensitivity of complex internally coupled systems.

AIAA J1990;28(1):153–60.

[137]Haftka R.Structural optimization with aeroelastic constraints:a survey of U.S.applications.Int J Vehicle Des1986;7:381–92.

[138]Bowman K,Grandhi R,Eastep F.Aerodynamic design optimization using sensitivity analysis and computational?uid dynamics.J Struct Optim 1989;1:153–61.

[139]Friedmann P.Helicopter vibration reduction using structural optimization with aeroelastic/multidisciplinary constraints:a survey.AIAA J1991;28:8–21. [140]Barthelemy J-F,Wrenn G,Dovi A,Hall L.Supersonic transport wing minimum design integrating aerodynamics and structures.J Aircraft1994;31:330–8. [141]Giunta A,Sobieszczanski-Sobieski J.Progress toward using sensitivity derivatives in a high-?delity aeroelastic analysis of a supersonic transport.

In:AIAA paper series,Paper98-4763;1998.

[142]Ghattas O,Li X.Domain decomposition methods for the sensitivity analysis of

a nonlinear aeroelastic problem.Int J Comput Fluid Dyn1998;11:113–30. [143]Moller H,Lund E.Shape sensitivity analysis of strongly coupled?uid–structure interaction problems.In:AIAA paper series,Paper2000-4823;

2000.

[144]Maute K,Lesoinne N,Farhat C.Optimization of aeroelastic systems using coupled analytical sensitivities.In:AIAA paper series,Paper2000-0560;

2000.

[145]Hou G,Satyanarayana A.Analytical sensitivity analysis of a static aeroelastic wing.In:AIAA paper series,Paper2000-4824;2000.

[146]Maute K,Nikbay N,Farhat C.Coupled analytical sensitivity analysis and optimization of three-dimensional nonlinear aeroelastic systems.AIAA J 2001;39(11):2051–61.

[147]Martins J,Alonso J,Reuther J.High-?delity aero-structural design optimization of a supersonic business jet.In:AIAA paper series,Paper 2002-1483;2002.

[148]Farhat C,Lesoinne M,Maman N.Mixed explicit/implicit time integration of coupled aeroelastic problems:three-?eld formulation,geometric conservation and distributed solution.Int J Numer Methods Fluids 1995;21:807–35.

[149]Jameson A.Re-engineering the design process through computation.J Aircraft1999;36(1):36–50.

[150]Reuther J,Jameson A,Remlinger M,Saunders D.Constrained multi-point aerodynamic shape optimization using an adjoint formulation and parallel computers,Part I.J Aircraft1999;36(1):51–60.

[151]Reuther J,Jameson A,Remlinger M,Saunders D.Constrained multi-point aerodynamic shape optimization using an adjoint formulation and parallel computers,Part II.J Aircraft1999;36(1):61–74.

[152]Leoviriyakit K,Jameson A.Aero-structural wing planform optimization.In: AIAA paper series,Paper2004-0029;2004.

[153]Leoviriyakit K,Jameson A.Multi-point wing planform optimization via control theory.In:AIAA paper series,Paper2005-0450;2005.

[154]Jameson A,Leoviriyakit K,Shankaran S.Multi-point aero-structural optimization of wings including planform variations.In:AIAA paper series, Paper2007-764;2007.

[155]Campobasso M,Duta M,G MB.Adjoint calculation of sensitivities of turbomachinery objective functions.AIAA J Propulsion Power2003;19(4): 693–703.

[156]Duta M,Giles M,Campobasso M.The harmonic adjoint approach to unsteady turbomachinery design.Int J Numer Methods Fluids2002;40(3–4):323–32. [157]Gerhold T,Galle M,Friedrich O,Evans J.Calculation of complex3d con?gurations employing the DLR TAU-Code.In:AIAA paper series,Paper 97-0167;1997.

[158]Mouton S,Laurenceau J,Carrier G.Aerodynamic and structural optimization of powerplant integration under the wing of a transonic transport aircraft.In: Proceedings of42th AAAF symposium on applied aerodynamics,Nice;2007. [159]Salah El Din I,Carrier G,Mouton S.Discrete adjoint method in elsA(Part2): application to aerodynamic design optimisation.In:Proceedings of the ONERA-DLR aerospace symposium(ODAS),Toulouse;2006.

J.E.V.Peter,R.P.Dwight/Computers&Fluids39(2010)373–391391

本文来源:https://www.bwwdw.com/article/qxje.html

Top