Radiceetal. RESEARCH Implicit large eddy simulations of anisotropic weakly compressible turbulence with application to core-collapse supernovae 5 1 David Radice, Sean M Couch and Christian D Ott 0 2 g u Abstract A Intheimplicitlargeeddysimulation(ILES)paradigm,thedissipativenatureofhigh-resolution 0 shock-capturingschemesisexploitedtoprovideanimplicitmodelofturbulence.TheILESapproachhas 2 beenappliedtodifferentcontexts,withvaryingdegreesofsuccess.Itisthede-factostandardinmany astrophysicalsimulationsandinparticularinstudiesofcore-collapsesupernovae(CCSN).Recent3D ] E simulationssuggestthatturbulencemightplayacrucialroleincore-collapsesupernovaexplosions, H howeverthefidelitywithwhichturbulenceissimulatedinthesestudiesisunclear.Especiallyconsidering thattheaccuracyofILESfortheregimeofinterestinCCSN,weaklycompressibleandstrongly . h anisotropic,hasnotbeensystematicallyassessedbefore.Anisotropy,inparticular,couldimpactthe p dissipativepropertiesoftheflowandenhancetheturbulentpressureintheradialdirection,favouringthe - o explosion.InthispaperweassesstheaccuracyofILESusingnumericalmethodsmostcommonly r employedincomputationalastrophysicsbymeansofanumberoflocalsimulationsofdriven,weakly t s compressible,anisotropicturbulence.Oursimulationsemployseveraldifferentmethodsandspanawide a [ rangeofresolutions.Wereportadetailedanalysisofthewayinwhichtheturbulentcascadeis influencedbythenumerics.OurresultssuggestthatanisotropyandcompressibilityinCCSNturbulence 2 havelittleeffectontheturbulentkineticenergyspectrumandaKolmogorovk−5/3 scalingisobtainedin v theinertialrange.Wefindthat,ontheonehand,thekineticenergydissipationrateatlargescalesis 9 6 correctlycapturedevenatlowresolutions,suggestingthatveryhigh“effectiveReynoldsnumber”canbe 1 achievedatthelargestscalesofthesimulation.Ontheotherhand,thedynamicsatintermediatescales 3 appearstobecompletelydominatedbytheso-calledbottleneckeffect,i.e.,thepileupofkineticenergy 0 closetothedissipationrangeduetothepartialsuppressionoftheenergycascadebynumerical . 1 viscosity.Aninertialrangeisnotrecovereduntilthepointwherehighresolution∼5123,whichwouldbe 0 difficulttorealizeinglobalsimulations,isreached.WediscusstheconsequencesforCCSNsimulations. 5 1 Keywords: turbulence;methods:numerical;supernovae : v i X 1 Introduction effectswherealloftheknownforces:gravity,electromag- r a Despitedecadesofstudiesandcompellingevidencethata netism, weak and strong interactions, are important. The significantfraction[1]ofstarswithinitialmassesinexcess task of modeling these systems is made particularly chal- of ∼8 solar masses explode as core-collapse supernovae lenging by the fact that the generation of the asymptotic (CCSN)attheendoftheirevolution,theexactdetailsofthe explosionenergies,althoughenormous(∼1044J),requires explosionmechanismarestilluncertain[2,3,4,5].Current arathersubtle,percent-levelimbalancebetweennon-linear state-of-theart3Dsimulationseitherfailtoexplodeorhave processesovermanydynamicaltimes. explosion energies that fall short of the observed energies Theflowofplasmainthecoreofastargoingsupernova by factors of a few for most of the progenitor mass range is known to be unstable to convection [7, 8, 9, 10] and/or [6,4,5]. toanotherlargescaleinstabilityknownasstandingaccre- tion shock instability [11, 12]. In any case, given the very Thedynamicsatthecenterofastarundergoingcorecol- large Reynolds numbers, as large as ∼ 1017 in the region lapse is shaped by a delicate balance between competing of interest [13] (the so-called gain region, where neutrino Fulllistofauthorinformationisavailableattheendofthearticle heating dominates over neutrino cooling), it is expected Radiceetal. Page2of17 that the resulting flow will be fully turbulent. It has been The aim of this work is to fill the gap between existing suggested [14, 15] recently that turbulence and, in partic- theoreticalstudiesandtheparticularapplicationsofourin- ular, turbulent pressure could tip the balance of the forces terest.Tothisendweuseapubliclyavailablecode,FLASH in favor of explosion. In this respect, anisotropy is of key [33, 34, 35], which is widely used in the computational importance, because it results in an effective radial pres- astrophysics community, and perform a series of simula- sure support with adiabatic index γ = 2, much larger tions of turbulence in a regime relevant for core-collapse turb thanthatofthermal(radiation)pressure(γ (cid:39) 4/3).This supernovae: driven at large scale, with large anisotropies th means that turbulent kinetic energy is a much more valu- and mildly compressible. We use five different numerical ablesourceofradialpressuresupportthanthermalenergy setupsand, foreach,several resolutionsinthe rangefrom (seeAppendixA). 643 to 5123 in a periodic domain. We study in detail the All of the current numerical simulations employ the wayinwhichtheenergycascadeacrossdifferentscalesis implicit large eddy simulation (ILES) paradigm [16, 17] represented by our ILES and we discuss the use of local (alsoknownasmonotoneintegratedLES(MILES))ofex- orlowerdimensionaldiagnosticsthatcanbeusedtoassess ploitingthedissipativenatureofhighresolutionshockcap- the quality of a global simulation in a complex geometry turing (HRSC) methods as an implicit turbulence model. where3Dspectraarenotreadilyavailable. However, the combination of the use of rather dissipative The rest of this paper is organized as follows. First, in schemes and the relatively low spatial resolution that can Section 2, we discuss the exact setup of our simulations andthediagnosticquantitiesusedinouranalysis.Then,in be achieved in global simulations is such that the fidelity Section 3, we discuss the basic characteristics of the flow withwhichturbulenceiscapturedisquestionable[13]. realizedinoursimulations.InSection4,wepresentade- TobeusefulinthecontextofCCSNsimulations,anILES tailed analysis of the way in which the energy cascade is should,attheveryleast,accountfortherightrateofdecay captured by the different schemes at different scales. In of the kinetic energy at the largest scales while avoiding particular, we quantify the accuracy with which different unphysical pile up of energy at smaller scales. Unfortu- methods capture the decay rate of energy from the largest nately, all of the current simulations seem to be strongly scales and the way in which energy is distributed across dominated by the so-called bottleneck effect [13], which scales. We discuss the role of anisotropies in the context corresponds to an inefficient energy transfer across inter- ofthe4/5−law,afundamentalexactrelationforisotropic mediatescalesduetotheviscoussuppressionofnon-linear andincompressibleturbulencerelatingthestatisticsofve- interactionwithsmallerscales[18,19,20,21,22].Current locityfluctuationswiththeenergydissipationrate(seeSec- globalsimulationsachieveresolutions,intheturbulentre- tion2.3),inSection5.Weexploretheuseofthe2D,trans- gion,comparabletothoseof303−703 latticesinperiodic verse,energyspectrumasadiagnosticfor3Dsimulations domains[23,15,13].Attheseresolutions,almostallofthe in Section 6. Finally, we present a brief summary of our dynamical range of the simulations can be expected to be mainfindings,aswellasadiscussionoftheirimplications directly affected by numerical viscosity [24]. The fidelity forCCSNsimulationsinSection7.AppendixA,contains withwhichturbulenceiscapturedinthesesimulationswill somesupplementalbackgroundmaterialontheroleoftur- thendependonthedegreewithwhichthenumericaltrun- bulenceintheexplosionmechanismofCCSN. cationerrorapproximatesanLESclosure. In this respect, it has been shown by [25] and [26] that 2 Methods many HRSC methods can be too dissipative to yield a 2.1 Numericalmethods faithfuldescriptionofturbulenceatlowresolutions.These We consider a compressible fluid with a prescribed accel- studies, however, considered a different regime, decaying eration,a,inaunit-boxwithperiodicboundaryconditions. isotropicturbulence,whileturbulenceinacore-collapsesu- The code that we employ for these simulations, FLASH, pernova,aswellasinmanyotherastrophysicalsettings,is solvesthegas-dynamicsequationsinconservationform.In oftenstronglyanisotropic[27,14,15]asrotationalinvari- particularweevolvethecontinuityequation ance is broken by gravity. [25] and [26] also considered different numerical schemes with respect to those used in ∂ ρ+∇·(ρv)=0 (1) t supernovasimulations.Bothoftheseaspectscan,inprinci- ple,beimportant.Firstofall,stronganisotropiescouldpo- andthemomentumequation tentially influence the turbulence dynamics at the level of the energy cascade and of the dissipation [28]. Secondly, ∂ (ρv)+∇·(ρv⊗v+pI)=ρa. (2) t some of the schemes used in computational astrophysics, suchasthepiecewiseparabolicmethod(PPM)[29]aswell These equations are closed with a simple isentropic equa- as some of the MUSCL [30] schemes, have been shown, tionofstate, differently from some of the methods considered by [25] and[26],tobewellsuitedforILES[31,32]. p=ρ4/3, (3) Radiceetal. Page3of17 thatcanbeconsideredasaroughdescriptionofagasdom- 2.2 Energytransferequations inatedbyradiationpressure.Sincetheequationofstateen- Inordertostudythecascadeofthespecifickineticenergy suresanadiabaticevolutionwedonotneedtosolvetheen- (whichwewillrefertosimplyas“kineticenergy”or“en- ergy equation as equations (1), (2) and (3) suffice to fully ergy”inthefollowing),|v|2/2,wewillconsideranenergy describetheflow. budget equation across different scales, analogous to the Equations (1) & (2) are solved using the directionally- one commonly employed in the study of incompressible, unsplit hydrodynamics solver of the open-source FLASH isotropic turbulence, e.g., [50]. In particular, we consider simulation framework. FLASH implements the corner themomentumequation(2)innon-conservationform, transport upwind method [36] for fully directionally- unsplit evolution of the Euler equations [37, 38]. FLASH ∂ v+(v·∇)v=−V∇p+a, (5) t includesseveraloptionsfortheorderofspatialreconstruc- tion [35], including 2nd-order TVD [30], 3rd-order PPM whereV =1/ρisthespecificvolumeofthegas. [29],and5th-orderWENOZ[39].Fluxesarecomputedat Wecanuseequation(5)toderiveanevolutionequation 2nd orderaccuracyusingoneofanumberofapproximate fortheFouriertransformofthevelocity Riemann solvers included in FLASH, such as HLLE [40] andHLLC[41].Second-orderaccuracyintimeisachieved vˆ(k)= e−2πik·xv(x)d3x. (6) viaacharacteristictracingevolutionoftheRiemannsolver (cid:90)R3 inputstatestothetimestepmidpoint[29].Weremarkthat, inaccordancewiththeILES,paradigm,wedonotinclude Transformingbothsidesofequation(5)weobtain any additional sub-grid scale model, but relied on the im- plicit turbulent closure built in the numerical schemes we ∂tvˆ+vˆ∗2πik⊗vˆ =−Vˆ∗2πikpˆ+aˆ, (7) usefortheintegrationofthehydrodynamicsequation. All of our simulations start with the fluid at rest ρ = 1, where∗denotestheconvolutionoperator,i.e., v = 0. Turbulence is driven using the stirring module of FLASH.ThismoduleusestheOrnstein-Uhlenbeckprocess [f ∗g](k)= f(q)g(k−q)d3q. (8) [42]togeneratestirringmodesinFourierspace.Thisyields (cid:90)R3 anaccelerationfieldwhichsmoothlydecorrelates[43]over Ifwemultiplybothsidesofequation(7)byvˆ∗andtakethe a timescale T . The FLASH implementation permits the s realpart,weobtainanequationforthe3Denergyspectrum use of any arbitrary combination of solenoidal and com- pressive modes [44]. For our runs, we set T = 0.5, we s use only solenoidal forcing and we restrict the accelerat- ∂ E(k)=T(k)+C(k)+(cid:15)(k), (9) ingfieldtobenonzeroonlyinthefirstfourFouriermodes. t This forcing is designed to mimic the influence of some where largerscaleweaklycompressibleflowand,forthisreason, it does not include any compressible component. This is 1 areasonableapproximationforlowMachnumberconvec- E(k)= vˆ·vˆ∗, (10) 2 tion which is well described by the anelastic approxima- T(k)=−2π(cid:60)(vˆ∗ik⊗vˆ)·vˆ∗, (11) tion,e.g.,[45].IntheCCSNcontext,simulationsshowthat the turbulence is highly anisotropic, being roughly twice C(k)=−2π(cid:60)(Vˆ∗ikp)·vˆ∗, (12) as strong in the radial direction as either tangential direc- (cid:15)(k)=(cid:60)aˆ·vˆ∗. (13) tion [46, 14, 47, 15] since it is driven by buoyancy due to anegativeradialentropygradient.Inordertoemulatethis Here E is the energy spectrum (the velocity power spec- behavior, the accelerating field in the x−direction (which tral density (PSD)) and T is the same transfer term as in is going to play the role of the radial direction) is scaled the classical incompressible equations and (cid:15) is the energy byaconstantfactor(beforethesolenoidalprojectionofthe injection rate. The C term vanishes in the incompressible accelerationfield)suchthatRxx (cid:39)2Ryy (cid:39)2Rzz,where limit and represents the interaction between kinetic and acoustic modes. In practice, in our models, C is found to R =(cid:104)ρv v (cid:105), (4) ij i j be at least one order of magnitude smaller than T at all scalesanditisthusnegligible.Inanycase,weretainC in is the Reynolds stress tensor (to simplify the notation we theanalysisbelow. considered a frame in which (cid:104)ρv(cid:105) = 0) and (cid:104)·(cid:105) denotes Foreachofthespectralquantities,S,beingE,T,C or(cid:15), an ensemble average. Finally, the overall strength of the wedefinetheintegratedspectrum,S(k),as stirringistunedtoachieveaRMSMachnumberof(cid:39)0.35, which is typically observed in realistic CCSN simulations [48,49]. S(k)= S(k)δ(|k|−k)dk, (14) (cid:90)R3 Radiceetal. Page4of17 δ(·)beingtheDiracdeltafunction. 2.3 Structurefunctions Integrating equation (9), we obtain the following one- Theenergyspectrumanditssources/fluxesgiveacompre- dimensionalenergybalanceequation hensive picture of the energy cascade and can be used to assess the level of convergence of the simulation. Unfor- ∂ E(k)=T(k)+C(k)+(cid:15)(k). (15) tunately, 3D energy spectra and fluxes are not easily ac- t cessible in calculations in complex domains and/or with Thiscanalsobewrittenintermsoftheenergyfluxacross inhomogeneousturbulence.Inthesecases,localquantities scales, in the physical domain are more easily extracted and ana- lyzed.Hence,oneofthegoalsofthisworkistovalidatethe k use of indirect measures of convergence of ILES. Among Π(k)=− T(ξ)dξ, (16) thesequantities,thestructurefunctionsofthevelocityap- (cid:90)0 peartobenaturalcandidatesforstudy. as Wedefinethevelocityincrements r ∂ E(k)+∂ Π(k)=C(k)+(cid:15)(k). (17) δv(x,r)= v(x+r)−v(x) · (22) t k r (cid:2) (cid:3) Noticethatwedidnotassumeisotropyinanyoftheabove. andstudythequantities Equation(15)isderivedintheinviscidlimit.Inpractice, ourevolutionmethodintroducesdissipationintheformof S (r)=(cid:104)δvp(cid:105) , (23) p j=0 “numerical viscosity”. This can be quantified in terms of theresidual where, (cid:104)·(cid:105)j=0 denotes an ensemble average as well as a meanoveralloftheanglesbetweenvandr(inotherwords R(k)=∂ E(k)−T(k)−C(k)−(cid:15)(k). (18) we are looking at the j = 0 component of the SO(3) de- t compositionofthestructurefunctions[54]).Inthecaseof This can be used to define a wave number dependent nu- homogeneous turbulence Sp does not depend on x and is mericalviscosity: thusafunctionofonlytheseparationr. Themostimportantrelationinvolvingthestructurefunc- 1 R(k) tionsistheso-called4/5−law,whichrelatesthethirdorder ν(k)=− . (19) structurefunction,S (r),withthemeanenergydissipation 2k2E(k) 3 rate, Weremarkthatνdoesnot,ingeneral,correspondtoaclas- ∞ sicalshearorbulkviscosity,butcanneverthelessbeinter- (cid:104)(cid:15)(cid:105)= (cid:15)(k)dk, (24) pretedasarelativemeasureofthedissipationactingatdif- (cid:90)0 ferentwavenumbers(see,e.g.,[51,52,53]foralternative and states that, for incompressible, homogeneous and approaches). isotropicturbulence[50]: In practice, since we will be working in the station- arycase,afterhavingtakentheappropriatetimeaverages, 4 S (r)=− (cid:104)(cid:15)(cid:105)r. (25) R(k)reducesto 3 5 R(k)=−T(k)−C(k)−(cid:15)(k). (20) Equation(25)canbederivedfromtheNavier-Stokesequa- tion for fully-developed, incompressible, homogeneous and isotropic turbulence and it is one of the few exact re- Finally,sinceweareworkinginaperiodicdomain,which lations in the theory of turbulence [50]. In the anisotropic we take of size L = L = L = 1, all of the spectra x y z orcompressiblecase,however,equation(25)isnotstrictly are quantized and non-trivial only for k ,k and k inte- x y z validandcouldbeviolatedinthedata.AsweshowinSec- gers.Furthermore,alloftheintegralsinwavenumberspace tion5,wefindequation(25)tobeverywellsatisfiedbyour reduce to summations. Integrals over spherical shells are data,suggestingthatthe3rdorderstructurefunctioncanbe transformedtoweightedsumsfollowing[43]: averyusefuldiagnosticinglobalsimulations. 4πk2 E(k)= E(k), (21) 2.4 Transverseenergyspectrum N k k−1/2<|k|≤k+1/2 Another alternative to the analysis of 3D spectra, which (cid:88) has been adopted by several authors in the core-collapse where N is the number of discrete wave-numbers in the supernova context [55, 23, 47, 13], is the use of 2D spec- k shellk−1/2<|k|≤k+1/2. tracomputedusingasphericalharmonicsexpansionofthe Radiceetal. Page5of17 velocity field tangential to one or more spherical shells in asinarealCCSNsimulationwherestrongshocksneedto thesimulation.Analogously,weemulatethisbylookingat behandledinsomepartsofthedomain. quantitiesinthey−zplaneandwedefinethe2Dspectra For each group of simulations we run four resolutions: 643,1283,2563 and 5123. The RMS velocity in all of the 1 runs is v (cid:39) 0.4, giving an eddy turnover time τ = E (k )= v˜ ·v˜∗ δ k2+k2−k dk dk , rms ⊥ ⊥ 2 ⊥ ⊥ y z ⊥ y z 1/v (cid:39) 2.5. All of the simulations are run until time (cid:90)R2 (cid:16)(cid:113) (cid:17) (26) t =rm1s00 ((cid:39) 40 eddy turnover times). The time evolu- tionofafewrelevantdiagnosticsisshowninFigure1for where v is the projection of the velocity perpendicular ourfiducialgroupofruns(PPM HLLC)atdifferentresolu- ⊥ to the x−direction and we introduced the partial Fourier tions.Wecanseehowtheflowisacceleratedfromrestand transformofv : quicklyreachesasteady,fullyturbulent,state.Inallcases, ⊥ steady state is reached after t (cid:38) 3 (∼ 1 turnover time) 1 Lx/2 and the diagnostics are insensitive to the resolution. The v˜ (k ,k )= lim dx ⊥ y z Lx→+∞Lx (cid:90)−Lx/2 (27) resultsfortheotherruns(notshown)areverysimilartothe ones of PPM HLLC as they all achieve very similar RMS e−2πi(kyy+kzz)v⊥(x,y,z)dydz. Mach numbers and Reynolds stresses. All of the analysis (cid:90)R2 shownintherestofthepaperareperformedusing3803D InthelimitofinfiniteReynoldsnumber/resolution,the2D snapshots(evenlyspacedintime)ofthedataintheinterval spectrumisexpectedtohavethesameasymptoticbehavior 5≤t≤100. asthe3Dspectrum,howeveritisa-prioriunclearifE is A first, qualitative, comparison between the different ⊥ agoodproxyforE atfiniteresolution.Forthisreasonwe methods can be done by looking at their visualizations. finditusefultoinvestigatethishere. In particular, in Figure 2, we show a visualization of the Aswasthecaseforthe3Dspectra,alsoherethespectrum magnitudeofthevorticityinthex–z planeforfourofthe isnon-trivialonlyforintegerk andk ,whenperiodicityis fiveschemes(excludingPPM HLLC CFL0.8)atthehigh- y z takenintoaccount.Theintegralinequation(26)istreated est resolution (5123). The data is taken at the final time analogouslytotheintegralintheequation(14)forthe3D (t=100).Asitcanbeseenfromthefigure,allofthesimu- case,whiletheaverageinthex−directioninequation(27) lationsshowthepresenceofthin,elongated,regionsofhigh isconvertedtoanaverageoverthex−extentofthesimula- vorticity, as typically seen in direct numerical simulations tionbox. (DNS) of homogeneous turbulent flows [57, 58]. How- ever, the width and the intensity of the vorticity at these 3 Basic flow properties smaller scales depend crucially on the numerical scheme. We employ the finite-volume HRSC (Godunov) approach Methods with small intrinsic numerical viscosity, such as in which physical states are reconstructed at inter-cell PPM HLLC and WENOZ HLLC, present smaller structures boundariesandlocalRiemannproblemsaresolvedtocom- andmoreintermittentvorticityfieldswithrespecttomore pute the physical inter cell fluxes. In particular, we per- dissipativemethods,suchasPPM HLLEandTVD HLLE. form five groups of simulations using different numerical methods. Each group is labeled using the name of the re- 4 The energy cascade construction algorithm and of the Riemann solver. For in- Inthissectionwefocusouranalysisontheaccuracywith stanceTVD HLLE,denotesagroupofsimulationsdoneus- which the energy cascade is captured by our ILES runs. ingTVDreconstructionandHLLERiemannsolver.Single First,wefocusonthelargestscalesofthesimulationwith simulations are labeled using their resolution so that, for the goal of quantifying the accuracy in the decay rate of instance, TVD HLLE N128, denotes the TVD HLLE run the energy as a function of the resolution for the different done using a 1283 grid. For all of the runs the timestep methods. Next, we will look at the energy distribution at ischosentohaveaCFL,i.e.,c∆t/∆x,of0.4,cbeingthe smaller scales where, in resolved simulations, the inertial maximum characteristics speed, with the exception of the range starts. Finally, we will look at the dynamics in the PPM HLLC CFL0.8 runs where we set the CFL to 0.8. dissipationregionandsummarize. For the TVD runs we use the monotonized central (MC) slopelimiter[30].TherunswithPPMusetheoriginalflat- teningandartificialviscosityprescriptionsfrom[29].The 4.1 Energydecayrate artificialviscositycoefficientis0.1.Weremarkthattheuse InthelimitofverylargeReynoldsnumberitisassumed, oftheartificialviscosityforPPMisnotreallynecessaryin instandardturbulencephenomenology[50],thatthereex- thisregime[56],howeverourgoalisnottoperformastudy istsarangeofwavenumbers(theinertialrange)whereen- oftheturbulentdynamics,buttoassesshoweachnumeri- ergyinjectionanddissipationcanbeneglectedinequation calmethodperformswhenusedunderthesamecondition (17). In this range we can write (compressible effects are Radiceetal. Page6of17 0.4 3 3 y z ms 0.3 Ry 2 Rz 2 Mr /x /x x x R R 0.2 PPMHLLCN64 1 1 PPMHLLCN128 PPMHLLCN256 PPMHLLCN512 0.1 0 0 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 t t t Figure1 TimeevolutionofthediagnosticquantitiesforthefiducialsetofrunsPPMHLLCwithdifferentresolutions.Theleftpanel showstherootmeansquare(RMS)Machnumber,whilethemiddleandleftpanelsshow,respectively,theratiosRxx/Ryyand Rxx/Rzz,RbeingtheReynoldsstresstensor(equation4).Sincethex directionistheanisotropicdirection(itwouldplaytheroleof − theradialdirectioninaCCSN)theratiosRxx/RyyandRxx/Rzz,offeraglobalmeasureoftheanisotropyoftheflowatthelargest scale.Allofthequantitiesappeartohavereachedstationarityaftertimet(cid:38)3andoscillatearoundtheirtargetvalues.Allresolutions producethesamequalitativebehavior. negligibleinoursimulations): 1.0 ∂ E(k)+∂ Π(k)(cid:39)0, (28) t k so that stationarity requires Π(k) (cid:39) const. In particular, /ǫhi0.9 since energy is conserved, one finds Π(k) (cid:39) (cid:104)(cid:15)(cid:105). This max meansthat,inthelimitoflargeReynoldsnumbers,theen- Π PPMHLLC ergydecayratedependsonlyonthemacroscopicproperties 0.8 PPMHLLCCFL0.8 PPMHLLE of the flow (and in particular not on the nature of the vis- TVDHLLE cosity),afactthathasalsobeenverifiednumerically[59]. WENOZHLLC 0.7 Thesignificanceofthispropertyanditsimportanceforthe 64 128 256 512 modelingofturbulencecannotbeoverstated. N In the context of CCSN simulations this means that the Figure4 Dissipationrateoftheenergyatthelargestscales largescalekineticenergy,acrucialquantityforthedynam- duetotheturbulentcascade(notincludingdirectdissipationby ics of the explosion [15], can be faithfully captured even thenumericalviscosity)asafunctionofresolutionandforallof theschemes.Thedissipationrateisnormalizedsoastobe1 withsimulationsachievingmodestReynoldsnumbers. inthelimitoflargeReynoldsnumbers/resolution.At1283 For an ILES, a basic requirement, then, is that a suffi- pointsalloftheschemesshowanerroroflessthan10%,with cientlyhighresolutionshouldbeachievedtocorrectlyrep- theHLLCschemesalreadyclosetothe2%level. resenttheenergycascadeatthelargestscales.Whatqual- ifies as a sufficiently high resolution is of course depen- dentonthedetailsoftheclosurebuiltintothescheme(and We can make a more quantitative statement concerning ontheaccuracyrequiredfortheparticularapplication).To theenergydecayratebylookingatthepeakoftheenergy quantifythis,wecanestimatethelevelofaccuracythatcan flux as a function of resolution, as shown in Figure 4. We be reached at any given resolution, using our local simu- canseethatat1283pointsallofthesimulationshaveade- lations.Inparticular,wecanstudydirectlytheenergyflux viationfromtheasymptoticenergydecayrateoflessthan across scales, defined by equation (16). This is shown in 10%. The least dissipative methods already have an error Figure3forallofthedifferentruns,. close to the 2% level. A comparison between PPM HLLE Asdiscussedbefore,weexpectthatΠ(k) (cid:39) (cid:104)(cid:15)(cid:105)overan andPPM HLLCrevealstheprofoundimpactthatthechoice extendedregioninFourierspaceshouldbeadirectindica- of the Riemann solver has even at relatively large scale tion that a simulation has been able to recover an inertial (moreonthedissipativepropertiesofthedifferentschemes range.Perhapsnotsurprisingly,inlightofpreviousresults inSection4.3). [24],wefindthatregionswhereΠ (cid:39) (cid:104)(cid:15)(cid:105)aswideasafew wavenumbers4(cid:46)k (cid:46)10onlyappearatthehighestreso- 4.2 Energyspectra lutions(wewilldiscusstheinertialrangeinmoredetailin Obviously, not all of the dynamics of turbulence can be Section4.2).However,theamountofenergydecayingfrom reducedtotherateatwhichkineticenergydecaysfromthe thelargestscalesreachesanasymptoticvaluemuchquicker injection scale. The internal dynamics of the energy cas- thanthatimplyingthatthetotalkineticenergybudgetatthe cade,farfromtheinjectionscaleandfarfromthedissipa- largestscalesiswellresolvedevenatmodestresolutions. tion range, can also play an important role in many appli- Radiceetal. Page7of17 PPM HLLE N512 PPM HLLC N512 TVD HLLE N512 WENOZ HLLC N512 x x   Figure2 Squarerootofthemagnitudeofthevorticity,(cid:112) v,forfourofthesimulationswith5123resolutioninaslicethroughthe |∇× | middleofthex–zplaneatthefinaltimeofthesimulations(t=100).ThepanelsshowsimulationsusingPPMHLLEN512, PPMHLLCN512,WENOZHLLCN512,andTVDHLLCN512clockwisefromthetopleft.Thedirectionoftheanisotropicdrivingisupin thesefigures.Thecolorcodegoeslinearlyfrom0(novorticity;darkcolors)to15(lightcolors)anditisthesameforallpanels. cations.ToanalyzethisaspectweconsiderinFigure5the grid point, 512k∆x = 256 corresponds to a wavelength energy spectrum of the velocity defined by equation (10). oftwogridpointsandsoon. The spectra are compensated by k5/3 to highlight regions Looking at any of the groups of runs in Figure 5, one withKolmogorovscaling,whichmightbeexpectedinthe can immediately notice that the spectra obtained at differ- ent resolutions do not collapse into a single curve in the inertial range. Since we want to focus on quantities that dissipationregion,aswouldberequiredbyKolmogorov’s donotdepend(ordependweakly)onthenatureoftheen- firstsimilarityhypothesis[50](cf.,[60]).Thislackofcon- ergy injectionat largescale, weshow allof the spectraas vergenceinthedissipationregioncouldbeduetothenon- afunctionofadimensionlesswavenumber,512k∆x.The linearviscosityofHRSCschemes.This,inturn,couldre- rationale behind this normalization is that, first of all, we sult in an anomalous scaling of η with the grid spacing. assume the Kolmogorov scale η to be proportional to the SuchscalinghasbeenreportedinthepastforILES,butit gridspacing.Secondly,the512factorisintroducedtohave is notvery wellunderstood [52].The good agreementbe- the dimensionless k, 512k∆x coincide with the dimen- tweenthethreedifferentgroupsofsimulationsemploying sionaloneforthehighestresolutionruns.Withthischoice, theHLLCRiemannsolverseemstosupportthishypothesis 512k∆x = 512 corresponds to a wavelength of a single andsuggeststhatthenonlinearviscosityintroducedbythe Radiceetal. Page8of17 1.0 643 1283 0.8 2563 k)/ǫhi0.6 643 643 5123 Π( 0.4 1283 1283 2563 2563 0.2 5123 5123 PPMHLLC PPMHLLCCFL0.8 PPMHLLE 0.0 1.0 643 N =5123 1283 0.8 2563 k)/ǫhi0.6 5123 643 PPMHLLC Π( 0.4 1283 WENOZHLLC 2563 PPMHLLE 0.2 5123 TVDHLLE TVDHLLE WENOZHLLC PPMHLLCCFL0.8 0.0 100 101 102 100 101 102 100 101 102 k k k Figure3 Energyflux,asdefinedbyequation(16),obtainedwithdifferentnumericalmethodsandresolutions.Theenergyfluxis shownnormalizedtotheaveragedissipationrategivenbyequation(24).Fromlefttorightandfromtoptobottomweshowtheresults obtainedwithPPMHLLC,PPMHLLCCFL0.8,PPMHLLE,TVDHLLEandWENOZHLLC.Thebottomrightpanelshowacomparisonof allofthemethodsat5123.Alloftheschemesshowagoodlevelofaccuracyintheenergyfluxfromthelargestscales,witherrors smallerthanafew%alreadyatlowresolutions.Thedifferencesbetweentheschemesbecomemoremarkedatlargewavenumbers wherethenumericaldissipationstartstointerferewiththeenergycascade. Riemann solver is an important ingredient in setting this correspondtoanextremelyhighresolutioninglobalCCSN scaling. simulationsthattypicallyhave∼30−70zonesacrossthe Convergence appears to be recovered at larger scales turbulentregion[13]. (cid:38) 8∆x (512k∆x (cid:46) 64), but the spectra appear to be The overall behavior of the spectra, as obtained by all dominated by the bottleneck effect. This manifests itself schemes,isconsistentwithKolmogorov’stheoryofturbu- as a bump in the compensated spectra extending from the lence.Theanisotropiccontributionstotheangle-integrated dissipation range until the end of the inertial range, for spectraaretoosmalltobedetectedinourdata. the simulations that show one (e.g., until 512k∆x = 10 for the HLLC runs), or until the energy injection scale 4.3 Numericalviscosity (512k∆x = 4), for the simulations that show no or lit- At very small scales (∼ several grid points) the dynam- tle inertial range (TVD HLLE). The bottleneck effect is a ics is dominated by the numerical viscosity. This can be viscous phenomenon which is also observed in direct nu- estimated from the residual of the energy equation (17) mericalsimulations.However,inthepresentcontextwhere or, equivalently, by the effective numerical viscosity ν(k) viscosityisofnumericalorigin,itisattheveryleastques- (equation 19). The latter is shown in Figure 6 for all tionableifapronouncedbottleneckisadesirablefeatureof schemesandresolutions. the modeling. In astrophysical flows, where the Reynolds The first thing to notice is that the numerical viscosity numbers are typically very large, this pile up of energy at providedbyallnumericalschemesisnotconstant,butdif- largescalesisunphysicalandcouldaffectthequantitative fers by roughly an order of magnitude between low and andqualitativeoutcomeofasimulation[13].Aquantifica- tionofthe bottleneckeffectintermsof theenergybudget highk.Havingawavenumberdependentviscosityisade- isdiscussedinSection4.4. sirablefeatureexpectedinanyLESmodel(explicitoroth- Atevenlargerscales,aninertialrange(E ∼ k−5/3 and erwise).Nevertheless,thismakesthedefinitionandcalcu- Π ∼ const, see Figure 3) seems to be recovered by the lationoftheeffectiveReynoldsnumberachievedinasimu- leastdissipativeschemes(PPMandWENOZwithHLLC) lationambiguous.MeaningfulwaystoestimateitforILES in the region 4 (cid:46) k (cid:46) 10. PPM HLLE and TVD HLLE have been proposed [53] and they can be used to ease the have a more limited region, a few wave numbers at most, comparisonbetweendifferentsimulationsandassesstheir thatcouldbeinterpretedasbeinganinertialrange.Wenote quality. However, one has to be very careful while using that this resolution is not particularly high in comparison anyquoted“Reynoldsnumber”fromanILES,toestimate with state of the art DNS [59, 61], but it would already things like the dynamical range achieved by a simulation, Radiceetal. Page9of17 PPMHLLC PPMHLLCCFL0.8 PPMHLLE 10−1 /3 5 k E(k) 10−2 643 2563 643 2563 643 2563 10−3 1283 5123 1283 5123 1283 5123 TVDHLLE WENOZHLLC N =5123 10−1 /3 5 k E(k) 10−2 PPPPMMHHLLLLCCCFL0.8 PPMHLLE 643 2563 643 2563 TVDHLLE 10−3 1283 5123 1283 5123 WENOZHLLC 100 101 102 100 101 102 100 101 102 512k∆x 512k∆x 512k∆x Figure5 Energyspectra(equation10)obtainedwithdifferentnumericalmethodsandresolutions.Theenergyspectraare compensatedbyak5/3spectrum,sothatanyregionwithKolmogorovscalingshouldappearroughlyflat.Furthermore,thespectra areallplottedasafunctionofthedimensionlesswavenumber512k∆x(the512factorisintroducedtohavethedimensionlesswave numbercoincidewiththedimensionaloneforthe5123runs).ThefirstfivepanelsshowthePPMHLLC(upperleft),PPMHLLCCFL0.8 (uppercenter),PPMHLLE(upperright),TVDHLLE(lowerleft)andWENOZHLLC(lowercenter)groupofruns.Thelastpanel(lower right)showsacomparisonofallofthemethodsatthehighestresolution(5123).Aninertialrangeseemstoberecoveredonlyatthe highestresolutions(perhapswiththeexceptionofTVDHLLEwherenoinertialrangeisvisible).AllschemesemployingtheHLLC Riemannsolverareinverygoodagreement. becausethedissipativepropertiesofILESdifferconsider- well captured by the ILES (Section 4.1), and with the en- ablyfromtheonesofthetrueNavier-Stokesequations. ergytransferintheinertialrange,whichwehaveseentobe Twootherfeaturescanbeobservedinmostofthenumer- describedaccuratelyonlyatmuchhigherresolutions(Sec- icalviscosityprofiles.First,manyofthemexhibitasudden tion 4.2). In a turbulent flow both of these aspects are im- reversalathighwavenumbers.Thisisduetothefactthat portant and a good ILES should display a distribution of the numerical viscosity does not behave like a shear vis- energy across vortical structures at different scales that is cositysothat,althoughthenumericaldiffusionisstrongat ascloseaspossibletotheasymptoticone.Obviously,there thosescales,thenumericalviscosityappearssmallbecause is a limit to the accuracy that any ILES can achieve at a of a partial decoupling between vorticity and dissipation. fixedresolution.Here,wemakethisstatementmorequan- Second,athighresolutionandatthelargestscales,thenu- titativebyconsideringtheamountofkineticenergythatis mericalviscosityisclosetozeroorevenslightlynegative. wellresolvedbyeachsimulationatagivenresolution. The reason is that the residual of equation (9) oscillates aroundzeroanditistoosmalltobereliablyextractedfrom We introduce the cumulative energy spectrum, the inte- ourdata:amuchlongerintegrationtimewouldbeneeded graloftheenergyspectrum: toaccumulateenoughstatisticsforit. Finally, a comparison between the numerical viscos- k ity reveals two interesting effects. First, by comparing E(k)= E(ξ)dξ. (29) PPM HLLC and PPM HLLE, we see that the choice of the (cid:90)0 Riemannsolveraffectstheviscosityatbasicallyallscales. Second, if we compare PPM HLLC, PPM HLLC CFL0.8 ThisquantityisplottedinFigure7,whereitisnormalized and WENOZ HLLC, we see that doubling the timestep ap- by pears to have an effect comparable to the difference be- tweenthePPMandWENOZreconstructionsatintermedi- v2 +∞ atescales(40(cid:46)k (cid:46)100). rms = E(k)dk (30) 2 (cid:90)0 4.4 Theenergydistribution So far we have been concerned with the energy decay toobtainthecumulativedistributionfunctionofthekinetic rate from the largest scales, which we have shown to be energy.Asareference,wealsoshowthecumulativeenergy Radiceetal. Page10of17 10−1 643 2563 643 2563 643 2563 1283 5123 1283 5123 1283 5123 10−2 k) ( ν 10−3 PPMHLLC PPMHLLCCFL0.8 PPMHLLE 10−4 TVDHLLE 643 2563 PPMHLLC 1283 5123 PPMHLLCCFL0.8 10−2 PPMHLLE TVDHLLE k) WENOZHLLC ( ν 10−3 643 2563 1283 5123 WENOZHLLC N =5123 10−4100 101 102 100 101 102 100 101 102 k k k Figure6 Numericalviscosityasafunctionofthewavenumbermeasuredforallschemesandresolutions.Thenumericalviscosityis estimatedusingtheprocedureoutlinedinSection2anditisdefinedbyequation(19).Thedifferentpanelsare,fromlefttorightand fromtoptobottom,theresultsobtainedwithPPMHLLC,PPMHLLCCFL0.8,PPMHLLE,TVDHLLEandWENOZHLLC.Thebottom rightpanelshowacomparisonofallofthemethodsat5123.Thenumericalviscosityshowslargevariationsacrossthewavenumber space.ThechoiceoftheRiemannsolverplaysarolethatisatleastasimportantasthechoiceofthereconstructionmethodin affectingthenumericalviscositythroughouttheentirethespectrum. spectrumestimatedfromKolmogorov’stheory: case, however, anisotropic contributions cannot be ex- cluded [64], although they are known to be subdominant k in some important cases [65, 66, 67]. In this section we E (k)= E (ξ)dξ, (31) K41 K41 show that equation (25) is consistent with our data over a (cid:90)0 widerangeofscales. E (k)= K41 We compute the 3rd order structure functions of the ve- EPPMHLLCN512(k), ifk ≤4, (32) locity,definedbyequation(23),inarathersimplewayus- (cid:40)EPPMHLLCN512(4) k4 −5/3, k >4. ing a random sample of 20,000 points in each of the 380 3Ddatadumpsofoursimulations.Ateachtime,wecom- We find that as the resolution(cid:0)in(cid:1)creases, all schemes ap- putethe3rd powerofthevelocityincrementsforeachpair pear to be converging to the predictions of Kolmogorov’s of points and accumulate and average in time the results in bins of size ∆x. The resulting structure functions are theory.Theresultsatfiniteresolution,however,arenoten- couraging:at643 only∼50%orlessofthekineticenergy shown in Figure 8, compensated by −54r−1(cid:104)(cid:15)(cid:105)−1, so that isinwellresolvedstructures,whiletheother∼ 50%have theresultingquantityshouldbeequaltooneifthe4/5−law issatisfiedinourdata.Aswasthecasefortheenergyspec- piled up at rather large scale, with a cumulative excess of tra, we assume η ∼ ∆x and plot the structure functions ∼ 20%atthegridscale,mostlybecauseofthebottleneck versusr/∆x. effect.Athigherresolutions,theamountofkineticenergy wellcapturedbytheILESincreases,butat5123thisisstill The degree with which the 4/5−law is satisfied in our data is very good. We see that anisotropic contributions onlyabout80%oftheenergyandthereisstillacumulative onlyplayaminorroleintheangle-integratedformulation excessofover∼5%atthegridscale((cid:96)∼∆x). ofthe4/5−law.Thisisinagreementwiththeincompress- ible DNS of [67] and has been known to be true also for 5 The 4/5−law Rayleigh-Be´nardconvectioninmostregimes[68].Ourre- The 4/5−law (equation 25) is not a-priori valid in the sultsprovideanimportantnewexamplewherethisappears regime of turbulence we are considering. However, the to hold true. Secondly, for all of our simulations at 5123, 4/5−law has been numerically verified to hold also in wefind somesituationsoutsidethedomainofvalidityofitsderiva- tion. For instance, for isotropic mildly compressible de- 5 caying [62] and driven [63] turbulence. In the anisotropic mrax −4r−1(cid:104)(cid:15)(cid:105)−1(cid:104)δv3(cid:105)j=0 (33) (cid:26) (cid:27)