Implicit large eddy simulations of anisotropic weakly compressible turbulence with application to core-collapse supernovae

1 MB·English

by  David Radice

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)