Data-driven modelling of the Reynolds stress tensor using random forests with invariance

2020

Document Version

Final published version

Published in

Computers and Fluids

Citation (APA)

Kaandorp, M. L. A., & Dwight, R. P. (2020). Data-driven modelling of the Reynolds stress tensor using

random forests with invariance. Computers and Fluids, 202, [104497].

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

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.

Data-driven

modelling

of

the

Reynolds

stress

tensor

using

random

forests

with

invariance

Mikael

L.A.

Kaandorp

a,b,∗

_,

_Richard

_P.

_Dwight

a

a Aerodynamics, Faculty of Aerospace, Delft University of Technology, 2629 HS, Delft, the Netherlands

b Institute for Marine and Atmospheric research Utrecht (IMAU), Utrecht University, 3584 CC, Utrecht, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 9 August 2019 Revised 31 December 2019 Accepted 4 March 2020 Available online 5 March 2020

Keywords:

Turbulence modelling Machine-learning Random forests

Reynolds anisotropy tensor Non-linear eddy-viscosity closures

a

b

s

t

r

a

c

t

Anovelmachinelearningalgorithmispresented,servingasadata-driventurbulencemodelingtoolfor ReynoldsAveragedNavier-Stokes(RANS)simulations.Thismachinelearningalgorithm,calledtheTensor BasisRandomForest(TBRF),isusedtopredicttheReynolds-stressanisotropytensor,whileguaranteeing Galileaninvariancebymakinguseofatensorbasis.By modifyingarandomforestalgorithmtoaccept suchatensorbasis,arobust,easytoimplement,andeasytotrainalgorithmiscreated.Thealgorithm istrainedonseveralflowcasesusingDNS/LESdata,andusedtopredicttheReynoldsstressanisotropy tensorfornew,unseenflows.Theresultingpredictionsofturbulenceanisotropyareusedasaturbulence modelwithinacustomRANSsolver.Stabilizationofthissolverisnecessary,andisachievedbya contin-uationmethodandamodifiedk-equation.ResultsarecomparedtotheneuralnetworkapproachofLing etal.[29].ResultsshowthattheTBRFalgorithmisabletoaccuratelypredicttheanisotropytensorfor variousflowcases,withrealizablepredictionsclosetotheDNS/LESreferencedata.Correspondingmean flowsforasquareductflowcaseandabackwardfacingstepflowcaseshowgoodagreementwithDNS andexperimentaldata-sets.Overall,theseresultsareseenasanextsteptowardsimproveddata-driven modellingofturbulence.Thiscreatesanopportunitytogeneratecustomturbulenceclosuresforspecific classesofflows,limitedonlybytheavailabilityofLES/DNSdata.

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

1. Introduction

The last few years have seen the introduction of supervised machine learning(ML)algorithms astoolstoexploit dataforthe purpose ofmodeling turbulence.RANS models are currently, and areexpectedtoremainthenormforsimulatingturbulentflowsin most industrial applications [46],because of their computational tractability, but suffer from poor accuracy and predictive power in a variety of important flows [9]. While a variety of nonlin-eareddy-viscositymodels(NLEVMs)andReynolds-stresstransport (RST) models have been developed using traditional techniques, it is the simplest lineareddy viscositymodels such asthe k₋



modelandk−

ω

modelsthatremainbyfarthemostwidelyused. This has motivated some to advocate alternative modelling ap-proaches that utilizeavailable referencedatamore fully,andrely lessonphysicalreasoning[10].Supervisedmachine-learning meth-odsdevelopedinotherfields,are– inthebestcase– abletodistill

∗ _{Corresponding author at: Institute for Marine and Atmospheric research Utrecht}

(IMAU), Utrecht University, 3584 CC, Utrecht, the Netherlands.

E-mail address: m.l.a.kaandorp@uu.nl (M.L.A. Kaandorp).

largedata-setsintosimplefunctionalrelationships.Thisoffersthe hopeofsubstantiallyimprovingcurrentRANSmodels,bybuilding closurescustomizedforaparticularclassofflowsbasedon appro-priatereferenceLESorDNSdata-sets[29,51].Howeverthereexist significantobstaclestorealizingthesemodelsinpractice.

Firstly,ahighsensitivityofthemean-flowtothedetailed char-acter of the turbulence hasbeen reported – even inthe case of channelflows[49].Thisplacesstringentaccuracyrequirementson anydata-derived closuremodel. Secondly,there exists nounique map from local mean-flow quantities to the needed turbulence statistics, due to non-local and non-equilibrium physics. Thirdly, any closure model should incorporate basic physical constraints, such asGalilean invariance: readilyachievable inanalytically de-rivedmodels,butdifficultwhen employingML procedures which generate arbitrary functions. Fourthly, a ML model should pro-ducesmoothfields(theymustbeincorporatedintoaPDEsolver), yet able to captureflows withrapid spatial variations, asin e.g. boundary-layers.Finally,RANSsolversarenotoriouslyunstableand difficulttodrivetoconvergence,evenwithstandardclosure mod-els. Stabilization of a data-derived model is therefore necessary. These challenges are in addition to the standard ML challenges. https://doi.org/10.1016/j.compfluid.2020.104497

1. Flow chart of the machine learning framework.

Fig. 2. Barycentric map (transformation of Lumley triangle), for a turbulent square duct at Re = 3500 . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Turbulence anisotropy state in the square duct (upper-right quadrant), visualized with the RGB colormap ( Fig. 2 (a)). TBNN = Tensor-Basis Neural-Network; TBRF = Tensor-Basis Random-Forest. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Thinkforexampletrainingagainstlargevolumesofdatawithhigh dimensionality[6],andavoidingoverfitting,whichcanbedoneby usingforexamplecross-validationmethods,ormethodsspecificto thealgorithmused,e.g. randomlydropping connectionsinneural networks[48],orsimplifyingdecisiontreesbypruningthemback

[14].

In this paper, a newapproach is presented that addresses all thesechallengestosome extent,resultinginaclosuremodelthat significantly outperforms a baseline RANS model for a specified class of flows. The model is closely related to nonlinear eddy-viscosity models (NLEVMs), of e.g. [39], in which the normal-izedReynolds-stressanisotropy tensorispredictedfromthe local mean-flow. We integrate a tensor basis forthe anisotropy intoa modifiedrandom-forestmethod,resultingintheTensorBasis Ran-dom Forest (TBRF), analogously to the Tensor Basis Neural Net-work(TBNN) ofLing etal.[29].Galilean invariancecan therefore be guaranteed; and in comparison to artificial neural networks, ourrandom forestsare easiertoimplement,easier totrain, have

fewer hyperparameters, andhave naturaltechniquesfor avoiding overfitting [16]. We introduce a methodfor stabilizing the RANS solver, using a combination of a modified k-equation, and relax-ation against a Boussinesq stress-tensor. With thissolver we can converge mean-flows to a steady state, with our TBRF closure model.

Manytypesof approachcan be identified inthe literaturefor improvementsof RANS closuremodelswithdata, weonly give a selection here.Broadly speakingthereare parametricapproaches, whichcalibrateorslightlymodifyexistingmodelsasin[7,12] (of-ten with statistical inference); and structural approaches, which attempt to relax fundamental modelling assumptions, especially Boussinesq as in [29,52]. In the latter case, uncertainty quantifi-cationhasbeenusedtodeveloppredictivemodelsintheabsence ofdata,by incorporatingstochastictermsintendedtocapturethe effectofmodellingassumptions,see[13,50,57].

Withdataavailable,machine-learninghasbeenusedtocapture potentially complex relationships between mean-flow and

mod-Fig. 4. Flow cases used in this work; the mean flow from the reference DNS/LES solutions is shown here. The color represents the velocity magnitude, and streamlines are plotted. Clockwise rotating regions of separation are indicated by dashed lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Curved backward-facing step, [ b ] i j components from LES, RANS, TBRF, and TBNN. (For interpretation of the references to color in this figure legend, the reader is

Fig. 6. Stress type for the curved backward-facing step, visualized with the RGB colormap ( Fig. 2 (a)): (a) LES from Bentaleb et al. [3] ; (b) RANS k −ω, (c) TBRF, and (d)

TBNN. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

ellingterms. These approaches can largely be divided into those which solve inverse problems to identify local correction terms needed to match data, as in [11,36,45]; and those which apply machine-learningto directly predict modelterms based on local mean-flowonly,see[30].Themachine-learningmethodsusedare various,Lingetal.[29] useneuralnetworks;LingandCo-authors

[28,51] userandomforests;andWeatherittandSandberg[53] uses gene-expression programming to obtain a concise form of their model.Despitethepopularityof random-forests,existingauthors havenotincorporatedframe-invariance,orsolverstabilizationthat weintroducehere.Thesedevelopmentsarecriticalforthe robust-nessandpracticalityofthemethod.

This paper isstructured asfollows.Section 2 will discussthe methodologyforthe TBRF framework, which includesunderlying theory,theimplementation ofthe algorithm itself, andthe data-setsusedintheframework.Section3 willdiscusstheresultsfrom theframework, whichincludepredictions forthe anisotropy ten-sor, flow fields obtained after propagating these predictions into theflowfield,andsomeproposalsforfutureworkwithregardsto quantifyinguncertaintyofthemodeloutput.Section4 willpresent theconclusions.

2. Methods

When performing Reynolds-averaging of the incompressible Navier-Stokesequations,thecompleteeffectoftheunresolved tur-bulentfluctuationsonthemean-flowiscontainedinasingleterm,

∇

·

τ

(where [

τ

]i j:=u_iu_j is the Reynolds-stress tensor), which modifiestheNavier-Stokesmomentum equation [42]. Toobtaina turbulence closure in this setting,it is sufficient to specify

τ

in termsof mean-flow quantities u¯ etc. The methodological goal of thisworkistouseplentifulDNS/LESdatatoestimateanonlinear eddy-viscositymodel(NLEVM)ofthegeneralform

τ

 h

(

∇

¯u,...

)

,

whereh

(

·

)

isa functionofmean-flow quantitiesonly. This prob-lem can be cast as a standard supervised machine learning task

[35]. However we demand additionallythat h isinvariant tothe choiceofGalileancoordinateframe,andthatthetrainingprocess is computationally cheap and robust. The first is achieved with an integrity basis (Section 2.3), andthe second by using a mod-ifiedrandom forest (Section 2.5). The model for

τ

thus obtained isusedwithin a RANS solverto predictthe mean-flow;this pro-cedure requiressolver stabilization (describedin Section 2.7). Fi-nally,trainingandvalidationcaseswithDNS/LESdataarelistedin

Section3.1.

2.1. Outlineofframeworkfordata-driventurbulencemodeling

Our framework consists of a training phase and a prediction phase, see Fig. 1. In the training phase, high-fidelity DNS/LES data is collected for a number of test-cases. The data should ideally contain full first- and second-order single-point statistics highlyresolved inthespatialdomain,fromwhichthenormalized anisotropytensor,b,iscomputed,seeSection2.2.Thesesame test-casesaresimulatedwithRANSusingthek₋

ω

closuremodel,and input features are derived from the mean-flow solution at every nodeofthespatialmesh.Themachinelearningalgorithmisthen trained to approximate the DNS ground-truth anisotropy tensor, fromtheRANS inputsover thefull solutionfields ofall training-casessimultaneously.

As we use a greedy algorithm to train the decision-trees, the training cost for a data-set of size N is O

(

Nlog2N

)

, so there is nopracticalrestrictiononthenumberofcaseswhichcanbeused intraining (indeeda much more severelimitation hasproven to be theavailabilityofhigh-qualityDNS/LESdata). However,we do notexpectthemapfromthemean-flowtounresolvedturbulence statistics tobe unique, even foravery limitedclass offlows. So, inordertocapturethisnon-uniquenessandtopreventoverfitting, multiple decision-trees are trained across random subsets of the data.

Inthepredictionphase,forapreviously unseenflow case,the anisotropytensorisestimatedusingthemeanoftherandom for-est predictions, with input from a baseline RANS mean-field. An updated mean-field is obtained by solving a modified RANS sys-temwiththemomentumequationsupplementedbythepredicted anisotropy.

Wewanttostressthatthisisacorrectiveapproach,witha sin-gle ML prediction providing an updated solution, similar to that practicedin [28,55].Inother words,itis nota replacementfora ’standard’turbulencemodel,wheretheReynoldsstressisassumed tobeaspecificfunctionofthemeanflowproperties,whichisthen iteratively solved until convergenceis achieved.Such an iterative approachwouldintheoryalsobeconceivable,inwhichML predic-tionsofbandmodifiedRANS formaclosedloop.Inthiscasethe MLshould betrainedonthe DNSmean-field, nottheRANS field. However such ambitious approaches are currently untested, it is unclear underwhatconditionsthe coupledsystemwill converge, and whetherthe convergedsolution will resembleground-truth. Thisishoweveranimportanttopicforfurtherresearch.

Adding to this point, it should be stressed that the approach hereis onlytested forsteadyRANS flow cases.It is untestedfor

Fig. 7. Barycentric map data at three sections (see Fig. 6 a) for the curved backward-facing step. Comparing LES [3] , RANS k −ωsimulation, TBRF, and TBNN.

Fig. 8. Barycentric map data at five sections for the backward-facing step. Comparing LES [25] , RANS k −ω, TBRF, and TBNN.

unsteady flows such aspresentin low-pressure turbines,see e.g.

[1].

Thegroundtruthforbcomingfromhigh-qualityDNS/LESdata will still contain some error due to e.g. dicretization errors, er-rorsduetoafiniteaveragingtime,andpossiblymodelerrors.The errors in the ground truth are assumed to be small here (com-paredtothemodelerrorsofthemachinelearningalgorithms),and thereforenottakeninaccount.Supervisedmachinelearning meth-odsexistwhichtheoreticallycouldbeabletocaptureuncertainties inthegroundtruth(seee.g.Gaussianprocessregression[41])

2.2. Reynoldsstressandrealizabilityconstraints

First we briefly review some basicproperties ofthe Reynolds stresses,relevant forconstructinga meaningfulML model.Firstly

τ

canbeintrinsicallydividedintoisotropicandanisotropic (devia-toric)parts:

τ

=2

3kI +a , (1)

where a is the Reynolds stress anisotropy tensor, k:=1 2trace

(

τ

)

is the turbulent kinetic energy, and I is the identity. It is the anisotropic part of the Reynolds stresses which is important: only this part is effective in transporting momentum, while the isotropicstressescanbesimplyabsorbedintoamodifiedpressure term[40]. Thenon-dimensional Reynolds stress anisotropy tensor,

b_,isdefinedas:

b :=₂

τ

_k−1₃I , (2)

and this is the quantity which we attempt to predict with ma-chinelearning.Intheremainder ofthispaper“anisotropy tensor” willrefertob.Tomodelb,eddy-viscosityclosuremodelstypically make the intrinsic assumption that b is function of local mean-flow quantities only. Linear eddy-viscosity models such as k−



andk−

ω

[23,54],goontomakethespecific,Boussinesq assump-tionthat b 1

2

ν

t

(

∇

¯u+

∇

¯uT

)

=

ν

tSˆ forsome scalar eddy viscosity

ν

t

(

x

)

, which remains tobe specified. Boththe intrinsicand spe-cificassumptions introduce modelling error. We aim to estimate

Fig. 9. Square duct, [ b ] i j components from DNS, RANS, TBRF, and TBNN. (For interpretation of the references to color in this figure legend, the reader is referred to the web

version of this article.)

andreducetheerrorinthelatterwithLESdatabasesand machine-learning.

Thepropertiesof

τ

andbleadtophysicalconstraintson mod-elsandmeansofvisualization.AmatrixAispositivesemi-definite if(andonlyif)xT_{Ax≥ 0,}

_∀

_x∈RN_{.Sincetheouter productofany} vectoru withitself(u u _{)ispositivesemi-definite;andsincethe}

Reynoldsstress isanarithmetic averageofsuchtensors,it isalso positivesemi-definite.Astrivialconsequences,alleigenvaluesof

τ

arerealandpositive,and

τ

αα≥ 0

∀

α

∈

{

1,2,3

}

, det

(

τ

)

≥ 0,

τ

2

αβ≤

τ

αα

τ

ββ

∀

α

=

β

.

(3) Thesepropertiesof

τ

haveimplications forb.Lettheeigenvalues of

τ

be

φ

_i,andthoseofbbe

λ

_i,then

λ

i=

φ

i 2k−

1

3, (4)

andboththeeigenvaluesanddiagonalcomponentsofbareinthe interval [₋1

3,23]. Furthermore using the Cauchy-Schwarz

inequal-ityin (3) the off-diagonal componentsof bare in[−1

2,

1 2].Since

trace

(

b

)

=0 only two independent invariants of the anisotropy tensor exist, e.g.: II:=[b]i j[b]ji and III:=[b]i j[b]in[b]jn. Therefore in the II-III plane all realizable states of turbulence anisotropy can be plotted, which are further restricted to a triangular do-main corresponding to the constraints onb justmentioned. This leads to the well-known “Lumley triangle” of Lumley and Co-authors[31,32] whichcapturestheanisotropicstateofturbulence. Thelesserknownbarycentricmapwasintroduced in[2],andisa transformationoftheLumleytriangleintobarycentriccoordinates, andwillbeusedforthepurposesofvisualizationandcomparison inthispaper,seeFig.2.

These triangles highlight three limiting states of turbulence anisotropy: 1-component turbulence (restricted to a line, one eigenvalue of b is non-zero), 2-component turbulence (restricted toa plane, twoeigenvalues ofbare non-zero),and3-component turbulence (isotropic turbulence,three eigenvalues are non-zero).

Fig. 2 shows these, along with invariants fora square-duct flow simulated with DNS from Pinelli et al. [38], and a k−

ω

RANS simulation.Forany2dlineareddy-viscosityRANS simulation,the predicted anisotropy invariants will lie entirely along the line

Fig. 10. In-plane mean velocity profiles at two sections of the square duct. Comparing: DNS, and the mean velocity fields obtained by propagating the DNS anistropy, and the TBRF-predicted anisotropy (labeled b ij,DNS and b ij,TBRF respectively). Also shown are two non-linear eddy-viscosity models.

indicated as“plane strain”. This illustrates the inabilityof linear eddy-viscositymodelstoadequatelyrepresentanisotropy.

One further method of visualization we will use is the Red-Green-Blue(RGB)map,inwhicheachanisotropicstateisassigned a color and the flow domain is colored accordingly, in a kind ofgeneralized contourplot[50].Fig.2(a) presentsthiscolormap, and in Fig. 3 the colormap is applied to the square-duct with DNS and RANS data (the same data used for Fig. 2(b-c)). The DNS data shows 1-component turbulence close to the wall, and 3-component nearthe centreline ofthe duct, whereas the RANS simulationonlyrepresentsturbulencenearthe3-componentlimit. Alsointhisfigure(c)and(d),machine-learningpredictionsofthe sameinvarientsareplotted,tobediscussedlaterinSection3.

2.3. Invarianceoftensor-valuedfunctions

The Navier-Stokes equations are Galilean invariant, i.e. un-changed by choice of inertialframe. It is a physical requirement that any model for the anisotropy tensor also be frame invari-ant,andtherebysatisfythesimplerequirementthatthefunctional model should notdepend on the choiceof coordinatesystem. In factthisprovestobeacriticalrequirementforthesuccessofour machine-learningstrategy,seeSection3.LetQ ∈R3×3 _bean

arbi-traryrealorthogonaltransformation;thenascalar-valuedfunction

f:R3×3_{→R,withtensorargumentS∈R}3×3 _{isframeinvariant if}

andonlyif:

f

(

S

)

=f

(

Q SQ T

)

,

∀

Q ,S . (5) Similarlyatensor-valuedfunctionh:R3×3_→R3×3_isframe

invari-antifandonlyif(e.g.[47]):

Q h

(

S

)

Q T= h

(

Q SQ T

)

∀

Q ,S . (6) One means of finding a h satisfying(6),is to start witha scalar function h:R→R, and specify that h is a natural generaliza-tion of h. See Higham[17] forseveral standard definitions for h

given h (e.g. in terms of power-series). Under all these defini-tions, notonly dowe haveh

(

XSX−1

)

=Xh

(

S

)

X−1 forany invert-ible matrix X (implying frame invariance by setting X≡ Q); but also other properties such as h

(

ST

)

=h

(

S

)

T _and

_λ

_=eigval

_{

_S

_}

_⇒

h

(

λ

)

=eigval

{

h

(

S

)

}

.

In addition we assume that h has a power-series representa-tion h

(

S

)

= ∞  i=0 a(i)

₍

_λ

₎

_S i_,

_λ

_=eigval

_{

_S

_}

_{, a}(i)_:R3_→R,

forsome scalar-valuedfunctionsa(i) _{whoseargumentsarethe}

in-variantsofS.Wereducethisinfinitesumtoafinitesumwiththe followingtrick:bytheCayley–Hamiltontheorem,everymatrix sat-isfiesitsowncharacteristicequationq

(

S

)

=0.Howeverqisa poly-nomialofdegree3(in3d),whosecoefficientsarefunctionsonlyof theinvariantsofS.Henceusingqwecanrecurrsivelyexpress pow-ersS3 andhigherinterms ofI, SandS2. Asa resulttheremust existanequivalentexpressionforh:

h

(

S

)

= 2  i=0 ˜ a(i)

(

λ

)

S i,

forsomedifferentscalar-valuedfunctionsa˜(i)_:R3_→R.

Inourapplicationweseeka functionfrommultipletensors to theanisotropy tensorb_,meaning ageneralization oftheabove is required.Theresultremains– undertheaboveassumptions– that themostgeneralhcanbewrittenintermsofafinitetensor-basis knownastheintegritybasis.

Inparticularwhenderivingnonlineareddy-viscositymodels,it issometimesassumedthattheReynoldsstressesareafunctionof the local, normalized,mean rates of strain S and rotation R. I.e.

b:=b

(

S,R

)

,where S =1 2 k





∇

u +

∇

u T



_{, R =}1 2 k





∇

u −

∇

u T



_{. (7)}

Inthiscase[39] thereare10integrity basistensorsT(m),making themostgeneralexpressionforb:

b = h

(

S ,R

)

=

10

 m=1

Fig. 11. Comparison of the streamlines for the backward facing step as given by ( a ) the RANS k −ωsimulation, ( b ) the propagated flow field using b ij,TBRF , and ( c ) DNS,

adapted from Le et al. [26] . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 12. Skin-friction coefficient for the backward facing step. Experimental data from Jovic and Driver [18] ; DNS data from Le and Moin [25] .

whereg(m) _{arescalarfunctionsoftheinvariants}

θ

i.Thebasis ten-sorsderivedfromSandRare[39]:

T (1)=S T (6)=R2_S+SR2₋2 3I· trace

(

SR 2

₎

T (2)=SR− RS T (7)=RSR2− R2 SR T (3)=S2₋1 3I· trace

(

S 2

₎

_T (8)=SRS2_{− S}2_RS T (4)=R2−1 3I· trace

(

R 2

₎

T (9)=R2S2+S2R2−2 3I· trace

(

S 2 R2

)

T (5)=RS2− S2 R T (10)=RS2R2− R2 S2R withinvariants

θ

1 =trace

(

S 2

)

θ

2=trace

(

R 2

)

θ

3=trace

(

S 3

)

θ

4 =trace

(

R 2S

)

θ

5=trace

(

R 2S 2

)

.

In [51] thisapproach is extendedto derive aset of47 invari-antsbasedon

∇

u¯,

∇

k,and

∇

p.Thisisthesystemweuseinthe following;thefullfeature-setwillbeshowninSection2.6.

2.4. Tensorbasisneuralnetwork(TBNN)

In [28] an artificial neural network was used to represent h. By careful choice of the network topology the idea of the ten-sorbasis isencodedintothenetwork.The networkcontainsa 5-node input layer receiving the5 scalar invariantsderived from S

and R. These inputs then go through a densely connected feed-forwardnetworkwith10hiddenlayers.Additionally,thenetwork contains an extra input layer consisting of the 10 base tensors

T(m). The output of the feed-forward network (representing the

g(m) _{functions)ismerged withthisextra inputlayer, reproducing}

(8).Therebythetensor-basis formofhandGalilean invarianceis achieved.

In this work TBNN is used as a competing method for com-parison. The implementation was obtained from the authors of

build-stopping was used, which terminates trainingwhen the training errorreduces,butthe validationerrorstartsconsistently increas-ing. Since the validation error as a function of the epoch has a noisybehaviour,amovingaverageoffivesampleswastakento de-termine when early-stoppingshould be activated. Initial network weights were randomly chosen, and the TBNN was trained five timesfromwhichthenetworkwasselectedwhichperformedbest onthevalidationset.

2.5. Tensorbasisrandomforest(TBRF)

Decisiontreesbasetheirpredictionsonaseriesofif-thentests ontheinput.Randomforestsconsistofcollectionsofdecisiontrees with some randomized component differentiating trees. Multiple decision tree algorithms exist, of which the CART (Classification AndRegressionTree)algorithmservesasthestartingpointforthe TensorBasisDecisionTree(TBDT),whichisusedintheTensor Ba-sisRandomForest (TBRF).AbriefoverviewoftheTBRFalgorithm ispresentedhere, foramoretechnicaloverviewthe readeris re-ferredtotheappendix.

In the training phase of the CART decision tree, the feature spaceisrecursivelysplitintotwobins.Ineachbinaconstantvalue isused toapproximatethetrainingdata.Givenpfeatures letthe trainingdataconsist ofX∈Rp×N _{point locationsinfeature-space,}

andcorrespondingy∈RN_{scalaroutputvalues.Eachsplitisaligned} with an input feature, and therefore the location of the split is completelyspecifiedbyasplittingfeatureindex j∈

{

1,...,p

}

⊂ N,

andvalue s. The twobins in whichthe datais splitare denoted

RL⊂ R(left)andRR⊂ R(right),andaregivenby

RL

(

j,s

)

=

{

X

|

[X ]j≤ s

}

RR

(

j,s

)

=

{

X

|

[X ]j>s

}

. (9) For the TBDT, constant values are chosen for the tensor ba-sis coefficients g(m) _{(see (8)) in each bin, which will be denoted}

by g(_Lm) and g_R(m) for RL and RR respectively. The values are cho-sensuchthatthemismatchwithrespecttothereferenceDNS/LES anisotropytensorbisminimized.Thecostfunctioncanbedefined as: J=  xi∈RL(j,s)







10  m=1 T (_im)g(_Lm)− bi







2 F +  xi∈RR(j,s)







10  m=1 T (_im)g(_Rm)− bi







2 F , (10) wherethe Frobeniusnorm isusedto calculatethedifference be-tween the reference DNS/LES anisotropy tensor, and the tensor resulting from the tensor basis. It can be shown by setting the derivative of Jwithrespect to g(m) _{ineach bin to zero,the}

opti-mumvalueforthe10tensorbasiscoefficientsineachbinisfound using g =



_N  i=1 ˆ T T_iT ˆi



−1



_N  i=1 ˆ T T_ib ˆi



. (11)

Inother words,during thetraining ofthe TBDT, each split in the treeis made by solving two least squaresproblems for g(m)_,

foreach j,andeachvalue s,andselectingthecombinationwhich minimizesJ.Theouter minimizationissolvedby brute-forceover features j, andone-dimensional optimizationis used over s. The Brent1d-optimizationalgorithm isused,offering agoodtrade-off betweenspeedandrobustness(byfallingbackonthegolden sec-tionsearch inthe worst case)[4].When the numberof samples ina bin fallsbelowa threshold, we switch to brute-forcesearch fors.Thisthreshold wassetto 150samplesto limit thetraining time of thedecision trees, a sensitivitystudy showedno signifi-cantinfluenceofthethresholdontheperformanceoftheTBRF al-gorithm.ThesplittingoftheTBDTbranchesisterminatedateither aspecified maximumbranchingdepth,orata minimumnumber ofsamplesper bin.Duetotheredundancyofthetensor-basisfor anygiven samplei∈

{

1,...,N

}

, (11) can become ill-posed, espe-ciallytowardstheleavesofthetree,whenonlyafewsamples re-maininabin.ThereforesomeL2_{-regularizationisaddedtoJwith}

coefficient



∈R+_{,c.f.ridgeregression,seetheappendixformore}

details.

In the tensorbasis random forest, multiple tensor basis deci-siontreestrainedonacollectionofbaggeddatasetsareaveraged. Baggingimpliesdataisrepeatedlyrandomlydrawn(with replace-ment) from the full data-set. Bagging is expected to work well whencombining a numberofhigh-variance, low-bias estimators, suchasthe decisiontree.Byaveraging manynoisy,butunbiased models, the variance of the prediction is reduced[16]. The vari-anceofthepredictions willbereducedmosteffectivelyifthe er-rorsinthecomponentmodelsareasfaraspossibleuncorrelated. Thisisencouragedbyintroducingsome additionalrandomness in theindividual trees:ateachsplit,notallfeatures,butarandomly selectedsubset of theavailable features is usedforsplitting.The specific parameters used in our computations will be stated in

Section3.

Instead of taking the mean over all the tensor basis decision trees,it proved to be more successfulto take the medianof the trees in the random forest (as forexample investigated in [43]), since this removed sensitivity to outliers in the predictions, see theappendixforacomparison.Sincetherandomforestisa piece-wiseconstantapproximationofb,andderivativesofbareneeded intheN-Sequation, thepredictionsfromthe TBRFare smoothed spatiallywithaGaussianfilter,beforetheyarepropagatedthrough thesolverto obtainaflow field (seeSection 2.7).TheTBRF algo-rithm has no explicit spatialcorrelation in the predictions since theseare basedonlocalfeatures oftheflow,so filteringthe pre-dictionswillintroducesomespatialcorrelation.Thefilterstandard deviationwassetto3celllengths,asthissufficientlysmoothsthe predictions while maintaining all important features of the pre-dictedtensor(seetheappendixformoredetails).Thisfilterwidth isanadhocchoice,andcanpossiblybeadjustedmorespecifically fornumerical stabilityin futurework by looking at e.g. required conditionnumbersforthesolver(seee.g.[56]).

Severalbenefits canbe identifiedinusingtheTBRF algorithm. Firstofall,theavailable datacanbe dividedintoatraining data-setandvalidationdata-setinanaturalmanner.Sincerandom

sam-  S  time scale  ∂p ∂xi ∂p ∂xi 1 2ρ∂∂xku¯ 2

k Ratio of pressure normal stresses to shear

stresses ¯

ui_∂∂xki † | u



j u k S jk| Ratio of convection to production of TKE

 u 

i u j k Ratio of total to normal Reynolds stresses

 u¯i ¯u j_∂∂ux¯ij



† u¯l ¯u l ¯u i_∂∂xu¯jiu¯k

∂u¯k

∂xj Non-orthogonality between velocity and its

gradient

plesaretakenfromtheavailabledatawithreplacementuntilthe originalsize ofthedata-setisreached, anumberofsamples will notbepresentinthetrainingdata-setforeachdecisiontree,called out-of-bag(OoB)samples.TheseOoBsamplescanthenbeusedto givea validationerror,or OoBerror.Duringtrainingofthe trees, thisOoB error can be observed to determine when trainingcan be stopped. Using the OoB erroris similar to performing N-fold crossvalidation[16].UsingtheOoBerrorallowsustooptimizethe numberoftrees intheforest duringtraining. While hyperparam-etersoftheTBRF weretuned(see Section3),thealgorithmis ro-busttothechoiceofhyperparameters.Itwillworkout-of-the-box withoutmuchtuningquitewell,seetheappendixformoredetails. Comparedtoneuralnetworks,therandomforestalgorithmis fur-thermoreeasytotrain,sinceonedoesnothavetoworryabout se-lectingtheappropriate optimizationalgorithmwhich hasitsown setofhyperparameters,anditsconvergence.

The TBRF algorithm presented here was implemented in python,forthesourcecodesee[20].

2.6.Choiceofinputfeatures

Under themodellingassumption thatthe Reynoldsstress ten-sorcanbewell approximatedusingonlythemeanstressand ro-tationtensors,S andR_, [39], the5 invariantsofthe tensorbasis

(8),namely

θ

=

(

θ

1,...,

θ

5

)

are sufficienttodescribe every

possi-bletensorfunction.Inthecontextofmachine-learning,thischoice ofinput features wasmadeine.g. [29].However, ifwerelax this assumption,thenitisreasonabletoalsouseotherquantities avail-ablein themean-flow asinputs, provided they are appropriately normalizedandGalileaninvariant.Inparticular,inthecaseofthe squareduct(seeSection3.1)itwasobservedthatduetothe sym-metryofthecasethereareonlytwodistinct“basisfunctions” de-finedby

θ

,andthesearenotsufficient toaccuratelydescribe the DNSanisotropytensorforthiscase.

Therefore here we will use the full set of invariants derived fromS,R,and

∇

kfromWangetal.[51].Inordertousethe turbu-lentkineticenergygradientitisfirstnormalizedusing√k_/



_,and thentransformedtoanantisymmetrictensor:

A k=−I×

∇

k. (13)

Furthermorenine extrascalarfeatures which aremore physically interpretable,suchasthewall-distancebasedReynoldsnumberare used,whichwereobtainedfromWuetal.[55] (whichwereinturn

basedonthosepresentedin[30]).Allfeatures whichareusedare presentedinTable1,wherefeatureset1(FS1)isbasedonS and

Ronly,featureset2(FS2)additionallyAk,andfeatureset3(FS3) areadditionallythe featuresfromWu etal.[55].Forthe features inFS3annormalizationfactorisincluded,whereasthetensorsin FS1 andFS2are normalizedusing kand



. Notethat all features inFS3arerotationallyinvariant, butsome (ortheirnormalization factors)arenotGalileaninvariantastheyincludetermsdepending onthevelocityoftheflow– thedistinctionismarkedwitha†in thetable.

2.7. Propagationofthepredictedanisotropytensor

TheopensourceCFDtoolboxOpenFOAMwasusedtocalculate RANSflowfieldsinthiswork.Thek−

ω

turbulenceclosuremodel wasused,together withthesecond-orderaccurate SIMPLE (Semi-ImplicitMethodforPressureLinkedEquation)scheme.

Simply setting the prediction of the anisotropy tensor bML in

themomentum equation adversely affects thenumericalstability ofthesolver.Asalreadyshownin[56],treatingtheReynoldsstress asanexplicitsourcetermintheRANSequationscanleadtoan ill-conditionedmodel.Twomainstrategiesareusedheretoimprove stability: (a) under-relaxingbML against the Boussinesq bB:=

ν

tSˆ with a relaxation parameter

γ

∈ [0, 1], and (b) simultaneously solving a modifiedk-equation to obtain a turbulence kinetic en-ergycorrespondingtothemodifiedanisotropy.

In detail, the incompressible Reynolds-averaged Navier-Stokes equationsare

∂

¯u

∂

t +¯u·

∇

¯u=

∇

·



− ¯p+

ν

S ˆ−

τ



(14)

where

ν

is the molecular viscosity. The prediction bML is

intro-ducedintothemomentumequation,bymodelling

τ

as

τ



τ

ML

(

γ

)

:=

2

3kI +2k[

(

1−

γ

)

b B+

γ

b ML]. (15) Theblendingparameter

γ

startsat0andisgraduallyincreased duringthesimulation,i.e.a continuationmethod,see e.g.[22].A linearrampbasedontheiterationcountnisused:

γ

n=

γ

maxmin



1, n nmax



,

where

γ

max≥ 0.8isachievedinalltest-casespresentedhere,and

itera-P=−

τ

:

∇

¯u, (16) isapproximatedinthek−

ω

modelbyreplacing

τ

withits Boussi-nesqapproximation[54].Hereweusethemodel

τ

MLfrom(15)

in-stead,includingtheblendingwith

γ

.

With these modifications, the solver converges for b-tensors originatingfromDNS,TBNNandTBRF.

3. Results

We comparepredictions oftheTBRF algorithm justdescribed, with baseline RANS (k−

ω

), DNS/LES references (withheld refer-ence data), andthe TBNN algorithm with the same feature sets as TBRF. Data for training and predicting comes from five flow caseswhichwillbediscussedbrieflyinSection3.1.Predictionsof theanisotropytensoritselfwillbepresentedinSection3.2; corre-spondingmean-flowpredictionsarepresentedinSection3.3.

3.1. Flowcases

Fiveflowcasesareusedintheframeworktotrainandtestthe machinelearningalgorithms.ForallflowcasesDNSdataorhighly resolved LESdatawasavailable.Themeanflowsforthereference DNS/LESsolutionsarepresentedinFig.4.Theflowcasesare:

(a) Periodichills(PH): FiveReynoldsnumbersare availablein the DNS/LES data-set from Breuer et al. [5], ranging from

Re=700toRe=10595basedonthebulkvelocityatthe in-letandthehillheight.

(b) Converging-diverging channel (CD): The DNS data for comes fromLaval andMarquillie [24] at Re=12600 based onthechannelhalf-heightandthemaximumvelocityatthe inlet.

(c)Curved backward-facing step (CBFS): The Reynolds num-beravailableisRe=13700basedonthestepheightandthe center-channel inlet velocity, withhighlyresolved LESdata fromBentalebetal.[3].

(d) Backward-facing step (BFS): Re=5100 based on the step heightandfreestreamvelocity.ThecorrespondingDNS sim-ulationcanbefoundin[26].

(e) Square duct (SD): Data-setsatmultiple Reynoldsnumbers are available from Pinelli et al. [38], with a total of six-teenrangingfromRe₌1100toRe₌3500basedontheduct semi-heightandthebulkvelocity.

Thefirstfouraforementionedcasesfeatureflowseparationand subsequentreattachment.Recirculationbubbles,non-parallelshear layers,andmean-flowcurvatureareall knowntoposechallenges for RANS based turbulencemodels. The square duct flow caseis symmetric;Fig.4(e)onlypresentstheupperrightquadrantofthe duct,wheretheflowintheductmovesout-of-plane.Prandtl’s sec-ondarymotionofthesecondkindisvisible,drivenbyturbulence anistropy.Assuchitisnotcapturedatallbylineareddy-viscosity models,whichmakes themidealforisolatingeffectsofnonlinear modelling[37].Forallcasesmeshindependencestudieswere per-formed for theRANS simulations, andmesheswere chosen such

caseshavebeenabbreviated,andthenumberbehindthe abbrevi-ation indicates theReynolds number. The table alsopresents the numberof samples usedfor training, N_sample (randomly sampled fromthetotaldata-set),andthenumberofusablefeaturespresent inthe trainingsets, Nfeature. From theavailable features the ones

withlowvariance(<1× 10−4_{)werediscarded,astheseeitherdid}

not contain anyinformation at all, orwere largely spurious. The startingfeature sets (FS) used are thosespecified in Table 1.For allcasesthek−

ω

turbulencemodelwasusedfortheRANS simu-lations.HyperparametersoftheTBRFweretunedforcasesC1,C2, andC4 using a validation set consistingof PH2800 andSD3200. A total of 100 TBDT’s were used to make the predictions. From the17availablefeatures11wererandomlyselectedforcalculating eachoptimalsplitintheTBDT.Theleafnodeswere settohavea minimumof9samples,andtheregularizationfactor



wassetto 1_{× 10}−12.ForcaseC3thesamesettingwereused,exceptthatthe treeswerefullygrown(i.e.eachleafnodeconsistsofonesample), andallfeatureswereusedtocalculatetheoptimalsplit.

First prediction for the curved backward facing step will be presented (C1), which is relatively similar to the training cases for which reliable data was available (periodic hills and the converging-divergingchannel).Next,resultsforthebackward fac-ingstepwillbepresented(C2),whichfeaturesstrongerseparation thantothetrainingcasesandwillthereforefeaturemore extrapo-lation.Lastly,resultsforthesquareductwillbepresented.A com-parisonwillbe madeforthecaseusingonlyfeatures basedonS

andRaswasalsodonein[29] (C3),andacasewhereallavailable featuresareused(C4).

3.2.1. Curvedbackwardfacingstep

For the curved backward facing step (case C1 in Table 2),

Fig.5 presentsthefournon-zerouniquecomponentsofb,asgiven bythe LESdata,the (k₋

ω

),andthe TBRF andTBNNalgorithms. TakingLESasa reference,RANS onlygivesacceptable predictions for the [b]12 component. By the Boussinesq assumption, results

willonlybeacceptable wherethe turbulenceanisotropy tensoris alignedwiththe meanrateof strain tensor, whichis empirically notagoodassumptioninthemajorityofthedomain.Incontrast, both machine learning algorithms give reasonable predictions of the tensorfield in all components, andpredictions are relatively smooth.Inparticular,featuresontopofthestepthe[b]11

compo-nentiscapturedqualitativelybytheTBRFandTBNNalgorithms,as wellasthe[b]33componentafterthestep.

Morerevealingisa plotofthestress typesin thedomain us-ingtheRGBmap,Fig.6.TheLESdatashows1-component turbu-lence ontop of the step,which is transportedinto theflow, be-yond the location atwhich the shear layerseparates. In [3] it is notedthatproductionofthestreamwisefluctuationisstrongly in-creasedattheshear-layerseparationlocation,leadingtoadditional 1-componentturbulence. As this shear layer gainsdistance from thewall,theturbulencerapidlyevolvestowardsthe3-component statedueto theredistributionprocess. Onthecurvedpartofthe step and the bottom wall, 2-component turbulence can be ob-served.In theremainder ofthe domain,3-componentturbulence

pletelymissed.Incontrast,theTBRFalgorithmaccuratelycaptures theturbulentstateasgivenbytheLESdata:1-component turbu-lence can be seen on top of the hill andat the separation loca-tion,it accurately predicts the 2-componentstate on the curved partof the walls and on the bottom wall after the step, and 3-componentturbulencecanbeobservedintherecirculationregion. Some noise is visiblehowever, mostnotably around x/h=0.0 to

x/h=1.0 away from the wall. The TBNN algorithm captures the differenttypesofturbulenceaccuratelyaswell. Closeto thewall ontop ofthe step it captures the1-component turbulencea bit lesswellthantheTBRFalgorithm,andsomespuriouspatternsare visibleabovethestepandclosetothelowerwallaroundx/h=6.0. To better quantify the accuracy of reconstruction, three sec-tionsthroughtheflowdomainareplottedinthebarycentricmap. Thesesectionsarelocatedatx/h=2,x/h=3,andx/h=4,which whichareatthefront,middle,andaftpartoftheseparatedregion. Thefirstsectionatx/h=2rangesfromy/h=0.5toy/h=1.5,the other twosections rangefromy/h=0.0 toy/h=1.5. Results are presentedin Fig. 7. As can be seen both machine learning algo-rithmsreproduce quitecloselytheLESreferencedata.They accu-ratelycapturethe 2-componentturbulence closetothewall, and move towards the3-component cornerwhen moving away from the wall. Discrepancies can be seen when moving upwards past tothewaketothechannelcenter,wheretheLESdataindicatesa movetowardsaxisymmetricexpansion,whichislessstrongly rep-resentedbytheMLalgorithms.

3.2.2. Backwardfacingstep

WeconsidercaseC2(c.f.Table2).FromLeandMoin[25] DNS dataisavailableforfivedifferentsectionsatspecifiedx/hlocations fortheReynoldsstressesandvelocities(histhestep-height). Loca-tionsonthebarycentric mapforthesefivesectionsareplottedin

Fig.8.Thesectionsrangefromy_/h=−1 (bottomwall)toy_/h=0 (thelocationofthestep).

Results are similartothose oftheCBFS:the machine-learning algorithms are able to give a qualitatively accurate prediction of theturbulencecharacter,withsomequantitativediscrepancies.For

x/h=4andx/h=6predictionsclosetothewallaremoreaccurate forTBRFthanTBNN. Thesituationisreversed forx/h=10,15,19,

whereTBNNslightlyoutperforms,atthecostofsomeunrealizable predictionsclosest to thewall. In all our studies, we have never observed unrealizable predictions from TBRF, despite no explicit realizabilityconstraintbeingimposedonthemethod.

Moving away from the wall into the shear layer TBRF erro-neouslyheadstoofarbacktowards thetwo-componentboundary atthesectionsclosesttothestep.Thereasonforthisisunclear,at similar(shear-layer)locationsinthetrainingflows,theturbulence doesnot exhibitsuch behaviour.FurthermoreTBNNisreasonably accuratehere. Diagnostictoolsareneeded, andwill bea focusof futureresearch.Nonetheless atthe section furtherfromthe step, bothMLmethodsperformwell.

C4 0.0521 0.0681

Table 4

Backward facing step, reattachment point locations. Model xreattach [x/h] RANS 5.45 RANS+ b ij,TBRF 6.32 DNS [26] 6.28 Experiment [18] 6.0 ± 0.15 3.2.3. Squareduct

Thelocalstresstypeforthesquare ductwasalreadyshownin

Fig. 3; individual components of the anisotropy tensorshown in

Fig.9.InbothcasesanistropyisvisualizedforDNS, RANS(k−

ω

), TBRFandTBNNpredictions.InFig.3 MLresultsareonlyshownfor caseC4(17features);inFig.3 additionallycaseC3isshown.Note thatthesearechallengingcasesduetothesubstantialdifferences betweenthetrainingandpredictionflows.

Asexpected, theBoussinesq modelyields non-zeropredictions onlyfor[b]12and[b]13– thoughthesearerelativelywellpredicted.

Anistropy is confined to the 3-component corner, on the plane-strainline.Examiningthepredictions ofML,itcangenerallyseen thattheintroductionofextrafeatureshassignificantlymoreeffect thanthechoiceofneural-networksversusrandom-forests.For ex-ample,lookingat[b]11,theanisotropyoftheReynoldsstressisnot

capturedclosetothewallsforC3,whereasitispresentinC4.The magnitudeof [b]12 and[b]13 is underpredicted,independently of

theMLmethod,butimprovedincaseC4comparedtoC3.Similarly in all cases the magnitudeof [b]23 is over-predicted by ML, but

themagnitudeof theover-prediction islessincaseC4. To quan-tifytheseobservations,the rootmeansquare error(RMSE)ofthe anisotropytensorwithrespecttotheDNSdataisgiveninTable3. TheRMSE’sarelowerwhenintroducingmorefeatures,forboth al-gorithms.Thiscanlargelybeexplainedbyvisualizingtheshapesof theinputfeaturesincaseC3.Doingthisitcanbeobservedthat,of the5features,3are approximatelyscaled versionsofthe other2 – effectivelyreducingtheinputspacetotwo-dimensions.This ex-plainsthedifficultynonlineareddy-viscositymodelsbasedononly

SandR,haveinreproducingthemagnitudeofthesecondaryflow inthesquareduct.

3.3. Anisotropytensorpropagation

This section presents flow fields obtained by propagating the predictedanisotropytensorsforthesquareductflowandthe back-wardfacingstep.Predictionsfortheanisotropytensorwere prop-agatedusingthestabilizedsolverpresentedinSection2.7.

Fig. 13. Illustration of the robustness of the TBRF algorithm.

3.3.1. SquareDuct

Twosectionsinthesquare ductwill beanalyzedwithrespect tothe in-planemeanvelocitymagnitude(i.e.indicatingthe mag-nitude ofthe secondary flow). Fig.10 presentsthe velocity mag-nitude forsections located at y/h=0.5and y/h=0.8.DNS from Pinelli et al. [38] is used as a reference. In order to verify the propagationmethodinisolation(withoutpredictinganistropy),the anisotropy tensor obtained directly from DNS (b_ij,DNS) is propa-gated, seethegraylines.Thisisa“bestcasescenario” wherethe anisotropytensorisassumedtobeperfect(uptostatistical conver-genceoftheDNS).Themean-flowfieldasobtainedbypropagating thepredictions fromtheTBRF algorithm(bij,TBRF,seecolumn5 of

Fig.9)isindicatedby theredlines.Furthermore,resultsfromthe quadratic eddy viscosity modelof Shih etal. [44], andthe cubic eddy viscosity model of Lien et al.[27] are presented. Since the lineareddy-viscosity modeldoesnot yieldanysecondary flow at all,thisresultisomitted.

WhenexaminingtheresultsofFig.10,thestoryisbroadlythe samefory/h=0.5andy/h=0.8.Mean-velocityfieldsobtained us-ing bij,DNSbroadlyreproduce theDNSmeanvelocity, both in

am-plitudeandlocationofkeyfeatures,withthebestfitnearthewall (z=1), andtheworst nearthechannel centerline (z=0). Subse-quently approximating bij,DNS by bij,TBRF causes additional errors,

but theses errors are of similar magnitude to the errors already made in the propagation. In particular, key features are correct, andamplitudesareappropriate.Whatisalsoclearhowever,isthat predictions arestill farmoreaccuratethan bothnon-linear eddy-viscositymodels.ThemodelfromShihetal.[44] isabletopredict thelocationofthepeaksofthein-planeflowmagnitudequite ac-curately,butsignificantlyunderpredictstheoverallmagnitude. Pre-dictionsbythecubiceddy-viscositymodelfromLienetal.[27] are faroff overall.

3.3.2. Backwardfacingstep

Fig.11 presentsstreamlinesinthe flowfield forthebackward facing step, withresultsfrom k−

ω

RANS, thepropagated veloc-ityfieldusingthepredictedanisotropytensorfromtheTBRF algo-rithm (bij,TBRF),and DNSdata fromLe etal.[26].The size ofthe

recirculationregionismorecorrectlypredictedforthepropagated velocityfieldusingb_ij,TBRFcomparedtothebaselineRANS simula-tion.Furtherawayfromthewallthesolver usingbij,TBRFdoesnot

introduce spurious effectsand results are similar to the baseline RANSsimulation.Thereattachmentpointlocations(x_reattach)forall threecasesarepresentedinTable4.Asignificantimprovementis shownforthepropagated velocityfield comparedto thebaseline RANSsimulation.

TheskinfrictioncoefficientsfromtheRANSsimulationandthe propagated flow field usingbij,TBRF arecompared to experimental

datafromJovicandDriver[18] inFig.12.Thepropagatedflowfield showsaveryclosematchtotheexperimentaldata,andthe major-ityofresultsfallwithintheerrorboundsgivenbytheexperiment

(_{± 0.0}005c_f).Thereattachmentpointofthepropagatedflowfield (6.32) comparesfavourably tothe experimentally found reattach-mentpoint(6.0± 0.15).

4. Conclusions

In this work, a novel random forest algorithm was intro-ducedforRANSturbulencemodeling,tomodeltheReynoldsstress anisotropy tensor. The algorithm was trained using invariant in-putfeaturesfromseveralRANS (k−

ω

)flowfields,andthe corre-spondingresponses fortheanisotropytensorfromDNSor highly-resolvedLES data.Galileaninvariance ofthepredictedanisotropy tensorisensuredbymakinguseofatensorbasis,derivedin[39]. The newrandom forest algorithmis calledthe Tensor-Basis Ran-domForest (TBRF) algorithm,similarlyto theTensor-Basis Neural NetworkfromLing etal.[29] fromwhich itwasinspired.Robust predictions ofthe Reynolds-Stressanisotropy tensor are obtained bytakingthemedianoftheTensor-BasisDecisionTree(TBDT) pre-dictionsinsidetheTBRF.

PredictionsfortheReynolds-stressanisotropy tensorwere pre-sentedforthesquareductflowcase,curvedbackward-facingstep, andbackward-facing step.Improvement is observedwith respect tothebaselinek−

ω

simulations,andtheTBRFalgorithmperforms onparwiththeTBNNalgorithm.ComparedtoTBNN,theTBRF al-gorithm is relatively easy to implement and train (one doesnot have to think about matters such as the optimization algorithm usedtotunetheneuralnetworkweightsanditsconvergence);the out-of-bagsamplesfromthedecisiontreesallowforanaturalway toquantifythevalidationerrorduringtrainingandthus selecting theamountoftrees tobeusedintherandomforest.The few re-maining hyperparameters are quite robust: the TBRF works well out-of-the-boxevenwhenusingstandardhyperparametersettings (fullygrowntrees,usingallavailable featuresforcreatingthe de-cisiontreesplits).

Acustomsolver forpropagatingtheanisotropy tensorwas in-troduced, which blends the predictions for theanisotropy tensor witha k−

ω

turbulencemodel. This solver greatly increases nu-merical stability of the propagation. Propagations for the square duct flow case and backward facing step are presented, which show a closematch with respect to corresponding DNSand ex-perimentaldata-sets.

ApossibilityforfutureworkmightbeusingtheTBRFfor quan-tifyinguncertaintyofthepredictionsaswell.Foreverylocationin theflowdomainthetreesintherandomforestcanbeanalyzedfor theirvariance,whichwouldmakeitpossibletouseananisotropy eigenvalueperturbationmethodologytoquantifyuncertaintysuch asproposedin[13].Inordertoachievemeaningfulboundsfor un-certaintyofthepredictions, onecould lookate.g. Bayesian Addi-tiveRegressionTrees[8],orjackknife/infinitesimaljackknife meth-ods[15],ormodifytherandomforestalgorithmitselffor meaning-fuluncertaintybounds,seee.g.[33,34].Itwasobserved,thatoften

influencetheworkreportedinthispaper.

The authors declare the following financial interests/personal relationshipswhichmaybeconsidered aspotentialcompeting in-terests:

CRediTauthorshipcontributionstatement

MikaelL.A. Kaandorp: Conceptualization, Methodology, Soft-ware, Investigation, Writing - original draft. RichardP. Dwight: Conceptualization, Methodology, Supervision, Writing - original draft.

Acknowledgments

The authors wouldlike to thank Javier Fatou Gómezfor sup-plying the OpenFOAM continuation solver used to propagate the Reynoldsstressesintothe flowfield;alsoJulia Ling forproviding herimplementationofTBNNandsupportforthiswork.

Appendix. TBRFimplementationdetails

Decisiontreesbasetheirpredictionsonaseriesofif-thentests ontheinput.Randomforestsconsistofcollectionsofdecisiontrees withsome randomized component differentiating trees. Multiple decision tree algorithms exist, of which the CART (Classification AndRegressionTree)algorithmisusedasastartingpointhere.

As mentioned in Section 2.5, the feature space is recursively splitintotwobins,RR andRL.Thesplittingisperformedgreedily, witheachsplitselectedtominimizethemismatchbetweenthe re-sponsey,andthebestconstant approximationoftheresponsein bothbins.Specificallyforeachsplitwesolve[16]:

min j,s



min cL∈R  xi∈RL(j,s)

(

y i− cL

)

2+min cR∈R  xi∈RR(j,s)

(

y i− cR

)

2



, (17)

whereyi denotes theresponse at Xi. Finding constantscL,cR∈R amountsto averaging y within RL andRR respectively, effectively minimizingthevariance inbothbins. Startingfromthefull data-setthe samemethodis then appliedto RL andRR in a recursive fashion. The procedure is terminated either at a specified maxi-mumbranchingdepth,oraminimumnumberofsamplesperbin. ThenewTBDTalgorithmiscomparablewiththeCARTdecision tree algorithm, but instead of approximating the response with constantvaluescL andcR,thealgorithmfindsaconstantvaluefor eachofthetensorbasiscoefficientsg(m)_{in(8),chosentominimize}

themismatchbetween thisexpression andthe anisotropy tensor fromDNS.Specificallywesolve

min j,s

⎡

⎣

min g(m) L ∈R10  xi∈RL(j,s)







10  m=1 T (_im)g(_Lm)− bi







2 F

optimizationiteration ofs.As forCART, (18) is repeatedforeach bin,untilastoppingcriterionisreached.

Explicitly,by flatteningthetensorateach point,anddefining:

ˆ T i=

⎡

⎢

⎢

⎢

⎣

[T (_i1)]11 [T i(2)]11 · · · [T (i10)]11 [T (_i1)]12 [T i(2)]12 · · · [T (i10)]12 . . . ... ... ... [T (_i1)]33 [T i(2)]33 · · · [T (i10)]33

⎤

⎥

⎥

⎥

⎦

, b ˆi=

⎡

⎢

⎢

⎣

[b i]11 [b i]12 . . . [b i]33

⎤

⎥

⎥

⎦

, (20) eachoftheminimizationproblemsovergbecomesmin_gJ where

J= N  i=1



_T ˆ_{ig − ˆb i}

_

2_{, (21)} withsolution g =



N  i=1 ˆ T TiT ˆi



−1



_N  i=1 ˆ T Tib ˆi



. (22)

which canbe solved separately forRL andRR toobtain g(Lm) and

g(_Rm).Theoverallcostofthisalgorithm(asforCART)isdomainated bysortingthedata-valueswithrespecttocoordinatej.Thiscostis

O

(

NlogN

)

inthenumberofdata-values,leadingtoanoverallcost ofO

(

Nlog2N

)

. Unliketrainingneuralnetworks, thisprocedure is fast, robust, easy to implement, andindependent of anystarting guess.

Dueto theredundancyof thetensor-basis foranygiven sam-plei∈

{

1,...,N

}

,thisproblemcanbecomeill-posed,especially to-wardstheleavesofthetree,whenonlyafewsamplesremainina bin.ThereforesomeL2_{-regularizationisaddedtoJwithcoefficient}



∈R+_{,c.f.ridgeregression.Tosummarize,(22) ismodifiedto}

g =



N  i=1 ˆ T TiT ˆi+



I



−1



_N  i=1 ˆ T Tib ˆi



. (23)

Inpracticeitwasobservedthatby takingthemedianofall deci-siontreesinsteadofthemean(see Section2.5),thisalready pro-videsalotofrobustness,andthat



canbesettoaverylowvalue ingeneral. The value of



inthe randomforest wastuned using validationdata-sets,seeSection3.2.

OneimportantdifferenceoftheTBRFoverthestandardrandom forestisthefactthatiteffectivelydoesnotdirectlypredictthe fi-naloutcome,butthecoefficientsg,whichmultiplythebasis ten-sorsT.Unlikethestandardrandomforest,thismeansthatthe val-uesforthefinalpredictionsdonothavetoliein-betweenthe val-uesofthepointsusedfortraining.Decisiontrees arewellknown to be sensitiveto smallchanges inthe data,andthis manifested during testingas highlyirregular and inconsistent predictions in smallregionsofthespatialdomain.ForasampleTBDTprediction