Intrusive generalized polynomial chaos with asynchronous time integration for the solution of the unsteady Navier–Stokes equations

ContentslistsavailableatScienceDirect

Computers and Fluids

journalhomepage:www.elsevier.com/locate/compfluid

Intrusive generalized polynomial chaos with asynchronous time integration for the solution of the unsteady Navier–Stokes equations

P. Bonnaire

^a

, P. Pettersson

^b,∗

, C.F. Silva

^a

aGroup of Thermo-Fluid Dynamics, Technical University of Munich, Germany

bNORCE Norwegian Research Centre, Norway

a rt i c l e i nf o

Article history:

Received 16 April 2020 Revised 5 March 2021 Accepted 29 March 2021 Available online 2 April 2021 Keywords:

Generalized polynomial chaos Stochastic Galerkin projection Navier–Stokes equations Kármán vortex street

a b s t r a c t

Generalizedpolynomialchaosprovidesareliableframeworkformanyproblemsofuncertaintyquantifi- cationincomputationalfluiddynamics.However,itfailsforlong-timeintegrationofunsteadyproblems withstochasticfrequency.Inthiswork,theasynchronoustimeintegrationtechnique,introducedinpre- viousworkstoremedythisissueforsystemsofODEs,isappliedtotheKármánvortexstreetproblem.

Forthispurpose,wemake useofastochasticclock speedthat providesthephaseshift betweenthe realizationsandenablesthesimulationofanin-phasebehavior.Resultsoftheproposedmethodareval- idatedagainstMonteCarlosimulationsandshow goodresultsforstatisticfieldsandpoint-wisevalues suchasphaseportraits,as wellasPDFsofthelimitcycle.We demonstratethatlow-orderexpansions aresufficienttomeetthedemands forsomestatisticmeasures.Therefore,computationalcostsarestill competitivewith thoseof thestandard form ofintrusivegeneralized polynomial chaos(igPC)and its non-intrusivecounterpart(NigPC).

© 2021PublishedbyElsevierLtd.

1. Introduction

Generalized polynomial chaos (gPC) [1–4] is an uncertainty quantification (UQ) approach devoted to the propagation of un- certainties frominputstooutputs ofanumericalmodel,e.g.,dis- cretized partial differential equations. In this method, a series expansion in orthogonal basis functions or modes is applied to the model, i.e., each stochastic variable is described as a linear combination of stochastic modes. The modes are generally or- thogonal polynomial functions of standardized random variables withknownstatisticalproperties.TheQuantitiesofInterest(QoIs), which aredirectly relatedtothe output(s)ofthe model,are also describedbyalinearcombinationofwell-knownstochasticmodes.

As aresult,theUQproblemreducestothefindingofthegPCco- efficients describing the QoIs. There are two classes of methods for their evaluation: Non-intrusive generalized polynomial chaos (NigPC)andintrusivegeneralizedpolynomialchaos(igPC).

NigPC has consolidated in the past years as an approach to performUQonsystemscharacterizedbycomplexnumericalmod- elsthatdescribesteadyQoIs[5–9].Non-intrusivemethodsrelyon random ordeterministic sampling. The formercase includesran- domdiscreteL²projection[10],sparsegPCbasedonleastanglere-

∗Corresponding author.

E-mail address: pepe@norceresearch.no (P. Pettersson).

gression[11],andcompressivesamplingforidentificationofsparse gPCapproximations[12].Inthecaseofdeterministicsamplingin- stead,theQoIisevaluatedatcarefullychosenpointsinparameter space,e.g.,Gausspoints.Oncetheinputrandomvariablesaredis- cretized– following one ofthetwo aforementioned methods– a collection of representative samples is generated in the stochas- tic space.Subsequently, associated deterministic numerical simu- lations are performed, where each simulation can be understood asarealization, i.e.,asampleofa stochasticsimulation.Provided thatthe QoIisasufficientlysmooth functionoftheuncertainin- putparameters,thenumberofdeterministicsimulationsneededis considerablysmallerthanin othersamplingbasedmethods,such asMonte-Carlo simulation orLatin-Hypercube methods. For very largenumbers ofrandomvariables, thesmoothnessrequirements andcomputational cost become prohibitive (the curse of dimen- sionality).Thiscan tosome extent be remediedbymodel reduc- tiontolimitthenumberofinputrandomvariablesandthenum- berofbasisfunctions,inparticularforellipticproblems[13].Fur- thermethodsusedtoimprovetherateofconvergencewithrespect tothenumberoftotaltermsintheseriesexpansionforlinear,el- lipticPartialDifferentialEquations(PDEs)withdifferentfunctional uncertaintiesinthediffusioncoefficientarepresentedinBachmayr etal.[14,15],whicharebasedontheresultsfrom[16].Anisotropic Smolyaktype polynomialspacescanalleviatethecurseofdimen- sionality underproper smoothnessassumptions ofthe parameter randomfields[17,18].

https://doi.org/10.1016/j.compfluid.2021.104952 0045-7930/© 2021 Published by Elsevier Ltd.

certaintiesintheSmagorinskycoefficient characterizedbyareac- tive LargeEddySimulation induce uncertaintiesinthemean and rms values ofvelocity, temperature andmixture fraction at spe- cific locations [9]. In contrast, few workshave been carried out, where the QoIs are steady fieldquantities. For example, Onorato etal.[21],whileusingthe2DTurbulent Navier–Stokesequations, investigated howuncertaintiesinthe angleofattack ofan airfoil induced uncertaintiesinthefieldofMachnumber.Congedoetal.

[6] investigated how uncertain inflow conditions induced uncer- tainties in the field ofmean velocity ofa turbulent swirled flow modeled by means ofboth Reynolds averaged Navier Stokes and LargeEddySimulation.Notethatinthosecasesthecomputational costs ofNigPCarehigh,asthegPCcoefficientsassociatedtoeach pointinthefieldofinterestmustbecomputed.

The intrusivecounterpart of gPC rose considerableinterest in the first decade of 2000’s [22–27]. Contrary to NigPC, igPC by means of stochastic Galerkin projection is not a samplingbased method: onlyone numericalsimulation ofan expandedmodelis neededtocalculatethegPC coefficients,whetherthesearescalars orfieldquantities.ThepracticaluseofigPCmaybelimiteddueto thefactthatstandardnumericalsolvers,oftentheresultofyearsof investigation, cannot beused withoutprior modification.Compu- tationalaspectsofigPCformulationsoftheincompressibleNavier–

Stokes equationsinthe laminarflow regimewere treatedinKnio andLeMaître[28],andboundaryconditionsforawell-posedfor- mulation anda stable numerical scheme were presented in Pet- terssonetal.[29].

Awell-knownlimitationofgPCisitsinaccuracyinuncertainty propagationofunsteadyproblems,wherethephasebetweenreal- izationschangesovertime,e.g.,turbulentflowsorflowsexhibiting periodic behavior. Intheory,unsteadystochastic problemscan be accuratelydescribedbyconsideringaninfinitenumberofstochas- tic modes[30–32].In practice,accountingforan infinite number of stochastic modesisinfeasible. Theenergy inperiodic systems, in contrast to turbulent flows, is usually distributed in a hand- ful of frequencies. As a result, the phase shift between realiza- tions of the stochastic periodicproblem canbe tracked. Twosit- uationsarise.(i)Deterministicfrequencies:caseswherethephase shift betweenrealizationisnegligible,so thatregular gPC isreli- ableandcanbeusedwithoutmodificationtostudyperiodicflows.

This fact is sometimes lost in claims in the scientific literature:

“gPC...cannotdealwiththeNavier–Stokesequationsforunsteady noisyflows,suchasflowpastastationarycylinder.Fortheseprob- lems gPC failsto convergeaftera shorttime”[33];“gPCtendsto break down forlong-time integration” [34]; “Awell-known diffi- cultyisthesimulationofuncertaintime-dependentproblemsover long times” [35]; “note also that polynomial chaos could intro- duce convergenceissues when dealingwith unsteady flow prob- lems”[6].Wewouldliketostressthattheprevioussentencesmay bemisleadingiftheentirecontextoftheassociatedstudiesisnot carefullyread.Therearenumerousproblemsregardingperiodicos- cillations, where regular gPC could be used to study uncertainty propagation. However, onlyfew studies,suchasLucorandKarni- adakis[24],havefocusedonthismatter.(ii)Stochasticfrequencies:

equation(ODE)that describesthe adjustmentofthe clock.Schick etal.[36]proposedamethodfortheadjustmentofthelocaltime in cases where the phase shift between realizations is irregular.

Schicketal.[36] arguethattheATImethodisnotrobust enough forthatpurpose,andthatitispreferabletousearobustoptimiza- tionmethodtodeterminetheclock adjustment,despiteincreased computational cost. Mai and Sudret [37] proposed time warping forphaseadjustmentforNigPC.Giraldietal.[38]alsousedNigPC andatimeshifttosynchronizewavearrivaltimeswithinanearth- quakemodel.Conceptuallyrelatedisthetechniqueofintroducing aphysicalcoordinatetransformationtoaligndiscontinuitiesinpa- rameterspace,appliedtomulti-layeredsedimentarybasins[39].A different approach wasproposed by Gerritsma etal. [34], where time-dependentgPCwasintroducedwithnewstochasticvariables andorthogonal polynomialsbeingconstructedastimeprogresses.

Consequently, the number of stochastic modes remains small as theassociatedbasisisnearly optimalfora sufficientlysmalltime interval.

Due to the constraints mentioned above, standard NigPC and igPC should be avoided for problems where the QoIs are peri- odicwithstochasticfrequencies,orturbulent.Accurateuncertainty propagation forturbulent or periodic (stochastic frequency) QoIs would bring significant advantages in many engineeringapplica- tions.We wantto stress,however,that a considerableamountof work(mostlyinturbulentQoIs)stillremainstoattainthatgoal.

Inthiswork,weinvestigatetheKármánvortexstreetgenerated by a flow around a rectangularcylinder by meansof the incom- pressibleNavier–Stokesequations.Thefollowingarethemainob- jectivesofthepresentstudy:

-Toshowthatstandard gPCissuitableforunsteadyproblems, wheretheQoIsareperiodic andthefrequencyremains de- terministic.

-To show that standard gPC failswhen investigating periodic QoIs that exhibit stochastic frequencies. We want to warn thatstandardigPCandNigPCcannot beusedtostudysuch cases,ormorecomplexcasesregardingturbulentQoIs.

-Toshow thatATI-igPCworks predominantlywell inperiodic flows that exhibit stochastic frequencies. This extends the work ofLe Maître etal. [35],who introduced ATI-igPCfor ordinarydifferentialequations.

-ToshowadvantagesofigPCwithrespecttoNigPCintermsof computationalcosts,whentheQoIsofinterestareunsteady fields.

Thepaperisstructured asfollows: InSection2the gPCseries expansion is put in a wider context, and it is shown how ATI- igPCcanbederivedfromageneralseriesexpansion.Section3ex- plainsthe fundamentalsoftheATImethod.InSection 4,the gPC expanded set of equations for the incompressible Navier–Stokes equationsisderivedusingtheATImethod.Section5givesabrief overviewofthestatisticalmeasuresusedinthiswork.Thenumer- icalresults,includingsomeinformationonthecomputationalcosts ofthesolver,arepresentedinSection6andconclusionsaredrawn inSection7.

2. Seriesexpansionsofrandomfields

Consider asuitableprobabilityspacedescribingallrelevantin- putuncertaintiesofthenonlinearspace-timedependentmodelof interest. The elementary events ofthe space are denoted

ω

∈, where ^{isthe samplespace.Torepresentandpropagate uncer-} taintythroughthemodel,considerthegeneralseriesexpansion

v ₍

x,

τ

^,

ω )

⁼

v ₍

x,

τ )

^+∞

i=1

i

( τ

^,

ω ) v

i

(

^x,

τ )

^{, (1)}

where

v

₍x,

τ

) ^{is the expected value function,} i(

τ

,

ω

) ^{are zero-} mean stochastic processes, and

v

i(^x,

τ

) ^{are linearly} independent fieldsinphysicalspace,with

{

i

}

,

{ v

i

}

^possiblyorthonormalbases inrandomandphysicalspace,respectively.Thescaledtime

τ

^isa

function ofphysical time that we may alsoallow to be stochas- tic, i.e.,

τ

=

τ

(^t,

ω

)^{. Eq. (1) is a} generalization of the classical Karhunen-Loèveexpansionwithtimedependentstochasticandde- terministic components. Severalmethods can be derived by sub- stituting the expansion (1) into the nonlinear model of interest.

Inallsubsequentexampleswewillassume thatthestochasticdi- mensionality can be restricted to n dominant dimensionsof un- certainty.The temporalevolution ofthebasis

{ v

i

}

^{isrestrictedto}

the instantaneous normal direction, i.e., any new contribution is orthogonal to the currentbasis. By projecting the governingPDE ontothebases

{

i

}

^and

{ v

i

}

^whilekeeping

τ

=tfixed,oneobtains the DynamicallyOrthogonal (DO) field equations[40,41].The DO method wasusedtosolve thestochastic Navier–Stokesequations inSapsisandLermusiaux[42].Employingalinearscalingof

τ

^with

respecttotimetinconnectionwiththeDualDOequations(where the stochastic insteadofthe deterministicbasis iskept orthonor- mal)ledtoimprovementintheperformance bymeansofsmaller effectivenumberoftermsin(1)[43].Byfixingthestochasticbasis withrespect totime andlettingi(

ω

)=i(

ξ

(

θ

))^{be orthogonal} polynomials withrespectto some choice ofrandom vector

ξ

(

θ

), theDOfieldequationsreducetothestochasticGalerkinprojection ofthegPC expansion.Controllingthedynamicsof

τ

(

ξ

,t)^withan ODEresultsintheATImethod.Next,wewilldescribegPCandATI insomemoredetail.

2.1. Generalizedpolynomialchaos

gPCoffers aframeworkbyspectralexpansionsinrandomvari- ables via orthogonal, stochastic polynomials. Assume that all in- put random variables and fields of relevance to the application of interest can be parameterized to sufficient accuracy by an n- dimensional random vector

ξ

=(

ξ

¹,. . .,

ξ

ⁿ), where for brevity of notation we have omitted the dependence on , i.e.,

ξ

=

ξ

(

ω

)^. The random variables are assumed independent with joint PDF w(

ξ

)=w₁(

ξ

1)...wn(

ξ

ⁿ)^.

Let

{

i(

ξ

)

}

^{beasetofnormalized}multivariatepolynomialsor- thogonalwithrespecttow,wherei=(ⁱ1,...,in)∈Nⁿ₀.Themulti- variatepolynomialsareproducts ofunivariateWiener-Askey poly- nomials,i.e.,i(

ξ

)=

ψ

i₁(

ξ

1). . .

ψ

in(

ξ

n)^.Theorthogonalityofpoly- nomialsimpliesthat

E_ξ

i

( ξ )

j

( ξ )

= n

k=1

δ

ikjk, ik,jk∈N0, (2)

where

δ

i j istheKroneckerdeltaandE_ξ(·)^denotestheexpectation operator,

E_ξ

(

^f

)

^{= f}

( ξ )

^w

( ξ )

^d

ξ

^{. (3)}

Anyfinite-variancerandomfieldx(

ξ

)^{canthenbeexpressedbyits} gPCexpansion,

x

( ξ )

=

i∈N0

xi

i

( ξ )

. (4)

Forpractical use,thegPC expansionmustbetruncatedtoafinite numberoftermsby restrictingtheorderof thepolynomials.The truncationschemeisdefinedbyselectingafiniteindexsetI⊂Rⁿ₀ givenby

I=

{

ⁱ∈Rⁿ₀:

ⁱ

q≤p

}

, (5) wheretheq-norm(orquasi-normifq<1)isdefinedby

ⁱ

q=

_n

j=1

i^q_j

1/q

. (6)

Settingq=1leadstothetotal-ordertruncationincludingallpow- ers up to p, which is used in this work. Choosing q<1 yields the hyperbolic index sets [11], favoring the influence of univari- ate polynomialsatthe expense of high-order mixedpolynomials including several random dimensions. Other options include the hyperbolic cross truncation strategy where compound nonlinear effects from multiple random variables are truncatedin favor of higher-ordernonlineareffectsinsingleandsmallnumbersofjoint variables [44].The totalnumber of terms (^P+1) ^{in an isotropic} total-orderexpansionofmaximumpolynomial order pandnun- certaininputvariablesisgivenby

P+1=

(

^n+p

)

^!

n!p! . (7)

Usingsingle-indexnotation,i.e.,mappingthevector-valuedmulti- index set to the non-negative integers, the resulting scalar in- dices correspond to different multi-dimensional basis functions.

ThetruncatedgPCexpansionthenreads x

( ξ )

^≈

P

i=0

x_i

i

( ξ )

^{. (8)}

Oncethe truncatedgPC coefficientsin (4)are calculated, theex- pected value E_ξ(^x) ^{and the} approximate variance V_ξ(^x) ^{of the} stochasticsolutionaregivenby

E_ξ

(

^x

)

=x₀, V_ξ

(

^x

)

=

P

i=1

x²_i.

3. Asynchronoustimeintegration

Generalized polynomial chaos fails to represent periodic QoIs thatarecharacterized byastochastic frequency.Toovercomethis problem, Le Maître etal. introduced asynchronous time integra- tion(ATI)[35].Wewillbrieflydescribe themethodinthefollow- ing.Foramoredetailedexplanation,thereaderisreferredto[35]. ConsiderasetofstochasticODEs,writtenas

d

dty

( ξ

,t

)

=F

(

^y

( ξ

,t

))

, (9)

wherey(

ξ

,t)^{isarandomvectorand}Fisadiscreteoperator.Now choosea certain realizationof

ξ

^asadeterministic referenceand callthissolutiony^(r)thatobeysEq. (9):

d

dty^(r)

(

^t

)

^=F(r)

(

^y(r)

(

^t

))

^{. (10)}

Eqs.(9)and(10)canbe seenassemi-discretizedPDEs,wherethe vectorywilllaterbe assignedtothegPCexpandedvelocitycom- ponents.Buildingonthis,Fcanbeinterpretedasanoperatorthat describestheconvection,diffusionandthepressuretermincluded intheNavier–Stokesequations.

ThefocusontheincompressibleNavier–Stokesequationsmakes ourworkdifferfrom[35],wheresmallsystemsofODEswerecon- sidered. Foralmost all

ξ

,Eq. (9)issatisfied bysolutions y(

ξ

,t),

Fig. 1. The behavior of the stochastic solution y (ξ, t) and the in-phase solution z (ξ, t) in the t- and τ-time space.

Fig. 2. Phase portrait of a realization of the stochastic solutions z (ξ, t) (orange) and the reference y ^(r)(blue), where y ^(r), z : R ⁿ→ R ². The vectors represent the input for = d (ξ, t)·F ^(r)(y ^(r)(t)). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

each associatedwith a particular oscillationfrequency f(

ξ

)^{. This} frequency is slightlyshifted from a referencefrequency f^(r) that correspondstoareferencesolutiony^(r).Wenowstretchandcom- press the solutions y(

ξ

,t) ^{so that all signals have the same fre-} quency f^(r). We call these uniformlyoscillating solutions z(

ξ

,t), which canbe generatedby introducinga stochastic timevariable

τ

(

ξ

,t)^{thatisafunctionofthephysicaltimet andoftherandom} realization

ξ

^{(i.e.,the time rescaling is differentfor each realiza-}

tion):

z

( ξ

^,t

)

^=y

( ξ

^,

τ ( ξ

^,t

))

^{. (11)}

Inordertogetaformulationfortheuniformlyoscillatingsolu- tions z(

ξ

,t) ^{onthe basisofthe operator}F,the chainruleisap- pliedtoEq. (11):

d

dtz

( ξ

^,t

)

⁼ _dt^dy

( ξ

^,t

)

t=τ (ξ,t) d

τ

dt =F

(

^y

( ξ

^,

τ ( ξ

^,t

))) τ

^{˙ =F}

(

^z

( ξ

^,t

)) τ

^˙.

(12) The termF(^z(

ξ

,t))

τ

^{˙ is a scaled version of} F(^z(

ξ

,t)), which guaranteesthatthesolutionsz(

ξ

,t)^areinphase.

Fig. 1illustrates theprocess of thetime transformation.From left to right: (1) The reference solution anda solution y(

ξ

,t) ^of a selectedrealizationof

ξ

^{are shown,which areout ofphase. (2)}

Usingatimescaling,itispossibletoalignbothsignalsbychoosing z(

ξ

,t)^{sothattheyareinphaseinthet-time space.(3)Sincethe} referencesolutionhasdoublethefrequencycomparedtoy(

ξ

,t)ⁱⁿ thisexample

τ

(^t)=2t.(4)Toretrievey(

ξ

,t),weplotz(

ξ

,t)^over

τ

(^t)^{.Thesignalsappeartobeinphase,butthe}

τ

^{-axisprogresses}

twiceasfast.Therefore,theoriginalfrequencyofy(

ξ

,t)^isrestored inthe

τ

^-timespace.

Thetimescalingcanbesummarizedas

z

( ξ

,t

)

=y

( ξ

,

τ ( ξ

,t

))

(push-forward),

y

( ξ

,t

)

=z

( ξ

,

τ

⁻¹

(

t,

ξ ))

(pull-backward), (13)

wherethepush-forwardoperationdescribestheprocessinwhich allsolutionsarebroughtintophase.Thepull-backwardtimescal- ingisusedtorecovertheoriginaloscillatorybehaviorofthesolu- tion.Thetasknowistofindawayofdescribing

τ

^{thatguarantees}

in-phasebehaviorbetweenthesolutionsz(

ξ

,t) ^{andthereference} solutiony^(r).

Thedistancebetweenthestochastic solutionandthereference solutioninthephaseportraitcanbeexpressedas

d

( ξ

^,t

)

^=z

( ξ

^,t

)

^−y(r)

(

^t

)

^{. (14)}

Whenthedistancevectorandthevelocity vectorofthereference trajectoryareorthogonal,thesolutionsoscillate inphase.The fol- lowingdotproductisthereforeameasureoftheirphaseshift:

=d

( ξ

^,t

)

·F^(r)

(

^y(r)

(

^t

))

. (15) ThisrelationisgraphicallyshowninFig.2,wherez(

ξ

,t)^andy(r) are two-dimensional vectors. When the trajectories are in phase (middle),then=0andtheclockspeeddoesnotneedtobead- justed.If<0(left),the localclockspeed mustbe increasedso thatthestochastictrajectorycancatchupwiththereferencesolu- tionandviceversa(right).

LeMaîtreetal.[35]proposedan ODEforthedeterminationof

τ

˙ ^{thatenforcesthe}=0criterion:

d

dt

τ

˙

( ξ

^,t

)

=−

α

⁰

τ

^˙

( ξ

^,t

)( ξ

^,t

)

+

α

¹

(

¹−

τ

˙

( ξ

^,t

))

. (16) In this expression

α

0 determines how fast the local clock speed respondstochangesin(

ξ

,t),

α

1 controlstheasymptoticbehav- ior of

τ

˙, and(¹−

τ

˙)^{ensures that theclock speed remains close} to1. Additionally,theseparametersare chosen ina waythat the fixedpoint

τ

˙^∗,i.e., thevalueof

τ

˙ ^{that makesthe}right-hand-side ofEq. (16)equalto zeroorequivalently _dt^d

τ

˙ =0,should behave asanattractorforalargerangeof^.Consequently,thedynamical adjustmentoftheclockspeedasillustratedinFig.2isassured.We summarizethesystemofequationsthatfollowtheATImethod:

d

dty^(r)

(

^t

)

=F^(r)

(

^y(r)

(

^t

))

, (17)

d

dtz

( ξ

^,t

)

=

τ

˙^F

(

^z

( ξ

^,t

))

, (18)

d

dt

τ

˙

( ξ

^,t

)

=−

α

⁰

τ

^˙

( ξ

^,t

)( ξ

^,t

)

+

α

¹

(

¹−

τ

˙

( ξ

^,t

))

, (19) withinitialconditions:

y^(r)

(

^t=0

)

^{=y0, (20)}

z

( ξ

,t=0

)

=y⁰, (21)

τ

˙

( ξ

,t=0

)

=1. (22)

Inallnumericalexperimentsweset

α

0=8e−5and

α

1=10. Notethat noexactknowledgeofthefunctionaldependenceof t and

τ

^{isrequiredwhenfollowingtheATImethodasthedynam-}

icsof

τ

˙ ^{iscontrolledviaanODE.}Accordingly,weconsidertheATI approachtobemoregeneralifcomparedtoothermethods,asthe onefollowedinMusharbashandNobile[43],whereanexplicital- gebraicrelation betweent and

τ

^{isrequired.Next,we willapply}

theATIframeworktotheincompressibleNavier–Stokesequations.

4. ApplicationtotheincompressibleNavier–Stokesequations 4.1. Stochasticflowformulation

The Navier–Stokes equations are obtained from conservation lawsofmassandmomentum.Incaseofincompressibleflows,the densityisassumedtobeconstantinspaceandtime.Conservation ofmassandmomentumreduceto:

∇

·u=0, (23)

∂

^u

∂

^{t +}

⁽

^u·

∇ )

^u=−

ρ

¹

∇

^p+

ν∇

^{2u, (24)}

withsuitableboundaryconditionstobedescribedlater,andwhere u, p,

ρ

^and

ν

^{represent the velocity, pressure, density, and the}

kinematic viscosity,respectively. In thisstudywe assume that u, p and

ν

^are stochastic. Therefore, these quantities are expanded asgPCseries(8)andinsertedintothegoverningequations,which yields

P

i=0

i

∇

^·ui=0, (25)

P

i=0

i

∂

^ui

∂

^{t +}

P

i=0

P

j=0

i

j

(

^ui·

∇ )

^uj=− P

i=0

i

1

ρ ∇

^pi

+ P

i=0

P

j=0

i

j

ν

i

∇

^2uj (26)

The system Eqs. (25) and (26) is projected onto the stochastic basis

{

k

}

, which is composed by orthonormal polynomials(see Eq. (2)).ThisresultsinP+1equationsformassconservationand 3(^P+1)^{equationsformomentum}conservation:

∇

·uk=0, (27)

∂

^uk

∂

^{t +}

P

i=0

P

j=0

E_ξ

i

j

k

(

^ui·

∇ )

^uj=−1

ρ ∇

^pk

+ P

i=0

P j=0

E_ξ

i

j

k

ν

i

∇

^2uj, (28)

where k=0,...,P. The expectations of basis function triples E_ξ

ijk

can becomputedexactlyfororthogonalpolynomials usingGaussianquadraturerules.Alternatively,closed-formexpres- sionsforsomefamiliesofexpectationsoftripleproductsbasedon polynomialrecurrencerelationsareavailableintheliterature[45]. Thenumericalmethodusedtosolvethestochastic,incompressible Navier–StokesequationsispresentedinAppendixA.

4.2. AsynchronoustimeintegrationofincompressibleNavier–Stokes

WedefineFasadiscreteoperatorassociatedwiththestochas- ticmomentumequation

F(^u(

ξ

,t))=

−

u(

ξ

,t)·

∇

u(

ξ

,t)− 1

ρ ∇

^p(

ξ

,t)+

ν

(

ξ

)

∇

^2u(

ξ

,t)

, (29) wherethebrackets‘[]’symbolizethespatialdiscretizationopera- tion.Inasimilarway,F^(r)isdefinedasadiscreteoperatorassoci- atedwiththedeterministicmomentumequation

F^(r)

(

^u(r)

(

^t

))

=

−

u^(r)·

∇

u^(r)−1

ρ ∇

^p(r)+

ν

^(r)

∇

^2u(r)

. (30)

Eq. (18)isrewrittenresultingin d

dtu

( ξ

,t

)

=

τ

˙

( ξ )

F

(

^u

( ξ

,t

))

. (31) Applying gPC expansion to Eq. (31) andprojecting the resulting expressiononthestochasticpolynomialbasisleadsto

d

dtu_k=F_k^gPC, fork=0,...,P, (32) where

τ

˙(

ξ

)^hasbeenincorporatedintheoperatorF_k^gPC:

F_k^gPC =

− P

i=0

P

j=0

P

m=0

E_ξ

i

j

m

k

τ

˙i

(

^uj·

∇ )

^um (33)

− P

i=0

P

j=0

E_ξ

i

j

k

τ

˙i

1

ρ ∇

^pj

+

P i=0

P

j=0

P

m=0

E_ξ

i

j

m

k

τ

˙ⁱ

ν

^j

∇

^2um

. (33)

Inasimilarway,Eq. (19)isexpandedleadingtoasystemofODEs forthecoefficient

τ

˙k,asshowninAppendix B.Notethat theex- pectationsofthequadrupleproductE_ξ

ij^mk

arecalculated usingGaussianquadraturerules.

5. Statisticalquantitiesofinterest

Numerical simulationof unsteady periodicflow permits ade- tailedinvestigationofthesystemasaresultofthelargeamountof information:temporalsignalsforeachflowvariableateverynode ofthecomputationaldomain.Aside fromthe instantaneousfields obtainedbythesimulation,twostatisticalquantitiesaregenerally ofmajorinterest instudiesfeaturing instability,vibrationsand/or limitcycles:thetemporalexpectedvalueEt(^y(^t)),andthetempo- ral varianceVt(^y(^t))^.Aschematicrepresentationofasignal y(^t), associatedwithagivenpointinphysicalspaceofthefieldy(^t),is giveninFig.3(graycurve).

Forillustrationpurposes,wediscretize thefield y(

ξ

,t) ^inran- dom space by N differentrealizations, as depictedin Fig. 3 (red curves), along with the PDFs of Et

y(

ξ

,t)

and Vt

y(

ξ

,t)

. To quantifythebehaviorofagivenperiodic,stochasticfieldy(

ξ

,t)^it isconvenienttodefinefourmeasures:

EE≡E_ξ

(

^Et

(

^y

) )

^=Et

E_ξ

(

^y

)

=y₀, (34)

Fig. 3. Illustration of signal y (ξ, t)for Nrealizations. The signal is characterized by four scalars EE , VE , EV and VV .

EV≡E_ξ

(

Vt

(

^y

) )

= P

i=0

Ety²_i−

(

Ety_i

)

²= P

i=0

y²_i−y²_i, (35)

VE≡V_ξ

(

^Et

(

^y

) )

⁼ P

i=0

y²_i, (36)

VV≡V_ξ

(

^Vt

(

^y

) )

⁼

2P

m=1

P i=0

P

j=0

y_iy_j−y_iy_j

E_ξ

(

i

j

k

)

2

= P

i=0

P

j=0

y²_iy²_j−y²_i y²_j =Vt

V_ξ

(

^y

)

, (37)

wherethebarnotationyhasbeentemporarilyintroducedtofacil- itatetodistinguishbetweentemporalandstochasticmeanvalues.

These expressions are the expected value in random space of thetemporalmean(Eq. (34))andvariance(Eq. (35)),respectively, as well asthe variance ofthe sametemporal meanEq. (36)and variancesEq.(37).Thelattertwocanbeusedtodefineerrorbars forthetwo formerexpectations.In theexampleofFig.3,we ob- serve in the histogramplot atthe bottom left that thetemporal mean value of each signal varies considerably along realizations

(9<Et(^y(

ξ

,t))<11). We alsoobservevariabilityin thetemporal variance(3<Vt(^y(

ξ

,t))<5).

The quantitiesEE, EV, VEandVV are usefulwhen assessing howthetemporalmeanandtemporalvarianceofaperiodicflowis affectedbyinputuncertainties.Ifonlythosemeasuresareofinter- est,NigPCshouldbethemethodconsidered,asitiseasytoimple- mentandperformssatisfactorily. Notethat theshapeoflimit cy- clesignals,describedbykcharacteristicamplitudes(period-klimit cycles)andfrequencies, arenot capturedbyEE,EV,VEandVV. Suchinformationcanbe,nevertheless,ofhighinterestincases,e.g.

whenthefrequencycontentofthesignalsandmaximaofthecor- responding crests are required. Accordingly, in the study of flow limit cycles, it isgenerally appropriate to analyze the signalsdi- rectly, e.g. by buildingphase portraitsof correlated quantities. It might also be of interest to evaluate the probability distribution thatcharacterizes aperiodicflowvariable.Forsuchrequirements, thequantitiesEE,EV,VEandVVarenotsuited.Inthiswork,we show that an adequate surrogate model for y(

ξ

,t) ^{is feasible by} meansofigPC(standardorATI)foroscillatinglaminarflows.

6. Numericalresults

Theimplementationofthemethod(Section 4)inourin-house solver is validated in the supplementary material using a back-

Table 1

Numerical values of characteristic length scales.

H/h c B/h c L u/h c L d/h c b c/h c Grid-size

8 16 3.5 10.5 2 80 × 160

Table 2

Coordinates of the probes.

P 1 P 2 P 3 P 4 P 5 P 6 x/h c 3.5 5.3 7.1 3.5 5.3 7.1

y/h c 0 0 0 2.5 2.5 2.5

ward facing step. The resulting tool is named SUN-S (Stochastic Unsteady Navier–Stokes) andcan take into account uncertainties inthefluidviscosityandinboundaryconditionsforvelocityuand pressure p. It wasdeveloped for MATLAB 2017b based on finite volumes.Thenumericaltestcasesareone-dimensionalinstochas- tic space. The extension to multiple dimensions is conceptually straightforward butleadsto increasedcomputationalcost, asdis- cussedinmoredetailinSection6.4.

6.1. Flowaroundarectangularcylinder

We investigate the Kármán vortex street generated by a rect- angularcylinder. Thevortexsheddingdownstreamthe interfering body isa well-knownobservationinfluidmechanics.Vorticesare formed attheback ofthecylinderandperiodicallyseparate from the upperandlower sideofthebody.Althoughthephenomenon isassociatedwithKármán’sname,thefirstexperimental observa- tions werereportedby Mallock[46] andBénard[47]in1907and 1908.The maincontribution ofKármán wasthe stabilityanalysis ofthevortexformationandtheexplanationofthemechanismsof wakedragin1911[48].TheKármánvortexstreetisaninteresting problem because we observeunsteady quantities while the flow patterncanstill beassignedtothelaminarregime.Therefore,we cansolvetheoriginalNavier–Stokesequationswithoutfurthertur- bulencemodeling.

ThegeometryandtherelevantdimensionsareshowninFig.4. Table 1presentstheassociatednumericalvalues.Thecoordinates of theprobes are listed in Table2. Theinlet velocity describesa flatprofilewiththebulkvelocityu_in.Thepressuregradientatthe inletiszero.Onthetopandbottomboundariesofthedomainwe set a zero gradient conditionforboth velocity and pressure.The

Fig. 4. Geometry, length scales and location of the probes of a 2D flow with a rect- angular bluff body.

outflowboundaryconditionsareasfollows:

∂

^u

∂

^{x =0,}

∂ v

∂

^{y=0 and p=1atm. (38)}

The rectangular cylinderis placed on the line ofsymmetry with respecttotheheightHanditssurfacesaremodeledassolid,non- slipwalls. TheReynolds numberforthisspecificcaseis givenby Re=u_in·hc

ν

^{, (39)}

where

ν

^{isthekinematicviscosityandhcistheheightofthecylin-}

der.

We study two cases and consider an artificial fluid. Case A refers to variations of kinematic viscosity, and case B refers to variations of inlet velocity. Both cases correspond to the same range in Re∈ [95,105]. Fig. 5 (left) shows the characteristic fre- quencyofthewakeoscillationasafunctionoftheStrouhalnum- berSr=f·hc/u_in.Probe5 ^{isusedtomeasurethe}oscillations.Re- sults,whichareinagreementwith[49–52]arescatteredalongone singlelinearoundSr≈0.135.

Fig.5(right)showstheoscillationfrequencyforthesamecases.

Changes in Re lead to stronger frequency variations in Case B.

CaseA exhibitsalmost no changein frequencyin therange Re∈ [70,130].

Fig. 5. Strouhal number as a function of Reynolds number for b c/h c= 2 (left). Results are shown in the range Re ∈ [100 , 250] for comparison with other studies (Okajima exp. [49] , Okajima num. [50] , Sohankar num. [51] and Islam num. [52] ). Frequency of the oscillations in the wake of the cylinder as a function of Re (right). Results are shown in the range Re ∈ [70 , 130] .

Fig. 6. Fields EE and VE for case A.

Fig. 7. Fields EV and VV for case A.

6.2. CaseA:uncertainviscosity

Inthissection,weconsiderthekinematicviscositytobeuncer- tain:

ν ( ξ )

=

ν

⁰⁺

ν

¹

ξ

^{, (40)}

where

ξ

∼U[−1,1] isa uniformdistribution and

ν

0=1 m² s⁻¹, and

ν

1=0.052m² s⁻¹.Theinletvelocityisdeterministicandcon- stantwithavalue100ms⁻¹.Consequently,Re∈[95,105].Wecon- sidernormalizedLegendrepolynomialsforgPC.

The simulationis carriedout usingthestandard intrusivegPC technique(notATI)andtheresultsarevalidatedinthreedifferent ways.

i) The fields associated with EE, VE, EV VV, introduced in Section5,arecomparedtothefieldsobtainedfromtheMonte Carlosimulation.

ii)Three realizations of

ξ

^{are selected and}corresponding values ofvelocity are computed fromthegPC model(8).Results are compared ina phase portrait to velocity time series obtained fromthreedifferentdeterministicsimulations.

iii) AMonteCarlosimulationisperformedwith500samples.The correspondingPDFattwolocationsandthreetimeinstantsare comparedtotheonesobtainedfromthegPCsimulation.

Figs. 6 and 7 show the fields of EE, VE, EV and VV calcu- latedfrom established limit cycles.We observe that EE captures the recirculationzone (dark blue) dowsntream of the bluff-body.

Thefield ofEV(Fig.7) suggeststhat thisrecirculationzonedoes not exhibit oscillations in the region immediately downstrem of the rectangular cylinder, as indicated by the dark blue zone of EV. As expected, a largeregion ofconstant inlet velocity (no os- cillations)isencounteredupstreamofthebluff-body.ThefieldEE alsoshowsregionswherethetemporalmeanvelocityishigh(dark red).Notethatonlyaportionofthisregionexhibitsoscillationsof smallamplitude. The field EV showszoneswhere theamplitude oftheoscillationsishigh,asobservedinthetwosymmetricdark redzones.Wealsoobservetwoelongated,symmetricareasdown- stream of the bluff-body. These regions characterize the increas- ing/decreasingamplitudeofthe wake downstreamof therectan- gularcylinder.

ThefieldVEsuggeststhatthetemporalmeanvalueoftheos- cillations(Et) isnotstronglyinfluencedbyuncertaintiesof

ν

^over

largeregions.Theovalshapeareadepictedinredistheonlypre- dominantregion,wherethetemporalmeanofvelocityoscillations isaffectedbyuncertaintiesin

ν

^.ThefieldVVshowstheinfluence of

ν

uncertaintiesontheamplitudeofoscillationsofvelocity.The variabilityisconfinedwithinthered-coloredregionofEV,which

Fig. 8. Limit cycle oscillations for probe 2 (left, period-two limit cycles) and probe 5 (right, period-one limit cycles) of the reference and the stochastic simulation with P = 4 for different realizations of ξ.

Fig. 9. Limit cycle oscillations for probe 1 (left) and probe 6 (right) of the reference and the corresponding stochastic simulation with P = 4 for different realizations of ξ. Note that other two probes are here considered to show additional phase portraits. The quality of agreement, as shown in this Figure, is the same among all probes considered in this study (not shown).

implies thatlarge amplitudeoscillations are affected by

ν

^uncer-

taintiesinastrongermannerthanthelowamplitudecounterparts.

From Figs. 6 and7, it is also observed that igPC agrees very well withMonte-Carlo simulations. Althoughnot shownin here, wehighlightthatanexpansionofP=1issufficientforaccurately capturingthefieldsofEE,VE,EVandVVincaseA.

Fig.8depictsaphaseportraitofthetwovelocitycomponentsu and

v

.Thecontrolpointscorrespondtoprobe2 ^{andprobe5 .The} computation isdone withan order ofP=4. Convergence results areachievedandnofurtherimprovementintheresultsisseenfor P>4. ThegPC trajectories,both period-two andperiod-onelimit cycles,agreeverywell withthe deterministiccounterpart,asalso showninFig.9.

We comparethestochastic signalsofthehorizontalvelocityu in Fig. 10,which are extracted atprobe 3 ^{and5 .Whereasthe} signal measured at probe 3 ^{exhibits two} characteristic frequen- cies (aperiod-two limit cycle), the signal of probe 5 ^{is related} toonlyonecharacteristicfrequency(aperiod-onelimitcycle).The period-twolimitcycleenclosestheinformationofthevortexshed- ding fromtheupperandlowerside ofthebluff-body.Duetothe modestcomputationalcostwegenerated4000samplesofgPC.By comparing the differences betweenestimators, it wasconsidered sufficient to use 500 samples for the much more expensive MC simulations. It is observedthat igPCnot only capturesthe shape of the two type ofsignals, butalsothe crests.Only smalldiffer- ences are encountered inthe spreadof therealizations. Thiscan be attributedto a very smallfrequencyshift,that is capturedby MCbutnotbyigPC.

Thesesignalsare thebasis forthePDFsshown inFig.11.The pointsintimechosen forevaluating thePDFarehighlightedwith bluelinesinFig.10andcorrespondto[0,0.125,0.75]forprobe3 and[0,0.25,0.5]forprobe5 ^{.TheshapeofthePDFsisfairlywell} recoveredby igPC,wherethemodeispredictedapproximatelyin themiddleoftherespectivePDF.ThesupportofthePDFsdiffers, whichmayarise fromsmallfrequencydependenciesofthedeter- ministicsolution.Toconclude:standardigPCcanaccuratelypredict thecharacteristicfrequencyandtheamplitudesforallrealizations, whereasthesupportofthePDFsisnotperfectlyreproduced.

6.3. CaseB:uncertaininletvelocity

Considerthecaseofuncertaininletvelocity

u

( ξ )

^=u0+u₁

ξ

^,

ξ

^{∈U[−1,1], (41)}

withu₀=100ms⁻¹,u₁=5ms⁻¹whiletheviscosity

ν

=1m²s⁻¹ remains constant and deterministic. The range ofRe investigated remains unchanged (Re∈[95,105]) with respect to case A. For the simulations using igPC-ATI, we set

α

0=8e−5 and

α

1=10. A suitable set of parameters can be found by evaluating (Eq. (15))fromprevious simulations andadjusting

α

0 and

α

1 so thatEq. (19)hasafixedpointforalargerangeof^.

WeinvestigatethefieldsEE,VE,EVandVVinFigs.12and13. Itisobservedthat EEandEVremainsunchangedwithrespectto caseA.Itmeansthat theoriginofuncertainties(eitherfrom

ν

^or

uinlet) do not play arole in thetemporal meanvalue andvari-

Fig. 10. Signals of the normalized u -velocity for probe 3 and 5 and different realizations of ξcompared to Monte Carlo. The results for igPC correspond to a simulation with P = 4 .

Fig. 11. PDFs of the normalized u -velocity for probe 3 and 5 and different snapshots of one period T refof the limit cycle for P = 4 .

anceofthevelocityoscillations.Incontrast,thefieldofVEisvery differentwithrespecttothecounterpartsofcaseA,whichimplies thatthepropagationofuncertaintiesis,therefore,verydifferentin case Aand incaseB. We observe inVEof caseB that the dark blue regionupstream ofthebluff-body doesnot exist(as incase A), anditisnow representedby alight blue-sky color.Thisis in agreementwiththefactthattheinletvelocityisconstantbutun- certain. We also observe that the regions of highvelocity inEE (dark red)are alsotheregions moresusceptibleforuncertainties, asillustratedinVE(darkred). ThefieldofVV,contrarytoVE,is very similarincaseAandincaseB,yetthemagnitudesarevery different:thelargestvaluesofVVincaseAarethreetimessmaller than in caseB. Itmeans that the amplitudeof thevelocity limit

cyclesaremoresusceptibletouncertaintiesintheinletvelocityu thantouncertaintiesin

ν

^.

Figs.12and13alsoshowthe fieldsEE,VE,EVandVVcom- putedbyigPCandigPC-ATI.Bothtechniquesareabletoaccurately capture the fields EE, EV. We also observe that the uncertainty in Et and Vt, which is represented by VE and VV, is not well captured by igPC. This is due to the fact that case B exhibits a stochasticfrequency.Accordingly,thephaseshiftbetweenrealiza- tions cannotbe modelledby thestandard igPCtechniqueforlow expansionorders(suchasP=4choseninthiswork).Incontrast, igPC-ATI accurately captures the fields VE and VV. Even an ex- pansion withP=1 is sufficient to reproduce all fields well (not shown).Sincethehighermodesof

τ

^{andudecaytoasimilarex-}

Fig. 12. Fields EE and VE for case B.

Fig. 13. Fields EV and VV for case B.

tent,itwasconsideredsufficienttotruncatetheexpansionofthose quantities to the same polynomial order, whichis also P=4 for igPC-ATI.

Wecomparethephaseportraitsobtainedbythreedeterminis- tic simulations withtheonesobtained byigPC andigPC-ATI. We consider thesamelocationsasincaseA.Fig.14showsa period-

two limit cycle, whose amplitude slightly increases for increas- ing valuesofRe.Thisisin correspondencewithcaseA, although in that casethe trajectoriesare closer to each other (see Fig. 8).

Fig.14alsoshowsthatigPC-standardfailswhenreproducingthese trajectories,whereasigPC-ATIdoesagoodjob.Thephaseportraits illustratedinFig.15,takenattwodifferentprobesfordifferentre-

Fig. 14. Limit cycle oscillations for an uncertain inlet velocity for probe 2 (left) and probe 5 (right) with P = 4 for different realizations of ξ.

Fig. 15. Limit cycle oscillations for an uncertain inlet velocity for probe 1 (left) and probe 6 (right) of the reference and the corresponding stochastic simulation with P = 4 for different realizations of ξ.

alizations of

ξ

, show a very good agreement with the reference simulation.

When looking atthe signals of the horizontal velocity u that were measured for probe 3 and 5 , the situation is similar to whatwasalreadyindicatedinthephaseportraits.InFig.16,where all realizationsor samplesare shown, we can observethat stan- dardigPCpredictsspuriousoscillationsatprobe3 ^andcannotre- producethecorrectamplitudes.Inaddition,anotherfactbecomes clearatthispoint:standardigPCisnotabletoreconstructamono- tonicfrequencyshiftbetweentherealizations.Thesignalofstan- dardigPCatprobe5 ^{exhibitsonlyafew}predominantfrequencies.

These frequencies can be detected followingthe red lines,which correspondto

ξ

∈[−1,−0.9]and

ξ

∈[0.9,1].Contrarytostandard igPC,eachrealizationcreatedbyigPC-ATIcanbeassigneditsown timescalebasedontheinformationavailablefrom

τ

(

ξ

,t)^.Thefre- quencyshiftpredictedbyMonteCarlocanthereforebereproduced wellbyasynchronoustimeintegration.

Based on the signals ofFig. 16, PDFs are created for MC and igPC-ATI (Fig. 17). In contrast to case A, in which all realiza- tionshavethesamefrequency,therealizationsareslightlyshifted against each other for case B. Therefore, the PDFs at the corre- spondingsnapshothavealargersupport,whatisadequatelyrepre- sentedbyigPC-ATI.Insummary,igPC-ATIcanreproducethephase shiftandpredictstheamplitudesofthelimitcycleverywell,while standardigPCfails.

6.4. Computationalcosts

Aspartofthedevelopmentofthiswork,simulationswerealso carried out withNigPC (standardandwithtime warping aspro-

posed in Mai and Sudret [37]), by means of which the stochas- ticquantities could be reproduced atgiventimes andspatial lo- cations. The quality of the results iscomparable to that of stan- dardigPCandigPC-ATI.Acomputationalcostcomparisonbetween NigPC (standardandwith time warping) andigPC (standardand withATI)ispronetosub-optimalimplementationsandothersoft- wareconstraints,suchasmemoryallocation.Weonlycommenton conceptualaspects that mayreveal advantagesanddisadvantages oftherespectivemethods.

The igPC method (standard or with ATI) may be well suited forthestudyofunsteady fields.Inthecasesinvestigated,weob- serve that the limit cycles are well resolved in time and space.

Once the gPC coefficients – which are actually unsteady, spatial fields – are obtained, it is simple to carry out thorough uncer- tainty quantification,wherethe evolution ofuncertaintiescanbe accuratelytracked intimeandspace.Weobservedthat theigPC- ATImethod ismore expensivethan the igPCcounterpart, mainly duetothecomputational costinvolvedinsolving themomentum equation. The latter is associated with a four-dimensional tensor (E_ξ(ijkm)^{)inigPC-ATIandwitha}three-dimensionaltensor (E_ξ(ijk)^)inigPC.

Ifonlytheconvective partP i=0

P j=0E_ξ

ijk

(^ui·

∇

)^uj of theexpanded Navier–Stokesequations isconsidered asan exam- ple, we observe that the double sum above consists of asmany termsasE_ξ

ijk

containsnon-zeroentriesforafixedk.Ifwe assume that the total computational cost for evaluating the mo- mentum equation is dominated by the number of terms in the summations,we concludethat thenumberofnon-zeroentriesof the aforementioned tensors is an indicator of the computational timeneededforsolvingthemomentumequation.