A high-order discontinuous Galerkin solver for the incompressible RANS equations coupled to the k−ϵ turbulence model

A high-order discontinuous Galerkin solver for the incompressible RANS equations

coupled to the k− turbulence model

Tiberga, Marco; Hennink, Aldo; Kloosterman, Jan Leen; Lathouwers, Danny

DOI

10.1016/j.compfluid.2020.104710

Publication date

2020

Document Version

Final published version

Published in

Computers and Fluids

Citation (APA)

Tiberga, M., Hennink, A., Kloosterman, J. L., & Lathouwers, D. (2020). A high-order discontinuous Galerkin

solver for the incompressible RANS equations coupled to the k− turbulence model. Computers and Fluids,

212, [104710]. https://doi.org/10.1016/j.compfluid.2020.104710

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

A

high-order

discontinuous

Galerkin

solver

for

the

incompressible

RANS

equations

coupled

to

the

k

−



turbulence

model

Marco

Tiberga

∗

,

Aldo

Hennink

,

Jan

Leen

Kloosterman

,

Danny

Lathouwers

Delft University of Technology, Department of Radiation Science and Technology, Mekelweg 15, 2629 JB, Delft, The Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 15 July 2020 Accepted 21 August 2020 Available online 2 September 2020 Keywords:

Discontinuous Galerkin FEM Incompressible RANS k −turbulence model Pressure correction Wall function

Symmetric interior penalty method

a

b

s

t

r

a

c

t

AccuratemethodstosolvetheReynolds-AveragedNavier-Stokes(RANS)equationscoupledtoturbulence modelsarestillofgreatinterest,asthisisoftentheonlycomputationallyfeasibleapproachtosimulate complexturbulentflowsinlargeengineeringapplications.Inthiswork,wepresentanoveldiscontinuous Galerkin(DG)solverfortheRANSequationscoupledtothek−model(inlogarithmicform,toensure positivityoftheturbulencequantities).Weinvestigatethepossibilityofmodelingwallswithawall func-tionapproachincombinationwithDG.Thesolverfeaturesanalgebraicpressure correctionschemeto solvethecoupledRANSsystem, implicitbackwarddifferentiationformulaefortimediscretization,and adoptstheSymmetricInteriorPenaltymethodandtheLax-Friedrichsfluxtodiscretizediffusiveand con-vectivetermsrespectively.Wepayspecialattentiontothechoiceofpolynomialorderforanytransported scalarquantityand showithastobethesameasthe pressureordertoavoidnumericalinstability.A manufacturedsolutionisusedtoverifythatthesolutionconvergeswiththeexpectedorderofaccuracy inspaceandtime.Wethensimulateastationaryflowoverabackward-facingstepandaVonKármán vortexstreetinthewakeofasquarecylindertovalidateourapproach.

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

1. Introduction

Engineering applications often requirethe simulation of com-plex turbulent flows via accurate Computational Fluid Dynamics (CFD)methods.DirectNumericalSimulation(DNS)andLargeEddy Simulation (LES)constitute superiorapproachesin thisregard, as they are able to resolve the stochastic fluctuations (though only the large-scale onesin case of LES) of anyturbulent flow quan-tity[1].However,nowadaystheyarestillverycomputationally ex-pensiveandoftenunaffordable forlargeengineeringapplications. Forthisreason,theReynolds-AveragedNavier-Stokes(RANS) equa-tions coupledtoturbulenceclosuremodelsoftenremainthe pre-ferred approach, even if it only allows for the modeling of the time-averagedflowquantities[2].Accurateandefficientnumerical methods to solve the RANS equations are therefore still of great relevance.

In this perspective, discontinuous Galerkin Finite Element methods(DG-FEM)areparticularlyattractive,duetotheir flexibil-ity,highaccuracy,androbustness.ThecharacteristicfeatureofDG isthatunknownquantitiesareapproximatedwithpolynomial

ba-∗ _{Corresponding author.}

E-mail address: M.Tiberga@tudelft.nl (M. Tiberga).

sisfunctionsdiscontinuousacrossthemeshelements.Thisrequires theuseofnumericalfluxestodealwiththediscretizationofface integrals,asinFiniteVolumemethods.However,thankstothelack ofcontinuityconstraints,conservationlawsaresatisfiedlocallyon eachelement,andtheresultingalgorithmstenciliscompact, mak-ingthe methodsuitable forefficientparallelization. Asin contin-uousGalerkinFEM,theaccuracy ofthesolutioncanbe improved byincreasingtheorderofthepolynomialdiscretization.Moreover, DGmethods caneasily handlecomplex geometries,structured or unstructuredmeshes,andnon-conformallocalmeshand/or order refinement.

Thisclass ofFiniteElementmethods wasoriginally developed inthe early ‘70s to solve radiationtransport problems [3]. How-ever,ithasbecomeincreasinglypopularforCFDapplicationsonly inthe past three decades, afterthe development ofeffective DG schemesforhyperbolic andellipticproblems.We refertothe re-viewsbyCockburnandShu[4] andArnoldetal.[5] foracomplete overviewofthesemethods.

Nowadays,substantial experiencehas been gainedin the DG-FEMdiscretization ofthe incompressibleNavier-Stokes equations, andavarietyofapproachescanbefoundinliterature.Anearly re-search effortinthe field is thework of Cockburn etal. [6], who proposed a locally conservative Local DG (LDG) method for the incompressible Oseen equations. The authors later extended the https://doi.org/10.1016/j.compfluid.2020.104710

approachto the full Navier-Stokes system[7],employing a post-processingoperator to obtain an exactly divergence-free convec-tive velocity field. Van der Vegt and Sudirham [8] presented an originalDGschemeforthesolutionoftheOseenequations,where discontinuousbasis functionsarealso usedtodiscretize thetime domain. The method, called space-time DG, is particularly suit-ableforsimulationsinvolvingdeformingormovingmeshes.More recently, Rhebergen et al. [9] extended the approach to the full Navier-Stokessystem.In[10],theincompressibilityconstraintis re-laxedby introducinganartificialcompressibilityterminthe con-tinuityequation. Thus, thenumericalflux forthe inviscidpart of the Navier-Stokes equations can be computed by solving a Rie-mannproblemjust likefor compressibleflows. The coupled sys-temof equations is solved by means of a Newton method with animplicitEulertime scheme.Higher-orderRunge-Kuttaschemes were instead employed by Bassi et al. [11]. Splitting methods basedonpressureorvelocitycorrectionapproacheshavealsobeen proposed in the context of DG solvers (e.g., [12–15]). In partic-ular, the solver presented in this paper is based on the work of Shahbaziet al. [16], where a second-order accurate pressure-correctionschemeisappliedatalgebraiclevel.TheSymmetric In-teriorPenalty(SIP)methodisusedtodiscretizethediffusive oper-ator,andtheLax-Friedrichsfluxischosenforthehyperbolicterm. Similarsolvers, butbasedon a SIMPLEalgorithm,were proposed byKleinetal.[17] and[18].

Regardingturbulentflowsimulations,DGsolvershavebeen de-velopedmostlyforthecompressibleRANSequationscoupledto ei-thertheSpalart-Allmaras(e.g.,[19–21])orthek−

ω

models(e.g.,

[22–25]).Forthelattermodel,itiscustomarytosolveforthe log-arithmof

ω

andtoimposearealizabilityconditiononit,toensure the positivityof the solution and enhance the robustness ofthe solver.Theshear-stresstransportmodelwasinsteadconsideredin

[26],where particularattentionwasgiventothe developmentof a stabilized continuous FE discretization of the eikonal equation forthecomputationofthedistancetothenearestwall. Farfewer applicationsofDGschemestoincompressibleturbulentflowscan be found in literature. The work of Crivelliniet al. [27] was the first. They extended the artificial compressibility flux method of

[10] to deal with the set ofincompressible RANS equations cou-pledto the SA model.Particular attentionwaspaid to the treat-mentofnegativevaluesoftheturbulentviscosity,thusincreasing therobustness of themethod. The approach was later testedon complexthree-dimensional flows [28] andextended tothe k−

ω

model[29] and tohigh-order Runge-Kuttatime schemes [30,31]. Morerecently,Kranketal.[32] presentedaDG solverforthe in-compressibleRANSequationscoupledtotheSAmodelbasedona semi-explicitvelocity-correction splittingschemeaugmentedwith anexplicitstepfortheturbulenceequation.

In this work, we extenda solver forlaminar flows presented in [33] to handle turbulence through the set of incompressible RANS equations coupled to the k−



model.We employ a pres-surecorrectionschemetosolvetheRANSequations.Thus,contrary tomost ofprevious literature, we do not rely on a free artificial compressibilityparameter whoseoptimalvalue isproblem-based. Moreover,our time schemeis fullyimplicit (asit relieson back-ward differentiationformulae) also for the turbulence equations. It constitutes the first step towards the development of a cou-pledCFD-neutronicssolver forlarge multi-physicsproblems,such astransientscenariosinmoltensaltnuclearreactors.Inthese sys-tems,theflow Reynoldsnumberisof theorderof105 _orhigher, and resolving the steep gradients that characterize the velocity profileclosetowallboundariesrequiresmassivecomputational re-sources.

To thebest ofour knowledge,all previous literatureon RANS DGsolversdealswithwall-resolvedturbulentflowsimulations.An approachfor efficientwall modeling in the DG context was

pro-posed by Krank et al. [32], who solved the RANS equations on coarse grids by enriching the polynomial function space for the velocity inthefirstlayer ofboundaryelements.In thisway,they couldimposeno-slipconditionsatthewallsandresolvethesharp solutiongradients.However,thek−



modeliswell-knownto per-formpoorlyinthevicinityofawall.Forthisreason,inthiswork weinvestigatethepossibilityofusinga morestandard wall func-tionapproach[34] tobridgethegapbetweentheviscouslayerand theloglayer.

Theremainderofthepaperisorganizedasfollows.InSection2, we presentthegoverningequationsandthe boundaryconditions that close them.We then describe ourspatial andtemporal dis-cretizationschemeinSections3 and4,andweprovidedetailson the solver in Section 5. In Section 6, we focus on the choice of polynomial orderforscalarquantities,which differsslightlyfrom what previously proposed in literature forsimilar DG solvers. In

Section 7, we verify our numericalscheme witha manufactured solutionandweassessthesoundnessofourapproachsimulating thesteadyturbulent flow overabackward-facing step anda Von Kármánvortexstreetinthewakeofasquarecylinder,finally draw-ingsomeconclusionsinSection8.

2. Governingequations

The RANS model equationsfor incompressibleflows read (we omitdependenciesforclarity)[2]

∂

u

∂

t +

∇

·

(

uu

)

+

∇

p−

∇

·

τ

=f, (1a)

∇

· u =0. (1b)

Here, u is the velocity andf represents a known momentum source.Thetotalstresstensoris

τ

=

(

ν

+

ν

t

)



∇

u+

(

∇

u

)

T



_{, (2)}

where

ν

and

ν

t are the molecular and the turbulent (or eddy) kinematic viscositiesrespectively. Finally, p= p˜

ρ+23k is a pseudo kinematic pressure,where p˜is the fluid pressure,

ρ

the density, andktheturbulentkineticenergy.

The system is closedby adopting thek₋



turbulence model

[34],where



indicatesthedissipationrateofk.Toensure positiv-ityoftheturbulencequantities(whichmightnotbepreservedby thenumericalscheme), themodelequationsare castinto a loga-rithmicform [35], whichis completely equivalentto the original model.Therefore,wesolvefor

K=ln

(

k

)

, E=ln

(



)

, (3) andthemodelequationsread(omittingdependencies)

∂

K

∂

t +

∇

·

(

uK

)

−

∇

·





ν

+

ν

t

σ

k



∇

K



=



ν

+

ν

t

σ

k



∇

K·

∇

K+

ν

te−KPk

(

u

)

− Cμ eK

ν

t + q_K, (4a)

∂

E

∂

t +

∇

·

(

uE

)

−

∇

·





ν

+

ν

t

σ





∇

E



=



ν

+

ν

t

σ





∇

E·

∇

E+C1

ν

te−KPk

(

u

)

− C2eE−K+qE, (4b) where P_k

(

u

)

₌

∇

u :



∇

u₊

(

∇

u

)

T



_{models the shear production} of turbulent kinetic energy. We used the standard values of the model constantsproposed by Launderand Spalding[34] and re-portedinTable1.Theeddyviscositycanbethencomputedas

ν

t=Cμ

k2

Table 1

Values of the closure parameters of the k −turbu- lence model.

Parameter Cμ σk σ C1 C2

Value 0.09 1.0 1.3 1.44 1.92

Clearly, thelatter expressionensures

ν

t is always positive. More-over,potentiallytroublesome divisionsbysmallvaluesof



(orby

kinsometermsoftheturbulenceequations)areavoided. Finally,q_K andq_E are extrasourcetermsthat,togetherwithf, are used inthe framework of themethod ofmanufactured solu-tions to verifythe numerical scheme(see Section 7.1). Normally, theyaresettozero.

2.1. Boundaryandinitialconditions

The coupled system of theRANS (1a-1b) andturbulence (4a-4b) equations can be solved only after imposition of proper boundary and initial conditions. The latter consist in imposing suitable values for the initial velocity, pressure, and turbulence fields,which we describe in detailfor each test casereportedin

Section7.Fortheformer,fourtypesofboundaryareconsideredin thiswork,andwedescribetheminthefollowing.

2.1.1. Inflowboundary

On inflow boundaries, Dirichlet conditions are prescribed for thefluidvelocityandtheturbulencequantities:

u=uD_{, K=K}D_=ln

₍

_kD

₎

_{, E=E}D_=ln

₍

_

D

₎

_{. (6)}

2.1.2. Outflowboundary

On outflowboundaries we prescribe a zero-traction condition for the momentum equation andhomogeneous Neumann condi-tionsfortheturbulenceequations:

[−pI+

τ

]· n=0, n·

∇

K=0, n·

∇

E=0. (7) Here,nis theoutwardunit normalvector atthe boundary,andI

representstheidentitytensor.

2.1.3. Wallboundary

The standardk−



modelcannot be usedtoresolve turbulent flows up to wall boundaries, in thevicinity of which viscous ef-fectsaredominantandturbulenceisdamped.Forthisreason,we usethewell-establishedwallfunctionapproachtodescribethe so-lutionintheseareas[34,36].

The computational boundary is shifted by a distance

δ

w from the actual physical wall, so that the first degrees of freedom of the solutionare placedin theregion offullyturbulent flow[37]. Thegapisbridgedbyanalyticalfunctionsthataresupposedtobe universal forallwall-boundedflows. Thisleadstoagreat gainin computationalefficiency,inparticularforhigh-Reynoldsflows, be-causemeshelementsarenot accumulatedintheviscous layerto resolvethesteepgradientsofthesolution.

Inthiswork,weadoptthetwo-velocityscalewallfunction de-scribedby Ignatetal.[38].Letutbe thevelocitycomponent tan-gential to the boundary and y the normal distance to the wall. Theyaremadenon-dimensionalasfollows:

u+= ut

u_τ, y

+₌

δ

wuk

ν

, (8)

where u_τ=√

τ

w/

ρ

is thefriction velocity (

τ

w istheshear stress atthewall),whereasukisasecondvelocityscalecomputedfrom

k:

uk=C_μ0.25



kw=C_μ0.25eKw/2, (9)

wherekwisthevalueoftheturbulentkineticenergyatthe bound-ary.Thewallfunctionisananalyticalrelationbetweenu+andy+:

u+= 1

κ

ln

Ey+

for y+>11.06, (10) where

κ

=0.41istheVonKármánconstant.Onlysmoothwallsare consideredinthiswork,sotheroughnessparameterisE=9.0.

Eqs.(9) and(10) are usedto prescribeahomogeneous Robin-like boundarycondition forthe momentum equation in the tan-gentialdirection.Infact,alinearrelationexistsbetweentheshear stressandthetangentialvelocitycomponent:

τ

· n− n

(

n·

τ

· n

)

=−ukuτt=−

uk

u+[u−

(

u· n

)

n], (11) where t is a unit vector that is tangential to the boundary and u₋

(

u_{· n}

)

n₌utt.Theboundarycondition(11) issupplementedby ano-penetrationconditioninthenormaldirection:

n· u=0. (12)

Regarding the turbulent quantities, homogeneous Neumann and Dirichletboundary conditions are prescribed for k and



respec-tively: n·

∇

K=0; E=ln



u3 k

κδ

w



. (13)

The value of

δ

w must be sufficiently large for the wall function (10) to be valid, that is,the first degreesoffreedom ofthe FEM solutionmust be in thelogarithmiclayer. This parameter canbe eitherestimatedasthe distancefromthe wall ofthe centroidof thefirstboundaryelementorfixedtoapredeterminedvalue[39]. Thisworkfeaturesstandardandpredictableflowconfigurations,so weadoptedthesecondapproach.

2.1.4. Symmetryboundary

On symmetry boundaries, we prescribe a homogeneous Neu-mannconditionforbothturbulencequantities:

n·

∇

K=0; n·

∇

E=0. (14)

Forthe momentum equation, we prescribevanishing normal ve-locityandtangentialstress:

n· u=0,

τ

· n− n

(

n·

τ

· n

)

=0. (15)

3. Spatialdiscretization

Inthissection,wedescribeindetailthespatialdiscretizationof thegoverningequationswiththediscontinuousGalerkinFinite El-ementmethod.Somenotationanddefinitionsmustbeintroduced first.

Let



be the computational domain and



its boundary. We indicate with



D_,



N_,



W_,and



S _{the non-overlapping} por-tions of boundary where Dirichlet (i.e., inflow), Neumann (i.e., outflow), wall, and symmetry conditions are imposed, such that



=



D_∪



N_∪



W_∪



S_.

Thedomainis meshedinto asetofnon-overlappingelements

Th.ThesetoftheirinteriorfacesisindicatedwithFi,whereasFD,

FN_,FW_,andFS _{indicatethesetsofDirichlet,Neumann,wall,and} symmetryboundaryfaces,whichdiscretizetherespectiveportions ofboundary. Given an element T∈Th, the set ofits faces is de-notedwithFT_{.TheneighborsoffaceFare groupedintoT}F

h.Any interiorfaceF∈Fi_{isequippedwithaunitnormalvectorn}F_with anarbitrarybutfixeddirection,whileforboundaryfaces nF coin-cideswithn,whichistheunitnormalvectoroutwardto



.

AnyscalarunknownisapproximatedonTh witha setofbasis functionsthatarecontinuouswithineachelement,but discontinu-ousattheelementinterfaces.Weuseahierarchicalsetof orthog-onalmodalbasisfunctions:thesolutionspacewithin anelement

isthespanofallpolynomialsuptoanorderP,withP≥ 1.Fora genericunknown

φ

,we denoteitsFEMapproximation by

φ

h and itspolynomialorderbyPφ,sothatthesolutionspacecanbe writ-tenas

Vh,φ:=

v

∈L2

()

|

v

|

T∈PPφ,

∀

T∈Th



, (16)

where_PPφ is the set of polynomials up to order _Pφ. The veloc-ityFEMapproximationliesinthecorrespondingvectorspaceVd

h,u, where d is the problem dimensionality. Despite not being a re-quirementofthe numerical method,Pφ isthe same on each el-ementandthe meshesareconforming forallthe test cases ana-lyzedinthiswork.

Thetraceof

φ

honafaceF∈Fiisnotunique.Forthisreason,it isnecessarytointroducetheaverage

{

φ

h

}

:= 12

φ

+

h +

φ

h−

andthe jump

φ

h:=

φ

_h+−

φ

_h−operators.Here,

φ

+_h and

φ

−_h arethefunction tracesdefinedas

φ

±

h

(

r∈F

)

=lim_ς↓0

φ

h

(

r∓

ς

nF

)

, (17) forapointronaninteriorfaceF∈Fi_{.Attheboundary,thejump} andaverageoperatorscoincidewiththeuniquetraceof

φ

h.

3.1.RANSDiscretization

The semi-discrete weak formulation of System (1a-1b) sub-jectedtothe boundaryconditions describedinSection 2.1 is ob-tained by substituting (u, p) with their DG-FEM approximations (uh,ph), by multiplying Eq. (1a) with the test function vh∈Vhd,u, bymultiplying Eq.(1b) withq_h∈Vh,p, andby subsequently inte-gratingovertheentiredomain:

Finduh∈Vhd,uand ph∈Vh,psuchthat

atime

₍

_u

h,vh

)

+aconv

(

uh,uh,vh

)

+adiff

(

uh,vh

)

+adiv

(

vh,ph

)

=l

(

uh,vh

)

∀

vh∈Vhd,u, (18a) adiv

₍

_u h,qh

)

=ldiv

(

qh

)

∀

qh∈Vh,p, (18b) where atime

₍

_u h,vh

)

=  T∈Th  T vh·

∂

uh

∂

t , (19)

and,lettingfh betheGalerkinprojectionoffontoVhd,u,

l

(

uh,vh

)

=lconv

(

uh,vh

)

+ldiff

(

vh

)

+  T∈Th  T vh· fh. (20)

Allothertermsaredescribedmoreindetailinthefollowing.

3.1.1. Convectiveterm

The convectivetermofthe momentumequation isdiscretized as aconv

₍

_β

_,w,v

₎

₌₋ T∈Th  T w·

(

β

·

∇

)

v+ F∈Fi  F vHF

_β

_,w

+ F∈FN  F

nF_·

β

_{w· v} + F∈FD  F max



0,nF_·

_β

D



w· v, (21) lconv

₍

β

_,v

₎

₌₋ F∈FD  F min



0,nF_·

β

D



uD_{· v, (22)}

where

β

istheconvectivefield(i.e.,thevelocity),andHF _isa nu-mericalflux functiondefinedontheinternal faceF.In thiswork, theLax-Friedrichsfluxisused[40]:

HF

₍

_β

_,w

₎

₌

α

F

(

β

)

2 w+n

F_·

_β

_w



_{, (23)}

where

α

F _{isevaluatedpoint-wisealongfaceFaccordingto}

α

F

₍

_β

₎

₌

_

_max



_nF_·

_β

+



_,



_nF_·

_β

−



_{, (24)} with



=2forthemomentumequation.

Notethatusually(nF_·

_β

D_{)<0onanyF∈F}D_{,andthereforethe} corresponding term in Eq.(21) drops out. However, we included thistermbecauseitisrelevantfortheTaylor-vortex-like manufac-turedsolutioninSection7.1,whereweimposeDirichletconditions ontheentireboundary,asitisusuallydoneforthisclassof prob-lems(e.g.,[16,41]).

3.1.2. Diffusivetermdiscretization

Todiscretizethediffusiveterm,we useageneralizationofthe SymmetricInterior Penalty(SIP)method [42, p.122]. Among the several methods available for elliptic problems [5], we optedfor theSIPbecauseitisconsistentandadjointconsistent,which guar-anteesoptimalconvergenceratesforanyorderofthe polynomial discretization.Moreover,itischaracterizedbycompactstencilsize, asthedegreesoffreedomoneachelementarecoupledonlywith thoseonthe firstneighbors,withpositive impactonmemory re-quirements,computationalcost,andparallelizationefficiency.

Thefollowingbilinearoperatorandright-handsidetermcanbe defined: adiff

₍

_w,v

₎

₌  T∈Th  T Ldiff

₍

_w

₎

_:

_∇

_v+  F∈Fi_∪FD  F

η

F_w·v −  F∈Fi_∪FD  F nF_·

_v·

_Ldiff

₍

_w

₎



₊

_Ldiff

₍

_v

₎



_·w

+adiff FW_∪FS

(

w,v

)

(25) and ldiff

₍

_v

₎

₌  F∈FD  F

η

F_uD_{· v− u}D_{· L}diff

₍

_v

₎

_{· n}F

_{, (26)} whereLdiff

₍

_w

₎

_{is a linearoperatorthat coincides withthe stress} tensor

τ

(Eq.(2))forthemomentumequation,thatis,Ldiff

₍

_u

₎

₌

_τ

_. Note that the effectiveviscosity (i.e.,

ν

₊

ν

t) is variable in space andso

∇

·

τ

=

(

ν

+

ν

t

)

∇

2u,which iswhythe regular SIPcannot beusedhere.

The penalty parameter,

η

F_{, must be sufficiently large to} en-surestabilitywithoutimpactingtoonegativelythecondition num-berofthesystem.We followtheoptimalvalueprescriptions pro-videdbyHillewaert[43] andDrossonandHillewaert[44] incase ofanisotropic meshesandhighlyvariableviscosity(as when tur-bulentflowsareresolved):

η

F ₌ D max T∈TF h



card

(

FT

₎

CP,T LT



, (27)

wherecard

(

_FT

₎

_{indicates thenumberoffacesofelement T;} C_P,T isafactordependingonthepolynomial orderofthesolutionand themeshelementtype,

C_P,T =



(

P+1

)

2_{, forquadrilaterals andhexahedra} (P+1)(P+d)

d , forsimplices.

(28) LT isalengthscaledefinedas

LT=f



T



leb



F



leb

, (29)

where



·



leb represents the Lebesgue measure of a geometrical entity(i.e.,length,area,orvolumedependingonthe dimensional-ity),andf=2forF∈Fi_{andf=1forboundaryfaces;finally,Dis} ascale dependingonthe diffusionparameterD(i.e.,theeffective viscosityforthe momentum equation),whichwe evaluate point-wisealongthefaceas

The contributionof wallandsymmetry faces muststill be speci-fied to complete the definitionof the SIPbilinear form (25). We extend what wasproposed by Cliffe etal. [40] to arbitrarily ori-entedsymmetry boundariesandtowalls, whereanon-zeroshear stress is specified inthe tangentialdirection, leadingto a Robin-likeboundaryconditionterm(see,forexample,[42,p.127]):

adiff FW_∪FS

(

w,v

)

=  F∈FW  F uk u+

w−

(

w· nF

₎

_nF

_·

_v−

₍

_{v· n}F

₎

_nF

+  F∈FW_∪FS  F

η

F_{w· n}F_{v· n}F_ −  F∈FW_∪FS  F nF_·

₍

_{v· n}F_

_Ldiff

₍

_w

₎

_{· n}F



+

Ldiff

₍

_v

₎

_{· n}F



_{w· n}F_

₎

_{, (31)} where u_k/u+ is the proportionality constant between tangential stressandvelocitycomponent(Eq.(11)).

3.1.3. Continuityterms

Thediscretizedcontinuityequationisconstitutedbythe follow-ingdiscretedivergenceoperatorandright-handsideterm[7,16]:

adiv

₍

_v,q

₎

₌₋ T∈Th  T q

∇

· v+  F∈Fi_∪FD_∪FW_∪FS  F

{

q

}

v· nF _{, (32)} ldiv

₍

_q

₎

₌  F∈FD  F q

uD_{· n}F

_{. (33)}

3.2. Discretizationofturbulenceequations

The spatial discretization of the turbulence Eq. (4a-4b), sub-jected to the boundary conditions described in Section 2.1, pro-ceedsalongthesamelineofwhatwasdescribedforthe momen-tumequation.Thesemi-discreteweakformulationsare

FindKh∈Vh,KandEh∈Vh,Esuchthat

atime

₍

_K h,

v

h

)

+aconv

(

uh,Kh,

v

h

)

+adiff

(

Kh,

v

h

)

=l

(

Kh,

v

h

)

∀

v

h∈Vh,K, (34a) atime

₍

_E h,

v

h

)

+aconv

(

uh,Eh,

v

h

)

+adiff

(

Eh,

v

h

)

=l

(

Eh,

v

h

)

∀

v

h∈Vh,E. (34b) Thetemporalandright-handsidetermsarestraightforward ex-tensiontothecaseofscalarunknownsofthosedescribedforthe momentum equation (Eqs. (19) and(20)). All termsonthe right-handside ofEqs. (4a) and(4b) are treatedexplicitlyintime and thereforeaddedtothesourceintegralsinl

(

Kh,

v

h

)

andl

(

Eh,

v

h

)

.

The convective terms are discretized as for the momentum equation, though the

α

F _{parameter in the Lax-Friedrichs flux} (Eq.(24))isevaluatedwith



=1.Diffusivetermsaretreatedwith theSIPmethodasinEqs.(25) and(26),wherethelinearoperator

Ldiff

₍

_v

₎

_reducestoLdiff

₍

_v

₎

_=D

_∇v

_,with

D=



ν

+νt

σk

, fortheKequation

ν

+νt

σ

, fortheE equation. (35)

These are also the diffusion coefficients used to evaluate the penalty parameter according to Eq. (27). Finally, wall faces con-tributetotheSIPbilinearformoftheE equationasotherDirichlet boundaries, whereasfortheK equationtheir contributionisnull, becausehomogeneous Neumannconditionsareimposedonthese boundaries. Likewise, for both equations symmetry faces have a nullcontributiontotheSIPbilinearform.

Table 2

Coefficients of the implemented BDF schemes. Order, M γ0 γ1 γ2

1 1 −1

2 3/2 −2 1/2

Table 3

Coefficients of the implemented extrapolation schemes.

Order, M ζ0 ζ1

1 1

2 2 −1

4. Temporaldiscretization

Time derivatives in weak forms (18a-18b) and (34a-34b) are discretized implicitly using backward differentiation formulae (BDF)oforder_M [14,16,18].Forthe genericunknown quantity

φ

andforaconstanttimestepsize



t,itis

∂φ

∂

t ≈

γ

0



t

φ

n+1₊M j=1

γ

j



t

φ

n+1− j_{, (36)}

wheren+1indicates thenewtimestepandtheBDFcoefficients

γ

arereportedinTable2.Weuseuptosecond-orderschemes.

The discretized equations are solved in a segregated way, so that the discrete solution vectors at the new time step,

u,p,K,E

n+1_{,arefoundwiththefollowingalgorithm:}

1.Obtain predictorsforallunknownsatthenewtime stepwith anMth_{-orderextrapolationfromtheprevioustimesteps:}

φ

n+1_≈

_φ

∗₌M−1 j=0

ζ

j

φ

n− j, (37)

wheretheweights(

ζ

j)arereportedinTable3.Apredictorfor the eddy viscosity,

ν

_t∗_, can be calculated from _K∗ and _E∗ ac-cordingtoEq.(5);

2.SolvetheRANSsystem(1)inasegregatedwaywitha second-orderaccuratealgebraicsplittingschemedescribedindetailin

Section 4.1. Use the predictors forthe velocity, the turbulent kineticenergy1_{,andtheeddyviscositytolinearizeallterms;}

3.Solvethe_KEq.(4a),usingun+1

h asconvectivefieldand

ν

t∗; 4.SolvetheEEq.(4b),usingu_hn+1asconvectivefieldand

ν

_t∗; 5.If necessary,update all predictorswiththenewsolutions and

iteratesteps2to4.

When theBDF2 scheme isused, the time extrapolation ofall quantities performed at the beginning of the time step guaran-teessecond-ordertimeconvergenceevenwhennon-linearitiesare unresolved (as proven in Section 7.1.1). During particularly vio-lentportions of transient(as it might happen atthe start-up of apseudo-transienttowardssteady-stateconditions),iterations be-tween RANSandturbulence equations(typically2or 3)are usu-allyperformedtobetterresolve non-linearities,thus avoidingtoo severerestrictionsonthetimestepsize.

4.1. Algebraicsplittingscheme

Thecoupledmomentumandcontinuityequationsaresolvedin asegregatedwaywithasecond-orderaccuratepressurecorrection method[45].Following[16],we performthesplittingatalgebraic level,whichdoes notrequirethe impositionofartificial pressure

1 Necessary to compute the u

boundaryconditions [46]. The fullydiscretized weak formofthe RANSsystem(18)canbecastintothefollowinglinearform:



_γ0 tM+N D T D 0



u p



n+1 =



bu bp



n+1 , (38)

where M is the mass matrix, N contains the implicit parts de-rivingfromthediscretization oftheconvective (21) anddiffusive

(25) terms, and D is the discrete divergence operator (32). The right-handsidevectorbu collectsall theknowntermsinthe mo-mentum equation (i.e., boundary conditions, external forces, ex-plicitterms fromthe discretizationof thetime derivative),while bp represents the fully-discrete right-hand side ofthe continuity Eq.(33). As explainedpreviously, predictorsforthe velocity field andtheeddyviscosityatthenewtimestepare usedtolinearize theconvectiveanddiffusivetermsinN.

The coupledSystem (38) issolved ina segregatedwayby re-placing(andapproximating)theleft-hand-sidematrixwithits in-completeblockLUfactorization[46],



_γ0 tM+N DT D 0



≈



_γ 0 tM+N 0 D −t γ0DM −1_DT



_

I −t γ0M −1_DT 0 I



, (39) where_Iistheidentitymatrix,andbyintroducing

δ

pn+1₌pn+1₋

pn_,suchthat



_γ0 tM+N 0 D −t γ0DM −1_DT



_

I −t γ0M −1_DT 0 I



_u

δ

p



n+1 =



_b u b_p



n+1 +



_−DT_p 0



n . (40)

Thediscretevelocityandpressurefieldsatthenewtime stepcan thereforebecalculatedwiththefollowingpredictor-corrector algo-rithm:

1.Findanapproximatevelocityfieldu˜n+1bysolving

 γ

0



tM+N



˜

un+1=bn_u+1− DT_pn _{; (41)}

2. Solve a Poissonequation to get thepressure atthe newtime step:





t

γ

0DM −1_DT



δ

pn+1=Du˜n+1− bp, (42) pn+1₌

_δ

_pn+1_+pn_{; (43)} 3. Correctthevelocityfield, sothatitsatisfiesthediscrete

conti-nuityequation

un+1_=u˜n+1₋



t

γ

0M

−1_DT

δ

_pn+1_{. (44)}

5. Solutionoflinearsystems

The numericalschemedescribedinSections3 and4 hasbeen implementedinthein-houseparallelsolver

DGFlows

.

We employ

METIS

to partition the mesh [47] and the MPI-basedsoftwarelibrary

PETSc

[48] to assemble andsolveall lin-earsystemswithiterativeKrylovmethods.Theimplicitconvection treatment leads to non-symmetric matrices for the momentum andturbulenceequations.Forthisreason,thoselinearsystemsare solvedwiththeGMRESmethodcombinedwithablockJacobi pre-conditioner,withoneblockperprocessandanincompleteLU fac-torizationon eachblock.Ontheother hand,thepressure-Poisson systemissymmetricandpositivedefinite.Hence,wesolveitwith

the conjugate gradient method and an additive Schwarz precon-ditioner, withone block per process andan incompleteCholesky decompositiononeachblock.

SolvingthepressurePoissonequationisthemost computation-ally expensivestep of our algorithm. This ispartially due to the factthatthepressurematrixinEq.(42) correspondstoaLaplacian operator discretized with the local DG method, as explained by Shahbazietal.[16],andsoitischaracterizedbylargerstencilsize than the SIPdiffusive operator. However, note that the matrix is thesameateachtimestep.Therefore,itispossibletoassembleit andcomputeitspreconditioneronlyonce,thuspartiallyalleviating thecomputational burden. Ontop ofthis, we initializeall Krylov solverswiththesolutionpredictorsdescribedinSection4 tospeed convergenceup.

6. Choiceofthesolutionpolynomialorder

Inthiswork, weusedamixed-orderdiscretizationforvelocity andpressure(i.e.,Pp=Pu− 1),whichsatisfies theinf-sup condi-tionandguarantees optimalerrorconvergenceratesandstability inthelow-



tlimit[33,49].

The polynomial orderof theturbulence quantities equals that ofthepressure,whichwebelieveisessentialforapressure-based DGsolver.Moreprecisely,thesolutionspaceofanarbitrary trans-portedscalarquantityshouldbea subsetofthesolutionspaceof thepressure.Wehavealreadymentionedthisin[33],butherewe corroborateitwithmorerigoroustheoreticalandnumerical argu-ments.Thischoiceofpolynomialorderisincontrastwithprevious literature on mixed-order DG schemes. For example, Klein et al.

[50] chosethesamesolutionspaceforthetemperatureasforthe componentsofthevelocityfield.Wesuspectthattheyfoundgood resultsbecausetheir testsweredone atalow Prandtlnumberof 0.7,whereastheproblemwiththesolutionspacesmanifestsitself whentheconvectivetermdominates,asshowninthefollowing.

Our choice of polynomial order for anyscalar solution stems fromthefactthatthecontinuityequationisweighedbythe pres-surebasisfunctions,sothatthenumericalvelocitysatisfiesthe in-compressibilityconstraintonlyinaweaksenseuptoorderPp(i.e., Du_h=0).Thismeansthattheconvectivediscretizationcanonlybe consistentuptoanorderPp.

Toshowthis,considertheconvectionofagenericpassivescalar quantity

φ

withthenumericalvelocityuh.Forsimplicity,assume that the domainis closed(i.e., n· u=0 on



), sothat substitut-ing

β

←uhandthecontinuoussolution(i.e.,w←

φ

)intoEqs.(21), (22),and(23) givestheconvectivediscretization

aconv

(

uh,

φ

,

v

)

=−  T∈Th  T

φ

uh·

∇v

+  F∈Fi  F

φ



v



{

uh

}

· nF, (45a) lconv

₍

_v

₎

_{=0, (45b)}

for a test function

v

∈Vh,φ. Here, we have used the fact that

φ

is continuous, so that 

φ

=0 and

{

φ}

=

φ

, and so

HF

₍

_u

h,

φ

)

=

φ{

uh

}

· nF.

Comparethistothediscretecontinuityequation(Eqs.(32) and

(33)), which, upon substituting v_←u_h, integrating by parts, and using the fact that quh=q

{

uh

}

+

{

q

}

uh on an interior face, canberewrittenas adiv

₍

_u h,q

)

=  T∈Th  T uh·

∇

q−  F∈Fi  F q

{

uh

}

· nF=0, (46)

fora test function q that lies inthe pressure solutionspace Vh,p [51].ItisclearthatEq.(45a) canbeconsistentonlyforatest func-tionv thatispartofthetest spaceofthecontinuityequation.In fact,considerthespecialcaseofaconstantsolution

φ

.The convec-tivetermshouldvanishforallv,butaconv

₍

_u

Fig. 1. Domain of the lid-driven cavity test case.

whichonlyvanishesifVh,φ isasubsetofVh,p,thatis,ifthe poly-nomialorderofthepressureisatleastashighasthatofthescalar. If Pφ>Pp, the constant solutioncannot be preservedin general duetonumericalerrorintroduced bythediscontinuityacrossthe elementsofthevelocityfield.

Wecorroboratethiswithasimplenumericalexamplebasedon a standard lid-driven cavity problem. The domain is depicted in

Fig.1.TheReynoldsnumberis40,basedonthelid-velocityUand the characteristiclength L,so theflow islaminar. Wefirst found the steady-state velocity (uh) and pressure fields by solving the

laminarNavier-Stokesequations.Then,wesolvedatransport equa-tion fora passive scalar

φ

with thenumerical velocity uh anda

diffusioncoefficientD:

∂φ

∂

t +

∇

·

(

uh

φ

)

=

∇

·

(

D

∇φ

)

, (47)

withhomogeneousNeumannboundaryconditionsfor

φ

.Weseta constantinitialcondition

φ

0=

φ

(

t=0

)

=1,andtookasingletime stepwiththeBDF1schemeandU



t/L=1,choosingarelative tol-eranceof10−12fortheGMRESandCGKrylovsolvers,tominimize anynumericalerrorinthesolution ofthelinearsystems.Clearly, the discrete solution

φ

h should remain unchangedafterthe time step.WethereforecalculatedtheL2_-normerrorof

_φ

hwithrespect to theexpectedexactsolution

φ

0.Tohighlightthe impactofthe velocitydivergenceandcontinuityerrors,werepeatedthetestfor variousspatialrefinements.Following[52],weused

ED=L  T∈T  T

|

∇

· u

|

 T∈T  T

||

u

||

and EC=  F∈Fi  F



u· nF



 F∈Fi  F

|

{

u

}

· nF

|

(48) toestimatethemassconservationerrorinastrongsense.

Table4 collectstheresultsfortwostructuredmeshesmadeof quadrilateralelements:(i)50by 50uniform(indicatedwithM1); and(ii)90by90progressivelyfinertowardstheboundaries(M2). On M2, we studied the effect of varying the velocity polynomial order.

Results supportour theoretical argument. When P_φ>Pp, the constantsolutionisnot preserved,duetothenumericalerror in-troducedbydiscontinuityofthevelocityfieldacrosstheelements. In fact, asexpected, the errordecreases withspatial refinement, thatis,withlowerdivergenceandcontinuityerrors.Theproblem ismorepronouncedinconvection-dominatedproblems,asincase ofturbulentflows.Infact,thebiggesterrorscorrespondtolow val-uesofthediffusioncoefficient.IncreasingD,asexpectedtheerror decreases,becausetheresidualbecomesprogressivelymore dom-inated by the elliptic term. Asymptotically, the error is inversely proportionalto D.Ontheother hand,theconstant solutionis al-wayspreserved(uptonumericalprecision)whenPφ=Pp. 7. Testcases

Inthissection,weverifyandbenchmarkournumericalscheme withthreetestcases.First,weemployamanufacturedsolutionfor theRANS and_{K− E equations}to verifythe spatialandtemporal convergenceratesof thesolver. Then, wesimulate thestationary turbulentflowoverabackward-facing stepandfinallya Von Kár-mánvortexstreet inthewakeofasquare cylinder.Inbothcases, wecompareourresultswiththoseobtainedfromexperimentsand othernumericalsimulationsreportedinliterature.

The meshes were generated with the open-source tool

Gmsh

[53],whichisalsousedtopost-processthesolutionfields.

7.1. Manufacturedsolution

Asfirsttest-case,weemploythemethodofmanufactured solu-tionstoverifythecorrectnessoftheimplementationofour space-timenumericalscheme.Startingfromthewell-knowTaylor-Green vortexsolutionofthelaminarNavier-Stokesequations(see,for ex-ample,[16]),wegeneralizedittotheRANSequationsandincluded anexpressionfortheturbulencequantities.Theanalyticalsolution is

uex_=exp

_−2t˜



− cos

(

x˜

)

sin

(

y˜

)

+sin

(

x˜

)

cos

(

y˜

)



, pex₌₋1 4 exp

−4t˜

(

cos

(

2x˜

)

+cos

(

2y˜

)

)

,

Kex_=exp

_−2t˜

_cos

₍

_x˜

₎

_cos

₍

_y˜

₎

_{− 4.5,}

Eex_=log

₍

_C μ

)

+Kex ,

ν

ex t =exp

(

Kex

)

, (49) where ˜ t:= 50

ν

t

(

L/

π

)

2 , x˜:= x L/

π

, y˜:= y L/

π

. (50)

Thissolution is definedon the domain

(

x_,y

)

_∈[_−L,L]2 _{andt≥ 0.} Itdoesnot satisfythecoupledRANS andK− E equationsexactly, sotheextrasourcetermsf,q_Kandq_E arenon-zero.Wecomputed

Table 4

Lid-driven cavity test: L 2 -norm error in the scalar quantity _{φwith respect to the exact solution φ}

0 = 1 , as a function of the diffusion coefficient ( D ) and for

various meshes and polynomial orders, which affect the divergence ( E D ) and continuity ( E C ) errors of the velocity field. Note that φ0 = 1 always lies in the

solution space. The scalar quantity is only conserved when P φ= P p .

M1 ( P u = 2 , P p = 1 ) M2 ( P u = 2 , P p = 1 ) M2 ( P u = 3 , P p = 2 )

ED = 5.1e-2, E C = 1.2e-3 ED = 1.1e-2, E C = 1.9e-4 ED = 9.0e-3, E C = 7.5e-5

P φ= P p P φ= P u P φ= P p P φ= P u P φ= P p P φ= P u

D = 1.0e-07 1.55E-14 2.53E-01 6.94E-14 8.74E-02 1.61E-13 2.74E-02 D = 1.0e-04 1.04E-14 1.98E-01 1.39E-14 3.01E-03 2.32E-14 1.07E-03 D = 1.0e-01 2.11E-13 8.36E-05 2.79E-13 2.86E-06 1.69E-13 1.06E-06 D = 1.0e+01 2.03E-11 8.37E-07 1.00E-12 2.86E-08 4.63E-12 1.06E-08 D = 1.0e+03 2.26E-13 8.37E-09 8.76E-13 2.86E-10 2.37E-12 1.07E-10

Fig. 2. L 2 -norm errors of the velocity, pressure, and turbulence quantities with respect to the exact solution _{(49) as a function of the time step size and for both first and}

second order BDF schemes. The correct convergence order are reached.

them by substituting (49) into Eqs. (1a-1b) and (4a-4b) symbol-ically with a Python script. We imposed Dirichlet conditions at allboundariesofthedomain.Togetherwiththeinitialconditions, theywereevaluatedfromtheanalyticalsolution.

Allrelative errorsreportedinthefollowingwereevaluated for

L₌1_,

ν

₌2_.5_{× 10}−4_, andat t₌4_,which corresponds to almost athreefold decayof theinitial condition. We seta relative toler-anceof10−12fortheGMRESandCGKrylovsolvers.Moreover,all integralswere evaluated with a polynomial accuracy of 20, con-siderablyhigherthanthe usual3Pu− 1.We didthistominimize thequadratureerrorintheintegration oftheextrasource terms, whicharecomplicatednon-polynomialfunctions.

7.1.1. Temporalconvergence

To verify the temporal convergenceof our numerical scheme, wesolvedEqs.(1a-1b) and(4a-4b) onastructuredmeshconsisting of100by 100quadrilateralelementsadoptingapolynomial order

Pu=5,inordertoensure thatthetime errorwasdominant.We carriedoutthestudyforboththefirst-orderandsecond-orderBDF schemes,progressivelyhalvingthetimestepsizefrom2−6to2−11.

Fig.2 showstherelative errorsintheL2_{-normofthevelocity,} pressure,andturbulencequantities withrespectto theanalytical solution(49) asa functionofthetimestepsize.Thedashed lines helpverifyingthat theerrorconvergencerateistheexpectedone forbothBDForders.Aslightdeviationfromthesecond-order con-vergencetrendcanbenoticedinthepressureandvelocityplotsat small



t,duetotheonsetofthespace-errorsaturation.

Given theavailability of an exact expressionfor theturbulent viscosity,we couldcompute an L2 _{erroralsoforthisquantity.As}

Fig. 3. L 2 -norm error of the eddy viscosity with respect to the exact solution

(49) as a function of the time step size and for both first and second order BDF schemes. The error exhibits the expected convergence rate.

shown in Fig. 3, thiserror also exhibits the correct convergence rates.

7.1.2. Spatialconvergence

The spatial convergence was verified by solving Eqs. (1a-1b) and (4a-4b) on increasingly refined structured quadrilat-eral meshes, progressively halving the element size from 2−1 to 2−5. The study was repeated for various polynomial orders (Pu∈

{

2,3,4

}

).Toensurethespatialerrorwasdominant,wechose theBDF2schemeandatimestepsize



t=2.5× 10−5_.

Fig. 4. L 2 -norm errors of the velocity, pressure, and turbulence quantities with respect to the exact solution _{(49) as a function of the element size and for several orders of}

the polynomial discretization. Scalar quantities exhibit the theoretical convergence rate, whereas that of the velocity is sub-optimal.

Fig. 5. L 2 -norm error of the eddy viscosity with respect to the exact solution

(49) as a function of the element size and for several orders of the polynomial discretization of K. The error exhibits the expected convergence order.

Fig.4 showstheL2_{-normerrorsoftheunknownquantitieswith} respecttotheanalyticalsolution(49) asafunctionoftheelement sizeandthepolynomialorder.The errorsofthepressureandthe two turbulencequantities convergewiththe expectedtheoretical order(highlightedbythedashedlines).Ontheotherhand,the ve-locityerror exhibitsasub-optimal convergencerate, inparticular forthesecond-orderandthird-orderdiscretizations.Thisisdueto thefactthattheresidualofthemomentumequationisdominated bytheerrorintheturbulentviscosity,whichconvergesatthe the-oreticalrateforthescalarquantities,asshowninFig.5.

Tocorroboratethisconclusion,we repeatedthetest,thistime imposingtheexactexpressionoftheeddyviscosity(fromEq.(49)) in the momentum equation, thus decoupling the RANS set from the turbulence equations. Results reported in Fig. 6 demonstrate that thetheoretical convergencerateofthe velocity errorcanbe recoveredinthiscase.

7.2.Flowoverabackward-facingstep

As second test case, we study the turbulent flow over a backward-facing step. The problem is a well-known benchmark forturbulencesolvers andmanyvariants havebeen investigated. Here,we consider thesame configurationasin theexperimental studyby [54]2_{, depictedin Fig. 7. The ratio ofthe inlet channel}

height (Hi) to step height (H) is 2 and Re=44737, based on the reference free-stream velocity U0, the step height and the kinematicviscosity

ν

.

Weimposed Dirichletboundaryconditionsattheinlet,placed ata distanceLi=4H upstreamofthestep.Weusedan analytical profileforthe stream-wisevelocity (ux)that matches the experi-mentaldataprovidedbyKimetal.[54] (reasonably)well,

uD x =U0f

(

y

)

, with f

(

y

)

=



1− e−18(y+0.0495)_{, 0≤ y≤ 1} 1− e−18(2+0.0495−y)_{, 1<y≤ 2}, (51)

2 Experimental measurements by Kim et al. _{[54] have been retrieved from the}

turbulent flow database described in [55] and available at https://turbmodels.larc. nasa.gov/bradshaw.html .

Fig. 6. L 2 -norm errors of the velocity and pressure when Eq. (1) is solved imposing the analytical expression for the eddy viscosity from _{Eq. (49) . The optimal convergence}

rate of the velocity is recovered in this case.

Fig. 7. Backward-facing step domain (not in scale).

Fig. 8. M3 structured mesh used for the turbulent backward-facing step test case. It consists of 55,500 elements.

andwe set a turbulence intensity of7.5% and a ratio

ν

t/

ν

=35, withwhichwecomputeduniformprofilesfor_KD_andED_.Outflow boundaryconditionswere imposed atadistanceLo=30H down-streamofthestep.Finally,weimposed wallfunctionsonthetop andbottom walls, witha distance from the physicalwall of the computationalboundarysetto

δ

w=0.04for



W3andto

δ

w=0.02 fortheothers.

We computed the steady-state solution on three structured meshes, finer towards the wall boundaries, consisting of 12,420 (mesh indicated with M1 in the following), 29,700 (M2), and 55,500(M3) quadrilaterals.Aportionofthelattermeshisshown inFig.8.Tostudytheinfluenceofthe polynomialorder,we per-formedcalculationswithbothPu=3andPu=4.

Fig.9 showsthesteady-statevelocity,turbulentkinetic energy, andeddyviscosity fieldsobtainedonthe finestmeshfor _Pu=4. Theyare in good qualitative agreement withthose, for example, reportedbyKuzminetal.[37].Foramorequantitativeassessment of our results, Table 5 compares the dimensionless recirculation length (LR) we obtained with the experimental measurement by Kimetal.[54] andothernumericalpredictions.Ourresultsdiffer

Table 5

Comparison of the dimensionless recirculation length ( L R ) obtained in the present

work (for different meshes and polynomial orders) with experimental and numeri- cal results (obtained with the k −model) reported in literature. The relative error is computed with respect to the experimental measurements by Kim et al. [54] .

Reference LR / H Relative error

Kim et al. [54] , experiment 7.0 ± 0.5

Zijlema [56] 5.90 −15.7%

Thangam & Speziale [57] 6.25 −10.7% Muller et al. [58] 6.33 −9.6% Kuzmin et al. [37] 5.40 −22.9% Shih et al. [59] 6.35 −9.3% Ilinca et al. [35] 6.20 −11.4% Present, P u = 3 , M1 6.42 −8.3% M2 6.54 −6.6% M3 6.56 −6.4% P u = 4 , M1 6.47 −7.6% M2 6.56 −6.2% M3 6.60 −5.7%

onlyafewpercentatmostfromthemostaccuratepredictions re-portedinliterature, andalways towards abetter agreementwith the experimental results. In all cases, the recirculation length is underestimated by around 6%, which is expected with the k₋



model.DifferencesinLRonM2andM3arelimitedto0.5%atmost, andtheimpactofaricherpolynomialrepresentationislimitedto lessthan1%, thus showingthemeshconvergence ofour calcula-tions.

Finally, Fig. 10 shows the vertical profiles of the stream-wise velocitycomponentatdifferentlocationsdownstream ofthestep, obtainedwithPu=4andfordifferentmeshes.ProfilesonM2and

M3 are barely distinguishable. Results compare reasonably well with the experimental measurements of Kim et al. [54] up to

x/H=8.0, then the discrepancies due to the limitations of the

k−



modelbecomeapparent, asalreadydocumentedin[35],for example.

7.3. Vortex-sheddinginthewakeofasquarecylinder

The last simulation is of turbulent flow past a cylinder with squarecrosssection,characterizedbytheappearanceofaVon Kár-mánvortexstreetinthewakeoftheobstacle.Eveniftheflow fea-tures stochastic three-dimensional fluctuations, whichcan be re-solved only with LES or DNS approaches (see, e.g., [60,61]), the problemisusedasatwo-dimensional test caseforRANS solvers. Theflowisinfactcharacterizedbyameanperiodicitythatcanbe capturedwithunsteadyRANSequations.

Fig. 9. Turbulent flow over a backward-facing step: steady state (a) velocity and streamlines, (b) turbulent kinetic energy, and (c) eddy viscosity fields obtained on the M3 mesh for P u = 4 .

Fig. 10. Turbulent flow over a backward-facing step: Comparison with the experimental measurements of Kim et al. [54] of the stream-wise velocity profiles obtained at different locations for P u = 4 and the three meshes employed. Results on the two most refined meshes are barely distinguishable.

Fig. 11. Geometry of the flow past a square cylinder test case (not in scale).

Fig. 11 shows the problem domain. The Reynolds number is

Re=22000,basedonthesquarecylinderside(D),thefar-field ve-locity(U_∞),andthekinematicviscosity

ν

.Experimental measure-mentsonthisproblemwerecarriedoutbyLynetal.[62].

Symmetry conditions were imposed on the top and bottom boundaries (distant H=30D), whereas an outlet was placed at

Lo=40Ddownstreamoftheobstacletoavoidanyinfluenceonthe flowinthecylinderregion.Forthesamereason,theinlet bound-arywasplacedatLi=20Dupstreamofthecylinder.Thefollowing uniformDirichletconditionforthevelocitywasimposed:

uD

(

t

)

=

(

U∞f

(

t

)

,0

)

T, with f

(

t

)

=



_t2

4

(

3− t

)

, 0≤ t≤ 2 1, t>2 .

(52) Moreover, following[63], wecalculated theinlet valuesofK and E tomatcha turbulenceintensitylevelof2%andaviscosityratio

ν

t/

ν

=10,whichcorrespondtotheexperimentsofLynetal.[62]. Finally,wallfunctionswere imposedontheedgesofthecylinder, settingthedistanceofthecomputationaldomainfromthephysical wallto

δ

w=0.015.

Wesampledthetimevariationsinthedragandliftcoefficients, definedas Cd= 2Fx DU2 ∞ and Cl = 2Fy DU2 ∞, (53) where Fx and Fy indicate the stream-wise and vertical compo-nents of the force acting on the cylinder. The mean periodicity

Fig. 12. M3 unstructured mesh used for the turbulent flow past a square cylinder test case: complete overview (left) and detail around the obstacle (right). It consists of 64,963 triangles, most of which are concentrated in the proximity and the wake of the cylinder.

Fig. 13. Flow past a square cylinder: Instantaneous velocity (top) and eddy viscosity (bottom) fields obtained on M3 for P u = 3 .

Table 6

Comparison of the force coefficients (in terms of average value, root mean square of the fluctuations, and frequency) obtained in the present work for three meshes and time step sizes with experimental and numerical results reported in literature.

Reference C¯d Cd,rms Cl,rms St

Lyn et al. [62] , experiment 2.100 − − 0.132

Trias et al. [61] , DNS 2.180 0.205 1.710 0.132

Bosch and Rodi [63] , k − 1.637 0.002 0.305 0.134

Rodi et al.; Voke [64,65] , RANS models 1.64–2.43 ≈ 0–0.27 0.31–1.49 0.134–0.159

Present, M1 t = 1 . 25 × 10 −3 _{1.687 0.008 0.460 0.129}

M2 t = 1 . 25 × 10 −3 _{1.775 0.014 0.632 0.131}

M3 t = 1 . 25 × 10 −3 _{1.788 0.015 0.657 0.132}

M3 t = 2 . 5 × 10 −3 _{1.788 0.015 0.657 0.132}

M3 t = 5 . 0 × 10 −3 _{1.788 0.015 0.657 0.132}

of the flow can be quantified in terms of the Strouhal number

St= f D/U∞,wherefistheoscillationfrequency,whichwe deter-minedfromthetimetrendofCl.

WeperformedourcalculationswiththeBDF2timeschemeand



t=1.25× 10−3_{.Tostudytheinfluenceofspaceerrorsonthe} re-sults,weusedthreeunstructuredmeshesconsistingof40,281 (in-dicated withM1),51,913 (M2), and64,963(M3,showninFig.12) triangles,progressivelymorerefined towardstheobstacle.Onthe finestmeshthen,wevariedthetime stepsize to



t=2.5× 10−3 and



t₌5_.0_{× 10}−3.

Fig.13 showsanexampleoftheinstantaneous velocity magni-tude andeddy viscosityfields obtainedfor



t=1.25× 10−3 _and

Pu=3on M3. The vortex sheddingis clearlyvisible. Table 6 re-portstheresultsweobtainedintermsofStrouhalnumber,average

Cd, androotmeansquare valuesofthe fluctuationsofCd andCl, comparingthemwithexperimentalmeasurements,DNSandother RANSresults.Weendedoursimulationsattend=152, correspond-ingtoalmost2.5flow-throughperiodswhichissufficienttoreach flowperiodicity, andevaluatedthe quantitiesof interestover the final8oscillationperiods.

ResultsonM3 andfor



t₌1_.25_{× 10}−3 arefullyconvergedin time (they do not change even quadrupling the



t) and nearly meshindependent. The Strouhal numberagrees withthe experi-mentalmeasurementsandtheDNSpredictions,anditdiffersonly 1.5%fromthevaluereportedbyBoschandRodi[63] forastandard

k−



model. The discrepancyis much more relevantwhen com-paringtheother quantitiesof interestthough. Thismightbe due tothe great sensitivityof theforce fluctuations tothe numerical