Customized data-driven RANS closures for bi-fidelity LES–RANS optimization

Publication date

2021

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Zhang, Y., Dwight, R. P., Schmelzer, M., Gómez, J. F., Han, Z. H., & Hickel, S. (2021). Customized

data-driven RANS closures for bi-fidelity LES–RANS optimization. Journal of Computational Physics, 432,

[110153]. https://doi.org/10.1016/j.jcp.2021.110153

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

Customized

data-driven

RANS

closures

for

bi-fidelity

LES–RANS

optimization

✩

Yu Zhang

a

,

Richard

P. Dwight

b

,

∗

,

Martin Schmelzer

b

,

Javier

F. Gómez

b

,

Zhong-hua Han

a

, Stefan Hickel

b

a_{InstituteofAerodynamicandMultidisciplinaryDesignOptimization,NationalKeyLaboratoryofScienceandTechnologyonAerodynamic} DesignandResearch,SchoolofAeronautics,NorthwesternPolytechnicUniversity,YouyiWestRoad127,710072Xi’an,People’sRepublicof China

b_{AerodynamicsGroup,FacultyofAerospaceEngineering,DelftUniversityofTechnology,Kluyverweg2,2629HSDelft,theNetherlands}

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory: Received6April2020

Receivedinrevisedform18September 2020

Accepted8December2020 Availableonline26January2021 Keywords:

Turbulencemodelling Algebraicstressmodel Reynolds-averagedNavier-Stokes Large-eddysimulation Multi-fidelityoptimization

Multi-fidelityoptimizationmethodspromiseahigh-fidelityoptimumatacostonlyslightly greater than a low-fidelity optimization. This promise is seldom achieved in practice, duetotherequirementthatlow- and high-fidelitymodels correlatewell.Inthisarticle, we propose an efficient bi-fidelity shape optimization method for turbulent fluid-flow applicationswithLarge-EddySimulation(LES)andReynolds-averagedNavier-Stokes(RANS) as the high- and low-fidelity models within a hierarchical-Krigingsurrogate modelling framework. Sincethe LES–RANScorrelationisoften poor,weusethefull LES flow-field atasingle pointinthe design spaceto deriveacustom-tailoredRANSclosuremodel that

reproduces the LES at that point. This is achieved with machine-learning techniques, specifically sparse regression to obtain high corrections of the turbulence anisotropy tensor and theproductionofturbulencekineticenergyasfunctions oftheRANS mean-flow. The LES–RANS correlation is dramatically improved throughout the design-space. Wedemonstratetheeffectivityand efficiencyofourmethodinaproof-of-conceptshape optimizationofthewell-knownperiodic-hillcase.StandardRANSmodelsperformpoorlyin thiscase,whereasourmethodconvergestotheLES-optimumwithonlytwoLESsamples.

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

1. Introduction

Numericalfluid-dynamicshape-optimizationis anincreasingly centraltool inengineeringpractice.Typically Reynolds-AveragedNavier-Stokes(RANS)isusedasthefluidmodel,duetoitsacceptablecomputationalcostandaccuracy.However inmanyimportantsituations,thephysicsdemandsscale-resolvingsimulationsofturbulencesuchaslarge-eddysimulation (LES).Thisincludes,forinstance,designswhereseparationoccurs;junctionsflowswithcomplexturbulencebehavior (Simp-son [34], Belligoliet al. [2]);andnovelboundarylayerflowcontrolapplications.LESisoftenappliedforafinalanalysisand validation,butitshighcomputationalcostprecludesitsusewithinthedesignloop.

✩ _{ThisdocumentistheresultsoftheresearchprojectfundedbyInnovationFoundationforDoctorDissertationofNorthwesternPolytechnicalUniversity} undergrantNo.CX201801andChinaScholarshipCouncilNo.201706290083.

*

Correspondingauthor.

E-mailaddresses:zhangyu91@mail.nwpu.edu.cn(Y. Zhang),r.p.dwight@tudelft.nl(R.P. Dwight),m.schmelzer@tudelft.nl(M. Schmelzer), jfatoug@gmail.com(J.F. Gómez),hanzh@nwpu.edu.cn(Z.-h. Han),s.hickel@tudelft.nl(S. Hickel).

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

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

However,thesuccessofalltheseapproachesreliesheavilyonthecorrelationbetweenthehigh- andlow-fidelity mod-els (Han et al. [12], Zhang et al. [39]).Since modelswithhigherfidelitytypicallyhavealsosignificantly highernumerical costs,inpracticethehigh-fidelitymodelischosenasthecheapestmodelthatpredictstherelevantphysics.Necessarily,the low-fidelity modelwilllacksome oftherelevantphysics,leading toapoorcorrelation betweenthe models,whichlimits thepracticalapplicationofMFOmethods.

In separate developments of the past few years, data-driven methods forturbulence modelling based on supervised machine learning have been introduced to improve RANS predictions when reference data is available (Xiao and Cin-nella [38], Duraisamy et al. [6], Durbin [7], Kaandorpand Dwight [22]). Notably, Parish andDuraisamy [29] introduced alocalmultiplicativetermtocorrecttheturbulenceproductioninthek-equationofthek

−

ω

model.Thistermwaschosen by solvingan inverseproblemto matchLES referencedata,andwas thenusedtotrain aGaussian processtocorrectthe baseline turbulencemodel.Linget al. [24] trainedadeepneuralnetworktopredict turbulenceanisotropy ai j andreplace thebaseline turbulencemodel.An alternativeistouserandom-forestsforthesametask (Wang et al. [35], Kaandorpand Dwight [22]). Even though these approaches generalize the linear eddy-viscosity concept, they generate complex black-box closure models. More promising are approaches that generatecompact explicit expressions formodels. Notable are gene-expression programming(GEP)(WeatherittandSandberg [36,37]),anddeterministic sparse regressionofturbulence anisotropy (SpaRTA)(Schmelzer et al. [33]).Thesemethods generatemodelsthat can berapidlyimplemented inexisting CFD codesandevaluated atevery iteration ofaRANS solution,andpotentially inspected toidentify thephysical mecha-nismsinfluencingtheflow. TheSpaRTAapproachhasshowntheabilitytoconsistentlyreproducespecifiedLESmean-flows inaRANSsolver,usingonlythreetofiveadditionalnon-linearclosureterms.Itisimportanttounderstandthat–attheir presentlevelofdevelopment –thesemethodsdonotgenerategeneral-purpose turbulencemodels,butratherdeliver cor-rectionsforflowssimilartotheflowtheyweretrainedon.Thus,forflowswithahighdegreeofsimilarity,e.g.modifications inthegeometry,themethodsgeneratemodelseffectivelycustomizedtotheflowathand.

Thekeyoriginalcontributionofthispaperisahighly-efficientbi-fidelityoptimizationprocedure,usingLESasthe high-fidelitymodel,anddata-enhanced RANSasthelow-fidelitymodel.Inoutline,we:

1. PerformasingleLESsimulationatthebaselinegeometry

ψ

0,i

=

0;

2. Usethefull-fieldLESdatatogenerateacustomizedRANSmodelmatchingtheLESat

ψ

i; 3. SamplethiscustomRANSmodelthroughoutthedesign-space(cheap);

4. Combinecustom-RANSandLESresultsinabi-fidelitysurrogate;

5. PerformanewLESatthepointofmaximumexpectedimprovement

ψ

i+1 (expensive); 6. GeneratenewcustomRANSmodel(s),basedonallLESdataavailablesofar;

7. i

←

i

+

1,goto3.

The key observation isthat thecustom RANS modelcorrelates well withLES ina region of the design-spaceforwhich the LEStraining dataisinformative, andthat thiscorrelation applies to mean-velocitiesaswell asthe Reynolds-stresses (originating fromthe LES training data). Providedthe physicalprocesses occurringdo not significantlychange, the data-driven RANS remains a viable model. Also, whereas typical multi-fidelity methods use only the objective function from the high-fidelity simulation,our approach uses much more information.There are manyvariations on the basic pattern described above,for example: Which LESsamples are used fortraining inStep 6?;What criteria are usedfor sampling RANSinStep3,andLESinStep5?.

In the followingwe demonstrate the methodon a proof-of-concept optimizationproblem: a generalizedperiodic-hill (P-H)geometry,inwhichthesteepnessofthehillsisvaried (Fatou Gómez [8]).Thebaselineperiodic-hillflow(asproposed byMellenet al. [25] basedonanexperimentalstudybyAlmeidaet al. [1])isnotoriousforpoorperformanceofessentially allstandardRANSmodels,duetoitssensitivityontheflowseparationlocationandlarge-scalelow-frequencydynamics.We useahierarchicalKrigingbi-fidelitysurrogatemodel,withanefficientglobaloptimization (EGO)samplingprocedureforthe LES (Joneset al. [20], Zhanget al. [39]).The customRANSmodelissampledonauniformgrid,andupdatedateachstep withthelatestLESdata.ToachievetheLESoptimalsolutionwerequireonlytwoLESsamples(andonefurthersampleto

verifythe result).Ourmethodthereforeoffers a promisingpath towardsLES-quality optimizationwithonlya handful of LESsamples.

The paperis structured asfollows: inSection 2 we describethe computational settingof thebaseline RANS andLES simulations. Section 3 briefly describes SpaRTA, themachine-learning methodused in thiswork to generatedata-driven RANS models.SpaRTAis demonstratedonthebaseline periodic-hillcase(Mellenet al. [25]).Ourproposed bi-fidelity op-timizationprocedure isdescribedindetail,andappliedtothe generalizedP-Hoptimizationcasewiththe hillwidthasa designvariableinSection4.

2. BaselineincompressibleLESandRANSsimulations

ThesolversusedhereforLESandRANSarebothfinite-volumemethods,butotherwisenumericallydistinct- wedescribe both briefly.Forthestandardperiodic-hill test-casewecompareourLES resultswiththeexperimentaldata ofRapp [32] and well-resolved LES referencedata from Breueret al. [3], and demonstrate excellent agreement. We show that RANS resultsobtainedwithastandardk-

ω

SSTmodelarepoor,whichisinagreementwithliterature.

2.1. BaselineReynolds-averagedNavier-Stokeswithk-

ω

SST

Assumingincompressibleflowwithconstantfluiddensityequaltounity,thesteadyRANSequationsare

∂

Ui

∂

xi

=

0

,

(1) Uj

∂

Ui

∂

xj

=

∂

∂

xj



−δ

i jP

+

ν

∂

Ui

∂

xj

−

τ

i j



+ δ

i1f

,

whereUi withi

,

j

∈ {

1

,

2

,

3

}

arethecomponentsofthemean-flowvelocity, P isthemeanpressure,and

ν

isthekinematic viscosity. Thevolumeforcing f servestodrivethe flowthroughthe doublyperiodicdomain.The effectofturbulenceon the momentum equation is confined to the Reynolds stress tensor

τ

i j, which must be modelled. In this paper, we use thepopulark-

ω

SST modelasa baselinemodel.The modeldefinestransport equationsfortheturbulencekinetic energy k

:=

1₂

τ

ii,andthespecificturbulencedissipationrate

ω

:

Uj

∂

k

∂

xj

=

τ

i j

∂

Ui

∂

xj

  

Pk

−β

∗k

ω

+

∂

∂

xj



(

ν

+

σ

k

ν

t

)

∂

k

∂

xj



,

(2) Uj

∂

ω

∂

xj

=

γ

ν

t

τ

i j

∂

Ui

∂

xj

− β

ω

2

+

∂

∂

xj



(

ν

+

σ

ω

ν

t

)

∂

ω

∂

xj



+

2

(

1

−

F1)

σ

ω2 1

ω

∂

k

∂

xj

∂

ω

∂

xj

.

(3)

Theeddyviscosityismodelledas

ν

t

:=

a1k

max

(

a1

ω

,

S F2)

,

andtheanisotropicpartofReynolds-stresstensorismodelledbytheBoussinesqassumption

ai j

:=

τ

i j

−

2 3

δ

i jk



a B i j

:= −

2

ν

tSi j

,

where Si j

:=

1 2



∂

Ui

∂

xj

+

∂

Uj

∂

xi

,

and the isotropic part is absorbed into the pressure. All remaining terms and coefficients are omitted for brevity, see (Menter [26])fordetails.

The computational mesh used for RANS simulation of the periodic-hill is two-dimensional, with 120

×

130 cells in stream-wise andwall-normal directionsrespectively, achieving y+

≤

1 onall walls. Periodic boundaryconditions are set at inlet and outlet, and no-slip conditions on the walls. Volume forcing is used to drive the flow, witha proportional– integral–derivative (PID)controllerusedtoachievethetargetvelocity.ThesimulationisperformedatReh

=

10595 basedon thehillheighth usingsecond-orderaccurateSIMPLEsolver.Thepressureissolvedusinggeometric-algebraicmultigridwith aGauss-Seidelsmoother,whileUi,k and

ω

aresolvedusingaDiagonal-Incomplete-LU-preconditionedBiCGmethod.

2.2. Large-eddysimulationfortheperiodic-hills

The high-fidelity and “ground-truth”model inthis work iswall-modelled LES withour in-housefinite-volumesolver INCA. We use the incompressible staggered-grid version of the solver with a block-Cartesian background mesh andthe conservative immersed boundary method of Meyer et al. [27] for representing the geometry. Our LES follow a holistic modellingapproach, wherethe numericaldiscretizationandthe subgrid-scale (SGS)turbulenceclosureare fullymerged: Theadaptivelocaldeconvolutionmethod(ALDM)ofHickelet al. [17], HickelandAdams [16] isanonlinear finite-volume

Fig. 1. Meanvelocityprofilesfor:referenceLES( ) (Breueret al. [3]);PIV(◦) (Rapp [32]),ourLES( )andRANSwithk-ωSSTturbulencemodelling ( ).(Forinterpretationofthecolorsinthefigure(s),thereaderisreferredtothewebversionofthisarticle.)

discretizationschemetailoredforLESofturbulent flows.Optimummodelformanddiscretizationparameterswere learnt by a physics-informed geneticalgorithm using turbulencedata, spectral analysis andconstraints fromturbulence theory (Hickelet al. [17]).

Thenear-wallturbulenceismodelledratherthan resolved (Hickelet al. [19]).Inparticularweusethemodelof(Chen et al. [5]),basedontheturbulenceboundary-layerequationsandamixing-lengtheddy-viscositymodel.Validationresults fortheperiodic-hills–amongst othercases–areavailable in(Chenet al. [5]).Atthewall ourgridsmaintain y+

<

20.A third-orderexplicitRunge-Kuttaschemeisusedfortimeintegrationandthepressure-correctionPoissonequationissolved using BiCGstab withincomplete LU andalgebraic-multigrid as preconditioners. Overall, ourLES solver is very similar to that usedby Hickel et al. [18] forwall-resolvedLES oftheperiodic-hillflow, andlargely identicaltothat usedby Meyer et al. [28] forwall-modelledLESofalow-speedaerospaceapplication.

Theperiodichilliscomputedwithadomainsizeof

(

9h

,

3

.

03h

,

4

.

5h

)

withh

=

thehillheight,thegeometryandboundary conditionsfollowexactlythestandarddefinitionsoftheERCOFTAC Testcase9.2 [25].Aparameterstudywasperformedon theeffectsofthespan-wiseextentofthedomain,themeshresolution,andthesolutionaveragingtime(Fatou Gómez [8]). Agoodcompromise ofaccuracyandcostwas achievedforaspanwiseextentof Lz

=

4

.

5h,about1

.

2 millioncells (giving

y+

<

20 everywhere),andaveraging over



55 flow-through timesstartedafteran initial transientof



89 flow-through times.

TheLESsimulationcodeisparallelizedby domaindecompositionbasedontheMessagePassingInterface(MPI)library, andwasrunon20Intel(R)Xeon(R)CPUE5-2670v2coresat2.5 GHz.Attheseconditionseachcasecostapproximatelytwo daysofwall-clocktime,withsteepergeometriesgenerallymorecostlythanshallowcases.

Withtheseparameters,thesimulationsareperformedontheTUdelftHPC,andconvergedstatisticsareobtainedinless than 2days ona 20coresIntel(R) Xeon(R) CPUE5-2670v2at2.5GHz and128GBofRAM. Fig. 1showsourLES results forthemeanvelocityprofilesatseveralstreamwiselocationsforthestandardperiodic-hillgeometryandthewell-resolved LESdataofBreueret al. [3] (13

.

1 millioncells)andthereferenceexperimentsofRapp [32].WeobservethatourINCALES resultsareinverygoodagreementwiththereferencedata.

AlsoshowninFig.1aretheRANSresultsforthebaselinek-

ω

SSTmodeldescribedinSection2.1.TheRANSsimulation significantlyover-predicts thelength oftheseparatedflowregion,whichdistortstheflowandleadstolargeerrors every-where inthe domain. Havingflow-featuresvery dissimilar to the high-fidelity LES, thebaseline RANS is therefore likely inappropriateasalow-fidelitymodelinabi-fidelityoptimizationorsurrogatemodellingprocedure.

3. EnhancingRANSwithdata-driventurbulencemodelling

Inthissection webriefly describetheSpaRTA proceduretogenerateimprovedRANS modelsfromLESdata–see[33] foradetaileddiscussion.Theprocedurehastwomainstages:

1. Weusethek-corrective-frozen-RANS approach(Schmelzeret al. [33])to solveforcorrectivefieldsfortheanisotropy tensor,andturbulenceproduction.Thesefields,whendirectlyinjectedintotheRANSsolver,leadthesolvertoreproduce theLESmean-flow.

2. We use deterministic symbolic regression to find a low-complexity approximate map from the resolved mean-flow quantitiestothesecorrectivefields.Thismaprepresentsourcorrectiontothebaselineturbulencemodel.

3.1. Identifyingmodel-formerrorwiththek-corrective-frozen-RANSapproach

It iswell known(XiaoandCinnella [38])that injectingLES/DNSReynolds stressesdirectlyintoRANS solversdoesnot always improve accuracy of the mean-flow prediction. This issue can be partially addressed by using a field inversion procedure to identify corrective fields that result in the correct mean-flow (Duraisamy et al. [6]). In our problem full-field LESdataisavailable,andso acostlyinversionprocedureisunnecessary.Insteadwedefine aniterationthat leadsto equivalentcorrectivefields.

LetLESstatisticsbedenotedwithasuperscript“



”,sob

i j

=

ai j

/

2kisthenormalizedLESanisotropytensor,kistheLES turbulencekineticenergy,andU theLESmeanvelocityfield. Wedefineacorrectiontothenormalizedanisotropy tensor b i j via: b_{i j}

=

b_{i j}B

+

b_{i j}

,

b_{i j}B

:= −

ν

ˆ

t kS  i j

,

(4)

wherethecorrection isdefinedwithrespecttotheBoussinesq assumptionevaluated usingtheLESstrain-ratetensor S

i j,

k _{andeddy-viscosity}

_ν

_ˆ

t.Parneix et al. [30] andWeatheritt andSandberg [37] model the specific dissipation rate

ω

by passivelyevaluatingthe

ω

-equation(3),usingLESdataforotherquantities.Wefurtherextendthisprocedurebyrequiring thatthek equationalsosatisfiestheLESdata,theso-calledk-corrective-frozen-RANSapproach.Tothisendweintroducea (spatiallyvarying)correctionR to

ˆ

thek- and

ω

-equationsgiving:

U_j

∂

k 

∂

xj

=

P_k

+ ˆ

R

− β

∗k

ω

ˆ

+

∂

∂

xj



(

ν

+

σ

k

ν

ˆ

t

)

∂

k

∂

xj



,

(5) U_j

∂

ω

ˆ

∂

xj

=

γ

ˆ

ν

t

P_k

+ ˆ

R



− β ˆ

ω

2

₊

∂

∂

xj



(

ν

+

σ

ω

ν

ˆ

t

)

∂

ω

ˆ

∂

xj



+

2

(

1

−

F1

)

σ

ω2 1

ˆ

ω

∂

k

∂

xj

∂

ω

ˆ

∂

xj

,

(6)

in which quantitieswith a “^”are to be solved for,andremaining quantitiesderive fromthe LES data.The LES datais interpolatedontoaRANSmesh,(5) and(6) aresolvedreadilyinalooselycoupledfashion,andtheresulting

ω

ˆ

isusedto compute b_{i j} from(4). Thus twocorrective fieldsresult, thetensor-field b_{i j} andthe scalar-field R.

ˆ

When thesefieldsare injectedintotheRANSsolver,theresultingmean-fieldsconsistentlyagreewellwithLES,seeFig.2.Thisistrueforboththe baselineperiodichillswiththewidth

ψ

=

1,andsteeperandshallowerhillswith0

.

25

≤ ψ ≤

4 - i.e.thefulldesignspace (Mellenet al. [25]).

3.2. Deterministicsymbolicregressionforcorrectivefieldmodelling

Havingidentifiedcorrectivefields,to makepredictionsitisnecessarytomodelthesefieldsintermsofknown (mean-flow)quantities.FollowingPope [31],weassumethatthenon-dimensionalstrain-rateandrotation-ratetensors

˜

Si j

:=

1

ω

Si j

,

˜

i j

:=

1 2

ω



∂

Ui

∂

xj

−

∂

Uj

∂

xi

aresufficienttodescribethecorrectivefields,leading toPope’swell-knownbasis-tensorseries.TheCayley-Hamilton theo-remdictatesthat themostgeneralformofthefunctionb_{i j}

( ˜

Si j

,

˜

i j

)

is(undertheassumptions ofanalyticityandGalilean invariance): b_{i j}

( ˜

Si j

, ˜

i j

)

=

10



n=1 T_{i j}(n)

α

n

(λ

1

, . . . , λ

5

),

(7)

where T(_{i j}n) are tenbasis tensors,

λ

m are thefiveinvariants of S and

˜

˜

,and

α

n

(

·)

arearbitraryscalar functionsofthese invariants.Inthispaper,weonlyconsiderthefirstthreebasetensors,andfirsttwoinvariants,whicharegivenby:

T_{i j}(1)

= ˜

Si j

λ

1

= ˜

SmnS

˜

nm

,

T_{i j}(2)

= ˜

Sik

˜

kj

− ˜

ikS

˜

kj

λ2

= ˜

mn

˜

nm

,

T_{i j}(3)

= ˜

SikS

˜

kj

−

1

3

δ

i jS

˜

mnS

˜

nm

.

TomodelR we

ˆ

assumethatittakestheform

ˆ

R

=

2k

∂

Ui

∂

xj bR_{i j}

,

Fig. 2. MeanvelocityprofilesforPHcases:referenceLES( )(INCA),FrozenRANS( ),SpaRTARANS( )andRANSwithk-ωSSTturbulence modelling( ).

forsome tensor-fieldbR_{i j},modelledsimilarlytob_{i j} by (7).It remainsto fitafunctionofthe form(7) totheLESdata,for whichweusesymbolicregression.

For details refer to Schmelzer et al. [33]. In short, we forma large library of candidate basis functions, consisting of powers andproducts of theinvariants

λ

m.We then regressthe data– consistingof

(λ

m

,

T(_{i j}n)

,

b_{i j}

)

triples ateach spatial

Fig. 3. Comparisonofexactcorrectivefield(basedonLES)andSpaRTAmodelprediction(left:a11 inReynoldsstress;Right: R residualofktransport equation).

meshpoint–againstthislibrarywithsparsity-promoting



1 _{regularization.Thisselectsasmallsetofbasis functionsthat} representthedatawell.Concretely,lettheregressingfunctionbe:

b_{i j}



4



n=1 T(_{i j}n)

_M



m=1

φ

m

(λ)

· θ

nm



=

L

(

T

, λ

; θ),

where

θ

arecoefficientsoftheM basisfunctionsinthelibrary

φ

m

(

·)

foreachbasistensor,andL

(

·)

istheimpliedoperator, whichislinearin

θ

.Thenwesolvetheelastic-netregularizedregressionproblem:

θ

fit

:=

arg max θ∈





L

(

T

, λ



; θ) −

b_{i j}



_

2 2

+

σρ

θ

1

+

σ

(

1

−

ρ

)

θ

2 2

,

(8)

wherethefirstnormisoverall meshpoints ofalltraining cases,andonceagain



denotesquantities evaluatedfromthe LESdata.Wesolveforanumberofdifferentvaluesof

σ

and

ρ

,toobtainanumberofcandidatesparsitypatterns,andwe discardthevaluesofthenonzerocoefficients

θ

.Foreach ofthesparsity patterns,weperformridgeregression,usingonly thebasisfunctionsidentifiedinthesparseregression.Let

θ

sbeavectorofnon-zerocoefficientsidentifiedby(8),thenwe solve

θ

s_fit

:=

arg max θ∈





L

(

T

, λ



; θ

s

)

−

b_{i j}



_

2 2

+

σ

r



θ

s



2 2

,

(9)

foreachsparsitypattern.Thefinalmodelisselectedbasedonacompromisebetweensparsity,goodness-of-fit,andideally performance in cross-validation if data is available (Schmelzer et al. [33]). This two step procedure allows us to inde-pendently controlthe levelof sparsity (with

ρ

,

σ

), andthemagnitudeof theresulting coefficients(with

σ

r), whichcan otherwisebecomeverylarge.

ApplyingthisproceduretoourLESdataforthestandardP-Hcase(

ψ

=

1),theresultingalgebraic modelsforb

i j andR

ˆ

are: b_{i j}



M_b i j

=

2

.

8T (2) i j

; ˆ

R



MR

=

0

.

4T (1) i j

·

2k

∂

Ui

∂

xj

.

(10)

Fig.3showstheoptimalai j

=

2kb_{i j} andR identified

ˆ

fortheP-Hcaseusingthek-frozenapproachoftheprevioussection. Thepredictionsofthemodel(10) arealsoshown.Theagreementisgenerallygood,withthenotableexceptionofthelarge

ˆ

R nearthewall,whichisnot reproducedby(10).Thisfailurecan beattributedtothenon-dimensionalizationof S,

by thespecificdissipationrate

ω

,whichbecomeslargeclosetothewall,pullingallbasis-tensors T(n)_{downtozero.Thiseffect}

can beovercomebyusingcoefficients

α

n thatgoto infinitynearthewall,asinKaandorp andDwight [22].However,we considerthisunnecessary,ascorrectionsclosetothewall,especiallyinb_{,haveverylittleeffectonthemean-flow.}

Note that the above regression procedure is independent of the CFD solver, andis performed on a mesh-point-wise basis withflat data-vectors.OncetheresultingSpaRTA modelis obtained,itisimplemented asan OpenFOAMturbulence modelmoduleinC++usingautomaticcodegeneration.Themoduleiscompiled,anddynamicallylinkedtothemainsolver, withouttheneedforrecompilingOpenFOAMitself.

Fig. 4. Flowchart of the adaptive variable-fidelity optimization.

HavingimplementedM_b i j and

ˆ

R intheRANSsolver,weverifythatthemean-flowofthetrainingcaseisreproduced,and wemakepredictionsfor

ψ

=

0

.

25 and

ψ

=

4

.

0 andcomparewithourLESreference,seeFig.2.Inallcasesthemean-flow issignificantlyimprovedoverthebaselinek-

ω

SST,butdonotalwaysreachtheaccuracyofthedirectlyinjectedcorrective fieldsshown(“frozenRANS”inFig.2)whichrepresentabestcasescenario.

4. Bi-fidelityoptimizationloopwithcustomRANSaslow-fidelity

WeaimtoefficientlysolvethediscretizedPDE-constrainedminimizationproblem min

ψ∈J

(ψ

;

U

)

subject to RLES(ψ

;

U

)

=

0

,

where

ψ

aregeometricdesign-variablesinthedesignspace



,U isthefullflow-state, J isthecost-functionandRLES isthe (highest-fidelity) discretizedPDE operator(which dependson

ψ

via boundary-conditions).Weassume that nosignificant modellingisrequiredtoevaluate J givenU ,whichisthecaseformostcost-functionsofinterest.Wealsohaveavailablea low-fidelitydiscretizedoperatorRRANS[M

](ψ;

U

)

,whichdependsonsomeclosuremodelM derivedusingtheprocedureof theprevioussection.

Theoutlineofourproposedbi-fidelityoptimizationisgiveninFig.4.Thestepsare:

1. SolveRLES

(ψ

0

,

U

)

=

0 forthebaselinegeometry

ψ

= ψ0

,toobtainturbulencestatisticsU,

τ

i j,k.Leti

=

0. 2. Identifythemodel-formerrortermsb

i j andR using

ˆ

thek-corrective-frozen-RANSapproachofSection3.1andtheLES datafromStep1.SubsequentlyfindthealgebraicmodelsM0

b i j

,M0_R usingSection3.2,toobtainRRANS[M0

](ψ;

U

)

. 3. Generate N

1 low-fidelitysamplepoints

ψ

₁

,

. . . ,

ψ

_N bydesignofexperiment(DoE)(e.g.Latin-hypercubesampling)in

the(possiblymulti-dimensional)designspace.SolveRRANS[Mi

](ψ

_j

;

U

)

=

0 andevaluate J ateachsample j

∈ {

1

,

. . . ,

N

}

. 4. Buildabi-fidelityKrigingsurrogateof J usingLESsamplesat

ψ

1

,

. . . ,

ψ

i,andlow-fidelitysamplesat

ψ

1

,

. . . ,

ψ

N. 5. Choosenewhigh-fidelitysample

ψ

i+1 bymaximumexpectedimprovementonthesurrogateofStep4.EvaluateLESat

thispoint.

6. Trainanewlow-fidelitymodelMi+1 _{basedonLESdataat}

_ψ

1

,

. . . ,

ψ

i+1. 7. i

←

i

+

1.Ifconvergenceisachieved,thenterminate;otherwisegotoStep3.

Note that in thiswork the low-fidelity modeldepends onthe result ofthe high-fidelity (andthe locations atwhich thisisevaluated),incontrasttostandardmulti-fidelitymodels.Thispresentsachallengefromthestandpointofstatistical modelling ofthe bi-fidelity system– forexamplein co-Kriging, wherea prior correlation relationship must be specified between the fidelities. This will be the subject of future work. We bypass the issue hereby using hierarchical Kriging. Theoptimizationalgorithmisconductedusinganin-housetool“SurroOpt” (Han [11]).Notethatourmethodologyis inde-pendentof theLESmesh-resolution,andcan becombinedwithanyother high-fidelity datasource,such asDESorDNS. Our surrogate based optimization framework already generalizes to high-dimensional parameter spaces, constraints, and multipleobjectives,seeAppendixA.

4.1. Variable-fidelitysurrogatemodelling:hierarchicalKriging

The hierarchicalKrigingmodel (HanandGörtz [13])(HK) isconstructed ina recursiveway. First,an ordinaryKriging model

ˆ

JLow is builtfor thelow-fidelity objective function. Then a surrogatefor thehigh-fidelity costfunction

ˆ

J is built

usingthescaledlow-fidelityKrigingmodel

ˆ

JLow asthemodeltrend.Inthiswaythehigh-fidelityresponses aretreatedas realizationsofarandomprocess:

ˆ

J

(ψ )

=

ρ

ˆ

JLow

+

Z

(ψ )

(11)

whereZ isastationaryGaussianprocesswithzeromeanandacovarianceCov

[

Z

(ψ),

Z

(ψ

∗

)

]

=

σ

2_r

_(ψ,

_ψ

∗

₎

_,where

_σ

2_isthe processvariance,andr

(

·,

·)

andthecorrelationfunction.Inthispaper,weuseacubicsplinefunction(HanandGörtz [13]), whichintroduceshyper-parameters

θ

.

Similar to the standard single-fidelity Kriging model, the corresponding predictor and mean-squared error (MSE) of Krigingforthehigh-fidelitycostfunctionare

ˆ

J

(ψ )

=

ρ

ˆ

JLow(ψ )

+

rTR−1

(

J

−

ρ

F

)

(12)

ε

(ψ )

2

=

σ

2

{

1

−

rTR−1r

+ [

rTR−1F

− ˆ

JLow(ψ )

]

2

/(

FTR−1F

)

−1

}

(13) where

F

= [ ˆ

JLow(ψ1), . . . , ˆJLow(ψN

)

] ∈ R

N

.

Here, J

∈ R

N _{is the vector of high-fidelity observations at sample sites}

_ψ

_{= [ψ1}

_,

_{. . . ,}

_ψ

N

]

. R is the correlation matrix betweentheobservedhigh-fidelitysamplesandr isthecorrelationvectorbetweentheobservedsamplesandthepredicting point.Therefore,theKrigingmodelforhighfidelitycostfunctioncanberegardedasasumofscaledlowerfidelityKriging predictorandthediscrepancybetweenthescaledlower-fidelityfunctionandhigher-fidelityfunction.

ModelfittingoftheHKmodelusesthemaximum-likelihoodestimate: max ρ,θ ,σ2L

(

ρ

, θ ,

σ

2

₎

₌

_max ρ,θ ,σ2 1



2

π σ

2

_|

_R

_|

exp



−

1 2

(

J

−

F

ρ

)

TR−1

(

J

−

F

ρ

)

σ

2



,

(14)

givenwhichanexplicitexpressionfor

ρ

is:

ρ

= (

FTR−1F

)

−1

(

FTR−1J

).

Itremainstospecifyinwhatmannernewhigh-fidelitysamplesarechosen.Weusethemaximumexpectedimprovement (EI) principle (Jones et al. [20], Zhang et al. [39]).Define the improvement with respect to a currentbest value Jmin as (formally):

I

(ψ )

=

max

(

Jmin

− ˆ

J

(ψ ),

0

).

SincethepredictionoftheHKmodelatanyuntriedsite

ψ

isthenormaldistribution

N ( ˆ

J

(ψ), ε(ψ)

2

₎

_{,wechoosethevalue} of

ψ

thatmaximizestheexpectedimprovement

E

I

(ψ)

.Thishasaclosedformexpression:

E

I

(ψ )

= (

Jmin

− ˆ

J

(ψ ))



Jmin

− ˆ

J

(ψ )

ε

(ψ )



+

ε

(ψ )φ



Jmin

− ˆ

J

(ψ )

ε

(ψ )



,

where

φ (

·)

and

(

·)

arerespectivelytheunitnormaldensityanddistributionfunctions. 4.2. Periodichilloptimization:costfunctionanddesign-spaceexploration

The baselineperiodic-hill geometry(hillwidth

ψ

=

1)isthe widelyknownERCOFTAC Testcase 9.2withthe geometry definition of Mellen et al. [25]. The design variable is the normalized hill width: we simply scale the geometry ofthe original hill as

(

x

,

y

,

z

)

−

> (ψ

∗

x

,

y

,

z

)

, asillustrated in Fig.5, whileholding constant the distancebetweensubsequent peaks.Forvery wide hills (e.g.

ψ

=

4) the hills overlap,andnoflat regionbetweenthe hills remains.The geometries of periodic-hillsatthelower- andupper-boundaryofthedesign-spaceareshowninFig.5.

We choosean costfunction depending onthetotaldrag on thehills andaveraged turbulencekinetic energyoverthe whole flowfield(at aconstantmass-flow). Theobjectiveistosearch forthegeometrywithmaximumturbulencekinetic energyandminimumdrag. Atypical engineeringapplicationwithsimilarobjectives isthedesign ofreactorsorheat ex-changers in process engineering, whereone tries to minimize pressureloss andmaximize turbulentmixing ofreactants orofhot- andcoldfluids.Forthepresentproof-of-conceptoptimization,weightsforthesetwotermssuchthattheglobal optimum isnot too closeto theupper andlower boundariesofthe design-space



= [

0

.

25

,

4

]

.This choiceallows usto demonstrateunconstrainedoptimizationinthisproof-of-conceptstudy.Thechosenfunctionis:

min

ψ∈

−(¯



k

−

2



.

5 f

)/

J0



J

Fig. 6. Mean-squared error of mean velocity obtained by custom RANS and baseline RANS.

wherek,

¯

istheaveragedturbulentkineticenergyovertheflowfield,and f ,thevolumeforcingtermin(1),actsasaproxy fordrag.Thecostfunction J isnormalizedbytheLESresultofthebaselinegeometry, J0

:= ¯

k

−

2

.

5 f.

Novalidationdataexistsforthesegeometries,soourconfidenceinthecorrectnessofourLESrestsonthevalidationwe performedforthestandardgeometryinSection 2.2.Notethatinthiswork,LESisusedasourground-truth:thetargetof ouroptimizationistheLESoptimum.As suchit’saccuracyisnot ofprimaryconcernforvalidatingourmethod,provided onlythatit’ssignificantlymoreaccuratethanbaselineRANS.

Before the optimizationwe explore thedesign spacewiththe referenceLES, andthebaseline- and customized RANS models.Weusethecustommodel(10) trainedat

ψ

=

1 fromSection3.2andcompareagainstourLESground-truth.

InFig.6,weseethatthemean-velocitypredictionissignificantlyimprovedcomparedwiththebaselineRANSsimulation, across themajorityofthedesign-space. Fig.7showsthat thepredictionoftheaveraged turbulencekineticenergyk and

¯

proxy-drag f arebothsignificantlyimproved- aconsequenceofmoreaccurateReynoldsstressesandmeanflowfieldsfor all test geometries. The custom RANS performs bestwhen

ψ

≥

1, forwhich the flow showsa medium-sized separation bubble.Forthesteepergeometries

ψ <

1,thelargerseparationbubblesareslightlyunder-estimatedbythecustomclosure model.InFig.8,weobservethatthecustomRANSmatchesthemaintrendoftheLESobjectivefunctionverywell.TheLES optimizationproblemhastwo localminima,near

ψ

=

0

.

75 and

ψ

=

2

.

0,whereas theobjectiveestimatedby thecustom RANSonlyshowsthefirst.Withadditionalhi-fidelityLESsamples,thecustomRANS modelbecomesmoreaccurateasthe optimizationprogresses. Thebaseline RANS model,however,doesnot reproducethemagnitudeoftheobjectivefunction, northegenerallocationoftheminima:aconsequenceofsignificantlyunder-predictingbothdragandReynoldsstress.

Despitethesimplegeometry,thisflowisgovernedbyhighlycomplexinteractionsbetweenmeanflowaccelerationand deceleration, separationon a curvedwall, turbulentmixing inthe separatedshearlayer, flow reattachment and redevel-opmentoftheboundarylayer.Therelationbetweenthedesignvariable

ψ

,turbulencekineticenergyandpressuredrag is thereforenonlinear andnon-monotonic, whichleads totwo minimawiththechosen objectivefunction.We don’thavea detailedphysicalexplanationforthisphenomena,northekinknear

ψ

=

0

.

7 (visibleinFig.9).Thiscouldbethesubjectof afurtherpaperlookingatthedetailedflow-physicsofthisparameterizedgeometry.

Fig. 7. Comparison of cost-function components obtained by LES, custom RANS and baseline RANS.

Fig. 8. Design-space exploration of the objective function J .

4.3. Optimizationconvergence

Ourexpectationisthat: byenhancingquality ofthe low-fidelity model,theconvergenceofthebi-fidelity optimization will be improved. The low-fidelity modelneed not be highly-accurate forthe whole design spacein orderto be useful. IndeedgivenasingleLEStrainingsample,weexpectthecustomized-RANSclosuretobemostaccurateclosetothetraining point,andleastaccurateatthelimitsofthedesignspace.Astheoptimizationprogresses,andmoreLESsamplesareadded, wethenexpecttobecomemoreaccuratewhereverthereareLESsamples.

Theoptimizationstarts,asshowninFig.9(a),withoneLESasthehigh-fidelitysample(redsolidcircles)and26custom RANS simulationsasthelow-fidelitysamples(redhollowcircles).TheinitialHKresponsesurface passesthroughthe LES, andsince there isonly asingle high-fidelity sample, the shape isdictated entirelyby thecustom RANS, andthescaling parameter

ρ

in HKmodelis found tobe 0

.

78.The expectedimprovementcriterion selects asecond LES sample at

ψ



0

.

93,Fig. 9(a). LES isperformed atthis point, anda custom closure isbuilt, which isfound to be b

i j



Mb i j

=

1

.

2T (2) i j , and R

ˆ



MR

=

0

.

45Ti j(1) ∂Ui

∂xj. Using this newclosure model, the low-fidelity samples everywhere in the design spaceare reevaluated. With the updated low-fidelity samples and two LES high-fidelity samples,the HK model is rebuilt, shown in Fig. 9(b). Again the closurepassesthrough the LES samplesexactly, and thehierarchical Kriging iscloser to thetrue cost functioninthe regionneartheoptimum.The scalingparameter

ρ

inthe HKmodelisfoundto be0

.

93,closerto 1 becauseofthebettermatchbetweenlow- andhigh-fidelity objectivefunctions.AthirditerationisshowninFig.9(c).The convergence oftheprocedure intermsof J isshowninFig.10 (redsolid line).InFig. 9,thereferenceforthe objective functionisobtainedbyperforming59LESsamplesinthedesignspace (Fatou Gómez [8]).

Fig. 9. Threeiterationsoftheproposedmethod.(referencebyLES( ) (Fatou Gómez [8]),HK(—),EI( ),hi-fidelityLESsamples( ),low-fidelity customRANS( )andglobaloptimum( )).

Note that in Fig. 9 the initial values of the EI peak are very low,and at later iterationshigher. This is fairly typical behaviorfortheEIinKriging-basedoptimization;aconsequenceofinitialpoorestimationoftheKriginghyper-parameters by maximum-likelihoodestimation(MLE),duetothevery lownumberofsamples.ThevarianceoftheKrigingmodel(as wellasthelength-scale)changessignificantlyaseachnewsampleisadded,leadingtolargechangesinboththeresponse, andespeciallytheEI.Thisproblemreducesasthenumberofsamplesincreases,anddisappearsinthelimit.Giventhatwe are usingvery fewsamples,an alternative toMLEisto specifyvariance andlength-scalesby handandkeep themfixed.

Fig. 10. Convergenceofbi-fidelityoptimization(customRANSaslow-fidelitymodel( , ),baselineRANSaslow-fidelitymodel( , )),andsingle fidelityoptimizationusingLESonly( , ).

Thishasthedisadvantageofaddingmoreparameters tothemethod,andpotentiallybiasingtheEIhowever,sowedonot useithere.

Asacomparison,wealsoperformedasingle-fidelityoptimizationusingonlythehigh-fidelityLESmodel,andabi-fidelity optimizationusingthebaselineRANSasthelow-fidelitymodel.WiththesameinitialLESsample,theconvergencehistories intermsof J areshowninFig.10(blackdash-dottedlineandbluedashedline,respectively),whereineachcaseweshow theobjective-functionvalueofthemostrecentlyobtainedsample.Theconvergencehistoryofbi-fidelityoptimizationusing the baselineRANS modelasthelo-fimodel,showsnoimprovementover thesingle-fidelity optimization–andisinfact significantly moreexpensivethanpure LESoptimization,dueto thelargenumberofadditionalRANS solvesneeded. This isatypicalexamplebi-fidelitymethodswithpoorlo-fitohi-ficorrectionresultinginnospeed-upcomparedtohi-fionly. Consequently,morethan10 LESsamplesarerequiredtolocatetheglobaloptimum.

Theconvergencecriterionusedintheseoptimizations,simplydetectswhen J achievestheglobaloptimumbasedonthe 59 reference LESsimulations.Thisisnotpracticalinrealapplications,wherewe wouldsuggest usingone ofthestandard EGO-basedconvergencecriteria.

Finally, the optimum solution is found at

ψ

=

0

.

75944 with the cost function value around

−

1

.

68. In summary, the optimizationwithourproposedmethodologyconvergedtowithin theaccuracyofthehigh-fidelitymodelwithonly2LES samples,withonefurthersampleusedtoverifytheresult.

5. Conclusions

In this work, a novel bi-fidelity fluid-dynamic shape optimizationmethodology is proposed, in which LES is used as thehigh-fidelity simulationtool andanautomatically customizedRANS turbulenceclosureasthe low-fidelitymodel.The full-fieldLESdataobtainedfromthehigh-fidelitysamplesareusedtotrainaRANSmodelfortheflowbeingoptimized.The SpaRTA methodisused asarobust andeffectivewayoftailoringRANS models withsimplealgebraic augmentation.The customizedlow-fidelityRANSmodelissuccessivelyimprovedasmoreLESdatabecomesavailablewithinthedesignspace. ThetwofidelitiesarecombinedinahierarchicalKrigingsurrogate,andtheEGOprocedureisusedtoaddnewhigh-fidelity samples.

Given thetechnical complexity of thisprocedure, we have demonstratedit on a proof-of-concept shape optimization of aperiodic-hill geometrywithvarying hill-width. In thiscase, the methodconverges tothe trueoptimum (withinthe tolerance oftheLESsimulations)intwo iterations.Shape optimizationforflow-problemswherethedetails ofturbulence areimportant,includingseparatedflows,junctionflows,andtransitionalflows,maybecomecomputationallyfeasiblewith thisscheme.OurmethodologyisindependentofLESresolution(e.g.wall-resolvedLEScouldbeusedwithoutmodification), andcouldbecombinedwithanyotherreferencemethod,URANS,DES,DNSetc.

Future work will target application to a three-dimensional optimization problem of engineering interest. Our goal is that thiswillserve asa demonstratorofLES-driven optimizationasaserious optionin engineeringoptimization. Sucha problemwilllikelyinvolvehigh-dimensionalparameterspaces,constraints,andpossiblymultipleobjectives.The surrogate-basedoptimizationweusehere,generalizeswelltoalltheserequirements,seeAppendixA,andwillnotbethebottleneck. ItremainstobeseenwhetherourcustomizedRANSmodelswillgeneralizewellwithinmoregeneralparameter-spaces,and howtheywillcopewithdesign-spacesincorperatingdramaticallydifferentphysicsatdifferentdesignpoints.

Fig. 11. Flowchart of a SBO-type optimizerSurroOpt.

Declarationofcompetinginterest

The authors declare that they haveno known competingfinancial interests or personal relationships that could have appearedtoinfluencetheworkreportedinthispaper.

Appendix A. Surrogatebasedoptimizationwith

SurroOpt

SurroOpt

(Han [11])isaresearchcodedevelopedforacademic- andengineeringdesigndrivenbyexpensivenumerical simulations. It is a surrogate-based optimization (SBO) code, whichcan be used to efficiently solve arbitrarysingle and multi-objective(Paretofront),unconstrainedandconstrainedoptimizationproblemsinhigh-dimensionalparameterspaces (Fig.11).

Generally,aSBOprocessconsistsofthefollowingsteps:

1. Initialsamplepointsarechosenbyadesignofexperiment (DoE)method,andtheobjective(s)andconstraintsare evalu-atedatthesamplesbyanexpensivenumericalanalysiscodesuchasCFDsolver.

2. Initialsurrogatemodelsfortheobjectivesandconstraintsarebuilt,basedonthesamplingdataset.

3. Asub-optimizationsolvesanoptimalinfill-samplingproblem,basedonthemodelsfromstep2.,whichinturngenerates newsamplepointstobeevaluatedbytheexpensiveanalysiscode.

4. The newly selected sample points, as well as the functional responses are added to the sample database and the surrogatemodelsareupdated.

5. Steps3.and4.arerepeateduntilaterminationconditionissatisfied.

ThemainingredientsofsuchaSBOprocessaretherefore:DoE,surrogatemodeling,infill-samplingandsub-optimization, andterminationconditions.

SurroOpt

hasbuilt-in modernDoE methodssuitedfordeterministiccomputer experiments, suchasLatinhypercubesampling(LHS),uniformdesign(UD)andMonteCarlodesign(MC).Avarietyofsurrogatemodels areincluded,suchasquadraticresponsesurfacemodels(PRSM),kriging,gradient-enhancedkriging(GEK),hierarchical krig-ing (HK),co-Kriging, radial-basisfunctions(RBFs),artificialneutralnetwork(ANN),support-vector regression(SVR),etc. A numberofinfill-samplingcriteriaandthededicatedconstrainthandlingmethodsareavailable,suchasminimizingsurrogate prediction(MSP),minimizinglowerconfidencebounding(LCB),maximizingexpectedimprovement(EI),maximizing proba-bilityofimprovement(PI),andmaximizingmeansquarederror(MSE).Somewell-knownoptimizationalgorithms,suchas Hooke&Jeevespatternsearch,BFGSquasi-Newtonmethod,sequentialquadraticprogramming(SQP),orsingle/multi objec-tiveGAs,areemployedtosolvethesub-optimizations,inwhichthecostfunction(s)andconstraintfunction(s)areevaluated

bythecheap-to-evaluatesurrogatemodels.Inturn,newsamplepointsaregeneratedandevaluatedbytheexpensive anal-ysiscode.Fourterminationconditions aredefinedin

SurroOpt

,whichset thresholdsforthemaximum numberofCFD evaluations, the minimum distance betweensamples, the minimum value of expected improvement or the accuracy of surrogate model.

SurroOpt

managescompute cores withMPI, giving theuser severaloptions forparallelization of the simulationcodeandsamplingstrategy.See (Han [11])fordetailsandreferences.

Thisframeworkwillallowustogeneralizethemethodologydescribedinthispapertohigh-dimensional,multi-objective, constrainedoptimizationproblemsofpracticalengineeringinterest.

References

[1] G.Almeida,D.Durao,M.Heitor,Wakeflowsbehindtwo-dimensionalmodelhills,Exp.Therm.FluidSci.7(1993)87–101,https://doi.org/10.1016/0894 -1777(93)90083-U.

[2] Z.Belligoli,A.J.Koers,R.P.Dwight,G.Eitelberg,Usingananti-fairingtoreducedragatwing/bodyjunctions,AIAAJ.57(2019)1468–1480,https:// doi.org/10.2514/1.J057481.

[3] M.Breuer,N.Peller,C.Rapp,M.Manhart,Flowoverperiodichills–numericalandexperimentalstudyinawiderangeofReynoldsnumbers,Comput. Fluids38(2009)433–457,https://doi.org/10.1016/j.compfluid.2008.05.002.

[4]K.J.Chang,R.T.Haftka,G.L.Giles,P.J.Kao,Sensitivity-basedscalingforapproximationstructuralresponse,J.Aircr.30(1993)283–288.

[5]Z.Chen,S.Hickel,A.Devesa,J.Berland,N.Adams,Wallmodelingforimplicitlarge-eddysimulationandimmersed-interfacemethods,Theor.Comput. FluidDyn.28(2014)1–21.

[6] K.Duraisamy,G.Iaccarino,H.Xiao,Turbulencemodelingintheageofdata,Annu.Rev.FluidMech.51(2019)357–377,https://doi.org/10.1146 /annurev-fluid-010518-040547.

[7] P.A.Durbin,Somerecentdevelopmentsinturbulenceclosuremodeling,Annu.Rev.FluidMech.50(2018)77–103,https://doi.org/10.1146 /annurev-fluid-122316-045020.

[8]J.FatouGómez,Multi-fidelityCo-KrigingoptimizationusinghybridinjectedRANSandLES,MScThesis,DelftUniversityofTechnology,2018. [9]S.E.Gano,J.E.Renaud,B.Sanders,HybridvariablefidelityoptimizationbyusingaKriging-basedscalingfunction,AIAAJ.43(2005)2422–2430. [10] R.Haftka,Combiningglobalandlocalapproximations,AIAAJ.29(1991)1523–1525,https://doi.org/10.2514/3.10768.

[11]Z.Han,SurroOptagenericsurrogate-basedoptimizationcodeforaerodynamicandmultidisciplinarydesign,in:30thCongressoftheInternational CounciloftheAeronauticalSciences,Daejeon,SouthKorea,2016.

[12] Z.Han,C.Xu,L.Zhang,Y.Zhang,K.Zhang,W.Song,Efficientaerodynamicshapeoptimizationusingvariable-fidelitysurrogatemodelsandmultilevel computationalgrids,Chin.J.Aeronaut.33(2020)31–473,https://doi.org/10.1016/j.cja.2019.05.001.

[13] Z.H.Han,S.Görtz,Hierarchicalkrigingmodelforvariable-fidelitysurrogatemodeling,AIAAJ.50(2012)1885–1896,https://doi.org/10.2514/1.J051354. [14]Z.H.Han,S.Görtz,R.Zimmermann, Improvingvariable-fidelitysurrogatemodelingviagradient-enhancedkrigingandageneralizedhybridbridge

function,Aerosp.Sci.Technol.25(2013)177–189.

[15] Z.H.Han,R.Zimmermann,S.Görtz,Alternativecokrigingmodelforvariable-fidelitysurrogatemodeling,AIAAJ.50(2012)1205–1210,https://doi.org/ 10.2514/1.J051243.

[16]S.Hickel,N.Adams,Efficientimplementationofnonlineardeconvolutionmethodsforimplicitlarge-eddysimulation,in:HighPerformanceComputing inScienceandEngineering,vol. 6,2007,pp. 293–306.

[17]S.Hickel,N.Adams,J.Domaradzki,Anadaptivelocaldeconvolution methodforimplicitLES,J.Comput.Phys.213(2006)413–436.

[18]S.Hickel,T.Kempe,N.A.Adams,Implicitlarge-eddysimulationappliedtoturbulentchannelflowwithperiodicconstrictions,Theor.Comput.Fluid Dyn.22(2008)227–242.

[19]S.Hickel,E.Touber,J.Bodart,J.Larsson,Aparametrizednon-equilibriumwall-modelforlarge-eddysimulations,in:Proceedingsofthe2012Summer Program,CenterforTurbulenceResearchatStanfordUniversity,2012,pp. 127–136.

[20]D.R.Jones,M.Schonlau,W.J.Welch,Efficientglobaloptimizationofexpensiveblack-boxfunctions,J.Glob.Optim.13(1998)455–492. [21]A.G.Journel,J.C.Huijbregts,MiningGeostatistics,AcademicPress,1978.

[22] M.L.Kaandorp,R.P.Dwight,Data-drivenmodelling oftheReynoldsstresstensorusingrandomforestswithinvariance,Comput.Fluids202(2020) 104497,https://doi.org/10.1016/j.compfluid.2020.104497.

[23] M.C.Kennedy,A.O’Hagan,Predictingtheoutputfromacomplexcomputercodewhenfastapproximationsareavailable,Biometrika87(2000)1–13, https://doi.org/10.1093/biomet/87.1.1.

[24] J.Ling,A.Kurzawski,J.Templeton,Reynoldsaveragedturbulencemodellingusingdeepneuralnetworkswithembeddedinvariance,J.FluidMech.807 (2016)155–166,https://doi.org/10.1017/jfm.2016.615.

[25]C.Mellen,J.Fröhlich,W.Rodi,Largeeddysimulationoftheflowoverperiodichills,in:Proceedings ofIMACS,WorldCongress,Lausanne,2000. [26] F.R.Menter,Two-equation eddy-viscosityturbulencemodelsfor engineeringapplications,AIAAJ.32(1994)1598–1605,https://doi.org/10.2514/3.

12149.

[27]M.Meyer,A.Devesa,S.Hickel,X.Y.Hu,N.A.Adams,Aconservativeimmersedinterfacemethodforlarge-eddysimulationofincompressibleflows,J. Comput.Phys.229(2010)6300–6317.

[28]M.Meyer,S.Hickel,C.Breitsamter,N.A.Adams,Wall-modelledimplicitlarge-eddy simulationoftheRA16SC1highliftconfiguration,in:31stAIAA AppliedAerodynamicsConference,AIAApaper2013-3037,AmericanInstituteofAeronauticsandAstronautics,SanDiego,2013.

[29] E.J.Parish,K.Duraisamy,Aparadigmfordata-drivenpredictivemodelingusingfieldinversionandmachinelearning,J.Comput.Phys.305(2016) 758–774,https://doi.org/10.1016/j.jcp.2015.11.012.

[30] S.Parneix,D.Laurence,P.A.Durbin,AprocedureforusingDNSdatabases,J.FluidsEng.120(1998)40–47,https://doi.org/10.1115/1.2819658. [31]S.B.Pope,Amoregeneraleffective-viscosityhypothesis,J.FluidMech.72(1975)331–340.

[32]C.Rapp,ExperimentellestudiederturbulentenStrömungüberperiodischeHügel,Ph.D.thesis,TechnischeUniversitätMünchen,2009.

[33] M.Schmelzer,R.P.Dwight,P.Cinnella,DiscoveryofalgebraicReynolds-stressmodelsusingsparsesymbolicregression,FlowTurbul.Combust.(2020), https://doi.org/10.1007/s10494-019-00089-x.

[34] R.L.Simpson,Junctionflows,Annu.Rev.FluidMech.33(2001)415–443,https://doi.org/10.1146/annurev.fluid.33.1.415.

[35]J.X.Wang,J.Wu,J.Ling,G.Iaccarino,H.Xiao,Physics-informedmachinelearningforpredictiveturbulencemodeling:towardacompleteframework, in:Proceedingsofthe2016SummerProgram,CenterforTurbulenceResearchatStanfordUniversity,2016,pp. 1–10.

[36] J.Weatheritt,R.Sandberg,AnovelevolutionaryalgorithmappliedtoalgebraicmodificationsoftheRANSstress–strainrelationship,J.Comput.Phys. 325(2016)22–37,https://doi.org/10.1016/j.jcp.2016.08.015.