Evaluation of scale-resolving simulations for a turbulent channel flow

Delft University of Technology

Evaluation of scale-resolving simulations for a turbulent channel flow

Klapwijk, M.; Lloyd, T.; Vaz, G.; van Terwisga, T.

DOI

10.1016/j.compfluid.2020.104636

Publication date

2020

Document Version

Final published version

Published in

Computers and Fluids

Citation (APA)

Klapwijk, M., Lloyd, T., Vaz, G., & van Terwisga, T. (2020). Evaluation of scale-resolving simulations for a

turbulent channel flow. Computers and Fluids, 209, [104636].

https://doi.org/10.1016/j.compfluid.2020.104636

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.

ContentslistsavailableatScienceDirect

Computers

and

Fluids

journalhomepage:www.elsevier.com/locate/compfluid

Evaluation

of

scale-resolving

simulations

for

a

turbulent

channel

flow

M.

Klapwijk

a,b,∗

_,

_T.

_Lloyd

b

_,

_G.

_Vaz

b,c

_,

_T.

_van

_Terwisga

a,b

a Delft University of Technology, Mekelweg 2, 2628 CD, Delft, The Netherlands

b Maritime Research Institute Netherlands, Haagsteeg 2, 6708 PM, Wageningen, The Netherlands c WavEC - Offshore Renewables, Edifício Diogo Cão, Doca de Alcântara Norte, Lisbon, 1350-352, Portugal

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 24 July 2019 Revised 19 March 2020 Accepted 19 June 2020 Available online 20 June 2020

Keywords:

Turbulent channel flow Verification DES LES PANS Synthetic turbulence

a

b

s

t

r

a

c

t

Differentvariableresolutionturbulencemodellingapproaches(Hybrid,BridgingmodelsandLES)are eval-uatedforturbulentchannelflowatRe τ=395,forcasesusingeitherstreamwiseperiodicboundary con-ditionsorasyntheticturbulencegenerator.Theeffectofiterative,statisticalanddiscretisationerrorsis investigated.ForLES,littledifferencebetweenthedifferentsub-filtermodellingapproachesisfoundon thefinergrids,whileoncoarsergridsILESdeviatesfromexplicitLESapproaches.TheresultsforHybrid modelsarestronglydependentontheirformulation,andthecorrespondingblendingbetweentheRANS and LES regions.The application ofPANS with differentratios ofmodelled-to-total kineticenergy, fk,

showsthatthereisnosmoothtransitionintheresultsbetweenRANSandDNS.Insteadacase-dependent thresholdwhichseparatestwosolutionregimesisobserved:fkvaluesbelow0.2yieldaproperturbulent

solution,similartoLESresults;higherfkvalues leadtoalaminarflowduetofilteringofthesmallest

scalesintheinverseenergycascade.Theapplicationofasyntheticturbulencegeneratorisobservedto yieldsimilarperformanceforallmodels.Thereducedcomputationalcostandincreasedflexibilitymakes itasuitableapproachtoenabletheusageofSRSforindustrialflowcaseswhichdependonthe devel-opmentofaturbulentboundarylayer.Itensuresthatsufficientlarge-scalestructuresdevelopoverthe fullboundarylayer height,therebynegating theproblemofrelyingonthe inverseenergycascadefor thedevelopmentofturbulence.BothLESandPANSwithturbulencegeneratoryieldabettermatchwith thereference datathanHybrid models;ofthesemethods PANSispreferable duetotheseparationof modellinganddiscretisationerrors.

© 2020TheAuthors.PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

The applicationofComputationalFluidDynamics (CFD)inthe maritime sectorismoving towards morecomplexproblems,such as massivelyseparated flows, blunt bodies, off-designconditions, cavitationandnoisepredictions.Forsuchcasestheassumptionsin traditional unsteady Reynolds Averaged Navier-Stokes (RANS) ap-proaches are too limiting, leading to an underprediction or ab-senceofcorrectunsteadyflowphenomena.Resolvingthefull spec-trumofturbulencebymeansofDirectNumericalSimulation(DNS) remains out of reach due to excessive computational cost, lead-ing to a focus on Scale-Resolving Simulations (SRS), where the larger scales are resolved, with the smaller scales modelled. The increase inavailablecomputationalpowermakes thispossiblefor highReynoldsnumberflows.InSRS,theaddedphysicalresolution

∗ _{Corresponding author.}

E-mail address: m.d.klapwijk@tudelft.nl (M. Klapwijk).

shouldleadtoamoreaccuratedescriptionoftheflowanda reduc-tionofthemodellingerroratamorereasonablecost.Thedecrease inmodellingerrorleadstoanincreaseinimportanceofnumerical errors,therebymakingerroranalysisevenmorerelevant.

SRScanbedividedintothreemaincategories:LargeEddy Sim-ulation (LES), for which an (implicit or explicit) filter is applied throughout the computational domain, resulting in scales larger than the filterbeing resolved and smaller scales being modelled using a ‘sub-filter’model. The need to resolve a substantial part oftheturbulencespectrumleadstostringentgridrequirementsin nearwall regions,andthereforeoftento excessivecomputational costforindustrialflow problems,whichtypicallyinvolvecomplex geometriesandReynoldsnumbers whichoftenexceed106 _in

hy-drodynamicapplications[1].Thecostcanbereducedthroughthe useofwallmodelledLES(WMLES),althoughthedefinitionofthe interface between wall model and resolved flow is difficult for complex geometries. This has lead to the rise of two alternative approaches:‘Hybrid’methods,whereLESiszonallycombinedwith https://doi.org/10.1016/j.compfluid.2020.104636

RANS,whichisapplied inregions wherethegrid cannotsupport LESresolution;and‘Bridging’methods,suchasPartiallyAveraged Navier-Stokes (PANS)[2]. Bridgingmethods consist of a blending ofRANSandDNS, but,incontrastto Hybridmethods,the blend-ingisnot location dependent.Insteadit dependsonuser-defined settings,therebyallowing asmoothtransitionbetweenthe turbu-lencemodellingapproaches.Thispreventscommutationerrors in thetransitionbetweenRANS andLESzones, asoccursforHybrid methods[3–5].InthecaseofthePANSmodeltheblenddepends onthemodelled-to-totalratiosofturbulencekineticenergy,f_k,and dissipation,f_.

Inallapproachestheratiobetweenresolvedandmodelled tur-bulencedependsonafilterlength.InLESandHybridmodelsthis filterlength isimplicitlydefinedbythegrid,i.e.refiningthegrid reduces the influence of the sub-filter model. Non-zonal Hybrid models have the additional issue that the blending function be-tween LES and RANS depends on the local grid density. Conse-quently,withgrid refinement notonly theeffect ofthe LES sub-filtermodelreducesinthe LESregion,butalsoalargerregion of theflowissolved usingLES.Themodellingerroristherefore en-tangled with, and forILES even dependent on, the discretisation error,leadingtoalargegridsensitivity.Thesepropertiesmake es-timating the discretisation error impossible, and grid refinement studiestoverifytheresultsdifficult.Bothoftheseareessentialto enableVerification andValidationprocesseswhich areneededto assessthecredibility ofindustrial CFDcalculations. Anadvantage ofBridgingmodelsisthatthefilterlengthissetexplicitly,thereby theoreticallydecouplingthediscretisationandmodellingerrors.Of course,whenperformingacomputationwithahighphysical reso-lutiononacoarse grid,thediscretisationerrorwillbecome dom-inant.Dueto theseproperties,Bridgingmodelsare becoming at-tractiveforindustrialCFD,whereextensivegridrefinementstudies areoftenunaffordablewhileanestimateofthenumericalerroris stillrequired[6,7].

Linked to the increased physics in SRS is the requirement of morephysicalinflowboundaryconditions.It hasbeenrepeatedly shown that the results of LES or DNS can be dependent on in-flow conditions, e.g. [8,9]. For SRS of turbulent flows which do notexhibitstrongseparatedvorticalstructures,itisnecessarythat theinflowcontainstime-varyingstochastic fluctuationswhich re-sembleturbulence.Ifthis isnot addressed,laminar solutionscan be obtained, and consequently integral quantities, such as mean forces,can beunderpredicted [10].Forattached, weaklyunsteady flows,unphysical laminarseparationmayoccur easily since tran-sitiontoturbulenceissignificantly delayed.Ironically,thisimplies thatforcomputingmeanforcesRANSmethodsoftenyieldthebest results forsuch flows. Hybrid methods also better predict mean forces,dueto theuse ofRANS inside the boundarylayer, result-ingfor exampleinsuperior prediction ofvelocity gradients close to thewall, when usingsteady inflow. However, atthe interface betweenRANSandLESregions notallthemodelledturbulenceis transferreddirectly intoresolved turbulence,leadingto anoverly laminarflowfield. Insuch cases,inflowturbulencemightstill be necessary. While this problemhas been known for a number of years,manypublicationsusingSRSstill donotapply aproper in-flow.Theseworksoftenfocusonturbulentstructuresand dynam-icswhichappearsatisfactory,yet theforcessimultaneously show alarge mismatchwithexperimental data.Thisdiscrepancyin in-tegralquantitiesleadstodifficultiesinthevalidationofnumerical results,anddeterioratesthecredibilityofSRSforpractical applica-tions.Currently,theneedforsyntheticturbulencegeneration ham-perstheusage ofSRS forindustrial cases,such ascavitation and noiseresearchforshippropellers.

Inflow turbulencecan be generatedeitherby precursor meth-ods(such asachannel flow), orby synthetic methods,which do notrelyonflow recycling.Precursormethods aregenerallynoted

for their accuracy, although they are more expensive to usedue totheneedtogenerateturbulenceina second,separate,domain. Theiruseisalsooftenlimitedtocanonicalflows.Synthetic meth-odsarecheapertouse,easiertotunetoadesiredsetofturbulent inflowstatistics,andmoregenerallyapplicable [11].These proper-tiesmakeasyntheticmethodpreferableforindustrialCFD,despite therequirementofadevelopmentlengthtoallowthe introduced fluctuationsintodevelopto‘real’turbulence.Inthecaseof cavita-tionandnoiseprediction,afurtherrequirementisthattheinflow velocityfluctuationsaredivergence-free,therebyavoidingspurious pressurefluctuationswhich canpollutethe entiredomainfor in-compressiblecomputations.

Itiswell knownthatthechoice ofSRSmodelaffectstheflow prediction around the object of interest, although a subject less well addressedintheliteratureis themodelinteractionwiththe propagation of inflow turbulence. Consequently, an evaluation of different SRS methods in combination with inflow turbulence is necessary. This paper attempts to provide a systematicoverview of the effect of SRS turbulence modellingapproaches, both with streamwise periodic boundary conditions and with a synthetic turbulencegenerator. Full comparisonsbetweendifferentSRS ap-proaches, including higher order moments, are rare in literature

[12]; often, different codes, grids and solver settings are used, whichmakes it difficultto assessthe modellingerrorof the tur-bulence modellingapproaches.Models includedinthispaper are Hybrid (DDES [13], IDDES [14]) and XLES [15]); Bridging (PANS

[2] withfk∈[0.05, 1.00]);andLES(Smagorinsky [16],Lilly[17,18],

WALE [19], KSGS [20] and ILES).Attention is paid to the assess-mentofiterative,discretisationandstatisticalerrors.

ThechosentestcaseisaturbulentchannelflowatRe_τ=395_,a canonicaltest caseforthestudyofwall-boundedturbulence,due tothe simplegeometryandabundance ofreferencedata. Experi-mentalresultsarefirstpublishedninety yearsago,andnumerical studieshavebeen performedusingLES andDNSsince the1980s.

Fig.1 showsanoverviewofnumericalresultsavailableintheopen literature.Itisobservedthata rangeintotalnumberofgridcells is applied for LES/DNS approaches atlower Re_τ values. DESand PANSresultshavebeenpublishedforsignificantly higherRe_τ,but oftenwithout refining the grid.Mostresults shownare obtained using a finite volume approach, with second-order accurate dis-cretisation schemes. Fig. 1 also shows an overviewof the usage ofsyntheticturbulence generationmethods available intheopen literature,sortedpertestcase.

In the paper, Section 2 describesthe different turbulence ap-proaches used, and Section 3 describes the test case, numerical setup and turbulence-generatingmethod. After an assessment of the numericalerrors in Section 4,Sections 5 and 6 compare the resultsforHybrid,BridgingandLESapproachesusingstreamwise periodic boundary conditions and a synthetic turbulence gener-ator respectively. Finally, Section 7 discusses the implications of theresultsforindustrialtestcases,followedbytheconclusionsin

Section8.

2. Mathematicalturbulenceapproaches

Inthisworkdifferentapproachesformodellingturbulenceare applied,includingHybridandBridgingmodelsandLES.Themain equationsof theseapproachesare describedinthis section,with fullmodeldetailsgiveninAppendixA.

Forall approachesconsidered theinstantaneous quantities,



, can be decomposed into a resolved,







, anda modelled (unre-solved) component,

φ

, according to



=







+

φ

[46]. Applying thisdecompositiontotheincompressible,single-phase,Newtonian Navier-Stokesequationsleadsto

∂



U i



∂

x i

Fig. 1. Literature overview with a selection of the available numerical results. Turbulent channel flow results (left) as function of the number of grid cells, N c , and the

turbulence modelling approach (left) [21–35] . Usage of synthetic turbulence generation (right) sorted by test case and method (Digital Filtering (DFM), Forward Stepwise (FSM), Fourier (FM) and Synthetic Eddy method (SEM) [8,9,36–45] .

∂



U i



∂

t +



U j



∂



U i



∂

x j =− 1

ρ

∂

∂



x Pi



+

∂

∂

x j



ν



∂



U i



∂

x j +

∂



U j



∂

x i



+1

ρ

∂

τ

(

u i, u j

)

∂

x j . (2)

IntheseequationsU_idenotesthevelocitycomponents,Pthe pres-sure,

ν

thekinematicviscosity,

ρ

thedensityand

τ

(ui,uj)the

sub-filterstresstensorwhichismodelledusingBoussinesq’s hypothe-sis,

τ

(

u i, u j

)

=



U iU j



−



U i



U j



=2

ν

t



S i j



−

2

3k

δ

i j, (3)

with

ν

t theeddy-viscosity, kthe modelledturbulencekinetic

en-ergy,

δ

ij theKroneckerdeltaand



Sij



theresolved strain-rate

ten-sor,definedas



S i j



= 1 2



∂



U i



∂

x j +

∂



U j



∂

x i



. (4)

The differencebetweentherespective approachesliesin the def-initionofthefilteringoperator



·



,whichistemporalinthecase ofRANSandspatialforLESandPANS.Consequently,thesub-filter stresstensorismodelleddifferentlybyemployingdifferent expres-sionsfor

ν

tandk.

2.1. RANS

In unsteady RANS, all turbulent time scales and motions are modelled, and a scale separation between deterministic and stochasticscalesisassumed.Inthiswork,thesub-filterstress ten-sorismodelledusingthek−

ω

ShearStressTransport(SST)model (2003version)[47].Thismodelformsthebasis oftheHybridand Bridgingmodels.Itisa blendofak−

ω

modelinthewallregion andak−



modelinthefarfield.Thetransportequationsare

Dk Dt =P k−

β

∗

ω

k +

∂

∂

x j



(

ν

+

ν

t

σ

k

)

∂

k

∂

x j



, (5) D

ω

Dt =

α

ν

t Pk−

βω

2+

_∂

∂

x j



(

ν

+

ν

t

σ

ω

)

_∂

∂

ω

x j



+2

(

1− F1

)

σ

_ω

ω2

∂

k

∂

x j

∂

ω

∂

x j. (6) Theotherrelevantexpressionsare

P k=min

(

ν

t

|

S

|

, 10

β

∗k

ω

)

and

ν

t=

a 1k

max

(

a 1

ω

,

|

S

|

F 2

)

.

(7)

Here



|S|



isthemagnitudeoftheresolved strain-ratetensor.The transport equations make use of two blending functions, F1 and

F2,andtwolimiters. Fordetailsofthesefunctionsandthemodel

constantsseeSectionA.1. 2.2.Hybridmodels 2.2.1. DDES

DelayedDetachedEddySimulation(DDES)[13] isanadaptation oftheDetachedEddySimulation(DES)[3] model.Theswitch be-tweenRANSandLESisbasedontheturbulentlengthscaleltand

ashieldingfunctionfd.Thektransportequationbecomes

Dk Dt = Pk− k 3/2 l t +

∂

∂

x j



(

ν

+

ν

t

σ

k

)

∂

k

∂

x j



, (8)

includingtheturbulentlengthscale

l t= l tRANS− fdmax



l RANS t − ltSRS, 0



. (9)

TheRANSandSRSlengthscalesaredefinedas

l RANS t =

√

k

β

∗

_ω

and l SRSt =C DDES



, (10)

inwhich



isthemaximumcelllength.Thecoefficientsand addi-tionalshieldingfunctionaregiveninSectionA.2.

2.2.2. IDDES

InDDESthetransitionfromRANStoLESinside awall bound-arylayer isprohibited,which leadstoa lackof resolvedvelocity fluctuationsclosetothewall.Thiscansuppressdynamicsinflow phenomenaclosetothewall,suchassheetcavitation.Adifferent HybridmodelisImprovedDelayedDetachedEddySimulation (ID-DES)[48], whichaims to useDESas aWall-Modelled LES (WM-LES),whileemployingRANS inthenearwallregion insteadofan analytical expression. Consequently IDDES is better able to han-dleseparatingflows,sinceLESisalsoallowedinsidetheboundary layer.A wall-normalresolutiony+=u_τy/

ν

<1atthewall isstill required.InIDDEStheblendingisachievedbyadifferentshielding functionf_dt.Notethatduetothedifferentformulation,themodel isprohibitedtoswitchtoRANSintheboundarylayerandfarfield. Thiscanthereforedeterioratetheresultsonacoarsegrid.The co-efficientsandauxiliaryfunctionsaregiveninSectionA.3.

2.2.3. XLES

Thefinal Hybridmodelconsideredinthisstudyis eXtra-Large EddySimulation (XLES)[15].It issimilartoDDES andIDDES,but withaswitchingfunctionwhichisnotdependentonthewall dis-tance(whichcanbecomputationallyexpensiveanddifficultto de-fineforacomplexgeometry).Furthermore,inLESmodeadifferent

sub-filterstress modeltoDESisused (the KSGSmodel).The dif-ferencebetweenRANSandLESmodesliesinthedefinitionofthe eddy-viscosityanddissipation.ForRANSthesearegivenas

ν

t = l



k and



=

β

k

k 32

l , (11)

withl=√k/

ω

,whilefortheLESsub-filtermodel

ν

t =C 1





k and



=C 2

k 32



, (12)

with



definedasthe maximumcell length. TheRANS model is closed witha modifiedequation for

ω

, based on the k−

ω

TNT model[49], D

ω

Dt = P ω−

β

ω

ω

2₊

σ

d

ω

max



∂

k

∂

x i

∂

ω

∂

x i, 0



+

∂

∂

x j



(

ν

+

σ

ω

ν

t

)

∂

ω

∂

x j



, (13)

withaproductionterm

P_ω=

α

ω



S



2_{. (14)}

TheswitchbetweenRANSandLESmodeismadeusinga compos-itelengthscalel˜definedas

˜

l =min

(

l, C 1

)

. (15)

ThecoefficientsandauxiliaryfunctionsaregiveninSectionA.4. 2.3.Bridgingmodel

ThecombinationofRANSwithLESinHybridmodelsimproves accuracybutmayleadtocommutationerrorsinthetransition be-tweentheRANSandLESzones.Anapproachwithoutcommutation errorsis the Bridgingfamilyof models, such asthe Partially Av-eragedNavier-Stokes (PANS)model[2].The PANSmodelis based onspatially filtering the Navier-Stokes equations,wherethe sub-filterstress tensor ismodelled using a setof reformulatedRANS equations including the modelled-to-total ratio of turbulence ki-neticenergyanddissipationrate,

f k= k K and f ω=

ω



= f _ f k . (16)

The PANS model in this work is based on the k−

ω

SST model

[47,50],withthetransportequations Dk Dt =P k−

β

∗

ω

k +

∂

∂

x j



ν

+

ν

t

σ

k f _ω f k



∂

k

∂

x j



, (17) D

ω

Dt =

α

ν

t P k−



P− P f _ω +

βω

f _ω



ω

+

∂

∂

x j



ν

+

ν

t

σ

ω f _ω f k



∂

ω

∂

x j



+2

σ

ω2

ω

f _ω f k

(

1− F1

)

∂

k

∂

x j

∂

ω

∂

x j (18) with P=

αβ

_ν

∗k t and

ν

t= a 1k max

(

a 1

ω

,



S



F 2

)

. (19) FortheauxiliaryfunctionsandconstantsseeSectionA.5.

In PANS,the closureisviscosity-based(the sub-filter viscosity is a function of the modelled flow field (k,

ω

)), whereas in LES theclosure is grid-based(sub-filter viscosity is a function ofthe cut-off lengthscale(



))[51].Consequently,incontrasttoLES,the cut-off lengthscale of the resolved flow isnot predetermined in PANS.The physicalresolutionis onlydetermined by the settings, whichleads toan overlapbetween thePANS resolved and mod-elledspectra[51].Sincethegridremainsfixed,computationswith

anfklargerthanwhatthegridallowsarecomparabletoan

explic-itly filtered LES (althoughbased on a different modelling frame-work). Computations with a very low fk value are effectively an

ImplicitLES(under-resolvedDNS).

The filteringdepends onlyon the value offk andf.fk

deter-minesthephysicalresolutionoftheflow,f_determinestheoverlap betweentheenergy-containingandthedissipationranges.The ef-fectofmodifyingf_isinvestigatedinKlapwijketal.[52],inwhich itwasconcludedthat f_ shouldbe keptequalto 1.0(alsoknown asthe‘highReynoldsnumber’approach),toavoidaddingexcessive dissipation. fk and f can eithervary inspaceand/or time,or be

keptconstantthroughoutthedomainandcomputation.The draw-backofthefirstapproachisthatthisentanglesthemodellingerror andnumericalerror, andthereby destroys oneof thekey advan-tagesofthePANSmodel.Althoughthisapproachisgaining popu-larityinliterature,howeverthereisnoconsensusonhowfkshould

be (apriorior dynamically)chosen [53]. Throughoutthiswork a constantvalueforfkisused.

2.4. LES

Large EddySimulation (LES) is founded on the principle that the larger scales of turbulence are resolved, while the smaller scalesaremodelled.Inordertoenablethis,thefilteringoperation isperformedspatially;eitherexplicitlyorimplicitly(on thegrid). Torelatethesub-filterstresstothefilteredstrainrate,anumberof modelscan beemployed. Inthiswork,a selectionismadebased onthemostcommonmodelsfoundintheopenliterature. 2.4.1. Smagorinsky

TheSmagorinskymodel[16] modelstheeddy-viscosityas

ν

t=

(

C s

)

2

|

S

|

, (20)

with

|

S

|

=



2



Si j



Si j



andamodelconstantCs.Themodel

con-stant depends on the flow and in literature values in the range 0.065− 0.23arefound.Inthispaperavalueof0.10isapplied.The filteringisdoneimplicitly,usingthefilterasdefinedby Smagorin-sky:



=

(



x



y



z

)

1

3. ₍₂₁₎

In the case of highly anisotropic cells, which occur often near walls, this filter is too optimistic leading to an underpredicted eddy-viscosity.

2.4.2. Lilly

Tocircumventdifficultiesinobtainingageneralconstantforthe Smagorinskymodelandtoimprovebehaviournearwalls,Germano et al. [17] suggested using a constant which varies in time and space, thereby adapting to the resolved scales. Lilly [18] applied aleast-squaresestimatetoobtaintheconstant,knownastheLilly Model.Inthisapproach,nexttothespatialfilterasecond,coarser, filterisapplied. Thisfilter,knownasthe ‘test’filter,indicated by ·,isusuallydefinedas



=2



.Thesub-testfilterstress,Tij,is

de-finedanalogouslytothesub-filterstress,thatis

T i j− 1 3

δ

i jT kk=2



C s





2_

|

S

|



S i j



. (22)

Themodelconstantisobtainedusing

C s= 1

2

L i jM i j

M i jM i j,

(23) inwhich theerrorbetweenthe resolvedscales ofmotion Lij and

thelocal closureM_ij isminimised.L_ij isdefinedasthe difference betweenthesub-filterstressonthenormalfilter(

τ

ij)andonthe

test-filterlevel,

andMijisdefinedas

M i j=



2

|

S

|



S i j



−



2

|

S

|

S i j



. (25)

Adownsideofthismodelisthatthemodelrequiresspatial averag-ingoftheconstanttoreducethevariabilityinspaceandtime,else the value of Cs can become either negativeorunphysically large

[54].Sincethisis oftennotpossible forindustrialflowcases due to theabsenceofflow homogeneity,theconstant isbounded be-tween0and10timestheupperlimitoftheSmagorinskyconstant asfoundinliterature,soCs∈[0,2.3].

2.4.3. WALE

The Wall Adaptive Local Eddy-viscosity (WALE) model by NicoudandDucros[19] wasproposedforLESincomplex geome-tries to account for the effects ofthe strain androtation rate of the smallest resolved velocity fluctuations. It should alsorecover proper near-wall scaling for the eddy-viscosity without dynamic procedures.Thesub-filterviscosityisdeterminedas

ν

t= L 2s



S d i jS di j



3/2





S i j



S i j





5/2 +



S d i jS di j



5/4, (26)

basedonthesquareofthevelocitygradienttensor

S d i j= 1 2





g i j



2+



g ji



2



−1₃

δ

i j



g kk



2 (27) with



g i j



=

∂



U i



∂

x j . (28) Thelengthscaleisgivenas

L s=min

(

κ

d, C s

)

(29)

inwhichdindicatesthewalldistanceand

κ

=0.41. 2.4.4. KSGS

Thefinalsub-filtermodelistheksub-gridStress(KSGS)model

[20],for whicha transport equation forthe sub-filterturbulence kinetic energy(Eq.5) issolved,withtheproductiontermdefined as

P k=

ν

t



S



2. (30)

Theeddy-viscosityanddissipationtermsaredefinedinEq.12. 2.4.5. ILES

An alternative approach to LES sub-filter modelling is known asImplicitLES (ILES)[55].Instead ofapplyingasub-filter model, it is assumedthat the added numerical diffusiondue to theuse of coarse(r) grids andlow order (upwinding) schemes acts as a sub-filter model. Whilst an attractive approach due to low com-putationalcost(nosub-filtermodelisrequired),cautionisneeded inemployingthisapproach,especiallyinresolvingstructuresnear thepresumedcut-off scale[56].Duetotheabsenceofasub-filter model, the only contributing factor in the ratio of modelled-to-totalkineticenergyisthegrid.Asaconsequencemakingaproper grid becomes even more important than usual; the reliability of ILES for industrial flow cases on highly non-uniform, anisotropic grids canleadto largeerrors.Thefilterlengthvariessignificantly inthedomain,andduetotheanisotropiccells,thenumerical dis-sipation due to the convection scheme varies in different direc-tions.Secondly,theapplicationofloworder(upwinding)methods usingcoarsegridscanleadtotoomuchaddeddissipation,thereby not capturing finer structures [56]. Finally, the results are even moregriddependentthanfornormalLES,whichmakesproper so-lutionverificationimpossible.Neverthelessitisawidelyemployed approach,andisthereforealsoincludedinthiswork.

2.5.Turbulencegeneratingmethods

Asshownintheliterature,theuseofscale-resolvingturbulence methods for a turbulent channel yields a so-called supercritical laminarsolutionforwhichmanyflow-throughtimesareneededto triggertransitiontothe turbulentregime [31].To thisendinthe currentwork twomethods are employed tospeed up the transi-tion.Forcaseswithstreamwiseperiodicboundaryconditions(see

Section 3)the methodsuggestedby Schoppa andHussain[57] is usedasreportedinKlapwijketal.[58].Thismethodisonly appli-cableto caseswithstreamwiseperiodicboundaryconditions,and istunedforaturbulent channelflow.Itisthereforenotageneral approach.

Secondly,a more general method,which doesnot depend on streamwiseperiodicboundaryconditions,isapplied.Synthetic tur-bulenceis generatedusinga modifiedversionofthe digitalfilter methodbyXieandCastro[36].Thismethodismoregeneral,and doesnotrequireflowrecycling.Inthecurrentimplementationthe methodisabletogenerateanisotropic,inhomogeneousturbulence whichissufficientformostindustrialapplicationssuchastheflow around wings, propellers, ship hulls, etc. In the method random numbers,rm,l,i,withzeromeanandunitvariance,aregeneratedon

aCartesiangridateach timestep.Herem,l indicatetheposition indicesandithevelocitycomponent.Thesenumbersarespatially correlatedusing

ψ

m,l,i= N  j=−N N  k=−N b jb kr m+j,l+k,i, (31)

afterwhichtheyaretemporallycorrelatedwiththenumbers gen-eratedattheprevioustimestepusing



i

(

t

)

=



i

(

t −



t

)

exp



−

π



t 2T



+

ψ

i

(

t

)



1− exp



−

π



t 2T



. (32) HereT=Ii/UiistheLagrangiantimescaleandbjandbkarefilter

coefficientsused togenerate spatialcorrelations, andare defined as b j= b _j



_Nj l=−Nj b 2 l



and b k= b _k



_N k m=−Nk b 2 m



, (33) with b _j=exp



−

π|

j

|

2n



and b _k=exp



−

π|

k

|

2n



. (34) The spatially andtemporallycorrelated numbers aretransformed tovelocityfluctuationsusing

U _i=a i j·



i (35)

in which aij indicates the Lund transformation matrix, which is

basedona Choleskydecompositionofthe Reynoldsstress tensor Rij[10], a i j=

⎡

⎣



R 11 0 0 R 11/a 11



R 22− a221 0 R 31/a 11

(

R 32− a21a 31

)

/a 22



R 33− a231− a232

⎤

⎦

. (36) In the current work a string of pseudo-random numbers is em-ployed, i.e. the same range of random numbers for every com-putation.For more details the readeris referred to Xie and Cas-tro [36] and [39]. In these works the velocity is modified di-rectly inside the first non-linear PISO loop, either at the inflow, orinthedomain.Inthecurrentworkthevelocityfluctuationsare

transformedtoabody-forceterminthemomentumequations, ex-plicitlyaddedontherighthandsideoftheequations.Thisisdone toimproveiterativeconvergenceandmassconservation.The trans-formationisachievedusing

Fb,i=



U _{i,in f low}+U _i− Ui



ρ

U i L tg · b (37) where U_i,inflow is the mean velocity as defined at the generator plane,U_icomesfromEq.35 andUiistheinstantaneousvelocityin

acell,obtainedfromthesolveratthecurrentnon-linearloop.Ltg

indicatesthe distanceintheflow directionover whichthe body-forcetermisapplied,andbisanarbitrarymultiplicationfactorto increasetheconvergenceofthevelocity towardsthedesired fluc-tuations. For now b is takenas 300. Notice that this body-force termgoestozerowhenthelocalvelocityequalsthedesiredmean plusfluctuation.

3. Numericalsetupandsolver

Computations are madeusinga rectangulardomain,withtwo no-slip walls oriented normal to the y-direction and periodic boundary conditions in spanwise direction. (see Fig. 2). The re-maining boundaries are either connected using periodic bound-aryconditionsinordertoapproximatean infinitechannel;or al-ternatively, an inflow and outflow boundary condition is speci-fiedifasyntheticturbulencegeneratorisapplied. Cartesiangrids withhyperbolic tangent clusteringtowards thewallsare used,as describedin Section 4.3. The non-dimensional time step forthe gridused for comparing turbulencemodelling approaches, G4, is



t∗=1/2u_τ



t/

δ

≈ 1× 10−3_{. This leads to}

_

_t+_=u2

τ



t/

ν

≈ 0.08

and a maximum Courant number below 0.2 (2000 time steps perflow-through time).To maintaintheproper bulk andfriction Reynolds numbers, Re_b=U_b2

δ

/

ν

and Re_τ=u_τ

δ

/

ν

respectively, a body-forceisappliedwhichisproportionaltothestreamwise pres-suregradientdp/dx=−

τ

w/

δ

,with

τ

w=

ρ

u2_τ [59].

The numerical solver used for all simulations in this work is ReFRESCO[60],amultiphaseunsteadyincompressibleviscousflow solver using RANS and Scale-Resolving Simulation models, com-plemented with cavitation models and volume-fraction transport equationsfor differentphases. Forthe simulations reported here worktime integrationisperformedusing asecond-order implicit threetimelevelscheme,andtheconvectiontermsinthe momen-tumequationsare discretisedusinga second-orderaccurate cen-tral differencing scheme (the Péclet number has a magnitude of O

(

10

)

).An investigationintothe effectoftheconvection scheme for the momentum equation can be found in Section 4.3. The turbulence equations are discretised using a first-order upwind

scheme.Duetothefinegridresolutionemployed,theeffectof dis-cretisationonthesub-filterstressesislimited.

4. Numericalerrors

As generallyaccepted,numerical errorscan be divide into in-put,round-off,iterative,discretisation,and,inthecaseofunsteady computations,statisticalerrors[61].Theinputerrorisassumedto benegligibleduetothewellcontrolledconditionsinthe computa-tions,andeithertheperiodicnatureofthesolution,orthe reliabil-ityoftheDNSresultsreportedintheopenliterature. The round-off errorisnegligibleduetotheuseofdoubleprecisionarithmetic. Therefore,inthiswork,onlytheiterative,discretisationand statis-ticalerrorareassessed.

Theresultsareobtainedalongalineperpendicular tothewall at the centreof the channel (see Section 4.2). Computations are performed usinggrid G4 (see Section 4.3) and using PANS (fk=

0_.10) unlessindicated otherwise.Only half the channel height is shown.

4.1. Iterativeerror

The iterative convergence is assessed based on the residuals, normalised by the diagonal element of the left-hand-side ma-trix of the linear system of equations. This is of particular in-terest since ifthe iterative errorwould be ofthe same order of magnitude as the turbulence fluctuations, the results would be strongly affected. Despite this, the influence of iterative error is rarely studied inthe open literature. Following the approach ad-vocatedbyEçaetal.[62],aPANScomputationwith fk=0.10was

performedusingdifferentiterativeconvergencecriteria(L2=10−3,

10−4,10−5,10−6,10−7and 10−8).The effect on the mean velocity (u+),Reynoldsstresses(Rei j=uiuj/u2τ) andturbulencekinetic

en-ergyspectra(Eu(f)aty+≈ 20)alongawall-normallineatthe

cen-tre ofthedomainis shownin Fig.3.ForvaluesofL2=10−3 and

10−4, the mean velocity showsan underprediction in the buffer layer(5<y+<30), whileforL2≤ 10−5theresultsare converged.

TheReynoldsstressesandspectraalsoshowalargemismatchwith the referencedatafor L2=10−3 and10−4. Themagnitude ofthe

peakvalueReuuandtheturbulencekineticenergyspectraconverge

for stricter convergence criteria. As a compromise between cost andaccuracy,thecriteriumL2=10−6 isusedintheremainderof

thiswork.ApplyingthiscriteriumleadstoaresidualofL_∞=10−5 in each time stepfor momentum; theresiduals for pressureand turbulenceequationsareatleastoneorderofmagnitudesmaller. In thismanner, theiterativeerror issmaller than theturbulence fluctuationsofinterest.

Fig. 2. Schematic overview of the domain and physical parameters. The dashed lines indicate the computational domain. The figure is based on the drawing of de Villiers

Fig. 3. Mean velocity ( u +_{) profiles, Reynolds stress ( Re}

uu and Re uv ) profiles and turbulence kinetic energy spectra ( E u,y+_≈20(f) ) using different iterative convergence criteria.

4.2. Statisticalerror

A potentiallydominating errorinunsteady andespecially tur-bulence resolving simulations is the statistical error. In order to remove the start-up effects and estimate the magnitude of this source of uncertainty, the Transient Scanning Technique (TST) is employed [63].Thistechniqueallows an estimationof the statis-ticaluncertaintybasedonasignaloffinitelength.Theuncertainty is expanded to obtain a 95% confidence interval. The TST is ap-pliedtothevelocitysignalsandReynoldsstressesatmeasurement points along the height of the channel. Based on the TST, it is found that the first 11 flow-through times(≈ 30000 time steps) mustberemovedtoeliminatethestart-upeffects.Thisconclusion is independent of the wall-normal distance of the measurement point. In the remainder of this work, the mean values are then computedbased onapproximately 45flow-through times, result-ing in a statistical uncertainty for the mean streamwise velocity below2%,andfortheReynoldsstresscomponentsbelow10%.

Theseuncertaintiesareinagreementwiththe samplingerrors as obtained from the engineering approaches suggested by Ries et al.[64]. The estimates ofthe sampling error efor u+_, u and u2

rms aregivenas:

e =



2

δ

I 2 U bt av, e =



δ

U bt av and e =



4

δ

U bt av. (38)

HeretavindicatestheaveragingtimeandItheturbulenceintensity

I=u/u. Theseestimates give samplingerrors of 1%,5% and10%

foru, u and u2 _{respectively when applied to theresults of the}

presentstudy.

Inordertoreducetheseerrorsmoreflow-throughtimesshould becomputed. Itmustbenotedthatinliteraturespatialaveraging isoftenappliedtotheresultssincetheflowisstatistically homo-geneous[64].Inthiswaya lowstatisticaluncertaintyisachieved using fewer flow-through times. However this implies that the flowisstatisticallyconvergedintheentiredomain,butnotat ev-erylocation.Inthisworkspatialaveragingisexplicitlynotapplied, toensureastatisticallyconvergedsolutionatall locations,andto properly compare with results obtained with a synthetic turbu-lencegenerator.

4.3.Discretisationerror

In orderto assess the discretisation error, four differentgrids (withrefinementratiosri=hi/h1=



ti/



t1=1.00,1.25, 1.57and

1.97)wereemployed.DetailsofthegridsaregiveninTable1.The griddesignatedG4 isequaltotheoneusedfortheDNSreference data[28].Notethat allthesegrids arewellwithinLESguidelines, interms of wall resolution, found in literature, andhave resolu-tionstypicalofDNS[56,65].Itiscommonlyassumedinliterature, whengridswithDNSresolution,inconjunctionwithsecond-order schemes,areusedforLES,thatdiscretisationerrorsare negligible anddonothavetobeassessed.However,discretisationerrors de-pendonboththenumberofcellsandtheaccuracyoftheschemes employed,andassessmentoftheseerrorsisstillnecessary.

Table 1

Details of used grids. r i ( = h i /h 1 = ti / t1 ) indicates the refinement ratio, N the number of

cells in different directions, y1 the initial wall-normal spacing, x + = u τx/ν, y + = u τy/νand

z+ = u _{τz/νthe average non-dimensional wall units in different directions and the maximum}

Courant number, Co max .

Grid ri Nx Ny Nz N /10 6 y1 · 10 4 x+ y+ z+ Comax G6 3.97 63 47 47 0.14 2.1 24 0.2 20 0.07 G5 2.63 95 71 71 0.48 1.4 16 0.13 13 0.12 G4 1.95 127 95 95 1.18 1.0 12 0.10 10 0.20 G3 1.56 159 119 119 2.30 0.8 9.6 0.08 8 0.21 G2 1.25 199 149 149 4.50 0.6 7.7 0.06 6.4 0.22 G1 1.0 249 187 187 8.84 0.5 6.1 0.05 5.1 0.25

Fig. 4. Mean velocity ( u +_{) profiles, Reynolds stress ( Re}

uu and Re uv ) profiles and turbulence kinetic energy spectra ( E u,y+_≈20(f) ) using different grids and PANS ( f _{k = 0 . 10 ).}

TheeffectofgridrefinementisshowninFig.4.Boththemean velocityandReynoldsstressesappearreasonablyinsensitivetogrid resolution, however Reuu and Reuv deviate slightly on the finest

grid.The maindifferencesare observedfortheturbulencekinetic energyspectra.Grid refinementleadsto a slightlyincreased cut-off frequency,since the smaller cellsallow for higher wavenum-bersto be resolved. This indicates that the employed f_k (0.10) is belowthegrid cut-off,i.e.thegrid andnumericalsettingsdonot resultinDNSresolution.Nevertheless,basedonthesimilarity be-tweentheresultsitisconcludedthatthecoarsestgridhasa suffi-cientresolution.Forthisreason,gridrefinementfortheLES mod-els was not pursued, since due to the fine grid resolution, little effectofthesub-filtermodelisexpected.Instead,to comparethe effectofgrid resolutiononthesub-filtermodelling,grid coarsen-ingwasperformed.Tothisend, twoadditionalgrids (with refine-mentratiosri=2.63and3.97)areemployed,incombinationwith

PANS (f_k=0_.10), LES KSGS andILES. Fig. 5 showsthat even the coarsestgridstill hassufficientresolutiontopredictthemean ve-locityprofiles well, in the casesof PANS andILES.ForReuu both

ILES andPANS overpredict the peak near the wall, especially on coarser grids. This overprediction is absent for LES KSGS,due to thesub-filtermodel.Onfinergridshowever,thepeakis underpre-dicted,indicating thattheSGSistoodissipative.Thepeak onthe coarsestgridisalsoshiftedawayfromthewall.ForReuv,againLES

KSGSpredictsthemagnitudebetteronacoarsegridthanILESand PANS. Thisdifference is absent onfiner grids. Interms of turbu-lencekinetic energyspectra,itisclearthat gridrefinementleads toan increaseinenergyathigherfrequencies,i.e.thecut-off fre-quencyincreases.Theeffectofthesub-filtermodelisalsoclearfor LESKSGS,thecut-off onallgridsbeingatalowerfrequencythan forILES andPANS.This effectis mostsignificant onthe coarsest grid.

Fig. 5. Mean velocity ( u +_{) profiles, Reynolds stress ( Re}

uu and Re uv ) profiles and turbulence kinetic energy spectra ( E u,y+_≈20(f) ) using different grids.

Toconclude, thegriddesignated G4hasa sufficientresolution tobeabletocomparethedifferentmodels.Gridcoarseningshows that,although theresultsobtainedbyILESappearreasonable,the absenceofasub-filtermodelcanleadtoan overpredictionof tur-bulentstressesoncoarsegrids.Inaddition,similarityinresults in-dicates that PANS with a low fk is comparable to LES without a

sub-filtermodel(ILES).

A second source ofdiscretisation error is thediscretisation of theconvection termsintheequations.Dueto theusageofa un-structuredFiniteVolumeCFDcode,wearelimitedtosecond-order accurateschemes.FortheReynoldsnumberstudiedhere,the dis-cretisation ofthemomentumequationsisdominant.Basaraetal.

[66] investigated blended upwind-CD schemes for Finite Volume LES andHybrid methods.Theystate that‘resultsobtainedwitha

blendingfactorlowerthan 0.98or0.96aretreatedassuspicious’, andthat the use ofa second-order accurate Central-Differencing schemeispreferred.HoweverinindustrialhighReynoldsnumber flows this is often not possible, meaning that (lower-order) up-winding, or blended schemes are used [66]. In order to investi-gatethis effect,theconvection schemeforthe momentum equa-tion is varied between First Order Upwind (FOU),Central Differ-encing(CD), a blendedupwind-CD witha blending factorof 0.5 (FOU-CD) and the approximately second-order QUICK (Quadratic UpstreamInterpolationforConvectiveKinematics)scheme.The in-vestigationis performed for ILES,LES KSGS and PANS with f_k= 0.10.Thecomputationsare deliberatelyperformedonthe coarser grid G6 to highlight the difference between different convection schemesanddifferentsub-filtermodels[66].Thesedifferencesare

Fig. 6. Comparison of normalised mean velocity ( u +), Reynolds stress ( Re uu ) profiles and turbulence kinetic energy spectra ( E u,y+_≈20(f) ) using different convection schemes for the momentum equations.

expectedto be smaller onfiner grids. The non-dimensional time step



t∗=u_τ



t/2

δ

≈ 2× 10−3 _leadsto

_

_t+_=u2

τ



t/

ν

≈ 0.5and amaximumCourantnumberof0.1.

Fig. 6 shows thenormalised meanvelocity (u+) and Reynolds stressprofiles(ReuuandReuv),togetherwiththeturbulencekinetic

energyspectra E(f). In terms of velocity profiles, the CD scheme yieldsthebestmatchwiththereferencedataforallturbulence ap-proaches.TheQUICK schemecapturesthetrend,butoverpredicts thevelocityintherange0.1_<y/

δ

_<0.4andunderpredictsthe ve-locityintheouter layer(y+>50), especiallyforthePANSmodel. FOU clearlyyields a laminar, parabolic, velocity profile. Boththe meanvelocityandReynoldsstressesindicatethatFOUandFOU-CD yieldnofluctuations,i.e.alaminarflow;CDandQUICKyielda tur-bulentsolution.GenerallythemagnitudeoftheReynoldsstresses

arelargerforQUICKthanforCD.Finallythespectraarecompared. DuetothelaminarflowpredictedusingFOUandFOU-CDthe en-ergycontainedinthespectrumismuchlowerforallmodels.Both CDandQUICKshowthecorrectshape,butforallmodelsCD con-tainsmoreenergyacrosstheentirefrequencyrange.Thereappears to be little difference between ILESand PANS. The spectrum for LES KSGSshowsa cut-off at alower frequency, sincepart ofthe turbulenceismodelledbythesub-filtermodel.

ToconcludeitisclearthatforallSRS,firstorderschemes add toomuchdiscretisationerrorandleadtoamismatchinflow pro-file.Fullysecond-orderschemesyieldthebestresults,whileusing QUICK(acommonlyused schemein industrialapplications), rea-sonableresultsare obtained.WithaQUICKschemelessenergyis resolved than by the CDscheme,which is inlinewith literature

Fig. 7. Normalised mean velocity and Reynolds stress profiles and turbulence kinetic energy spectra using LES and DNS [28] . From left to right, and top to bottom u + , Re

uu ,

Reuv , Re vv , Re ww and E u,y+_≈20(f) .

[35,67].Intheremainderofthisworktheconvectiontermsinthe momentumequationarediscretisedusingasecond-orderaccurate CDscheme. Forindustrialcases theapplicationofCD istypically not possible,duetohighlocalcell Pécletnumbers[68],although reasonableresultscanstillbeobtainedusingQUICK.

5. Comparisonofturbulenceapproacheswithstreamwise periodicboundaryconditions

5.1. LES

Fig.7 showsthemean velocityprofiles, Reynoldsstresses and turbulence kinetic energy spectra of the streamwise velocity at y+≈ 20. The mean velocity profiles for all models capture the trend of the DNS data. In terms of Reynolds stresses, it is clear that allmodels capturethe trend,butdeviate intermsof magni-tude.ThehighestReynoldsstresstermsaregenerallyobtainedfor the ILESandKSGS,followed by theSmagorinsky, Lilly andWALE model. The Lilly model performs adequately forReuu, but

under-predicts Reuv,Revv andReww.Forall modelsthemagnitudeofthe

Reynolds stress terms is generally underpredicted, with Reww an

exception. For Revv, the peak is shifted towards the right for all

modelsindicating thatthestrongestturbulencefluctuationsoccur

furtherfromthewall.Theturbulencekineticenergyspectraforthe LESmodelsarecomparable,althoughasshowninSection4.3 this isgriddependent.Fortheresolutionemployedhere,thereislittle differencebetweenthemodels. OnlyILESshowsaslightlyhigher cut-off frequencydue to the absence ofa sub-filter model. Gen-erallyitappearsthat ILESyieldsthe bestresults,whichhasbeen observedbeforeforachannelflowonafinegrid.Thisisrelatedto excessive diffusion andnon-monotonic grid convergence for LES, i.e.onfinegrids ILEScangive betterresultsthanLESwitha sub-filtermodel[69,70].

Turbulentstructures in theflow are visualisedinFig. 8,using isocontourplotsofQ,basedon theinstantaneous flow field.Q is defined as Q=1/2

(

|



|

−

|

S

|

)

, with







defined as the anti-symmetric part of

∇

u, representing local flow rotation [71]. The isocontoursandthe sidesofthedomainare colouredbythe nor-malisedstreamwise velocity(u∗=u/Ub).Based ona visual

obser-vation, it appears that the Lilly and WALE model predict larger structures than the other models. Forthe Lilly model, thisis re-latedtotheapplicationofthe‘test’filter,whichislargerthanthe gridsize; fortheWALE model,thisis aresultoftheinclusion of walldistanceinthedeterminationofthelengthscaleLs.

Toconclude,thecomparisonshowsthatthebestmatchforthe Reynoldsstresses is obtainedwithILES, although thisis strongly

Fig. 9. Normalised mean velocity and Reynolds stress profiles and turbulence kinetic energy spectra using Hybrid models and DNS [28] . From left to right, and top to bottom

u+ , Re uu , Re uv , Re vv , Re ww and E u,y+_≈20(f) .

griddependent. Smagorinsky, KSGS andWALE perform compara-bly,andtheLillymodelshowsthepoorestperformanceduetothe largertest filter.TheWALEmodelshouldperformbetterforflows involvingcomplexgeometries,meaningitsadvantagesdonotshow uphere.

5.2.Hybridmodels

Fig.9 showstheflowfield statisticsfortheHybridmodels.For themean velocity, XLES showsan underprediction inthe buffer-layer(5≤ y+_{≤ 30).DDES showsanoverprediction inthelog-law}

regionandan underpredictionnearthecentreline.IDDESmatches thereferencedata well. In termsof Reynoldsstresses, XLES cap-tures the trend but underpredicts the magnitude forall compo-nents.The peak in theRevv and Reww distribution issignificantly

shiftedaway fromthewall, anindicationof thehybridnature of themodel.TheDDESmodelhasamoreinterestingbehaviour.The shielding function of the model is formulated such that closeto thewall, insidethe boundarylayer, RANS should be used.In the farfieldLESshould beemployed, withthe RANSmodelactingas sub-filtermodel[32].Sinceaturbulentchannelflowconsistssolely ofaboundarylayer,withnofarfieldregion,onemightexpectthe solution to be fully RANS. This explains the good match for the averagedvelocity.HowevertheReynoldsstress componentsshow thatturbulenceisfullymodelledonlyintheregiony+<50,while closertothecentreturbulenceisresolved.Theemployedgrid res-olutionleads to lSRS

t <ltRANS, thereby forcing the switch to occur

insidetheboundarylayer.ThedivisionbetweenRANSandLES re-gions,isvisibleinFig.10.Notethateventhoughthisisthesame behaviour which lead to the developmentof DDES asa replace-mentof DES, thiscan still occur forDDES under certain circum-stances(i.e.combinationsoftestcaseandgrid density).An effect ofusingLESonlycloseto thecentreline isan underpredictionof theenergycontainedinthespectrumattheinvestigatedlocation (y+≈ 20). In IDDES and XLES, the RANS region is much smaller (seeFig.10), whichisreflectedinthe magnitudeofthe Reynolds stresses.Due to the use ofRANS closeto the wall an underpre-dictionoccursinthisregion.ForIDDESandXLES modelsthe tur-bulencekineticenergyspectrumalsomatcheswellwiththeother

Fig. 10. Instantaneous, spatially averaged, LES regions ( l t /l tRANS ) for the DDES and

IDDES model (right). 1 indicates RANS, 0 LES.

LES results.Note that the underpredictionfor DDES isrelated to thelocationoftheprobe(y+≈ 20),whichisinsidetheRANSlayer. ThefluctuationsatthispointareLESfluctuationswhichinfluence theRANSlayer.Theturbulencekineticenergyspectrummatchthe DNSdatabetterforspectraclosertothecentreline(notshownin thispaper).

The turbulent structures in the flow are visualised in Fig. 11. The effectofusingRANS nearthewallsinthe DDESformulation isobvious,onlylargerstructuresinthecentreofthechannelexist. ThestructuresintheIDDESandXLESmodelaresimilartotheLES Smagorinsky,KSGSandILESresults.

Theresultsindicatethatforcaseswheretheinstantaneousnear wallflowfieldisofimportance(forinstancesheetcavitation), Hy-bridmodelsarelesssuitablethanLESorPANS.DDESisnotableto properly resolve the boundary layer. In contrast, IDDES performs better for the mean velocity and Reynolds stresses butwith the exception ofRevv. XLES underpredicts the all components of the

Fig. 11. Hybrid models instantaneous turbulent flow fields ( Q = 0 . 3 ), coloured by u ∗_{= u/U b . From left to right DDES, IDDES and XLES.}

Fig. 12. Normalised mean velocity and Reynolds stress profiles and turbulence kinetic energy spectra using PANS and DNS [28] . From left to right, and top to bottom u + ,

Reuu , Re uv , Re vv , Re ww and E u,y+_≈20(f) .

5.3. PANS

ThePANSmodelisappliedwithfk=0.75,0.50,0.25,0.20,0.15,

0.10and0.05.Fig.12 clearlyshowstheeffectofreducingthefk

pa-rameteronthemeanvelocityprofile: fk=1.00yields aRANS

re-sult asexpected, while when moving from fk=0.75 to 0.20the

modelled turbulence kinetic energy is reduced and more turbu-lence kinetic energyshould be resolved. However, by reducing fk

theflowprofilebecomesmorelaminar.Ratherthanresolving tur-bulence the velocity perturbations are damped and the flow re-mainslaminar.Thehighestfkvaluesstillhaveareasonablematch

withthe DNS, sincetheseyield mostlyturbulentRANS solutions. Below f_k=0_.20thefluctuationsarenotdampedandafully turbu-lentflowdevelops.The resultsfor f_k=0.15,0.10and0.05match theDNSdatavelocityalmostperfectly.Thesamebehaviouroccurs ifPANSisappliedtoafullydevelopedturbulentflow(forinstance obtained froma LES computation); for higherfk values the

fluc-tuations are damped after between five and seven flow-through times.The normalisedReynolds stressprofiles andturbulence ki-neticenergyspectrumyieldadditionalinsightintothisbehaviour. It is clearthat forcomputations withfk inthe range 0.75− 0.20

alaminarsolutionisobtained;theReynoldsstresstensor compo-nentsequal zero. TheReynolds stress profiles and turbulence ki-netic energyspectra fk=0.15, 0.10and0.05are comparableand

obtainthepropertrendsandorderofmagnitude.Interestingly,the peaksfor f_k=0_.10 arehigherthanfor f_k=0_.15 and0.05. Thisis an indicationoftheneedforfinergrids anditerativeconvergence criteriaforlowerfkvalues.Forboth fk=0.10and0.05thematch

is better than for explicit LES, which is related to the increased turbulenceresolution.

TheturbulentstructuresintheflowarevisualisedinFig.13 for the lower fk values. There is little difference between the

differ-ent simulations.It appears that theresults forcases with low fk

are identicalto ILES. Thisis not true when looking atthe eddy-viscosity,however.Themaximumeddy-viscosityratio,

ν

t/

ν

,inthe

field for fk=0.15 has a magnitude of O

(

102

)

, whereas for fk=

0_.05thisis_O

(

10

)

_,andforILESitis zerobydefinition.For com-parison, for a turbulent RANS solution

ν

t/

ν

=O

(

105

)

. It is clear

that the magnitude of eddy-viscosity has little effect on the re-sults,providedthat f_k isbelowthethresholdtoallowa turbulent flow.

Theobservedstrongdependenceonfk,resultingindistinct

lam-inar andturbulent flow regimes,is related to the physics of the problem. In a wall-bounded turbulent flow, such as a turbulent channel or flat plate, the small scales near the wall move away fromthewall,andmergeintoincreasinglylargerscalesawayfrom thewall.Thelarge scalesthenbreakup intosmallscales andare dissipated.Thisprocessisknownasenergybackscatter,orthe in-verse energycascade [33]. A turbulent flow only develops if the filterlength issmallerthan the lengthscales of thesmallscales, otherwisethemechanismresponsibleforcreatingafullyturbulent flowisfilteredout.Thiscanoccuronacoarsegridforallmethods, or,inthecaseofPANS,whenusingalargef_kvalue.Thishypothesis isconfirmedbythe guidelinethat forSRStheeffective computa-tionalReynoldsnumber,

Re e= U

δ

ν

+

ν

modelled = U

δ

ν

+ f 2 k

ν

t , (39)

must exceed the critical transition Reynolds number needed for theonsetofinstability,Rec [72].Foraturbulentchannelflow,this

Fig. 13. PANS instantaneous turbulent flow fields ( Q = 0 . 3 ), coloured by u ∗_{= u/U b . From left to right f k = 0 . 15 , 0.10 and 0.05.}

guideline issatisfied for thecases withfk<0.20(Rec≈ 2300,

ob-tainedfrom experiments [73]). To enable turbulent results forfk

valuesintherange0_.20_{− 1.}00asyntheticturbulencegeneratoris needed,whichnot onlyfeedsthesmallscales, butalsofeeds en-ergydirectlyintothelargerscalestoallowaturbulencecascadeto develop.Notethat thesyntheticturbulenceis insertedacrossthe heightofthechannel,incontrasttothephysicswhereturbulence developsparalleltothewalls.

6. Comparisonofturbulenceapproacheswithsynthetic turbulencegenerator

The secondsetofresultsisobtainedbyapplyingthesynthetic turbulencegeneratordescribedinSection2.5.The turbulent fluc-tuationsareaddedtothemeanturbulentvelocityprofile,whichis prescribed at the inlet. The domain is initialised withthe mean turbulent flow profile, while the turbulence generator is located justbehindthe inlet.The body force terms are added overLtg=

0_.1

δ

which corresponds totwo cellsin the streamwisedirection. Inthiscaseperiodic boundaryconditionsare only appliedinthe spanwisedirection. Based on the TST (see Section 4.2), temporal averaging ofthe results is performedover 7 flow-through times, after removing two flow-through times. The resulting statistical uncertaintyforthemeanstreamwisevelocityisbelow3%,andfor theReynoldsstresscomponentsbelow12%.Themeanvelocityand Reynoldsstress profiles prescribedare takenfromthe DNS refer-encedata.Theprescribedintegrallengthscalesareanisotropicbut ahomogeneousapproximationismadebasedonthelengthscales givenby Kim et al. [39]; in the streamwise directionthe length scaleistakenas0.9

δ

,inthespanwise andwall perpendicular di-rectionas0.045

δ

.Itisnotedthat,duetothisassumption,the pre-scribedlengthscalesatthewallaretoolarge.

The applicationof a turbulence generator affects the iterative convergence;theL2 andL∞ normsareshowninTable2.Theratio

L_∞/L2 isO

(

101

)

forsyntheticcasescomparedto O

(

102

)

for

recy-clingresults,indicating that the residuals stagnate ina consider-ablepartofthedomain.Fig.15 showsthatthisoccursinthecells closeto the location at which the body force termsare applied, duetothelocal,explicit,additionoftheseterms,whichalsovary pertime step.Thecloseproximity oftheturbulencegeneratorto theinflow (where velocity is prescribed)also contributes to this situation.

Itisobservedthatastreamwisedevelopmentlengthisneeded toallowtheaddedperturbations todevelopintoaturbulentflow profile,andobtainareasonablematchwiththeinputvalues.With thecurrent implementationof the turbulencegenerator, after 6

δ

the results for u and Reuu are self similar, as shown in Fig. 14.

Table 2

Average residual norms for recycling and synthetic cases. Equation Recycling Synthetic

L∞ L2 L∞ L2

Momentum: 10 −5 ₁₀ −6 ₁₀ −2 ₁₀ −4

Pressure: 10 −6 ₁₀ −7 ₁₀ −4 ₁₀ −6

Turbulence: 10 −6 ₁₀ −7 ₁₀ −7 ₁₀ −9

Atthislocation,Reuvisstillunderpredicted,Reuvkeepsdeveloping

untilapproximately10

δ

.ThisiscomparabletoKimetal.[39] and XieandCastro[36],whobothgivearequireddevelopmentlength of10

δ

.Theresultsshownintheremainderofthissectionare ob-tainedatalocation6

δ

behindtheturbulencegenerator.

In order to compare the results obtained from recycling and syntheticcomputations,Fig.16 showstheresolvedturbulence in-tensity versus the channel height for all models. As shown in

Section 5, forthe recycling cases,the LES models show the best matchwiththereferencedata,togetherwiththePANSresultswith fk<0.2. The synthetic resultshowever,are similar forall models.

WiththeexceptionofPANS,Iisoverpredictedatthecentreofthe domain.BothDDESandPANSwithhigherf_kshowan underpredic-tionofIatthewall.ThePANSresultsshow anincrease inIwith decreasingfk,asmightbeexpected.Neverthelesstheresultsshow

that the applicationof a syntheticturbulence generator can pre-vent the occurrenceof laminar flow for cases where the critical transition Reynolds number exceeds the effective computational Reynoldsnumber(asisthecaseforf_{k≥ 0.20,}seeSection5.3). 6.1. LES

Fig. 17 shows the results obtained using different LES mod-els. Asexpected theshape ofthe meanvelocity profiles matches theDNSdatawell, howeveran underpredictionoccursacrossthe channel height. For Reuu, ILES and WALE have the correct peak

value, with the other models showing the correct shape but a lower peak.In contrastto therecyclingresults(Fig.7), whereall models underpredicted Reuv and Revv, theresults lie onor above

the DNS data for all models except the Lilly model. This model does accurately predict Reww however,which is overpredicted by

theothermodels.Fortheturbulencekineticenergyspectrum,the ILESandWALEmodelcontainthemostenergy,whilethespectrum fortheLillymodelshowslessenergythanfortherecyclingcases. The underpredictionby the Lillymodel is similar tothat already showninFig.7 andisagainrelatedtotheapplicationofthelarger ‘test’filter,whichfiltersoutsyntheticturbulence.

6.2. Hybridmodels

The resultsobtainedusingdifferentHybrid models areshown in Fig. 18. For the mean velocity, the same underprediction oc-curs as for the LES models. The results forthe DDES model are improved with respect to recycling cases (Fig. 9): the magni-tudeisstillsignificantlyunderpredictedbutthedistributionofthe ReynoldsstresscomponentsisclosertotheDNSdata,andthereis moreenergyinthespectrum.FortheIDDESmodel,theReynolds stressesarealsohighercomparedtotherecyclingresults.Thepeak ofReuu,ReuvandRevvisbettercaptured,butReww isoverpredicted.

The XLES results are similar to the IDDES results, with the ex-ception of an underprediction of Reuu near the wall which was

alsoobservedforrecyclingcases.Theenergyintheturbulence ki-netic energyspectrum islower for XLESthan forIDDES, exhibit-ing a magnitude similar to DDES. For all models Reuv is

under-predictedatthe centre,which wasnot thecasefortherecycling results.

Fig. 14. Normalised mean velocity and Reynolds stress profiles at different downstream locations using PANS with f k = 0 . 05 and a turbulence generator. From left to right,

u+ , Re uu and Re uv .

Fig. 15. Velocity (left) and pressure (right) residuals in the domain for a PANS computation with turbulence generator ( f k = 0 . 05 ). The turbulence generator is located at the

left of the images.

Fig. 16. Resolved turbulence intensity ( I = 1 / 3 u_{i u} _{i / u ) for recycling (top row) and turbulence generator (bottom row) for LES, Hybrid and PANS models. Results obtained at} x = 2 δand averaged in transverse direction.

6.3. PANS

Fig.19 showstheresultsobtainedusingPANSwithdifferentfk

values.Asseenfortheothersyntheticresults,themeanvelocityis underpredicted.ForallcomponentsoftheReynoldsstressesa sim-ilarbehaviourisobserved;withdecreasingfktheresultsconverge

towards the DNS data with decreasing fk. Reuu is overpredicted

at the centre with f_k=0_.05. The same occurs for Reww when

fk≤ 0.25, butthen across the full channel height. Theincrease in

fluctuations with decreasing fk can be related to the results

ob-tainedwithoutturbulencegenerator:theadditionofsynthetic tur-bulenceleadstoaturbulentprofile,yetcomputationswithhigher fkvaluesadddissipation,therebydampingtheresolvedturbulence.

Consequently, the resolved Reynolds stresses decrease in magni-tude.Thisismostvisibleintheturbulencekinetic energyspectra, wherethecut-off frequencyincreaseswithdecreasingfk.

Fig. 20 visualises the turbulent structures from PANS, based on Q. The most important difference with the results presented in Section 5.3 (Fig. 13) is that for fk>0.15 turbulent structures

Fig. 17. Normalised mean velocity and Reynolds stress profiles and turbulence kinetic energy spectra using LES and a turbulence generator. From left to right, and top to bottom u + , Re

uu , Re uv , Re vv , Re ww and E u,y+_≈20(f) . Results obtained at x = 2 δand averaged in transverse direction.

Fig. 18. Normalised mean velocity and Reynolds stress profiles and turbulence kinetic energy spectra using Hybrid models and a turbulence generator. From left to right, and top to bottom u + , Re uu , Re uv , Re vv , Re ww and E u,y+_≈20(f) . Results obtained at x = 2 δand averaged in transverse direction.

fluctuations are filtered by the model as they propagate down-stream of the turbulence generator, as can be most clearlyseen nearthe walls.The imagesaresimilar totherecyclingresultsfor thelowerfkvalues(fk≤ 0.15),howeverclosetotheturbulence

gen-eratormorestructurescanbeseeninthecentreofthechannel. 7.Discussion

Inthisworkturbulenceisgeneratedbothusingstreamwise pe-riodic boundary conditions and a pressure gradient, i.e. as pre-cursor (Section 5), and using a synthetic turbulence generator (Section6).ForindustrialCFDitisshownthat,independentofthe selectedturbulencesimulationapproach,thesyntheticmethodcan

produce a turbulent inflow at significantly lower computational cost.Onlytwoflow-throughtimesfromgeneratortoobjectof in-teresthavetobecomputedbeforestatisticscanbecollectedbased ontheTST,incontrasttotheprecursorcomputationsforwhich11 flow-throughtimeswererequired(Section4.2).Theabilitytomore easilytunethemethodtoobtainthedesiredReynoldsstressesand lengthscalesisalsoattractiveforindustrialapplicationsforwhich thesequantitiesmayalreadybeknown.Finally,itisnotedthatfor arange ofindustrialtest cases,suchasfoils,wingsorpropellers, only homogeneous inflow turbulence is required.For such cases, theadvantagesofasynthetic generatorclearlyshow;no flow re-cyclingisrequired,andthe gridonly needs tobe refined around theobjectofinterestandupstream.