Fully-turbulent adjoint method for the unsteady shape optimization of multi-row turbomachinery

Fully-turbulent adjoint method for the unsteady shape optimization of multi-row

turbomachinery

Rubino, Antonio; Vitale, Salvatore; Colonna, Piero; Pini, Matteo

DOI

10.1016/j.ast.2020.106132

Publication date

2020

Document Version

Final published version

Published in

Aerospace Science and Technology

Citation (APA)

Rubino, A., Vitale, S., Colonna, P., & Pini, M. (2020). Fully-turbulent adjoint method for the unsteady shape

optimization of multi-row turbomachinery. Aerospace Science and Technology, 106, [106132].

https://doi.org/10.1016/j.ast.2020.106132

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

Aerospace

Science

and

Technology

www.elsevier.com/locate/aescte

Fully-turbulent

adjoint

method

for

the

unsteady

shape

optimization

of

multi-row

turbomachinery

Antonio Rubino,

Salvatore Vitale,

Piero Colonna,

Matteo Pini

∗

Propulsion&Power,FacultyofAerospaceEngineering,DelftUniversityofTechnology,Kluywerveg1,2629HS,Delft,Netherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received31October2019 Receivedinrevisedform1July2020 Accepted5August2020

Availableonline12August2020 CommunicatedbyJérômeMorio

The possibility of taking into account unsteady flow effects if performing turbomachinery shape optimization is attractive to accurately address inherently time dependent design problems. The harmonic balance method is an efficient solution for computational fluid dynamics problems of turbomachinery characterized by quasi-periodic flows. If applied in combination with adjoint methods, it enables the possibility to deal with unsteady fluid-dynamic design in a cost effective manner, opening the way towards multi-disciplinary applications. This paper presents the development of a novel fully-turbulent discrete adjoint based on the time domain harmonic balance method and its application to the constrained fluid dynamic optimization of an axial turbine stage. As opposed to previous works, the proposed method does not require any assumption on the turbulent eddy viscosity and on the set of input frequencies. The results show that the method provides accurate gradients, if compared with second order finite differences, and significant deviation with respect to the sensitivity computed with the constant eddy viscosity approximation. The application of the method to the fluid-dynamic shape optimization of the exemplary stage leads to improve the total-to-static efficiency of 0.8%. The efficiency increase is found to be higher than that obtained by means of a steady state optimization method.

©2020 The Author(s). Published by Elsevier Masson SAS. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

1. Introduction

Adjoint-basedshapeoptimizationmethodsareincreasingly be-comingessential forautomated design. Dueto their efficiencyin obtaining design sensitivities irrespectively of the numberof de-signvariables, these methods have allowed the possibility to ef-fectivelytacklemulti-disciplinaryoptimizationproblems character-izedbyahighnumberofdesignvariablesanddiscretizedonlarge domains [1].

Althoughoriginallyformulatedforaircraftdesign [2,3], adjoint-basedmethodshavebeensuccessfullyextendedtoturbomachinery designproblems.However, themajorityofthemethods currently adopted is based on steady state flow computations, essentially becausethis enablesthe reduction of thecomputational cost for designoptimization [4–9].

Giventheinherentlyunsteadynatureofturbomachineryflows, theuseof unsteadydesign methodsis expectedtoprovide steps forwardinfluiddynamicperformanceascomparedtosteadystate methods.Furthermore,iftransient flow effectsare accountedfor, intrinsically unsteadymulti-disciplinary optimization problemsof

*

Correspondingauthor.

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

turbomachinerycanbeeffectivelyaddressed.Examplesincludethe minimization oftonal noiseintransonic fans [10], the minimiza-tion of structural excitations caused by dynamic fluid-structure interactionphenomena [11,12],andtheaero-thermalperformance improvementofcascadessubjecttounsteadyheattransfer mech-anisms [13–15]. Performing adjoint-based unsteady design for multi-rowturbomachineryproblemsishoweveraformidable chal-lengebecauseofi)thepresenceofmultipleinteractingbladerows, which results in very costly CFD computations ii) the inherent difficulty ofattaining sufficiently converged flow andadjoint so-lutions, which ultimately affects the accuracy of the calculated gradients.

Due to the high computational cost and memory storage re-quirementsassociated withunsteady adjoints [16], several meth-ods have been proposed to improve the efficiency of obtaining time-accuratedesignsensitivities.Thesemethodsmainlytargetthe reduction of memory storage requirements. The algorithm pro-posedhaveresultedinlessaccurategradientcomputationsbytime and space coarsening techniques [17] or in higher I/O overhead ifcheckpointingalgorithmsareadopted [18].Recently,a discrete-adjointmethodhasbeenappliedtotime-accurateturbomachinery optimizationincombinationwithtimecoarsening [19].

Reduced order models have been investigated as a possible effectivealternative to time-accuratesimulations, in orderto de-https://doi.org/10.1016/j.ast.2020.106132

crease the computational cost and storage requirements of the primalsolver.Theharmonicbalance(HB)method,basedon spec-traldiscretizationintimeoftheunsteadyflowequations,isacost effectiveoptionfornon-lineardynamic problemsdominatedbya knownsetoffrequencies.

Past workhasbeenconductedto obtainHB-basedadjoint de-signgradientsforturbomachineryapplications.Nevertheless,these studies are limited tothe design of a single blade row, constant eddyviscosityandtheinabilityto solveforspectral gaps,i.e.the inability to deal with frequencies that are not harmonically re-lated [20,21]. A design algorithm based on a fully-turbulent HB adjointhasbeenrecentlydevelopedandappliedto the optimiza-tion of problemscharacterized by quasi-periodic flows [22]. This algorithmisrestrictedtoasinglecomputationaldomain,thusonly suited forthe automated design of isolated turbomachinery cas-cades.

ThispaperdocumentstheextensionofthenovelHB-based de-sign method proposed in Ref. [22] to fully-turbulent multi-row simulations, enabling the solutionof quasi-periodic unsteady op-timization turbomachinery problems, without any restriction on the turbulent eddy viscosity and on the set of input frequen-cies to be resolved. The method is based on a duality-preserving

approach [23] and it is implemented in the open source code SU2 [24,25].

The design gradients obtainedfrom the HB adjoint equations are verified using second-order central finite differences and ap-pliedtotheconstrainedshapeoptimizationofagasturbinestage. Twoexpansion ratiosareconsideredfortheselectedstage, corre-spondingtosubsonic andtransonicflow conditions.Furthermore, thebaseline stage shape isoptimized by meansofboth the pro-posed HB-based unsteady method and of a steady state method based on the mixing plane (MP) row interface. The objective of thiscomparisonistoassesswhethertheHB-basedautomated de-signprovidesagain incomputedfluiddynamicperformanceover the MP-based one, and if theseare dependent on the operating conditions.Finally,thecomputationalperformance ofthemethod isanalyzedindetailintermsofcomputationalcost,memory,and storagerequirements.

2. Method

2.1. Flowsolver

Let

ρ

be the density, E the total specific energy, t time and

v thevelocity vector ina Cartesianframe ofreference, the semi-discreteformoftheNavier-Stokesequationscanbewrittenas



∂

U

∂

t

+

R

(

U

)

=

0

,

t

>

0

.

(1)

U

= (

ρ

, ρ

v1

, ρ

v2

, ρ

v3

, ρ

E

)

isthevector ofconservativevariables and

R

theresidual operatorapplied to thespacialintegration of the convective and viscous fluxes Fc and Fv. The application of anArbitraryLagrangian-Eulerian(ALE) formulationonthedomain



, moving with velocity u without deforming in time, and its boundary

∂



[26] resultsin

R

(

U

)

=

f

(

Fc

,

Fv

)

in

 ,

t

>

0

,

v

=

u on

∂

 ,

t

>

0

.

(2) Theconvectivefluxesare

Fc

=

⎛

⎝

ρ

v

× (

ρ

(

vv

−

−

uu

)

)

+

p¯I

ρ

E

(

v

−

u

)

+

pv

⎞

⎠ ,

(3)

andtheviscousfluxesare

Fv

=

⎛

⎝

τ

·

¯

¯

τ

·

v

+

κ

∇

T

⎞

⎠ .

(4)

Here, p andT arethestaticpressureandtemperature,

κ

the ther-malconductivity,

μ

thedynamicviscosityand

τ

¯

theviscousstress tensor.Moreingeneral,forRANSequations,thevectorofthe con-servativevariablesU canberedefinedas

U

:=



Ul Ut



,

R

(

U

)

=

R

(

Uf

,

Ut

)

:=



R

f

(

Ul

,

Ut

)

R

t

(

Ul

,

Ut

)



,

(5)

in whichUl

= (

ρ

, ρ

v1

, ρ

v2

, ρ

v3

, ρ

E

)

andUt isthe vector ofthe

conservativevariablesassociatedtotheselectedturbulencemodel. For example, in case of the Menter Shear Stress Transport (SST) model [27],Ut

= (

ρκ

, ρω)

with

κ

theturbulentkineticenergyand

ω

thespecificdissipation.

Using an implicit Euler scheme for time-discretization of (1) leadsto



D

t

(

Uq+1

)

+

R

(

Uq+1

)

=

0

,

(6)

where q is the physical time step index, and

D

t is the

time-derivative operator. After time-integration and linearizing the residualoperatoronecanobtainthefollowingexpressionapplying theharmonicbalancemethod [22] withdual-timestepping [28] of pseudo-time

τ





I



τ

+

J





Un

= − 

R

n

(

Uq

) ,

n

=

0

,

1

, ...,

N

−

1

,

(7)

with



U

=

Uq+1

−

Uq. N is the total number of resolved time instances, linked to the number of input frequencies K by N

=

2K

+

1.Theoperator

R



nisdefinedas



R

n

(

Uq

)

=

R

n

(

Uq

)

+ 

N−1

i=0 Hn,i



Ui

+ 

N

−1 i=0 Hn,iUq_i

,

(8) inwhich H

=

⎛

⎜

⎜

⎜

⎝

H1,1 H1,2

· · ·

H1,N H2,1 H2,2

· · ·

H2,N

..

.

..

.

. .

.

..

.

HN,1 HN,2

· · ·

HN,N

⎞

⎟

⎟

⎟

⎠

,

(9)

istheharmonicbalanceoperator,calculatedas

H

=

E−1DE

.

(10)

E and E−1 are the direct and inverse Fourier matrix, and D is

the diagonal matrixcontaining the K input frequencies, i.e., D

=

diag

(

0

,

i

ω

1

,

...,

i

ω

K

,

−

i

ω

−K

,

...,

−

i

ω

1

)

.Amoredetaileddescription isgiveninRef. [29,22].

Forasteady-statecalculation(7) reducesto





I



τ

+

J





U0

= −R

0

(

Uq

) .

(11) 2.2. Fully-turbulentdiscreteadjointmethod

The equations for the time-domain HB method formulated in Ref. [22] are here extended to account for multi-row fully-turbulentoptimizationofturbomachinery.Thegeneralformulation givenherecanbeappliedtoanydesignprobleminvolving multi-pletime-zonesandgeometricalzonesanditdoesnotrequireany restrictiveassumptionontheeddyviscosity.

Theexpressionof(7) writtenasafixed-pointiterationis

inwhichUz,n and

G

z,narethevectorofconservativevariablesand

theiterationoperator ofthepseudotime-steppingrelative tothe physicalzone z andfortimeinstancen.Eachphysicalzonez

cor-respondstoabladerow.

Usingthedefinitionoftheoperator

G

z,n givenin(12),the

op-timizationproblem for thegeneric objective function

J

can be writtenas minimize αz

J

(

U

(

α

z

),

X

(

α

z

))

subject to Uz,n

(

α

z

)

=

G

z,n

(

U

(

α

z

),

Xz,n

(

α

z

)),

n

=

0

,

1

, ...,

N

−

1 z

=

0

,

1

, ...,

Z

−

1 Xz,0

(

α

z

)

=

M

z

(

α

z

).

(13)

Foreachphysicalzone z andtimeinstancen,

α

z isthesetof

de-sign variable, Xz,n the physical grid, and

M

z is a differentiable

function representing the mesh deformation algorithm. The ob-jective function

J

is obtained asthe spectral average over the resolvedtimeinstancesandthephysicalzones

J

=

1 Z N∗ Z−1

z=0 N

∗−1 n=0

J

z∗,n

,

(14) with

J

z∗,n

=

E∗−1

(

J

z,nE

) ,

(15)

inwhich E∗−1 _{is theextended inversediscreteFourier transform} matrixof size N

×

N∗ calculated for N∗ time instances, whereas

E is the N

×

N discrete Fourier transform matrix computed for the N input time instances (N

<

N∗). Equation (15) allows one, bymeansofFourierinterpolationonuniformlyspacedsamples,to reconstructthetrendoftheobjectivefunctionintime.

The Lagrangian of the constrained optimization problem is givenby

L

=

J

+

Z−1

z=0 N

−1 n=0



(

G

z,n

(

U

(

α

z

),

Xz,n

(

α

z

))

−

Uz,n

(

α

z

))



λ

z,n

+ (M

z

(

α

z

)

−

Xz,0

(

α

z

))



μ

z

.

(16) ThedifferentialoftheLagrangianis

d

L

=

Z−1

z=0 N

−1 n=0

⎛

⎝ ∂

J

∂

Uz,n 

+

Z−1

i=0 N

−1 j=0

∂

G

i,j

∂

Uz,n 

λ

i,j

− λ

z,n

⎞

⎠dU

_z,n

₊

Z−1

z=0 N

−1 n=0



∂

J

∂

Xz,n 

+

∂

G

z,n

∂

Xz,n 

λ

z,n



d Xz,n

−

μ

zd Xz,0

+

Z−1

z=0

∂

M

z

∂

α

z 

μ

_zd

α

z

,

(17) fromwhichtheadjointequationscanbeobtainedas

∂

J

∂

Uz,n 

+

Z−1

i=0 N

−1 j=0

∂

G

i,j

∂

Uz,n 

λ

i,j

= λ

z,n

,

(18) and

∂

J

∂

Xz,n 

+

∂

G

z,n

∂

Xz,n 

λ

z,n

=

μ

z

.

(19) Table 1

Axialturbinestagemainsimulationparameters.

Parameter C1 C2 Unit

Stator inlet blade angle 0 0 [◦]

Total inlet reduced temperature 2.3 2.3 [-]

Expansion ratio 1.5 1.9 [-]

Isentropic work coefficient 0.9 2.2 [-] Inlet turbulence intensity 5% 5% [-] Turbulent viscosity ratio 100 100 [-] Reynolds number 5·106 _6·106 _[-]

Fig. 1. MeshconvergencestudyforthetransonictestcaseC2usingtheharmonic balancemethodwithtwoharmonics.

Equation (19) can be solved directly once the solutionof (18) is known.Similartotheflowsolver,(18) canbeseenasafixed-point iterationin

λ

z,n,namely

λ

qz+,n1

=

∂

N

∂

Uz,n

(

U∗_z,n

, λ

q

,

Xz,n

) ,

(20)

whereU∗_n isthenumericalsolutionfortheflowequation(12) and

N

istheshiftedLagrangiandefinedas

N

=

J

+

Z−1

z=0 N

−1 n=0

G

z,n

(

U

,

Xz,n

)λ

z,n

.

(21)

Finally,the gradientoftheobjective function

J

withrespect tothevectorofthedesignvariables

α

z canbecomputedfromthe

convergedflowandadjointsolutionsusing d

L

 d

α

z

=

d

J

 d

α

z

=

∂

M

z

(

α

z

)

∂

α

z

μ

_z z

=

0

,

1

, ...,

Z

−

1

.

(22) The righthandside of(20) is obtainedusingAlgorithmic Dif-ferentiation applied to the underlyingsource code, including the boundary conditions and the stator-rotor sliding-mesh interface. TheADtooladopted [30] makesuseoftheJacobitapingmethodin combinationwiththeExpressionTemplatesfeatureofC++,leading onlytoaruntimeoverheadintheorderof10

−

20% ascompared totheflowsolver.

3. Application

Thetwo-dimensionalaxial stagedepictedinFig.2waschosen inthisworkinordertotest theproposed method.The blade ge-ometriescorrespondtothemid-spanprofilesadaptedfromthe1.5 stageexperimentalsetupoftheInstituteofJetPropulsionand Tur-bomachinery atRWTHAachen [31]. Comparedtothe original ge-ometry,thestator-rotorbladecountratiohasbeenmodifiedfrom 36

:

41 to41

:

41.Inordertoresemble theflowcharacteristicsofa typicalgasturbinestage,thetestcaseissimulatedunderthe op-eratingconditionsgiveninTable1,whichcorrespondto subsonic (

C1

)andtransonic(

C2

)flowconditions.

Fig. 2. Relative Mach number contour plot forC1. The harmonic balance simulation results depicted in Fig.2b are relative to t=0.

Fig. 3. Verificationofthetotal-to-totalefficiencyasafunctionoftimeobtainedwiththeHBmethodforadifferentnumberofinputfrequenciesandthetime-accurateURANS results(TA).Theconstantvaluecorrespondstothesteadystatesimulationresultsadoptingthemixingplane(MP)stator-rotorinterface.

Thefluid dynamicsimulations arecarriedout usingthe open-sourcecodeSU2 [24,25],extendedinthisworktoallowfor multi-rowHBbased flowsolutions andunsteadyconstrained optimiza-tion using the method discussed in Sec. 2. Forboth

C1

and

C2

theRoescheme [32] isusedtodiscretizetheconvectivefluxesand secondorderaccuracyisobtainedbymeansoftheMUSCL [33] ap-proachwithgradientlimitationbasedontheVenkatakrishnan lim-iter.TheSSTturbulencemodel [27] isemployedforbothtestcases withahybridquad-triangularmeshofapproximately80

,

000 ele-ments,inordertoensureavalueofy+lowerthan1allalongthe blade surface. Non-reflective boundaryconditions are imposed at theinletandoutletoftheturbinecascadeaccordingtothe formu-lationdescribedinRefs. [9,34].Totalconditionsandflowdirection areimposed atinlet,whilethestaticback-pressureisimposed at outlet.Non-reflectivitytreatmentisalsoappliedtothestator-rotor interface.

The mesh convergence study for the

C2

case using the har-monicbalance methodwithtwo harmonicsisdisplayedinFig. 1. Ameshofabout80kelements isdeemed suitedforoptimization purposesanditisthenusedforbothsteadyandunsteady compu-tations.

The selectedobjective functionis the dimensionlessstage en-tropygenerationcalculatedas

sgen

=

s

s,out

 − s

s,in



v2₀

/

T0s,in

+

s

r,out

 − s

r,in



v2₀

/

T0s,in

,

(23)

in which



ss,in and



ss,out are the entropy values calculated as mixed-outaverage [35] overthestatorinletandoutlet.Thesame procedureisused toretrievetheaverageentropyattherotor in-let,i.e.



sr,in,andattherotoroutlet,i.e.



sr,out.T0s,in isthetotal

temperature at the statorinlet and v0 the spouting velocity de-fined as v0

=



2

(

h0,in

−

his,out

)

, whereh0,in is thetotalenthalpy at the inlet ofthe stage and his,out isthe isentropic stage outlet static enthalpy.The stage entropy generationis calculated asthe summation ofthestator androtorgeneration separately inorder to prevent spurious entropy drops across the interface resulting fromnumericalaccuracyandtruncationerrors.

3.1. Flowfieldanalysis

The results from the harmonic balance (HB) simulations are first verified by comparison with the results obtained using a second-orderaccurate intime(TA)simulationandthoseobtained using asteady statemixing plane(MP) modelatthe stator-rotor interface. The time-accurate simulations are based on the dual time stepping method [28], using 50time steps per period (cor-responding tothe bladepitch) and80inner iterationsfora total of10periods.

Fig.3showsthetotal-to-totalstageefficiency,

η

tt,asafunction of time obtained fromthe TA andfrom the HB simulations.The HB-basedunsteady

η

ttapproachestheTAresultsbyincreasingthe numberofresolvedharmonics,withthesubsonicconfiguration

C1

approximatingtheURANS resultsbyresolvinga lowernumberof frequencies if comparedto

C2

.This can be explainedby observ-ing that, incaseoftransonic simulations performedwiththe HB solver,floweffectsthatarenonlinearintimecausedbytheshock interactionbetweenstatorandrotorhavetobemodeled.

Table2summarizes

δ

MP−HB,i.e.,therelativedifferencebetween themaintime-averagedquantitiescharacterizingthestage perfor-mancecomputedfromtheHB simulations with2 harmonicsand the correspondingonesresultingfromtheMP-basedcalculations.

Table 2

Comparisonofmixingplane(MP),harmonicbalance(HB),and timeaccurate(TA) simulationsresults ofthe stageperformance.δMP−HBistherelativedifferencebetweenMPandtime-averagedHBresultsfortheselected performance.RMSETA−HBrepresentstherootmeansquareerrorbetweenthetimedependentTAandHBresults. Parameter MP HB TA δMP−HB[%] RMSETA−HB

C1 ηts[%] 83.15 83.19 83.13 −0.05 1.01e−4

ηtt[%] 95.45 95.29 95.19 +0.17 4.44e−5

C2 ηts[%] 84.24 84.26 84.26 −0.02 2.00e−3

ηtt[%] 96.08 95.61 95.64 +0.49 6.43e−4

Fig. 4. Dimensionless static pressure distribution over the stator and rotor blade surfaces forC1. The total inlet pressure P0is used as reference, for both stator and rotor.

Fig. 5. Dimensionless static pressure distribution over the stator and rotor blade surfaces forC2. The total inlet pressure P0is used as reference, for both stator and rotor.

Furthermore,Table2reportstherootmeansquareerrorofthe HB-basedandthe TA-basedstage performance asafunction of time (RMSETA−HB). Note the RMSE is in the order of 10−4, meaning that the two models provide resultswell inagreement. As a re-sult,twoharmonicsare usedforcomputingtheHB-basedadjoint sensitivities described in thefollowing. Forboth

C1

and

C2

,the total-to-staticstageefficiencygivenbythesteadystatesimulations is characterized by a low deviation compared to the HB time-averagedresults.Therelativedeviationisapproximately0

.

05% for

C1

and 0

.

02% for

C2

. However, the total-to-total efficiency ex-hibits a larger relative difference between the HB-based andthe MP-basedsimulationresults,withthetransonicconfiguration hav-ingarelativedifferenceof0

.

49% andthesubsonicconfigurationof 0

.

17%.

Finally,Fig.4andFig.5showthedimensionlessstaticpressure distribution along the blade profiles retrieved from the MP and HBsimulationresults.Forboth

C1

and

C2

,thesteadystateblade loading differs from the time-averaged harmonic balance blade loading. Furthermore, the shock wave intensity andthe location of the associated flow discontinuity computed by the MP-based simulationdeviatefromthetime-averagedHBresults,forthe tran-sonicconfiguration

C2

.Themainreasonsforthisdifferenceare:i) a steady-state model with the MP interface cannot simulate the

unsteady potential effects generated by the stator-rotor interac-tion;ii)thestatorwakeisnottransportedtotherotorwhenusing the mixingplane interface;iii) fortransonic calculations, theHB methodis abletomodelthe unsteadynon-linear effectsderiving from the shock waves appearing due to the imposed flow con-ditions and to the time-dependent mutual position of the blade rows.

3.2. Adjoint-baseddesignsensitivities

The designgradients ofthe objectivefunction withrespect to the design variables, as defined by (22), are calculated for both thestatorandtherotoroftheselectedtestcase. Tothispurpose, two freeformdeformation (FFD)boxes [36] containing thestator androtorblade profilesare employed.The designvariables(DVs) correspondtothecontrolpointsoftheFFDboxasshowninFig.6

for a setof twelve DVs per row. Therefore,there are in total 48 degreesoffreedomintheoptimizationproblem.

The gradients obtained using the HB-based adjoint equations are first verified using second order finite differences (FD). The results of this verification for the sensitivities in the y-direction are reportedinFig. 7a together withtheresults obtainedusinga steady-stateadjointsolver [9] withthemixing-planeMPinterface.

Fig. 6. Blade profiles and example of free-form deformation (FFD) box.

Fig. 7. ValidationoftheC1normalizedadjoint-baseddesigngradientsusingsecond-orderfinitedifferences,forbothMPandHB(Fig.7a).Rotordesigngradients( y-direction) obtainedwiththeMPandtheHBmethodforadifferentnumberofinputfrequencies(Fig.7b).ThenumberofthedesignvariablescorrespondstotherotorFFDboxgiven inFig.6b.

Forboth HB and MPresults there isa very good agreement be-tweenadjoint-basedandfinitedifferencesgradients,withaRMSE ofapproximately4e

−

3.Thesameresultsare foundforthe sensi-tivitiesintheaxialdirection.

Fig. 7b depictsthe normalized values of the design gradients foran increasingnumberofresolvedfrequencies aswell asthose computedwiththesteadystateMPapproach.Thedesignvariable numbering corresponds to the DV labels given in Fig. 6b. There aretwomainobservationsthatcanbedrawnbyanalyzingFig.7b: i)thevalueofthegradientscomputedwiththeHB-basedadjoint iscomparatively thesame formorethantwo frequencies;ii)the largestdiscrepancybetweenHBandMPsimulationresultsoccurs intheproximity ofthe rotorleading edge.Thiscanbe attributed to theeffects ofwake-rotor interaction, that are not capturedby thesteady-statemodel.All gradientsconvergeto thesamevalues towardstherotoroutlet.

3.3. Constanteddyviscosity(CEV)assumption

AsdiscussedinSec. 2,theproposedHB-basedadjoint method allows one to avoid the use of any restrictive condition on the turbulent eddy viscosity. Past work focused in adopting a con-stanteddyviscosity(CEV)approximationtoeasethedevelopment processoftheadjointsolverandtoimproveitscomputational effi-ciencybutatthecostofalowergradientaccuracy [37,21].Because ofthisconsideration,theimpactoftheCEVassumptiononthe de-signgradientsisassessed. Theaimoftheanalysisdescribed here is to quantify the importance ofadopting fully-turbulent adjoint methodsforunsteadyturbomachinerydesign.

Fig. 8 reports the design gradients of the entropy generation obtained by using a fully-turbulent SST model and the CEV as-sumption.Forboththe

C1

and

C2

operatingconditions,thelargest

Fig. 8. Fully-turbulent(basedontheSSTturbulencemodel)vsconstanteddy viscos-ity(CEV)designgradients,inthe y-direction.Thenumberofthedesignvariables correspondstothestatorandrotorFFDboxesgiveninFig.6.

differencesbetweenthetwocomputedgradientsarethoserelative to therotor design variables.In addition,the deviationsbetween CEVandSST-basedgradientsaremoremarkedforthe

C1

configu-ration,incasetheflowissubsonic.Thiscanbearguablyattributed

Fig. 9. Comparison of the design gradients based on the SST turbulence model and those obtained by using a constant eddy viscosity (CEV) assumption.

Fig. 10. Optimization history. to the larger share of viscous losses in the subsonic case, and

therefore to a higher sensitivity of the objective function to the turbulentflowquantities.Suchdependenceissomewhatnot prop-erlymodeledwhenusingafrozenturbulenceapproximationinthe solutionoftheadjointequations.

Inordertoquantifythesedifferences,therelativedeviation be-tween the sensitivities computed by the CEV and SST model is presented in Fig. 9. Relativedifferences in excess of 20% can be observed for the

C1

configuration whereas they are up to 12% forthe

C2

configuration.Differentlyfromwhathasbeenreported in Ref. [37], the results of this analysis show that the CEV ap-proximation has a relevant effect on the final design gradients ifcompared tothe fully-turbulentadjoint solution.This outcome confirms theresults ofa previous similar study, butlimitedto a comparison between the solutions obtained from a steady-state adjointmethod [38].

Thesefindingsshowthattheproposedadjointmethodcan pro-videgradients thataresignificantly moreaccuratethanthose ob-tained by using the CEV approximation with minimal additional runtime cost. This way, it is possible to obtain fast convergence totheactuallocaloptimumwithanygradient-basedoptimization algorithm.

3.4.Constrainedoptimization

The fully-turbulent design sensitivities, computed with the steady state MP and the unsteady HB method, are employed in agradient-baseddesignproceduretoperformaconstrainedshape optimizationof theselected turbomachinery test case. The mod-ifiedversion of the nonlinear least-squares method (SLSQP) [39] wasselectedasgradient-basedoptimizationalgorithm.

The constrained optimization problem of the turbine stage is formulatedas minimize α sgen

(

α

) ,

α

= {

α

1

,

α

2

}

subject to:P∗

=

P₀∗

,

δ

t,z

= δ

t0,z

,

z

=

1

,

2 n

=

1

,

2

, ...,

N Uz,n

=

G

z,n

,

Xz,n

=

Mz,n

,

(24)

inwhichtheobjectivefunctionisgivenbytheentropygeneration ofthestage,sgen,averagedoveralltheN resolvedtime instances. sgen is afunction of theensemble ofthe statorandrotor design variables,i.e.

α

1 and

α

2.Thedimensionlessnominalpoweroutput ofthestage P∗₀ aswellasthetrailingedgethicknessofbothblade rows

δ

t0 areimposedasconstraints.Thedimensionlesspower

out-putisdefinedasfollows

P∗

=

wm

˙

ρ

0,inypu3_b

,

(25)

with w the Euler work, m the

˙

2D mass flow ratebased on the bladepitch yp,

ρ

0,inthetotaldensityatthestageinlet,andub the bladespeed.

Five time instances are selected to perform the optimization, from the analysis of the spectrum of the objective function and fromthedesignsensitivitiesgivenbyan increasingnumberof in-putfrequencies(seee.g.Fig.7b).

Fig.10showstheevolutionoftheoptimizationforboth

C1

and

C2

.Theentropygenerationandthestagepoweroutputarescaled inordertobettervisualizethedeviationsbetweenthesteadyand theunsteadyresults.Inthecaseofthe

C1

configuration(Fig.10a)

Fig. 11. Shape optimization for theC1configuration.

Fig. 12. Shape optimization for theC2configuration.

Fig. 13. Total-to-total efficiency as a function of time for the baseline and the optimized blade profiles. The constant lines correspond to the average values. the computed stage entropy generation is reduced by about 7%

withthesteadystate optimization methodandby approximately 20% withtheHB-basedoptimization.Theconstraintonthepower outputis satisfiedwithin0

.

6% inbothcases.Fig. 10b depictsthe optimizationconvergenceforthe

C2

configuration.Inthiscase,the objectivefunctionisreducedby approximately11% forthesteady optimizationand14% fortheunsteadyone.Theequalityconstraint onthenon-dimensionalstagepowerismaintainedwithin0

.

8%.

Figs. 11 and 12 report the baseline and the optimized blade profiles.Forboth

C1

and

C2

thelargestdeformationsarelocated intheareaoftherotorleadingedge,asconsequenceofareduction oftherotorincidenceangle.Inotherwords,inbothcases,the opti-mizationprocesssucceededinbetteraligningtherotornosetothe incomingflowdirection.TheHB-basedoptimizedshapeis further-morecharacterized,for

C1

only,byasignificantshapedeformation onthestatorsuctionsideandontherotorreararea.Forthestator,

thisis duebythe higherimpact ofpotential stator-rotor interac-tion effects on the suction side pressure distribution, as can be observedinFig.4.Fromtheanalysisofthefinaldesign,thelargest differencesbetweenthesteadyandtheunsteadyoptimization re-sultsareassociatedwiththesubsonicoperativeconditions,i.e.

C1

. This occursdespite the fact that the discrepancy inperformance between thesteady andthe HB simulation results,computed on thebaselineprofile,islowerforthe

C1

configuration(Table2).

Additionally,atime-accuratesimulationisperformedusingthe optimizedshapesobtainedwithboththesteadyandtheunsteady optimization method. Figs. 13 and 14 display the total-to-total stageefficiencyandthenon-dimensionalstagepowerasafunction of time computedwith URANS simulations as well asthe time-averagedvalues.Table3reportsasummary ofthefinal optimiza-tionresultsbasedonthetotal-to-totalefficiency.Thetotal-to-total efficiency of the

C1

configuration is increased by approximately

Fig. 14. Dimensionless stage power as a function of time for the baseline and the optimized blade profiles. The constant lines correspond to the average values.

Table 3

Summaryofharmonicbalance(HB)andmixingplane(MP)optimizationresults. δHB andδHBaretherelativedifferencebetweenthebaselineandthe optimized resultsbasedonHBandMP,respectively.

Parameter Baseline HBopt MPopt δHB δMP C1 ηtt[%] 95.10 96.05 95.30 +1.00 +0.21 C2 ηtt[%] 95.61 95.95 95.76 +0.36 +0.17

one percentage point using the HB-based optimization. The in-creaseinefficiencyachievedbymeansofthesteadyMP optimiza-tionis of0

.

2 percentage points. Thus, the unsteadyoptimization methodresultsinahigherperformanceimprovement.

Furthermore, the unsteady-based optimization better satisfies the power constraint, as depicted in Fig. 14a. The final time-averagednon-dimensionalpowerfromtheHB-optimizedstage dif-fers by 0

.

3% from the baseline one whereas the MP-optimized solutiondiffersby2

.

0%.

The time-accurate simulation resultsfor the

C2

configuration indicatethatthe HB-basedoptimizationleadsto anefficiency in-creaseofabout 0

.

4% compared to the 0

.

2% obtained by the MP-basedone.Forboth finaldesignsolutionsthepowerconstraintis satisfiedwithin1%.

Finally,forbothoperatingconditionstheoptimizationleads to areducedamplitudeofthetime-dependentefficiencyandthe di-mensionlesspower.Thedecreaseinamplitudeofthetotal-to-total efficiency is 24

.

4% for

C1

and 14

.

1% for

C2

. This demonstrates thattheapplicationoftheproposedunsteadymethodintrinsically affects the variation in time of the objective function. As con-sequence, the variability can be optimized, if needed, by simply reformulating the objective function,i.e. including theamplitude ofthequantityofinterestintheoptimizationproblem.

4. Performanceassessment

The performance of theproposed HB-baseddesign method is assessedintermsofcomputational cost,memoryandstorage re-quirements.Thisanalysisisconductedforboththeprimalandthe adjointsolverona2Dandona3Dgeometry.

4.1.2Dstage

Fig.15showstheperformanceresultsoftheprimalflowsolver for the

C1

configuration. The results are given as a function of theresolvedinputfrequencies.Thetime-accuratesimulations(TA) are initialized from a converged steady simulation. This is done inordertodecreasethenumericaltransientsnecessarytoreacha convergedperiodicsolution.

Fig. 15a reports the computational cost for the steady state (MP),harmonicbalance (HB),andtime-accurate(TA)simulations.

The computational cost forthe HB simulations scales as 2N

+

1, with N the number of input frequencies. This can be explained byrecalling that,withtheproposedHBmethod,inordertosolve

N frequencies 2N

+

1 time instances are required as expressed, e.g., in (8)).However,becauseofthesemi-implicitHBformulation adopted in (7), fora number offrequencies higherthan 4 a de-teriorationoftheconvergencerateisobserved.Thisexplainswhy the computational cost increases at a higher rate if the number ofresolved frequenciesisgreater thanfour(Fig.15a). Inthecase of two input frequencies theHB simulation is approximately 6.5 fasterthanthe time-accurate(TA) simulation.TheTAandtheHB simulations featurenearly thesameCPUtime for10 frequencies. The computationaltime associated tothesteady statesimulation isapproximately3timeslowerthanthatoftheHBsimulation ob-tainedforoneinputfrequency.

Fig. 15b depictsthe memoryand storagerequirements foran increasing number of resolved frequencies. The results are nor-malized usingthevalues obtainedfromthe TAsimulations. Both storage and memory increase linearly at a rate of 2N

+

1. For 2 input frequencies the memory requirement of the TA simula-tions is about 4 times lower than the HB simulations, whereas thenecessarystorageis41timeshigher.Theseresultsoutlinethe performance advantage of usingHB-based over TA-based adjoint methods,forunsteadyturbomachinerydesign.

The CPUtimeandmemoryrequirements oftheadjoint solver are1

.

2 and4

.

5 timeshigherifcomparedtotheprimalflowsolver.

4.2. 3Dstage

The performance is analyzedfora 3D turbine stage operating at the same working conditions of

C1

. The goal of this analysis istoassessthecapabilitytoobtainadjoint-basedsensitivitiesfora three-dimensionalgeometryandtoevaluateitsscalabilityinterms of computational cost, memory requirements and storage. Given theobjectiveoftheanalysis,thecalculationsareconductedby as-suming shrouded blades with and free-slip boundary conditions areappliedtothehubandshroud.Basedonthesemodel assump-tions for the problemathand, a structured mesh of about300k elements was selectedafter amesh independencestudy.Forthis test casethe numericalschemesemployed arethose usedfor

C1

and

C2

.

A HB simulation based on 3 time instances is considered in order tocompute the design sensitivities.Fig. 16ashowsthe ge-ometry ofthe stage aswell asthemid-spancontours ofthe en-tropygenerationnormalizedwiththeinletconditions.Theresults are relativetothe timeinstancet

=

T

/

3.Fig.16bdepictstheHB adjoint-basedsensitivitycorrespondingtot

=

T

/

3.

Furthermore,theprimalandtheadjointsolveraretested con-sideringa varyingnumberof inputfrequencies toinvestigatethe

Fig. 15. Performanceassessmentoftheprimalflowsolver,asafunctionoftheinputfrequencies:(a)Non-dimensionalcomputationalcost;(b)memoryandstorage require-ments.Theresultsofthetimeaccurate(TA)simulationsareusedasreferencevalues.

Fig. 16. 3D axial stage simulation results: (a) normalized entropy generation contours; (b) normalized adjoint-based surface sensitivity.

Fig. 17. Performanceassessmentoftheprimalflowsolver,asafunctionoftheinputfrequencies:(a)Non-dimensionalcomputationalcost;(b)memoryandstorage require-ments.Theresultsofthetimeaccurate(TA)simulationsareusedasreferencevalues.

computational performance. Fig. 17 reports the performance re-sults of the primal flow solver obtained with a number of re-solvedharmonicsrangingfrom1 to5.Similarlytothe2Dtestcase thecomputationalcost,thememory,andstoragerequirements in-creaselinearlyatarateof2N

+

1.When2frequenciesareresolved, thecomputationalcostandthestoragerequiredbytheTA simula-tionareapproximately3

.

5 and42 timeshigherthantheHB-based simulation.ThememoryrequiredbytheHBsolverfor2 harmonics isabout4

.

7 timeslargerthanthatoftheTAsolver.

The computational cost ofthe adjoint solver isapproximately 1.2 times higher than the primal solver, whereas the required

memoryofadjointsolverisabout4.9timesthatoftheflowsolver. The memory requiredby the adjoint solver for the 250000 ele-mentsmeshwas32

.

2Gb.Thetotalsimulationtimefortheadjoint solver was of approximately 450 minutes using a 10 cores Intel XeonE5-2687Wv3 CPUwithhyper-threading.

5. Conclusions

Thisworkdocumentsthedevelopmentofafully-turbulent har-monic balance (HB) discrete adjoint method for multi-row

tur-bomachinery design.The method was applied to theconstrained shapeoptimizationofanexemplaryaxialturbinestage.

Thekeyfindingsofthisstudycanbesummarizedasfollows 1. Thedesignsensitivities can beaccurately calculatedwiththe

proposed HBdiscreteadjointmethod.Thesesensitivitieswere verified usingsecond orderfinitedifferences,withoutany as-sumptionontheturbulenteddyviscosity.

2. Forfluiddynamicdesignpurposesthenumberofrelevant in-putfrequenciestoberesolvedcanbelowerthanthose neces-sarytoaccuratelymodeltheflowbehavior.

3. The assumption of constant eddy viscosity (CEV) was found to significantly affectthe accuracyof thedesign sensitivities. Relative differencesin excessof 20% between theCEV-based andthefully-turbulentresultswerecalculated.

4. Computational cost, memory and storage requirements in-crease linearly ata rate proportional to the number of time instances.Foranumberofinputfrequencieshigherthan4the computational cost featureda slowerconvergenceduetothe semi-implicitformulationadoptedfortheHBflowsolver. 5. The HB-based simulations exhibited higher memory

require-ments but lower storage if compared to time-accurate (TA) RANS simulations. For the analyzed test case, if 2 input fre-quencies areconsidered, the memoryrequirements were ap-proximately 4 times larger than that of the TA simulations, whereasthestoragerequiredwasabout41timessmaller. 6. TheHB adjoint solver featured a computational cost

approx-imately 1.2higher whencomparedto theprimal flow solver. The ratio between the memory required by the adjoint and flowsolverwasapproximately4.5.

7. The optimization results achieved by the proposed HB ad-jointshowremarkabledifferenceswhencomparedwithsteady state optimization results.Differences inthe optimized total-to-total stageefficiencyup to0.8percentagepoints were ob-tainedfortheexemplary2Dtestcase.

8. Forthe analyzedtestcase, theuseoftheunsteady optimiza-tionmethodalwaysledtobetterfluiddynamicperformanceif comparedtothesteadystateoptimizationresults.

Thefocusofthepresentanalysiswas ontheunsteady adjoint-basedfluiddynamic optimizationofaturbinestage characterized byanequalnumberofbladecountperrow.Futuredevelopments aredevotedtoextendtheperiodicboundaryconditionoftheflow solverinordertosimulateasingleblade passagehavingunequal azimuthalbladepitch.Thiswouldenable theresolutionof multi-stageunsteadyproblemscharacterizedbyasetoffrequenciesthat arenotintegermultipleofonefundamentalharmonic.

Declarationofcompetinginterest

Theauthorsdeclarethattheyhavenoknowncompeting finan-cialinterestsorpersonalrelationshipsthatcouldhaveappearedto influencetheworkreportedinthispaper.

Acknowledgements

Thisresearch hasbeensupported by RobertBosch GmbHand theAppliedandEngineeringSciencesDomain(TTW)oftheDutch OrganizationforScientificResearch(NWO),TechnologyProgramof theMinistryofEconomicAffairs,grantnumber13385.

References

[1]K.W. Kenway Gaetan, R.R.A. Martins Joaquim, Multipoint high-fidelity aerostructuraloptimizationofatransportaircraftconfiguration,J.Aircr.51 (1) (2014)144–160.

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

[3]M.B.Giles,N.A.Pierce,Anintroductiontotheadjointapproachtodesign,Flow Turbul.Combust.65 (3–4)(2000)393–415.

[4]D.Wang,L.He,Adjointaerodynamicdesignoptimizationforbladesin mul-tistage turbomachines—part I: methodology and verification, J. Turbomach. 132 (2)(2010)021011.

[5]D.Wang,L.He, Y.Li,R.Wells,Adjointaerodynamicdesignoptimizationfor bladesinmultistageturbomachines—partII:validationandapplication,J. Tur-bomach.132 (2)(2010)021012.

[6]J.Luo,J.Xiong,F.Liu,I.McBean,Three-dimensionalaerodynamicdesign opti-mizationofaturbinebladebyusinganadjointmethod,J.Turbomach.133 (1) (2011)011026.

[7]B.Walther,S.Nadarajah,Optimumshapedesignformultirowturbomachinery configurations usingadiscreteadjoint approachand anefficientradial ba-sisfunctiondeformationschemeforcomplexmultiblockgrids,J.Turbomach. 137 (8)(2015)081006.

[8]M.Pini,G.Persico,D.Pasquale,S.Rebay,Adjointmethodforshape optimiza-tioninreal-gasflowapplications,J.Eng.GasTurbinesPower137 (3)(2015) 032604.

[9]S. Vitale,T.A.Albring,M. Pini,N.R. Gauger,P.Colonna, Fullyturbulent dis-creteadjointsolverfornon-idealcompressibleflowapplications,J.Glob.Power Propuls.Soc.1(2017)252–270.

[10]J.J. Defoe,Z.S. Spakovszky, Effectsofboundary-layeringestion on the aero-acousticsoftransonicfanrotors,J.Turbomach.135 (5)(2013)051013. [11]L.He,D.Wang,Concurrentbladeaerodynamic-aero-elasticdesignoptimization

usingadjointmethod,J.Turbomach.133 (1)(2011)011021.

[12] A.Engels-Putzka,J.Backhaus,C.Frey,Forcedresponsesensitivityanalysisusing anadjointharmonicbalancesolver,J.Turbomach.141 (3)(2019),https://doi. org/10.1115/1.4041700.

[13]V.Michelassi,F.Martelli,R.Dénos,T.Arts,C.Sieverding,Unsteadyheattransfer instator-rotorinteractionbytwoequationturbulencemodel,in:ASME1998 InternationalGasTurbineandAeroengineCongressandExhibition,American SocietyofMechanicalEngineers,1998,V001T01A058.

[14]A.A.Ameri,D.L.Rigby,E.Steinthorsson,J.Heidmann,J.C.Fabian,Unsteady anal-ysisofbladeandtipheattransferasinfluencedbytheupstreammomentum andthermalwakes,J.Turbomach.132 (4)(2010)041007.

[15]N.R.Atkins,S.J.Thorpe,R.W.Ainsworth,Unsteadyeffectsontransonicturbine blade-tipheattransfer,J.Turbomach.134 (6)(2012)061002.

[16]D.J.Mavriplis,Multigridsolutionofthediscreteadjointforoptimization prob-lemsonunstructuredmeshes,AIAAJ.44 (1)(2006)42–50.

[17]S.Nimmagadda,T.D.Economon,J.J.Alonso,C.Silva,B.Y.Zhou,T.Albring, Low-costunsteadydiscreteadjoints foraeroacousticoptimizationusingtemporal and spatialcoarsening techniques, in:2018 AIAA/ASCE/AHS/ASC Structures, StructuralDynamics,andMaterialsConference,2018,p. 1911.

[18]Q.Wang,P.Moin,G.Iaccarino,Minimalrepetitiondynamiccheckpointing al-gorithm forunsteadyadjointcalculation,SIAMJ.Sci.Comput.31 (4)(2009) 2549–2567.

[19]G. Ntanakas,M.Meyer, K.C. Giannakoglou,Employing thetime-domain un-steady discreteadjointmethod forshapeoptimization ofthree-dimensional multirowturbomachineryconfigurations,J.Turbomach.140 (8)(2018)081006. [20]H.Huang,K.Ekici,Adiscreteadjointharmonicbalancemethodfor

turboma-chineryshapeoptimization,Aerosp.Sci.Technol.39(2014)481–490. [21]C.Ma,X.Su,X.Yuan,Anefficientunsteadyadjointoptimizationsystemfor

multistageturbomachinery,J.Turbomach.139 (1)(2017)011003.

[22]A. Rubino,M. Pini, P.Colonna, T.Albring,S. Nimmagadda, T.Economon, J. Alonso,Adjoint-basedfluiddynamicdesignoptimizationinquasi-periodic un-steadyflowproblemsusingaharmonicbalancemethod,J.Comput.Phys.237 (2018)220–235.

[23]T.Albring,M.Sagebaum,N.R.Gauger,Efficientaerodynamicdesignusingthe discreteadjointmethodinSU2,in:17thAIAA/ISSMOMultidisciplinaryAnalysis andOptimizationConference,2016.

[24]F. Palacios, J. Alonso, K. Duraisamy, M. Colonno, J. Hicken, A. Aranake, A. Campos,S.Copeland,T.Economon,A.Lonkar,etal.,StanfordUniversity Un-structured (SU2):anopen-sourceintegrated computationalenvironmentfor multi-physicssimulationanddesign,in:51stAIAAAerospaceSciencesMeeting IncludingtheNewHorizonsForumandAerospaceExposition,2013,p. 287. [25]T.D.Economon,F.Palacios,S.R.Copeland,T.W.Lukaczyk,J.J.Alonso,SU2:an

open-sourcesuiteformultiphysicssimulationanddesign,AIAAJ.54 (3)(2015) 828–846.

[26]J. Donea, Antonio Huerta, J-Ph Ponthot, A. Rodriguez-Ferran, Arbitrary Lagrangian-EulerianMethods,Wiley,2004,Ch.14.

[27]F.R.Menter,Zonaltwoequationk-ω,turbulencemodelsforaerodynamicflows, in:FluidDynamicsandCo-locatedConferences,AIAA,1993,pp. 1993–2906. [28]A.Jameson,Timedependentcalculationsusingmultigrid,withapplicationsto

unsteadyflowspastairfoilsandwings,in:AIAApaper1991-1596,1991. [29]K.R.Naik,Fouriercollocationmethodsforunsteadyflows,Ph.D.thesis,Stanford

University,2017.

[30]M.Sagebaum,T.Albring,N.R.Gauger,High-performancederivative computa-tionsusingCoDiPack,arXivpreprint,arXiv:1709.07229,2017.

[31]B.Stephan,H.Gallus,R.Niehuis,Experimentalinvestigationsoftipclearance flowanditsinfluenceonsecondaryflowsina1-1/2stageaxialturbine,in: ASMETurboExpo2000:Power for Land,Sea,and Air,AmericanSocietyof MechanicalEngineers,2000,V001T03A099.

[32]P.L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes,J.Comput.Phys.43 (2)(1981)357–372.

[33]B.VanLeer,Towardstheultimateconservativedifferencescheme.V.A second-ordersequeltoGodunov’smethod,J.Comput.Phys.32 (1)(1979)101–136. [34]M.B.Giles,NonreflectingboundaryconditionsforEulerequationcalculations,

AIAAJ.28 (12)(1990)2050–2058.

[35]A.P.Saxer,Anumericalanalysisof3-Dinviscidstator/rotorinteractionsusing non-reflectingboundaryconditions,Tech.Rep.,GasTurbineLaboratory, Mas-sachusettsInstituteofTechnology,Cambridge,Mass.,1992.

[36]J.Samareh,Aerodynamicshapeoptimizationbasedonfree-formdeformation, in:10thAIAA/ISSMOMultidisciplinaryAnalysisandOptimizationConference, 2004,p. 4630.

[37]A.C.Marta,S.Shankaran, OnthehandlingofturbulenceequationsinRANS adjointsolvers,Comput.Fluids74(2013)102–113.

[38]S.Vitale,AdvancementsinautomateddesignmethodsforNICFD turbomachin-ery,Ph.D.thesis,DelftUniversityofTechnology,2018.

[39]D.Kraft,Asoftwarepackageforsequentialquadraticprogramming,Tech.Rep., DLRGermanAerospaceCenter—InstituteforFlightMechanics,Koln,Germany, 1998.