Adjoint-based fluid dynamic design optimization in quasi-periodic unsteady flow problems using a harmonic balance method

Adjoint-based fluid dynamic design optimization in quasi-periodic unsteady flow problems

using a harmonic balance method

Rubino, A.; Pini, M.; Colonna, P.; Albring, T.; Nimmagadda, S.; Economon, T.; Alonso, J.

DOI

10.1016/j.jcp.2018.06.023

Publication date

2018

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Rubino, A., Pini, M., Colonna, P., Albring, T., Nimmagadda, S., Economon, T., & Alonso, J. (2018).

Adjoint-based fluid dynamic design optimization in quasi-periodic unsteady flow problems using a harmonic balance

method. Journal of Computational Physics, 372, 220-235. https://doi.org/10.1016/j.jcp.2018.06.023

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

Adjoint-based

fluid

dynamic

design

optimization

in

quasi-periodic

unsteady

flow

problems

using

a

harmonic

balance

method

A. Rubino

a

,

M. Pini

a

,

∗

,

P. Colonna

a

,

T. Albring

b

,

S. Nimmagadda

c

,

T. Economon

d

,

J. Alonso

c

a_{Propulsion&Power,DelftUniversityofTechnology,Kluyverweg1,2629HS,Delft,theNetherlands} b_{ChairforScientificComputing,TUKaiserslautern,67663Kaiserslautern,Germany}

c_{DepartmentofAeronautics&Astronautics,StanfordUniversity,Stanford,CA94305,USA}

d_{MultiphysicsModelingandSimulation,BoschResearchandTechnologyCenter,Sunnyvale,CA94085,USA}

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received5December2017

Receivedinrevisedform23May2018 Accepted6June2018

Availableonline15June2018

Keywords:

Unsteadyoptimization RANSadjoint Harmonicbalance Quasi-periodic

Shapeoptimizationinunsteadyflowproblemsenablestheconsiderationofdynamiceffects ondesign.Theabilitytotreatunsteadyeffectsisattractive,asitcanprovideperformance gains when compared to steady-state design methods for a variety of applications in whichtime-varying flows are of paramount importance.This is the case, for example, in turbomachinery orrotorcraft design. Given the high computational cost involved in time-accuratedesign problems, adjoint-based shape optimizationis apromising option. However,efficientsensitivityanalysisshouldalsobeaccompaniedbyasignificantdecrease incomputational costfor the primalflow solution,as well.Reduced-order models,like those based on the harmonic balance concept, in combination with the calculation of gradientsvia adjointmethods, are proposed for the efficient solution ofacertainclass ofaerodynamicsoptimizationproblems.Theharmonicbalancemethodisapplicableifthe flowischaracterizedbydiscretefinitedominantflowfrequenciesthatdonotneedtobe integermultiplesofafundamentalharmonic.Afully-turbulentharmonicbalancediscrete adjointformulation basedonaduality-preserving approachis proposed.The methodis implementedbyleveraging algorithmic differentiationand is applied totwo test cases: theconstrainedshapeoptimizationofbothapitchingairfoilandaturbinecascade.Akey advantageofthecurrentapproachistheaccuratecomputationofgradientsascompared tosecondorderfinite differenceswithoutanyapproximation inthe linearizationofthe turbulentviscosity.Theshapeoptimizationresultsshowsignificantimprovementsforthe selectedtime-dependentobjectivefunctions,demonstratingthatdesignproblemsinvolving almost-periodicunsteadyflowscanbetackledwithmanageablecomputationaleffort.

©2018TheAuthor(s).PublishedbyElsevierInc.Thisisanopenaccessarticleunderthe CCBYlicense(http://creativecommons.org/licenses/by/4.0/).

1. Introduction

The advancement of computational resources has enabled the application of CFD-based design methods to complex

shapes,discretizedonlargedomains,oftenincombinationwithhighfidelitymodels [1,2].

*

Correspondingauthor.

E-mailaddress:m.pini@tudelft.nl(M. Pini).

https://doi.org/10.1016/j.jcp.2018.06.023

0021-9991/©2018TheAuthor(s).PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

Optimization methods fordesign purposes have significantly improved, offering the possibility to deal withcomplex problems at a reduced computational cost [3–5]. In particular, adjoint-based optimization methods [6,7] provide a very efficientapproachforcomputingdesignsensitivities irrespectiveofthenumberofdesignvariables.Inapplications where thenumberofdesignvariablesisconsiderablygreaterthan thenumberofobjectivefunctions, adjointmethodsallowthe computation ofoptimalsolutions inthemostcost-effectiveway,making themwell suitedforcomplexindustrial applica-tions [8–10].

Todate,mostoftheworkonadjointmethodshasbeenbasedontheassumptionofsteadyflow. Obtainingtheadjoint

solution of an unsteady flow problem can pose a challenge because of the associated large computational andmemory

requirements [11].However,accountingfortime-dependentflowphenomenaintheoptimizationprocessisoftenessential in applicationscharacterized by intrinsicunsteady effects, such asthe aerodynamic design ofrotorcraft, turbomachinery, open rotors, andwind turbines,toname afew.Furthermore, unsteadyadjoint-basedoptimization methodscan pavethe waytothesolutionofmultidisciplinaryproblemscharacterizedbytime-dependentphenomena suchasthoseencountered inaeroelasticityornoisereduction [12],forexample.

Giventhecomputationalcostofaccuratelyobtaining anunsteadysolutionandthusitsadjoint,sufficientlyaccurate re-ducedordermethods [13] mustbeused, iftheobjectiveistosolvedesignproblemsroutinely.Amongothers,theHarmonic Balancemethod(HB)isapromisingoptionforapplicationsinvolvingquasi-periodicflowscharacterizedbyafiniteamount ofdominantfrequenciesthatneednotbeharmonicsofeachother [14–16].

ConsiderabledevelopmenthasbeendedicatedtoHBmethodsforturbomachineryapplications,wherebytheflow spec-trumisdominatedbythebladepassingfrequencies [15].Onetechniquethatisasubsetoftheharmonicbalancemethod, sometimes referred toasthe TimeSpectral method [14,17,15], hasbeen originallyformulated forperiodic flow problems, andassociatedadjointshavebeenderived [18–21].

HB adjoint-basedoptimizationopens up thepossibility todeal withproblems featuringflow frequencies that are not harmonically related, thus without the restriction of resolving harmonics of a single fundamental flow frequency. Cur-rentlyavailable methodsin the literature are limitedto inviscid flow problemsor frozenturbulence assumptions during design [22,23], which can havea strongimpact on the final optimizationresult [24]. In thiswork, a fully-turbulentHB

adjoint-based shape optimization method is proposed and implemented within the SU2 open-source software

environ-ment [25,26].Thealgorithmisbasedontheduality-preservingapproach [27–30],whichallowstheadjointsolvertoinherit thesameconvergencepropertiesoftheprimalflowsolver.Thechosenalgorithmensuresrobustconvergenceofthe numer-icalsolutionoftheturbulentadjointequations.

Themethodisdescribedindetailfirst,followedbytheillustrationoftwoapplicationcases.Anairfoilpitchingatarate characterizedbytwofrequenciesthatarenotharmonicallyrelatedisfirstconsidered,followedbyaturbinecascadesubject tounsteadyinletconditions.

2. Method

2.1. Timediscretization

Thesemi-discreteformoftheNavier–Stokesequations,foragenericcellvolume

,

is



∂

U

∂t

+

R

(

U

)

=

0

,t

>

0

.

(1)

U

= (

ρ

,

ρ

v1

,

ρ

v2

,

ρ

v3

,

ρ

E

)

isthevectorofconservativevariables,with

ρ

thedensity,v thevelocityvector,andE thetotal specificenergy.



anditsboundary

∂

areassumedtovarytheirpositionintime,withvelocityu,withoutdeforming.

R

istheresidualoperatorforthespatial integrationoftheconvective andviscousfluxes, i.e., Fc and Fv.UsinganArbitrary Lagrangian–Eulerian(ALE)formulation [31],

R

(

U

)

=

f

(

Fc

,

Fv

)

in

 ,t

>

0 v

=

u on

∂

 ,t

>

0 (2) with Fc

=

⎛

⎝

ρ

v

× (

ρ

(

vv

−

−

uu

)

)

+

pI

¯

ρ

E(v

−

u

)

+

p v

⎞

⎠ ,

(3)

wherep isthestaticpressure.Theviscousfluxesaregivenby

Fv

=

⎛

⎝

μ

·

τ

¯

μ

τ

¯

·

v

+

κ∇

T

⎞

⎠ ,

(4)

¯

τ

= ∇

v

+ ∇

v

−

2

3I

¯

(

∇ ·

v

) .

(5)

The turbulencemodelingisconsidered, accordingto theBoussinesqhypothesis, bydefining

μ

=

μ

l

+

μ

t and

κ

=

κ

l

+

κ

t.

μ

l and

μ

t are the laminar andturbulent dynamic viscosities, whereas

κ

l and

κ

t are thelaminar andturbulent thermal

conductivities.

Theapplicationoftime-discretizationto(1),usinganimplicitEulerscheme,leadsto



D

t

(

Uq+1

)

+

R

(

Uq+1

)

=

0

,

(6)

whereq isthephysicaltimestepindex,and

D

t isthetime-derivativeoperator.Foradualtime steppingapproach [32,33]

withpseudo-time

τ

,onewouldobtainthefollowingdiscretization





U

q+1



τ

+ 

D

t

(

U

q+1

₎

₊

_R

₍

_Uq+1

₎

₌

₀

_,

₍₇₎

wheretheadditionaltermisusedtorelaxthesolutionateachphysicaltimestep.

2.2. Harmonicbalanceoperator

AsgiveninRef. [15],theFouriercoefficientsresultingfromtheapplicationofthediscreteFouriertransform(DFT)are

ˆ

uk

=

1 N N

_

−1 n=0 Une−iωktn

,

(8)

whereU

˜

= [

U0

,

U1

,

...,

UN−1

]

,isthevectoroftheconservativevariablesevaluatedatN timeinstances

t

= [

t0

,t

1

, ...,

tN−1

] .

(9)

HencethecorrespondingFouriercoefficientsare

ˆ

u

= [

u0

,

u1

, ...,

uK−1

] ,

(10)

with N

=

2K

+

1,

ω

k

=

2

π

kfk and K being the number of frequencies. An odd number of time instances is used in

this work in order to prevent numerical instabilities [34]. The ensemble of the K resolved frequencies is denoted by

ω

= [

0,

ω

1

,

...,

ω

K

,

ω

−K

,

...,

ω

−1

]

.BydefiningtheDFTmatrixas

Ek,n

=

1 Ne

−iωktn

_,

_n,k

∈ [

₀

_,

_N

] ,

₍₁₁₎

anditsinversediscreteFouriertransform(IDFT)

E−_n,1_k

=

eiωktn

_,

_n,k

∈ [

₀

_,

_N

] ,

₍₁₂₎

onecancalculatetheFouriercoefficientsas

ˆ

u

=

EU

˜

,

(13)

withthecorrespondingvectorofconservativevariables

˜

U

=

E−1u

ˆ

.

(14)

Ifthefrequencies fk (andhence

ω

k)arenot multiplesof f1,oncetheDFT matrixfrom(11) ortheIDFTmatrixfrom(12) areconstructed,itisnotpossibletoobtainananalyticalexpressionforthecorrespondinginversematrix.

The timeoperatorof(6) canbeapproximatedusingspectralinterpolation.Applyingthespectraloperatortothevector ofconservativevariablesU ,

˜

evaluatedatN timeinstances,yields

D

t

(

U

)

≈

D

t

( ˜

U

).

(15)

Giventhatu is

ˆ

independentfromtime,using(14) and(13),onecanwrite

D

t

( ˜

U

)

=

D

t

(

E−1u

ˆ

)

=

∂

E−1

∂t

u

ˆ

=

∂

E−1

∂t

EU

˜

.

(16) From(12)

∂

E−1

∂t

=

E −1_D

_,

₍₁₇₎

where

Dk,n

=

i

ω

k

δ

k,n

.

(18)

D isthediagonalmatrixgivenby

D

=

diag

(

0

,

i

ω

1

, ...,

i

ω

K

,

i

ω

−K

, ...,

i

ω

−1

) .

(19)

Onecancombine(16) with(17) anddefinethespectraloperatormatrixH as

H

=

E−1D E

.

(20)

Here,E−1isgivenanalyticallyby(12),andE iscomputedbyinvertingE−1usingGaussianelimination.Iteventuallyfollows that

D

t

( ˜

U

)

=

HU

˜

.

(21)

2.3. Time-domainharmonicbalance

ConsideringU ,

˜

thesetoftheconservativevariables evaluatedat N timeinstances,onecan write(7) forasingle time instancen as





U q+1 n



τ

+ 

D

t

(

U q+1 n

)

+

R

(

Uqn+1

)

=

0

.

(22)

Linearizationoftheresidualyields

R

(

Uqn+1

)

=

R

(

Uqn

)

+

∂

R

(

Uqn

)

∂

U_nq



Un

=

R

(

U

q

n

)

+

J



Un

.

(23)

From(21),

D

t isalinearoperator,sothefollowingmanipulationsarepossible:

D

t

(

Un

)

=

D

t

(

Unq+1

−

Uqn

)

=

D

t

(

Uqn+1

)

−

D

t

(

Uqn

) ,

(24) and

D

t

(

Uqn+1

)

=

N−1



k=0 Hn,k



Uk

+

N−1



k=0 Hn,kUq_k

.

(25)

Substitutionofthelinearizedexpressionstransforms(22) into





I



τ

+

J

+ 

Hn,n





Un

+

R

(

Uqn

)

= −

N

_

−1 k=0

(

1

− δ

n,k

)H

n,k



Uk

−

N

_

−1 k=0 Hn,kUqk

.

(26)

Inthiswork,asemi-implicitapproachisusedtosolve(26) as





I



τ

+

J





Un

= ˜

R

n

( ˆ

U q

) ,

(27) where

˜

R

n

(

Uq

)

= −

R

(

Uqn

)

−

N−1



k=0 Hn,k



Uk

−

N−1



k=0 Hn,kUq_k

.

(28)

Equation(27) is solvedforeachtimeinstanceinasegregatedmanner.Therefore,anunsteadyflowproblemcharacterized by Kfrequencies requiresthat the solutionof2K

+

1 nonlinearsystems ofequationsneed tobe computed. The current harmonicbalanceapproachisimplementedintheopen-sourceSU2code [25,26].

2.4. Governingequationsoftheadjointsolver

The adjoint equationsare now derived fortheproposed Harmonic Balance formulation described inSec. 2.3for both laminarandturbulentflows.

Equation (27) isreformulatedintermsofafixed-pointiterationforeachUn as

Uq_n+1

=

G

n

(

Uq

) ,

(29)

where

G

n isthe iteration operator of the pseudo time-stepping attime instance n. If

G

n is contractive, i.e.,

||

∂_∂G_Un_n

||

<

1,

accordingtotheBanachfixed-pointtheorem [35],(29) admitsauniquefixed-pointsolutionU_n∗suchthat

˜

R

n

(

U∗

)

=

0

⇐⇒

U∗n

=

G

n

(

U∗

) .

(30)

Theaerodynamicdesignproblem,includingapossibleexplicitdependenceoftheobjectivefunction

J

onthevectorofthe designvariables

α

,canbeexpressedas

minimize

α

J

(

U

(

α

),

X

(

α

))

subject to Un

(

α

)

=

G

n

(

U

(

α

),

Xn

(

α

)),

n

=

1

,

2

, ...N

Xn

(

α

)

=

Mn

(

α

).

(31)

Xn arethephysicalgridsconstructedforeachtimeinstance,andMn

(

α

)

isadifferentiablefunctionrepresentingthemesh

deformation algorithm.The objective function

J

isobtained asthe spectral average (using (42)) over theresolved time instances

J

=

f

(

J

(

U1

,

X1

),

J

(

U2

,

X2

), ...,

J

(

UN

,

XN

)) .

(32)

TheLagrangianoftheconstrainedoptimizationproblemisgivenby

L

=

J

+

N



n=1

λ

n

(

G

n

(

U

(

α

),

Xn

(

α

))

−

Un

(

α

))



+

μ

n

(M

n

(

α

)

−

Xn

(

α

))



,

(33)

with

λ

and

μ

beingtheadjointvariables.Giventheconstraintequationsin(31)

Un

(

α

)

−

G

n

(

U

(

α

),

Xn

(

α

))

=

0

,

n

=

0

,

1

, ...,

N

−

1

Xn

(

α

)

−

Mn

(

α

)

=

0

,

(34)

andomittinginthenotationtheexplicitdependencefromtheindependentvariables,onecanexpressthedifferentialofthe Lagrangianas d

L

=

N

_

−1 n=0



∂

J

∂

Un 

+

N

_

−1 k=0

λ

k

∂

G

k

∂

Un 

− λ

n



dUn

+

N

_

−1 n=0



∂

J

∂

Xn 

+ λ

n

∂

G

n

∂

Xn 

−

μ

n



d Xn

+

N

_

−1 n=0

μ

n

∂

Mn

∂

α

d

α

,

(35) fromwhichtheadjointequationscanbeobtainedas

∂

J

∂

Un 

+

N

_

−1 k=0

λ

k

∂

G

k

∂

Un 

= λ

n

,

(36) and

∂

J

∂

Xn 

+ λ

n

∂

G

n

∂

Xn 

=

μ

n

.

(37)

Equation (37) canbesolved directlyoncethe solutionof(36) isknown.Similarto theflowsolver, (36) canbe seenasa fixed-pointiterationin

λ

n,namely

λ

qn+1

=

∂

N

∂

Un

(

U_n∗

, λ

q

,

Xn

) ,

(38)

whereU

ˆ

∗n isthenumericalsolutionfortheflowequation(30) and

N

istheshiftedLagrangiandefinedas

N

=

J

+

N



n=1

Table 1

SimulationparametersofthepitchingNACA64A010airfoiltestcase.Thereducedfrequencyω∗

isbasedonthesemi-chordlength.

Symbol Value Units

Free-stream temperature T∞ 288.15 [K]

Mach number Ma∞ 0.78 [–]

Reduced frequencies [ω∗1,ω∗2] [0.197,0.512] [–] Since

G

n iscontractive, _∂∂_UN_n isalsocontractivebecause

_∂λ

∂

_n



_{∂ ˆ}

∂

_U

N

n



=

∂

G

n

∂ ˆ

Un 

=

∂

G

n

∂ ˆ

Un

<

1

.

(40)

Therefore,accordingtothefixed-pointtheorem,(38) willconvergeatthesamerateastheprimalflowsolver.

Therighthandsideofequation(38) isobtainedusingAlgorithmicDifferentiationappliedtotheunderlyingsourcecode oftheprogramthatcomputes

G

n.TheADtooladopted [36] makesuseoftheJacobitapingmethod incombinationwiththe

ExpressionTemplates featureofC++,leadingonlytoasmallruntimeoverhead.

Finally,thegradientoftheobjectivefunction

J

withrespecttothevectorofthedesignvariables

α

canbe computed fromtheconvergedflowandadjointsolutionsusing

d

L

 d

α

=

d

J

 d

α

=

μ

n

∂

Mn

(

α

)

∂

α

.

(41) 3. Results

TheproposedHB-basedadjointmethodhasbeenappliedtotwotestcases:thefluiddynamicshapeoptimizationofboth apitchingairfoilandaturbinecascade.Forbothcases,aharmonicbalanceflowsolutionhasbeenfirstobtainedandverified against a fully-unsteadysimulation using a second-order dual time stepping method [33]. For the spatial discretization, second-order accuracyhasbeenobtainedforthe convectivefluxesusinga centeredscheme ortheMUSCL approach [37] andfortheviscous fluxesusingacorrected average-of-gradientsmethod.Thereaderisreferredto Ref. [26] fordetailson thenumericalmethodsandtoolsavailableinSU2.

The resultsobtainedwiththeHBsolver at N timeinstances forageneric quantity ofinterest

are interpolatedtoa largertimevectort∗ oflengthN∗using

∗

=

E∗−1

(

E

) ,

(42)

whereE∗−1 isthelargerinterpolatedIDFTmatrixofsizeN∗

×

N givenby E∗−_n,1_k

=

eiωkt∗n_,andE−1 _istheN

×

_{N IDFTmatrix.}

Theconstrainedshapeoptimizationproblemissolved oncethegradientisobtainedfromtheadjointsolution,usinga modified version ofthenonlinearleast-squares method(SLSQP)[38]. Accurategradientvaluesare obtainedbyleveraging Algorithmic Differentiation (AD). The reversemode ofthe open-source AD tool CoDiPack [29,39] isused in thiswork to linearizetheprimalsolver.

3.1. NACA64A010pitchingairfoil

Atwo-dimensionalpitching NACA64A010airfoilininviscidflow isconsideredfirst.The rigidbody motionisimposed byassigningatime-varyingangleofattackcomputedas

α

=

1

.

01

[

sin

(

ω

t

)

+

sin

(

2

.

6

ω

t)

] .

(43)

The flow is transonic withshocks appearing along theupper andlower surfaces ofthe airfoilwhile pitching, ascan be observed in Figs. 5a, 5b and 5c. The convective fluxes are computed with the JSTscheme [32]. The mesh is composed entirely of triangular elements. It contains approximately 11 000 grid points with 200 points along the airfoil and 100 pointsonthefar-fieldboundary.ThemainsimulationparametersaresummarizedinTable1.

Areferencetime-accuratesimulationwasperformedwith160timestepspersmallestperiod(correspondingto2.6

ω

1) in order to get a well-resolved solution in time. Fig. 1a and 1b show the mesh close to the airfoiland the frequency spectrumofthedrag coefficient.Fig.1bshowsthattheoveralldominantfrequenciesare thepitching frequencies

ω

,2.6

ω

andtheirlinearcombinations.ThisanalysisisperformedtoselecttherelevantfrequenciesfortheHBcomputation. Thetestcaseisthensimulatedusingtheharmonicbalancemethodwithseveralchoicesforthenumberoftimeinstances withthetwoinput frequenciesaspitchingfrequencies.Fig. 2showstheconvergenceoftheharmonicbalancesolutionto the fully-unsteadysolution withan increase inthenumber ofinput frequencies.The convergedsolution fortheselected time instances ofboth thelift andthedrag coefficient is interpolatedusing(42). Accordingto thespectrum ofthe drag coefficient(Fig.1b),theselectedinputfrequenciesarereportedinTable2.

Fig. 1. NACA 64A010 airfoil mesh at t=0 and spectrum of the drag coefficient.

Fig. 2. Liftcoefficientclanddragcoefficientcd:HBsolutionobtainedfordifferenttimeinstancesvstimeaccuratesolution.Trefisthereferencetime intervalcorrespondingtothehighestfrequency.

Table 2

OptimaltimeperiodratioTopt/T0andinputfrequenciesforthepitchingairfoiltestcase,with ω2=2.6ω1.T0istheperiodcorrespondingtothelowestfrequencyvalue.

Time instances Topt/T0 Input frequencies 5 1.13 0,±ω1,±ω2 7 1.00 0,±ω1,±ω2,±2ω1

9 1.41 0,±ω1,±(ω2−ω1),±ω2,±(ω2+ω1) 11 1.38 0,±ω1,±(ω2−ω1),±2ω1,±ω2,±(ω2+ω1)

Inordertoensureconvergenceforanysetoffrequencies,anoptimaltimepseudo-periodisselectedusingthealgorithm proposed inRef. [40].Auniformtimesamplingwithinthepseudo-periodT isadopted,andthevaluesoftheoptimaltime periodarereportedinTable2.

Overall,theagreementoftheliftcoefficientisexcellent,resultingina root-mean-square-error(RMSE)of0.00275 with justtwoinputfrequencies,duetothefactthatthemovementoftheairfoilisprimarilyinthedirectionofthelift(Fig. 2a). Conversely, inthe caseofdrag, alarger numberofinput frequencies isrequiredto capturethefully-unsteadyresultdue to thedominantnon-linearitiesin theflow(shocks).Hence, morethanfivetime instancesare necessaryfortheaccurate determinationofthedragcoefficient,ascanbeseeninFig.2b.Sincetheselectedobjectivefunctionisthedragcoefficient,a suitablenumberofinputfrequenciesshouldbeidentified.Thetime-averageddragcoefficientis0.0027for5timeinstances, 0.0028for9,11and13timeinstancescomparedto0.0028forthetime-accurateunsteadysolution.Itcanbeinferredthat 9timeinstancesaresufficientforshapeoptimization.Fig.2reportsalsothesteady-statevaluesofCl andCd obtainedwith

α

=

0.

Aftertheanalysisoftheflowsolution,theadjoint-basedshapeoptimizationproblemisconsidered.Thetime-averaged drag coefficient cd isminimizedover thepseudo-period.Equalityconstraintsonthe liftcoefficient (cl)andthemaximum

Fig. 3. ComparisonbetweenthenormalizedgradientsoftheobjectivefunctioncalculatedwiththereversemodeoftheADandthesecondordercentral finitedifference (a);flowandadjointsolverconvergencehistoryofthedensityresidualrelativetothefirsttimeinstance (b).

Fig. 4. Shape optimization convergence history of the drag coefficient (a); baseline vs optimized airfoil (b). minimize α cd

(

Un

,

Xn

,

α

)

subject to: cl

=

cl0

,

δ

max

= δ

max0

,

Un

=

G

n

,

n

=

1

,

2

, ...N

Xn

=

Mn

.

(44)

50 Hicks–Hennebump function variables [41] distributed uniformlyover the upperandlower surfaces ofthe airfoilare chosenasthedesignvariables.

The gradients ofthe discrete adjoint are first comparedwith thegradients obtained witha second-order accurate fi-nite difference (FD) method.This verificationis performedovera set of8 Hicks–Henne bumpfunctions withthefirst 4 locatedonthesuction sideoftheairfoilandtheothersonthepressureside.Fig.3adepictsthecomparisonbetweenthe adjoint-basedgradients andthe FDgradients. Fig.3bshowstheconvergenceofboththe flowandtheadjointsolverasa functionofthenumberofiterations.Asdiscussed inSec. 2.4,theconvergencerateoftheadjointsolverisinheritedfrom theprimal flowsolver.The computationaltime ratiobetweentheadjointsolutionandtheflowsolution isabout1.3. The averagecomputationaltime forone iterationoftheprimalsolverwas approximately0.28 seconds,usingthe4coresIntel XeonE5-1620Processor withhyper-threading.

Fig.4a reportsthereduction indrag coefficient withthe numberofoptimizeriterations aswell asthe liftcoefficient constraint.The shapeoptimizationprocess (Fig.4)achievesadecreaseindragcoefficientofabout50%, whilemaintaining theliftcoefficientandthemaximumthicknesswithin1% oftheconstraintvalue.

TheMachcontoursofthebaselineairfoilcomparedwiththecontoursoftheoptimizedairfoilatdifferenttimeinstances arereportedinFig.5.The optimizedairfoilshapeleads toanattenuationofthestrongshockscharacterizingthebaseline configuration,resultingindragreduction.

Finally, the unsteady optimization results are compared with those obtained froma similar constrained steady state optimizationoftheairfoilat

α

=

0.The finalshapegivenbytheunsteadydesignmethoddiffersfromthesteadyone and

Fig. 5. PitchingairfoilMachnumbercontourscalculatedatthreedifferenttimeinstanceswiththeHBmethod,forboththebaseline(a),(b),(c)andthe optimized(d),(e),(f)profile (forinterpretationofthereferencestocolorinthisfigure,thereaderisreferredtothewebversionofthisarticle).

Fig. 6. Schematic geometry of the T106D-EIZ turbine cascade [42] and blade mesh.

itisshowninFig.4b.Furthermore,whenthesteady-optimizedairfoilundergoesthepitchingmotionprescribedbythetest case,thereductionofthetime-averageddragcoefficientis44%,whereasitwasreducedbymore,i.e.50%,fortheunsteady optimization.Thisshowsthattheinclusionofunsteadinessisworthwhileinthiscase.

3.2. T106D-EIZturbinecascade

The aim ofthis test case is to assess the capabilityof the HB-based design method for fully-turbulentflows. Inthe experimental setup [42], the unsteadiness in the wake ofan upstream blade rowis approximated by a moving bar, as depictedinFig.6a.Themoving barsare locatedatxb

/

l

=

0.7 upstreamofthecascadeinlet plane,havingavelocity vb

=

21.4 m/s paralleltotheinlet.AschematicrepresentationisshowninFig.6a,andthemainoperatingconditionsforthetest casearereportedinTable3.

The two-dimensionalflow domainis discretizedwithapproximately40 000elements (Fig. 6b), andtheRoe schemeis selectedfortheconvectivefluxdiscretization.Asuitablespacingofquadrilateralelements isusedto clusterthenearwall cells such that y+ islessthan1.Thistest caseisa benchmarkforthestudyoflaminar-to-turbulent transition.However,

Table 3

SimulationparametersoftheT106D-EIZtestcase,asdescribedin[42].Re2istheReynolds numberbasedontheexitvelocityand density;Ma2theexitMachnumber;T u1 theinlet turbulenceintensity.

Symbol Value Units

Exit Reynolds number Re2 200 000 [–]

Exit Mach number Ma2 0.593 [–]

Bars speed vb 21.4 [m/s]

Background turbulence level T u1 2.5 [%] Eddy viscosity ratio μt/μ 100 [–]

Fig. 7. Validation of steady state simulation results with experimental data [42].

since the present work aims to assess the methodology for design optimization only, the computations are performed assumingfully-turbulentconditions,employingtheSST turbulencemodel [43].

Inordertocalculatethecascadeperformance,thetotalpressurelosscoefficientisevaluatedas

ζ

P

=



Ptot,1



∂1

−



Ptot,2



∂2



Ptot,2



∂2

−

P2

∂

2

,

(45) where



Ptot,1



∂1 and



Ptot,2



∂2 aretheinletandoutlettotalpressureaveragedovertheircorrespondingboundary, respec-tively.

P2

∂2 is theaverage staticpressureatthe cascadeoutlet. Theaverages atthe boundariesarecalculated usinga mixed-outaveragingprocedure [44].

First,avalidationisperformedusingtheexperimentaldata [42] oftheturbinecascadeoperatingatsteadystate(Fig.7). Thesimulationresultsshowaverygoodagreementwiththeexperimentaldata.Asexpected,themaindeviationbetween CFDandexperimentsoccursataboutx

/

l

=

0.75,possiblyduetotransitionnotbeingaccuratelymodeled.

Forthistestcase,twoconfigurationsareconsideredfordesign:i)aspatiallynon-uniform,time-dependentinletboundary condition;ii)aspatially uniform, time-dependent inletboundarycondition. Theterminology OptC1and OptC2 isusedto refertothefirstandthesecondshapeoptimizationproblem,respectively.

3.2.1. OptC1 configuration

Inthisconfiguration,an inletboundaryconditionisimposedinordertoreproducethewakesgeneratedbythemoving bars.Theimposed valuesofthetotalpressure,temperature, andflowdirectionattheboundaryareinterpolatedfromthe results ofa steadystate simulation ofthe flow pastthe bars. Withthisboundary condition, onlymultiples of theblade passing frequencyareexpected. The fundamentalblade passingfrequency, givenaratiobetweenthe bladepitchandthe barpitchyp

/

yb

=

3,isdefinedas

f1

=

3 vb

yp

.

(46)

ToverifytheHBsolution,asecond-ordertime-accurateURANSsimulationusingthedualtimesteppingmethodis per-formedwithatime-step150x smallerthanthelowestperiod(1/f1).Thetotalpressurelosscoefficientfromthissimulation iscompared inFig.8a withthe HBsolutionobtainedwith3,5,and7 timeinstances. Theselected timeinstances corre-spondtothesolutionforthefrequencyvectors

ω

N3

= [

0,

±

ω

0

]

,

ω

N5

= [

0,

±

ω

0

,

±

2

ω

0

]

and

ω

N7

= [

0,

±

ω

0

,

±

2

ω

0

,

±

3

ω

0

]

. Theresolvedfrequencies are,therefore,multiplesofthefundamentalbladepassingfrequencyonly.Thetotalpressureloss coefficient,definedin (45),andshowninFig.8aasfunctionoftime,isobtainedbyspectralinterpolationoftheharmonic balanceresultsusing (42).TheRMSEofthetotalpressureloss coefficientforthesolutionobtainedwith5 timeinstances isequalto0.010.Theharmonicbalancesolutionobtainedwith5timeinstancesisabout9x fasterthanthetime-accurate

Fig. 8. Comparisonofthetotalpressurelosscoefficient

ζ

Pasfunctionoftime,obtainedwithtimeaccurateandHBsimulations(a);convergencehistoryof

thefirsttimeinstancedensityresidual(b).

Fig. 9. Verificationbetweenthenormalizedgradientsoftheobjectivefunction

ζ

pwithrespecttothedesignvariablesαcalculatedwiththeadjointmode

oftheADandsecondordercentraldifferencefinitedifference(FD).

solutioncalculatedoveratotalsimulationtimeoffiveperiods,whichincludestheinitialtransientbeforereaching conver-gence toaperiodic flowfield solution.5timeinstances areusedforshape optimization,asatrade-off betweenaccuracy andcomputationalcost.

InFigs.11a,11b,and11c,theMachnumbercontoursfromtheHBsimulationarereportedfor3differenttimeinstances withthesimulationperiodgivenby T

=

1/f1.Theresultsshowthebarwakesenteringthecascadeandaseparationarea occurringataboutx

/

l

=

0.7.

Next,theshapeoptimizationproblemofthecascadeconfigurationisconsidered.Itcanbeexpressedas

minimize

α

ζ

P

(

Un

,

Xn

,

α

)

subject to:

α

out

<

α

out0

+

4◦

,

δ

t

= δ

t0

,

Un

=

G

n

,

n

=

1

,

2

, ...N

Xn

=

Mn

,

(47)

wherethetime-averagedtotalpressurelosscoefficient

ζ

P,obtainedfrom (42),isselectedasobjectivefunction.Inequality

constraintsontheabsoluteexitflowangle(

α

out)andtrailingedgethickness(δt)areimposed.Theoptimizationisperformed

using an ensemble of 16 geometrical design parameters

α

based on a free-form deformation (FFD) approach [45]. The gradients ofthe objectivefunctionare againobtainedwiththeproposed adjointtechnique andcomparedwiththe same gradients obtained bysecond-order central finite differences(FD). Theresults ofthiscomparison arereported inFig. 9a, showingexcellent agreementbetweenADandFDgradients(RMSE

=

2

·

10−5_{).Theratiobetweenthecomputationaltime} oftheadjointsolutionandtheprimal flowsolutionisapproximately1.7.TheaverageprimalsolverCPUtimeperiteration wasabout1.41 secondsona4coresIntelXeonE5-1620Processor withhyper-threading.

Fig.10showsthattheconvergenceoftheoptimizationtotheminimumobjectiveisnearlyreachedafteronly7 evalua-tions,althoughsatisfyingtheconstraintrequiresmoreevaluations.Fig.10bhighlightsthattheperformanceoftheoptimized

Fig. 10. Shape optimization history of the total pressure loss coefficient and comparison between baseline and optimized blade profile (OptC1).

Fig. 11. MachnumbercontourscalculatedatthreedifferenttimeinstanceswiththeHBmethod,basedonthe OptC1testcase,forboththebaseline(a), (b),(c)andtheoptimized(d),(e),(f)bladeprofile.

bladeissignificantlyimproved,asthetotalpressurelosscoefficientisapproximately38% lower,whiletheconstraintonthe absoluteoutletflow angleissatisfied.The separationarea,asseen inFig.11,is considerablysmallerwiththeoptimized bladeshape.Theunsteady optimizationleads toadecreaseinthepeakofthetotalpressurelosscoefficient of44% anda reductionof54% of thesignalamplitude (Fig.12a). Furthermore,theobjectivefunction spectrumobtainedfromaURANS simulation of the optimized blade (Fig. 12a) doesnot contain additional frequencies when compared with the baseline configuration.

3.2.2. OptC2 configuration

InordertoinvestigatethecapabilitiesoftheHB-baseddesignmethodtodealwithproblemscharacterizedbyfrequencies that are not multiples of one fundamental harmonic, a second configuration of the T106D-EIZ cascade is considered. In analogywithpreviouswork [46],atimefluctuatinginlettotalpressureisprescribedas

˜

Fig. 12. Total pressure loss coefficient evolution in time calculated with URANS simulation for both the baseline and the optimized configuration.

Fig. 13. Comparisonofthetotalpressurelosscoefficient

ζ

Pasfunctionoftime,obtainedwithtime-accurateandHBsimulations(a);flowandadjointsolver

convergencehistoryofthefirsttimeinstancedensityresidual(b).

Table 4

Inputfrequenciescorrespondingtoadifferentnumberofresolvedtimeinstances, forthe OptC2(i.e.ω2=2.7ω1)configuration,andcorrespondingRMSEofthetotal pressurelosscoefficientbetweentime-accurateandHBsimulationresults.

N. time instances Input frequencies RMSE

5 0,±ω1,±ω2 0.041

7 0,±ω1,±(ω2−ω1),±ω2 0.019

Fig. 15. MachnumbercontourscalculatedatthreedifferenttimeinstanceswiththeHBmethod,basedonthe OptC2testcase,forboththebaseline(a), (b),(c)andtheoptimized(d),(e),(f)bladeprofile.

where A1

=

A2

=

0.04,

ω

1

=

2

π

f1,and

ω

2

=

ω

1

/2.7.

Asforthe OptC1configuration,atime-accuratesimulationisperformedtoverifytheHBsolutionandselecttherelevant inputfrequenciesby analyzingthespectrumofthetotalpressurelosscoefficient.Fig.13adepictstheevolutionintime of thetotalpressurelosscoefficientresultingfromboththeURANSandHBcomputations.Table4reportstheinputfrequency

vectors for the HB simulation with the associated RMSE between the HB and URANS values of the total pressure loss

coefficient.

Figs.15a,15b,and15cshow Machnumbercontourplots,obtainedfor3ofthe7resolvedtime instances.Inthiscase, asopposedtothe OptC1configuration,theinletflowfieldattheupstreamboundaryisuniforminspacewithnoincoming wakes,giventhattime-varying,butuniform-in-space,valuesoftotalpressureandtotaltemperatureareimposed.

Afteranalysisoftheflowsolution,theshapeoptimizationisperformedbasedon7timeinstances,andthecorresponding averaged total pressure loss coefficient is chosen as the objective function. The adjoint-based gradients needed for the optimizationarefirstcomparedwithsecond-orderfinitedifferences.TheresultsofthiscomparisonarereportedinFig.9b, andarecharacterizedbyaRMSElowerthan2

·

10−5.Furthermore,asseeninFig.13b,theadjointsolveragaininheritsthe convergenceratefromtheprimalsolveranditscomputationalcostisabout1.7higherthanthatassociatedwiththeprimal solver.Theaveragecomputationalcostoftheprimalsolver,forasingleiteration,wasabout1.98 secondsona4coresIntel XeonE5-1620Processor withhyper-threading.

The design problem considered inthis test caseis analogous to that forthe OptC1 configuration, and it is formally expressedby (47).Fig.14ashowsthatthe finaldesignisreachedafterapproximately14optimizerevaluations,satisfying theconstraintonthe outletflow angle.TheMach contourplotsfortheoptimizedshape areshowninFig.15.Again,the flow separationontherear suctionside issignificantly mitigatedforall oftheresolved time instances,duetothe lower camber angle of the optimizedblade profile (Fig. 14b). In this case, a reduction in the average pressure loss coefficient ofapproximately14% is achievedwithboththeamplitude andmaximumpeakreducedby44% and29%,respectively.The time-dependentsimulationoftheoptimizedbladeconfirmstheseresults(Fig.12b)andrevealthatnoadditionalfrequencies intheobjectivefunctionspectrumarepresentwhencomparedwiththoseappearinginthesolutionforthebaselineshape. 4. Conclusions

Afully-turbulentharmonicbalancediscreteadjointformulationhasbeendevelopedandappliedtotheshape optimiza-tionoftwo test casesinunsteadyflows:i) apitching airfoilwitha moving grid,andii)anaxial turbinecascadesubject tounsteady flowconditionsattheinlet. Theproposed methodisbasedona duality-preservingalgorithm,whichenables theadjoint solverto inheritthe convergencepropertiesfromtheprimal flow solver. Duetoits efficiency,the framework

enablesshapeoptimizationforquasi-periodicallyforcednonlinearfluidproblemscharacterizedbyasetoffrequenciesthat arenotnecessarilyintegermultipleofonefundamentalharmonic.

Theresultsofthetwotestcaseshaveclearlydemonstratedthatthemethodiscapableofprovidingaccurategradientsin theunsteadysetting,ascomparedtosensitivitiescomputedbysecond-orderfinitedifferences.Thegradient-basedunsteady optimization hasledto improvementsofpractical significance.The meanandthe amplitude ofthetime-varying aerody-namiclosseshavebeenminimizedwithrespectto thebaselineconfiguration.Thishasbeenaccomplishedbyconsidering the minimumnumber ofinput frequencies,according toa spectral analysisofthe flow field,which resultsinsignificant computationalcostsavings.

The development of a fully-turbulent adjoint optimizationframework basedon HB paves theway to the solution of additionalunsteadydesignproblemsthatareencounteredinnumerousadvancedapplications,suchasthemulti-disciplinary optimizationofturbomachinery, includingnovelconceptsofpropulsionsystemsbasedonboundarylayeringestion.Future effortswillbedevotedtothecomparisonbetweenfully-turbulent,time-accurateandHBadjoint-basedshapeoptimization methodsinordertoassessadvantagesanddrawbacks.

Acknowledgements

Thisresearch hasbeensupportedby RobertBoschGmbH, GermanyandtheAppliedandEngineering SciencesDomain

(TTW) ofthe DutchOrganization forScientific Research(NWO), Technology ProgramoftheMinistry ofEconomicAffairs, grant number13385. Theauthorsaregratefultotheir colleagueS. Vitaleforthevaluableinputsanddiscussions. Further-more,thisworkwouldnothavebeenpossiblewithoutthecontributionsfromtheSU2DevelopmentTeam.

References

[1]F.T.Johnson,E.N.Tinoco,N.J.Yu,ThirtyyearsofdevelopmentandapplicationofCFDatBoeingCommercialAirplanes,Seattle,Comput.Fluids34 (10) (2005)1115–1151.

[2]D.N.Mavris,O.J.Pinon,AnOverviewofDesignChallengesandMethodsinAerospaceEngineering,Springer,2012,pp. 1–25.

[3]J.Sobieszczanski-Sobieski,R.T.Haftka,Multidisciplinary aerospacedesignoptimization:surveyofrecentdevelopments, Struct.Multidiscip.Optim. 14 (1)(1997)1–23.

[4]Y.Lian,A.Oyama,M.-S.Liou,Progressindesignoptimizationusingevolutionaryalgorithmsforaerodynamicproblems,Prog.Aerosp.Sci.46 (5)(2010) 199–223.

[5]Z.Li,X.Zheng,Reviewofdesignoptimizationmethodsforturbomachineryaerodynamics,Prog.Aerosp.Sci.93(2017)1–23.

[6]O.Pironneau,Onoptimumdesigninfluidmechanics,J.FluidMech.64 (01)(1974)97–110.

[7]A.Jameson,Aerodynamicdesignviacontroltheory,J.Sci.Comput.3 (3)(1988)233–260.

[8]G.K.Kenway,J.R.Martins,Multipointhigh-fidelityaerostructuraloptimizationofatransportaircraftconfiguration,J.Aircr.51 (1)(2014)144–160.

[9]Z.Lyu,J.R.Martins,Aerodynamicdesignoptimizationstudiesofablended-wing-bodyaircraft,J.Aircr.51 (5)(2014)1604–1617.

[10]D.Wang,L.He,Y.Li,R.Wells,Adjointaerodynamicdesignoptimizationforbladesinmultistageturbomachines—partII:validationandapplication,J. Turbomach.132 (2)(2010)021012.

[11]D.J.Mavriplis,Solutionoftheunsteadydiscreteadjointforthree-dimensionalproblemsondynamicallydeformingunstructuredmeshes,in:Aerospace SciencesMeetings,AIAA,2008.

[12]T.D.Economon,F.Palacios, J.J.Alonso, Unsteadycontinuousadjointapproachforaerodynamicdesignondynamicmeshes, AIAAJ.53 (9)(2015) 2437–2453.

[13]D.J.Lucia,P.S.Beran,W.A.Silva,Reduced-ordermodeling:newapproachesforcomputationalphysics,Prog.Aerosp.Sci.40 (1)(2004)51–117.

[14]K.C.Hall,J.P.Thomas,W.S.Clark,Computationofunsteadynonlinearflowsincascadesusingaharmonicbalancetechnique,AIAAJ.40 (5)(2002) 879–886.

[15]A.Gopinath,E.VanDerWeide,J.Alonso,A.Jameson,K.Ekici,K.Hall,Three-dimensionalunsteadymulti-stageturbomachinerysimulationsusingthe harmonicbalancetechnique,in:45thAIAAAerospaceSciencesMeetingandExhibit,2007,p. 892.

[16]K.Ekici,K.C.Hall,E.H.Dowell,Computationallyfastharmonicbalancemethodsforunsteadyaerodynamicpredictionsofhelicopterrotors,J.Comput. Phys.227 (12)(2008)6206–6225.

[17]E.VanDerWeide,A.Gopinath,A.Jameson,TurbomachineryApplicationswiththeTimeSpectralMethod,AIAAPaper4905,2005.

[18]S.Choi,K.Lee,M.M.Potsdam,J.J.Alonso,Helicopterrotordesignusingatime-spectralandadjoint-basedmethod,J.Aircr51 (2)(2014)412–423.

[19]H.Huang,K.Ekici,Adiscreteadjointharmonicbalancemethodforturbomachineryshapeoptimization,Aerosp.Sci.Technol.39(2014)481–490.

[20]J.P.Thomas,K.C.Hall,E.H.Dowell,Discreteadjointapproachformodelingunsteadyaerodynamicdesignsensitivities,AIAAJ.43 (9)(2005)1931.

[21]S.K.Nadarajah,A.Jameson,Optimumshapedesignforunsteadythree-dimensionalviscousflowsusinganonlinearfrequency-domainmethod,J.Aircr. 44 (5)(2007)1513–1527.

[22]A.Engels-Putzka,C.Frey,Adjointharmonicbalancemethodforforcedresponseanalysisinturbomachinery,in:ProceedingsoftheVIInternational ConferenceonCoupledProblemsinScienceandEngineering,InternationalCenterfor NumericalMethodsinEngineering,Barcelona,Spain,2015, pp. 645–656.

[23]C.Ma,X.Su,X.Yuan,Anefficientunsteadyadjointoptimizationsystemformultistageturbomachinery,J.Turbomach.139 (1)(2017)011003.

[24]Z.Lyu,G.K.Kenway,C.Paige,J.Martins,AutomaticdifferentiationadjointoftheReynolds-averagedNavier–Stokesequationswithaturbulencemodel, in:21stAIAAComputationalFluidDynamicsConference,vol.10,SanDiego,CA,2013,pp. 6–2013.

[25]F.Palacios,J.Alonso,K.Duraisamy,M.Colonno,J.Hicken,A.Aranake,A.Campos,S.Copeland,T.Economon,A.Lonkar,etal.,Unstructured(SU2):an open-sourceintegratedcomputationalenvironmentformulti-physicssimulationanddesign,in:51stAIAAAerospaceSciencesMeetingIncludingthe NewHorizonsForumandAerospaceExposition,StanfordUniversity,2013,p. 287.

[26]T.D.Economon,F.Palacios,S.R.Copeland,T.W.Lukaczyk,J.J.Alonso,SU2:anopen-sourcesuiteformultiphysicssimulationanddesign,AIAAJ.54 (3) (2015)828–846.

[27]E.J.Nielsen,J.Lu,M.A.Park,D.L.Darmofal,Animplicit,exactdualadjointsolutionmethodforturbulentflowsonunstructuredgrids,Comput.Fluids 33 (9)(2004)1131–1155.

[28]D.J.Mavriplis,Multigridsolutionofthediscreteadjointforoptimizationproblemsonunstructuredmeshes,AIAAJ.44 (1)(2006)42–50.

[29]T.Albring,M.Sagebaum,N.R.Gauger,EfficientaerodynamicdesignusingthediscreteadjointmethodinSU2,in:17thAIAA/ISSMOMultidisciplinary AnalysisandOptimizationConference,2016.

[30]S.Vitale,T.A.Albring,M.Pini,N.R.Gauger,P.Colonna,Fullyturbulentdiscreteadjointsolverfornon-idealcompressibleflowapplications,J.Glob. PowerPropuls.Soc.1(2017)252–270.

[31]J.Donea,AntonioHuerta,J-PhPonthot,A.Rodriguez-Ferran,ArbitraryLagrangian–EulerianMethods,Wiley,2004,Ch.14.

[32]A.Jameson,W.Schmidt,E.Turkel,NumericalsolutionoftheEuler equationsbyfinitevolumemethodsusingRunge Kutta timesteppingschemes,in: 14thFluidandPlasmaDynamicsConference,1981,p. 1259.

[33]A.Jameson,TimeDependentCalculationsUsingMultigrid,withApplicationstoUnsteadyFlowsPastAirfoilsandWings,AIAAPaper1596,1991.

[34]A.Gopinath,A.Jameson,Applicationofthetimespectralmethodtoperiodicunsteadyvortexshedding,in:44thAIAAAerospaceSciencesMeetingand Exhibit,2006,p. 449.

[35]S.Banach,Surlesopérationsdanslesensemblesabstraitsetleurapplicationauxéquationsintégrales,Fundam.Math.3 (1)(1922)133–181.

[36]M.Sagebaum,T.Albring,N.R.Gauger,High-PerformancederivativecomputationsusingCoDiPack,arXiv:1709.07229.

[37]B.VanLeer,Towardstheultimateconservativedifferencescheme,V:asecond-ordersequeltoGodunov’smethod,J.Comput.Phys. 32 (1)(1979) 101–136.

[38]D.Kraft,ASoftwarePackageforSequentialQuadraticProgramming,Tech.Rep.,DLRGermanAerospaceCenter—InstituteforFlightMechanics,Köln, Germany,1998.

[39]M.Sagebaum,T.Albring,N.R.Gauger,High-performancederivativecomputationsusingCoDiPack,arXiv:1709.07229.

[40]S.Nimmagadda,T.D.Economon,J.J.Alonso,C.R.IlariodaSilva,Robustuniformtimesamplingapproachfortheharmonicbalancemethod,in:46th AIAAFluidDynamicsConference,2016,p. 3966.

[41]R.M.Hicks,P.A.Henne,Wingdesignbynumericaloptimization,J.Aircr.15 (7)(1978)407–412.

[42]P.Stadtmüller,L.Fottner,AtestcaseforthenumericalinvestigationofwakepassingeffectsonahighlyloadedLPturbinecascadeblade,in:ASME TurboExpo2001,PowerforLand,Sea,andAir,AmericanSocietyofMechanicalEngineers,2001.

[43]F.Menter,M.Kuntz,R.Langtry,TenyearsofindustrialexperiencewiththeSST turbulencemodel,in:Turbulence,HeatandMassTransfer,vol.4(1), 2003,pp. 625–632.

[44]A.Prasad,Calculationofthemixed-outstateinturbomachineflows,J.Turbomach.127 (3)(2004)564.

[45]J.Samareh,Aerodynamicshapeoptimizationbasedonfree-formdeformation,in:10thAIAA/ISSMOMultidisciplinaryAnalysisandOptimization Con-ference,2004.

[46]T.Guédeney,A.Gomar,F.Gallard,F.Sicot,G.Dufour,G.Puigt,Non-uniformtimesamplingformultiple-frequencyharmonicbalancecomputations,J. Comput.Phys.236(2013)317–345.