A localised subgrid scale model for fluid dynamical simulations in astrophysics II: Application to type Ia supernovae

0.38 MB·English

by  W. Schmidt

Astronomy&Astrophysicsmanuscriptno.3618 February5,2008
(DOI:willbeinsertedbyhandlater)
A localised subgrid scale model for fluid dynamical simulations
in astrophysics II: Application to type Ia supernovae
W.Schmidt1,2,J.C.Niemeyer1,W.Hillebrandt2,andF.K.Ro¨pke2
1 Lehrstuhlfu¨rAstronomie,Institutfu¨rTheoretischePhysikundAstrophysik,Universita¨tWu¨rzburg,AmHubland,D-97074
Wu¨rzburg,Germany
6
2 Max-Planck-Institutfu¨rAstrophysik,Karl-Schwarzschild-Str.1,
0
0 D-85741Garching,Germany
2
Received/Accepted
n
a
J Abstract. Thedynamicsoftheexplosiveburningprocessishighlysensitivetotheflamespeedmodelinnumericalsimulations
3 of type Ia supernovae. Based upon the hypothesis that the effective flame speed is determined by the unresolved turbulent
2 velocityfluctuations,weemployanewsubgridscalemodelwhichincludesalocalisedtreatmentoftheenergytransferthrough
theturbulencecascadeincombinationwithsemi-statisticalclosuresforthedissipationandnon-localtransportofturbulence

1 energy. Inaddition, subgrid scale buoyancy effectsareincluded. Inthelimitof negligible energytransfer and transport, the
v dynamicalmodelreducestotheSharp-Wheelerrelation.Accordingtoourfindings,theSharp-Wheelerrelationisinsuffcient
0 toaccountforthecomplicatedturbulentdynamicsofflamesinthermonuclearsupernovae.Theapplicationofaco-movinggrid
0
techniqueenablesustoachieveveryhighspatialresolutionintheburningregion.Turbulenceisproducedmostlyattheflame
5
surfaceandintheinteriorashregions.Consequently,thereisapronouncedanisotropyinthevicinityoftheflamefronts.The
1
localised subgrid scale model predictssignificantly enhanced energy generation and lessunburnt carbon and oxygen at low
0
velocitiescomparedtoearliersimulations.
6
0
/ Keywords.Stars:supernovae:general–Hydrodynamics–Turbulence–Convection–Methods:numerical
h
p
-
o 1. Introduction Niemeyer1999;Zingaleetal.2005).Intheaforementionednu-
r
merical models, a DDT is artifically triggeredat more or less
t
s ForsupernovaeoftypeIa,Hoyle&Fowler(1960)proposeda arbitraryinstantsoftime.
a
thermonuclear runaway initiated in C+O white dwarfs close
: As for the flame speed model, the controversyis whether
v to the Chandrasekhar limit as the cause of the explosion.
subgridscale (SGS) turbulenceis mostly drivenbyRayleigh-
i
X Since the original proposal, there has been vivid controversy Taylorinstabilitiesordominatedbytheenergytransferthrough
ofhowsuchanexplosionmightcomeaboutandwhattheex-
r theturbulencecascade.Theformerpointofviewholdsthatthe
a act physical mechanism could be. Today the computational
magnitude of SGS velocity fluctuations v is basically given
′
facilities to process three-dimensionallarge-eddy simulations
by the Sharp-Wheeler relation (Davies&Taylor 1950; Sharp
(LES)oftheexplosioneventareavailable.Remarkably,these
1984)
powerfulmeanshavenotaidedinarrivingata consensusyet.
v (l)=0.5 lg (1)
RT eff
The disagreementstems from some crucial questions. Firstly,
what is the appropriate flame speed model? Secondly, does where v (l) is the asymptotic risepvelocity of a perturbation
RT
the explosion completely proceed as a deflagration, or does of size l due to buoyancy. The effective gravity g is deter-
eff
a transition to a delayed detonation set in at some point? minedbythe densitycontrastat theinterfacebetweenburned
The deflagration to detonation transition (DDT) proposed by andunburnedmaterial.Settingtheturbulentflamespeedequal
Khokhlov(1991)andWoosley&Weaver(1994)appearstore- to v (∆), where ∆ is the resolution of the numerical grid,
RT
solve the drawbacks of the pure deflagration model. In par- has been used in some simulations of type Ia supernovae
ticular, the energy output obtained from simulations with ar- (Gamezoetal. 2003; Calderetal. 2003; Gamezoetal. 2004).
tificial DDT is closer to the observed one, and less car- However,simple scaling argumentsdisfavourthis proposition
bon and oxygen is left behind (Gamezoetal. 2004, 2005; (Niemeyer&Hillebrandt 1995; Niemeyer&Kerstein 1997).
Golombek&Niemeyer2005). Forthe theoreticalunderstand- Assuming that non-linear interactions between turbulent ed-
ing of thermonuclearsupernovae,however,the lack of a con- dies of different size l set up a Kolmogorov spectrum, the
vincingexplanationfortheinitiationofthetransitionisunsat- rootmeansquareturbulentvelocityfluctuationsobeythescal-
isfactory (Khokhlovetal. 1997; Niemeyer&Woosley 1997; ing law v(l) l1/3. Since the Sharp-Wheelerrelationimplies
′
∝
2 W.Schmidtetal.:AlocalisedsubgridscalemodelforfluiddynamicalsimulationsII
v (l) l1/2, we have v (l)/v(l) l1/6 0 towards de- fluctuations. The form of the Archimedian force will be pro-
RT RT ′
∝ ∝ →
creasinglengthscales.Consequently,Niemeyer&Hillebrandt posedinSect.2.1.Moreover,theEuleriantimederivativemust
(1995) adopted a SGS model based on the dynamical equa- account for the rescaling of the turbulence energy due to the
tion for the turbulence energy k , i.e. the kinetic energy of temporalshift of the cutoff length. We will denote the partial
sgs
unresolvedvelocityfluctuations(Schumann1975).Themajor derivative with respect to the rescaled quantities by ∂⋆. The
weaknessoftheirapproacharisesfromthefairlytentativeclo- Lagrangiantimederivativethusbecomes
sureswhichwereformulatedadhocforLESofstellarconvec-
D⋆ ∂⋆
tion(Clement1993). = +v ∇. (3)
Variousrefutationsof the scaling argumenthave been put Dt ∂t ·
forward. To begin with, the spectrum of turbulence energy ThecompletedynamicalequationfortheSGSturbulentve-
might be different from the Kolmogorov spectrum. However, locityis
recent direct numerical simulations support the hypothesis of
aKolmogorovspectruminbuoyancy-driventurbulentcombus- D⋆ 1
q ∇ ρℓ q ∇q ℓ ∇q 2
tion(Zingaleetal.2005).Amoreseriousconcernisthatthere Dt sgs−ρ · κ sgs sgs − κ| sgs|
mightbe notenoughtime to reach the state of developedtur- 1(cid:16) (cid:17) 7 q2 (4)
bouflaenscuepewrintohvaaKexoplmloosgioonr.ovThspisecqturuesmtioinntihsedtirffiancsuiletnttoscseetntlaeriao = √2CAgeff+ℓν|S∗|2− 30qsgsd− ℓsǫgs.
priori.Forthisreason,wetookanunbiasedpointofviewand The rate-of-strain scalar S is defined by S 2 = 2S S ,
accommodatedbuoyancyeffectsintheformofanArchimedian | ∗| | ∗| i∗j i∗j
whereS isthetrace-freepartofthesymmetrisedJacobianma-
forcetermintheSGSturbulenceenergymodel. i∗j
trixofthevelocityfield,S = 1(∂ v +∂v ),andd =S =∂v
In contrast to the previously used SGS turbulence energy ij 2 j i i j ii i i
isthedivergence.Thecharacteristiclengthscalesℓ ,ℓ andℓ
model,thenewlocalisedmodelwhichisthoroughlydiscussed κ ν ǫ
arerelatedtoSGSturbulenttransport,therateofenergytrans-
inpaperIneitherpresumesisotropynoracertainturbulenceen-
fer from resolved toward subgrid scales and the rate of vis-
ergyspectrumfunction.Thisispossiblebyvirtueofadynam-
cous dissipation. Each characteristic length can be expressed
icalprocedureforthedeterminationoftheSGSeddy-viscosity
in termsofthe effectivecutofflength∆ anda similaritypa-
ν = C ∆k1/2 which was adapted from Kimetal. (1999). eff
sgs ν sgs rameter:
Furthermore, we apply the co-expanding grid introduced by
Ro¨pke (2005). The growth of the cutoff length ∆ due to the C ∆ 2√2∆ C ∆
ℓ = ν eff, ℓ = eff, ℓ = κ eff. (5)
grid expansion poses a challenge for the SGS model because ν ǫ κ
√2 Cǫ √2
of the partitioning between resolved energy and SGS energy
changesintime.Wewillshowthatthisrescalingeffectcanbe For the determination of the closure parameters Cν, Cǫ and
takenintoaccountbyutilisingthedynamicalprocedureforthe Cκ see Sect. 4 of part I. The advection of qsgs by the re-
calculation of eddy-viscosity parameter C . The rescaling al- solvedflowiscomputedwiththepiece-wiseparabolicmethod
ν
gorithmaswellasthecomputationoftheArchimedianforceis (Colella&Woodward 1984). Due to the dissipative effects of
explainedinSect.2,followedbythediscussionofresultsfrom this numericalscheme on the smallest resolved length scales,
three-dimensionnumericalsimulationsinSect.3.Itisdemon- weset∆eff 1.6∆(Schmidtetal.2005b).Thediffusionterms
≈
stratedthatthenewlyproposedSGSmodelsubstantiallyalters on the left hand side of equation 4 is computed by means of
thepredictionsofthedeflagrationmodel.Inparticular,wewill fourth order centred differences, and for the source term on
analysethesignificanceofSGSbuoyancyaffects. theright-handsideasemi-implicitAdams-Moultonmethodis
used.Intheremainderofthissection,wedescribethecalcula-
tionoftheSGSArchimedianforceandtherescalingprocedure.
2. Theflamespeedmodel
For the relation between the turbulent flame speed s and the
t
2.1.Archimedianproduction
SGSturbulencevelocityq ,weadheretotheresultsfoundby
sgs
Pocheau(1994)fromatheoreticalanalysisandset There has been an ongoing debate whether the production of
SGS turbulence is dominated by the buoyancy of SGS per-
q 2
s = s 1+C sgs , (2) turbations in the interface between ash and fuel or by eddies
t lam t
s slam! producedthroughtheturbulencecascade.Inthefirstcase,the
whereC =4/3.Intheasymptoticregimeofturbulentburning, sourceofenergyisthegravitationalpotentialenergy,whereas
t
s 2q /√3whichisconsistentwithPeters(1999). non-linear transfer supplies kinetic energy from larger scales
t sgs
≃
Theevolutionofq isgivenbyanon-linearpartialdiffer- in the second case. In general, it is quite difficult to sepa-
sgs
entialequationwhichisobtainedbydividingtheequationfor rate the energy injection caused by gravity in wave number
the specific SGS turbulence energy k = 1q2 (see Sect. 3 spacebecausegravitationaleffectsaregenuinelynon-local.For
sgs 2 sgs
ofpartI)byq .ForturbulencedrivenbytheRayleigh-Taylor Rayleigh-Taylordriventurbulence,however,weknowthesim-
sgs
instability,anadditionalsourcetermstemsfrombuoyancyef- ple Sharp-Wheeler scaling relation (1) which provides an al-
fectsonsubgridscales.ThisArchimedianforce,whichispro- gebraic closure for the Archimedian force density Γ intro-
sgs
portionaltotheeffectivegravityduetothedensitycontrastbe- duced in Sect. 3 of part I. Therefore, we propose a novel ap-
tweennuclearfuelandash,directlyinducesturbulentvelocity proachwhichcombinesboththeproductionofSGSturbulence
W.Schmidtetal.:AlocalisedsubgridscalemodelforfluiddynamicalsimulationsII 3
throughthecascadeandtheRayleigh-Taylormechanisminthe Theclosure(8)doesnotincludealloftheintricateeffects
dynamicalequationforq . that contribute to SGS buoyancy. It merely captures what is
sgs
Because Γ has the dimension of an acceleration times presumablytheleadingordereffect.Infact,onewouldhaveto
sgs
velocity, we interpret this term as the product of a specific model the interaction between turbulent potential and kinetic
Archimedianforce and the SGS turbulencevelocityq . This energyfluctuationsonunresolvedscales. Unfortunately,there
sgs
means the following:In any finite-volume cell those portions existsnotheoreticalframeworkforthistaskyet.Moreover,the
ofthefluidwithdensitylessthanthesmootheddensityρwill conceptofaSGSArchimedianforceentailsaviolationofen-
experience buoyancy relative to the other portions of higher ergy conservation because the contributionof Γ to the pro-
sgs
density in the mean gravitational field. In subsonic turbulent duction of turbulence energy is not balanced in the resolved
flows, the small random density fluctuations caused by com- energybudget.Consequently,thetotalenergyofthesystemef-
pressionandrarefactionareexpectedtoproduceverylittlenet fectively increases. However, Γ is non-zero only in a small
sgs
buoyancy. However, a special situation is encountered in the volumefraction.Forthisreason,theresultingviolationofen-
cells intersected by the flame fronts. Since the flames are far ergyremainsnegligiblerelativetothetotalenergybudget.
from being completely resolved in numerical simulations of In order to determine the parameter C , we observe that
A
SNe Ia, the substructure of the front in combination with the the Sharp-Wheeler SGS model is obtained as an asymptotic
density gradientacross the frontwill produceSGS buoyancy. relation in the limiting case of neglectingthe non-localtrans-
Equivalently, one can think of of perturbations in the flame port D , the turbulent energy transfer Σ and the pressure-
sgs sgs
frontonlengthscalesλ .l.∆asbeingRayleigh-Taylorun- dilatationλ (seeSect.2ofpaperIfordefinitions).Dropping
fp sgs
stableandproducingSGSturbulence.Thefirepolishinglength thecorrespondingtermsinequation(4),oneobtains
λ then marks the lower threshold for perturbations to grow
fp
(see Zingaleetal. 2005). Once turbulence has developed, λfp Dq 1 C g q2sgs. (11)
can be identified with the Gibson length lG, i.e. the smallest Dt sgs ≃ √2 A eff − ℓǫ
lengthscaleonwhichtheflamepropagationisaffectedbytur-
for a fluid parcelin the vicinity of the flame front. In the sta-
bulent eddies. Perturbationsof size l & ∆, on the other hand,
tionaryregime,thisequationhasthefixedpointsolution:
setfluidintomotiononnumericallyresolvedlengthscales.But
the transfer of energythrough the turbulence cascade eventu-
2C ∆ g
allyproducesSGSturbulenceaswell.Thisproductionchannel q A eff eff =v (∆ ). (12)
sgs RT eff
correspondstothetermΣsgsintheequationfortheSGSturbu- ≃ s Cǫ
lenceenergy(seesect3ofpartI).
ConsistencywiththeSharp-Wheelerrelation(1)impliesC =
The Archimedian force generated by the density gradient A
C /8. Since C 0.5...1.0 for developed turbulence (see
acrossflamefrontsisgivenbytheeffectivegravity ǫ ǫ ≈
Schmidtetal.2005a),theestimateC 0.1isobtained.
A
≈
g =Atg, (6)
eff
wheretheAtwoodnumber 2.2.Rescalingofthesubgridscaleturbulenceenergy
ρ ρ
At= f − b (7) If a non-static, co-moving grid is used, the implicit filter
ρ +ρ
f b introducedin paper I, Sect. 3, becomes time-dependent.
eff
h i
is a measureforthe densitycontrastbetweenburnedmaterial Therefore,time-derivatesdonotcommutewiththeoperationof
and nuclearfuel. In the vicinity of the flame front,the lowest filteringandadditionaltermsariseinthedynamicalequations.
order estimate for the SGS buoyancyterm is Γsgs ρqsgsgeff, These terms are equivalent to the additional fluxes which are
providedthat∆>λfp.If∆becomessmallerthanλfp∼,SGSper- includedintheimplementationoftheRiemannsolverformov-
turbationsinthe∼flamefrontarenotsubjecttotheRT-instability inggrids(Mu¨ller1994;Ro¨pke2005).However,thereisasub-
andΓsgsvanishes.Inconclusion,weproposethefollowingclo- tletyrelatedtotheshiftingcutoffwhichseparatestheresolved
sure: andunresolvedscales.Asthegridexpandshomologouslywith
1
Γ = C ρg q . (8) the bulk of the white dwarf, the grid resolution ∆ gradually
sgs A eff sgs
√2 decreasesintimeandthegrowingcutofflengthentailsagrad-
Heretheeffectivegravityismorepreciselydefinedby ual rise of the SGS turbulenceenergy.This rise is inherentto
the grid geometry and immediately affects the decomposition
g =χ (G =0)θ(∆ λ )At(ρ)g, (9)
eff ±δ − fp oftheenergybudget.Thetwo-thirdslawfordevelopedturbu-
where χ n∆(G = 0) is the characteristic function of all cells lenceimplies qsgs ∆1/3 (seeFrisch1995, Sect.5).Thus,it
± h i∝
for which the distance from the flame front (represented by iseasytorescalethemeanvalueofq if∆changesbyasmall
sgs
G(x,t) = 0) is less than δ, and θ is the Heaviside step func- fractionδ∆/∆:
tion, i.e. θ(∆ λ ) = 1 for ∆ > λ and zero otherwise. The
fp fp
− 1δ∆
Atwood number is expressed as a function of the mean den- q(1+δ)∆) 1+ q∆ (13)
sity which is obtained from a fit to the numerical data from h sgs i∝ 3 ∆ h sgsi
!
Timmes&Woosley(1992):
Becausethisscalinglawisastatisticalruleonecannotex-
1 0.145 0.0100 pectthatitholdslocally.However,thedynamicalprocedurefor
At(ρ)= 0.0522+ (10)
2" √ρ9 − ρ9 # the calculationof the eddy-viscosityalso can beutilised fora
4 W.Schmidtetal.:AlocalisedsubgridscalemodelforfluiddynamicalsimulationsII
localisedrescalinglaw.Let∆ bethegridresolutionattimet, set method in the passive implementation (Osher&Sethian
t
and ∆ the slightly smaller resolution prior to the last time 1988; Reineckeetal. 1999). The implemented equation of
t δt
−
step.ApplyingthetestfilterintroducedinpaperI,Sect.4.1,at state for electron-degenerate matter is described in Reinecke
timet,theturbulentvelocityq associatedwithvelocityfluctu- (2001),Sect.3.2.Thermonuclearburningismodelledbysim-
T
ationsonlengthscalesgreaterthanβ∆ andsmallerthanγ β∆ ple representative reactions (Reineckeetal. 2002): 12C and
t T t
isobtained.Hereβ = ∆ /∆istheconstantratiooftheeffec- 16C is fused to 56Ni and α-particles at densities higher than
eff
tivecutofflengthtothesizeofthegridcells.SincetheRiemann 5.25 107gcm 3,whereas24Mgisproducedatlowerdensities
−
·
solverdoesnotaccountforthefractionalgrowthoftheturbu- in the late stage of the explosion. Finally, all reactions cease
lenceenergyduetotheshiftofthecutoff,theresultofadvanc- below 107gcm 3. This threshold presumably marks the tran-
−
ing the dynamicalequationsfrom t δt to t is q∆t δt. Now an sitionfromtheflamelettothebrokenreactionzoneregimeof
− sg−s
estimate of q∆t can be made locally via interpolation of the turbulentburning(seeNiemeyer&Kerstein1997).Thecorrect
sgs
turbulence energyin length scale space. Using the contracted treatmentoftheburningprocessinthisregimeisstillamatter
Germanoidentity (paperI, Sect. 4.1),the total turbulenceen- ofdebateand,forthisreason,notincludedinthepresentsim-
ergyassociated with the test filter lengthat time t is givenby ulations.
As initial model, we choose a white dwarf of mass M =
ρk∆t 1.4M composedof equalmass fractionsof carbonand oxy-
ktγuTr∆bt = h ρs(Tgs)iT +k(T). (14) genw⊙ithacentraldensityρc =2.0 109gcm−3andtemperature
·
T =7.55 108K.AssuggestedbyWunsch&Woosley(2004),
The unknown is the rescaled SGS turbulence energy ks∆gts. thceradialt·emperatureprofileisgivenbyaparabolawithacut-
Neglecting compressibility effects and variations on sub-test offatthethermalradiusΛ=7.35 107cm:
filter lengths, we set kγT∆t k∆t + k(T). Then linear interpo- ·
turb ≃ sgs r 2
lationbetween∆t−δt andγT∆t yields: T(r)=Tc 1− Λ θ(r−Λ)+T0θ(Λ−r), (19)
f " (cid:18) (cid:19) #
ks∆gts ≃ks∆gt−sδt + 1 fkT, (15) where θ denotes the Heaviside step function. The thermal ra-
− dius specifies the size of the convectivecore prior to the run-
wheretheinterpolatingfactor f isgivenby away. At larger radii, the matter is isothermal with T = 5
0
·
105K. In the centre, we set an axisymmetric initial burning
∆ ∆
f = t− t−δt . (16) regionwithsinusoidalperturbations(seeRo¨pke&Hillebrandt
γ ∆ ∆
T t− t−δt 2005). In order to achieve higher resolution, only single oc-
DuetothesmallnessofaCFLtimestep,thefractionalchanges tantssubjecttoreflectingboundaryconditionswereevolvedin
ofthecutofflengthwillbesmall.Hence, f 1. thesimulationsdiscussedhere.
≪
Theproblemwiththerescalinglaw(15)isthatitfailstoac- Moreover,we applied the co-expandinggrid technique of
countforthecorrectasymptoticbehaviourinthelimitoffully Ro¨pke(2005).Thereby,itispossibletomaintainanequidistant
developedturbulence.OnaccountoftheGermanoidentity(pa- grid geometry over the whole domain of turbulent burning at
perI,Sect.4.1),onewouldexpectthestatisticalrelation anystageoftheexplosion,evenwhentheejectahaveexpanded
byalargefactorcomparedtotheinitialsizeofthewhitedwarf.
(γ2/3 1) k = k(T) (17)
T − h sgsi h i Recently, Ro¨pkeetal. (2005) have combined this technique
withthegridgeometryusedbyReineckeetal.(2002)inorder
for regions of nearly homogeneous turbulence obeying
tocapturetheburningprocessintheinteriorwithoptimalres-
Kolmogorov scaling. This relation is asymptotically repro-
olution,whileusingacoarsergridwithexponentiallyincreas-
ducedbytherescalingmodifiedlaw
ingcellsoutside.Thehybridgrid,evenatmoderateresolution,
γ2/3(1 f) f enables us to resolve details of the turbulent flame dynamics
ks∆gts ≃ Tγ2/3 −f ks∆gt−sδt + γ2/3 fkT, (18) which used to be inaccessible for non-adaptive schemes. All
T − T − numericalsimulationspresentedinthisarticlefeatureahybrid
which results from the interpolation of the turbulence energy grid.
dividedbythe associated lengthto the power2/3.Thisis the The non-uniform grid geometry poses certain difficulties
rescalinglawwhichwasimplementedforthenumericalsimu- when applying the localised SGS model. In Sect. 2.2, we
lationsofthermonuclearsupernovaexplosionsdiscussedinthe showedthatitisrelativelyeasytoaccountforthevariationin
nextsection. thetimedomainduetheco-expansionofthegrid.Inthecaseof
non-uniformgrids,however,thefilteroperationdoesnotcom-
mutewithspatialderivativesinthedynamicalequations.Apart
3. Numericalsimulations
fromthat,theweighingofnodesforthediscretefilteringpro-
In the following, we will present results from several cedure becomes dependent on the location. This would lead
numerical simulations using the methodology outlined in to substantialcomplicationsin the numericalimplementation.
Ro¨pke&Hillebrandt (2005). In essence, the piece-wise Fortunately, we found a simple solution: Since the turbulent
parabolic method is used to solve the hydrodynamical equa- dynamics mostly takes place in the burning region which is
tions (Colella&Woodward 1984; Fryxelletal. 1989) and the containedwithintheuniformpartofthegrid,wecomputedthe
evolutionoftheflame-frontsiscomputedbymeansofthelevel eddy-viscosityparameterC onlyinthisregionandsettherate
ν
W.Schmidtetal.:AlocalisedsubgridscalemodelforfluiddynamicalsimulationsII 5
Fig.1.Timeevolutionofintegratedquantitiesforthreesimulationswithidenticalinitialconditionsandresolution2563.Inone
case,Clement’sSGSmodelwithwallproximityfunctions(WPF)wasused.Fortheothertwosimulationsweappliedthelocalised
SGSmodelwithArchimedianproduction,respectively,switchedoffandon.
ofturbulentenergytransferequaltozerointheexterior.Later dynamicalprocedurewhichimplieslessenergytransferiftur-
inthissection,wewilldemonstratethatneglectingtheenergy bulence is still developing but enhanced transfer in the fully
transferoutsidetheburningregioncanbejustifiedaposteriori. developedcase(Schmidtetal.2005a).
TheoriginalSGSturbulenceenergymodelimplementedby
Thisisclearlyreflectedinthetimeevolutionoftheturbu-
Niemeyer&Hillebrandt (1995) is based on statistical closure
lence energyin two simulationswhich differ only in the SGS
parameters.Clement(1993)suggestedtosetC =0.1W,where
ν model.Thegraphsofthemass-integratedSGSturbulenceen-
W isanempiricalwallproximityfunction(WPF),
ergy are plotted in the top panel on the left of Fig. 1. There
e are two variants of the localised SGS model, one including
W =min 100,max 0.1, 10 4 int . (20)
− Archimedian production as described in Sect. 2.1, whereas it
· k
" sgs!# isassumedthatenergytransferthroughthecascadeistheonly
Sincee c2,theratioe /q2 isbasicallyaninverseMach source of turbulence production in the alternative model. In
int ∼ s int sgs
numbersquared.Ifthereislittleturbulenceenergy,C iscon- contrast to Clement’s model, there is initially very little SGS
ν
siderably enhanced. On the other hand, C becomes smaller turbulenceenergyfollowedbyamuchsteeperriseinthesimu-
ν
than0.1if theSGSturbulencevelocityexceedsafewpercent lationswiththelocalisedSGSmodel.Therapidgrowthoftur-
of the speed of sound. This behaviour of the eddy-viscosity bulence energy can be attributed to the substantially stronger
parameter is qualitatively different from the prediction of the turbulenceproductionwithinthetimeintervalfrom0.3to1.0s
6 W.Schmidtetal.:AlocalisedsubgridscalemodelforfluiddynamicalsimulationsII
Fig.2.TimeevolutionofintegratedSGSquantitiesforaseriesofsimulationswithvaryingresolution.
(seerightbottompanelinFig.1).Inthesecondhalfofthecom- explosionbutturbulentenergytransferfromresolvedscalesis
bustion phase, the total SGS turbulence energy is almost one neverthelesstheprimarysourceofSGSturbulenceproduction.
order of a magnitude larger which enhances the flame propa-
Accordingtothescalingargumentmentionedintheintro-
gation speed accordingly. At later times, the discrepancy be-
duction,buoyancyshouldbecomeevenlessimportantrelative
comesevenmorepronouncedbecausetherescalingofk im-
sgs totheturbulentenergytransferwithincreasingresolution.This
plementedinthelocalisedmodelfeedskineticenergyintothe
is indeed observed in a series of simulations with the reso-
subgridscalesagainsttheactionofSGSdissipation.Thenetre-
lution varying between 1283 up to 3843 grid cells. The time
sultisasignificantenhancementoftheexplosionenergyanda
evolutionoftheintegratedrateofenergytransferandspecific
largeryieldofburningproducts(seebottompanelontheleftof
Archimedian force, respectively, is plotted in the top panels
Fig.1).AlsonotethattheadditionalproductionofSGSturbu-
of Fig. 2. In order to interpretthese graphs, it is importantto
lencebytheArchimedianforceterminequation(4)increases
note that the SGS energy transfer is expected to become sta-
the explosion energyeven further. Compared to the reference
tisticallyscaleinvariantinthecaseofKolmogorovturbulence.
simulation with Clement’s model, the final kinetic energy of
This follows from the scaling law v(l) l1/3 for the turbu-
0.472 1051erg is about25% greater. However,this does not ′ ∝
imply·thatArchimedianproductiondominatesovertheturbu- lent velocity fluctuations. Hence, ksgs ∝ ∆2/3. Since the char-
acteristic time scale of turbulent velocity fluctuations scales
lencecascade.Infact,theevolutionoftheSGSturbulenceen-
with l2/3, it follows that the average time derivative of k is
ergydiffers onlylittle, as one can see from the plotin Fig. 1. sgs
scale invariant. The rate of energy transfer per unit mass, on
In conclusion, SGS buoyancy effects appear to influence the theotherhand,isproportionalto∆k1/2S 2.Thisexpressionis
sgs| ∗|
W.Schmidtetal.:AlocalisedsubgridscalemodelforfluiddynamicalsimulationsII 7
also scale-invariantbecause the rate-of-strainscalar S mea-
∗
| |
sures the inverse time scale of the smallest resolved velocity
fluctuations,i.e. S 2 ∆ 4/3, while ∆k1/2 ∆4/3. The com-
| ∗| ∝ − sgs ∝
putedrateofenergytransferplottedinthetoppanelontheleft
ofFig.2exhibitpeakvalueswhicharewithinthesameorderof
magnitude,althoughthe energytransferseemsto be underes-
timatedforthelowestresolutions.Initially,theturbulenceen-
ergyrisesexponentiallyataratewhichchangesonlylittlewith
resolution(rightbottompanelinFig.2).Thisisreflectedinthe
nearlycoincidinggraphsoftherateofSGSenergytransferup
tot 0.3sshownintheleftpanelonthetopofFig.2.
≈
The SGS turbulence energy in the regime of fully devel-
opedturbulencedecreasesforhigherresolution.Thistrendcan
be discerned particularly in the post-burning phase, in which
no further energy is injected and the turbulent flow begins to
decay. From the initial production phase to the post-burning
phase,however,arathercomplicatedbehaviourbecomesman-
ifest asthe resultof the interplaybetweenturbulenceproduc-
tion by the strain of the resolved flow, the SGS Archimedian
forceandthegridexpansion.Theexplosionenergeticsplotted Fig.3.Timeevolutionofthetotalenergyforthesamesimula-
inFig.3showsthefollowingbehaviourdependingonthenu- tionsasinFig.2
mericalresolution (also see Table 1): Whereas the final value
ofthetotalenergyisaboutthesameforthelowerresolutions,
there is a significantly enhanced yield of energy in the case Table 1. Total release of nuclear and kinetic energyand total
N = 3843. However, this does not imply that the model fails massesofirongroup(“Ni”)andintermediatemass(“Mg”)el-
toconvergewithincreasingresolution.Turbulentflowregions ementscorrespondingtoFig.3.
are confined in a fraction of the numerical grid, whereas the
greaterpartofthe gridisoverheadrequiredformodellingthe
N E [1051erg] E [1051erg] M /M M /M
non-turbulentouterpartsofthe expandingstar andsome por- nuc kin Ni ⊙ Mg ⊙
1283 0.963 0.433 0.523 0.175
tion of the surroundingquasi-vacuum.In the case N = 2563,
1923 0.970 0.442 0.529 0.172
for example, significant SGS turbulence production occurs in
2563 1.000 0.472 0.548 0.172
theinner1003cellsatt=0.5s.Thisisthetimeofmaximalin-
3843 1.087 0.560 0.586 0.206
tegratedturbulenceproduction.However,inpaperIwedemon-
stratedthat100cellsineachspatialdimensionisdefinitelynot
sufficient to resolve developed turbulent flow sufficiently far
down toward the inertial subrange using PPM. Although this DetailsoftheSGSdynamicsareillustratedbycontourplots
is merely a crude estimate, it appears plausible that even the of two-dimensional spatial sections from the simulation with
supernovasimulationwithN =3843resolvestheturbulentdy- N =3843gridcells(Figs.4,5,and6).IneachFig.,thefollow-
namics only marginally. As a further indication, the plateau- ingdynamicaltermsofequation(4)areplotted:
likeflatteningoftherateofproductionanddissipation,respec-
1. Rateofproductioncausedbystrain,ℓ S 2(lefttoppanel).
tively, can be seen in the left panels of Fig. 2 for the highest ν| ∗|
2. SpecificArchimedianforce0.1g (righttoppanel).
resolution only. We interpret the flattening as a consequence eff
3. Rateofdissipation 7q d q2 /ℓ (leftbottompanel).
oflocalstatisticalequilibriumbetweenproductionanddissipa- −30 sgs − sgs ǫ
tion. 4. Rateofdiffusion ρ1∇· ρℓκqsgs∇qsgs −ℓκ|∇qsgs|2(rightbot-
tompanel).
Unfortunately, we were not able to perform a run of still (cid:16) (cid:17)
higher resolution due to the limitations of our computational The flame surface as given by the zero level set is indicated
resources. In any case, we expectthat N = 5123 grid cells in by the contours in white. Note that these quantities have the
oneoctantwouldbesufficientforanaccuratemodellingofthe dimensionof acceleration.Fig. 4 shows the typicalRayleigh-
turbulent dynamics, whereas simulations with less resolution Taylormushroomshapeswhichhaveformedoutoftheinitial
canbeutilisedtodiscerntrendsinparameterstudies.Asimilar sinusoidal perturbations at time t = 0.3s. Significant energy
conclusionwas drawnby Ro¨pke(2005) froma series oftwo- transferisconcentratedinsmallregionsandthereislittledis-
dimensional simulations, in which a pronouncedjump of the sipationyet.ComparingtoFig. 2, onecansee thatturbulence
totalenergywasfoundbetweentheN =2562andthe5122run, productionisjustabouttorise.Att=0.45s,therateofenergy
respectively,whilemoreorlessthesameenergywasobtained transferhasreacheditsmaximumandisspreadalloverthein-
for N 5122. Moreover, snapshots of the zero level set for terioroftheflames(seeFig.5).TheaccelerationofSGSfluid
≥
varyingresolutionsuggestedthatsecondaryKelvin-Helmholtz parcelssubjecttothelargeststrainexceeds106timesthegravi-
instabilitiesarebarelyornotatallresolvedwithN 2562. tationalaccelerationonEarthrelativetotheresolvedflow.The
≤
8 W.Schmidtetal.:AlocalisedsubgridscalemodelforfluiddynamicalsimulationsII
Fig.4.ContoursectionsshowingthecontributionstotheevolutionoftheSGSturbulencevelocityq givenbyequation(4)at
sgs
t = 0.3s.OnlytheinnerregionofthegridwithN = 3843 cellsisshown.Thewhitecontoursrepresentthesectionsthroughthe
flamesurface.
SGSbuoyancyistypicallybyanorderofamagnitudesmaller. neticenergythroughtheturbulencecascadeanddiffusionacts
Both dissipation and transportdue to SGS turbulentdiffusion toredistributethisenergyintoregionswithlittleproduction.
are comparable to the rate of energy transfer at this time. In
theunburnedmaterialoutside,ontheotherhand,thereisvirtu- Fig. 10 shows the evolution of the SGS turbulence veloc-
allynoSGSturbulence.Thereby,itisconfirmedthatswitching ityqsgs attheflamefrontsinthree-dimensionalvisualisations.
offtheenergytransfertermsinthenon-uniformgirdregionsat Thegridlinesroughlyindicatetheuniformpartofthenumer-
sufficientdistancefromtheflamefrontsisareasonablesimpli- icalgrid.Thecorrespondingabsolutescaleisindicatedbythe
fication.Obviously,theflowishighlyanisotropicinthevicin- sizeXuniandthecorrespondingnumberofcellsNuni.Inthefirst
ity of the flames which highlightsthe necessity of a localised threesnapshotsonecansee thegrowthoftheinitialperturba-
SGSmodel.InFig.6,onecanseethatturbulentenergytransfer tions. The axial symmetry is gradually broken by the forma-
is declining and becoming small relative to the Archimedian tionofsecondaryinstabilities.Att=0.45sthesmallerplumes
force near the flame front. However, this does not imply that originatingfrom these instabilities are highly turbulent. From
the amount of SGS turbulence and, thus, the turbulent flame t = 0.6sonwards,the system increasinglylooses its memory
speedisdominatedbythe SGSbuoyancybecausethe bulkof oftheinitialconditionandtheturbulenceintensityattheflame
SGS turbulence energy has been produced by transfer of ki- surfaceis abatingandlevelling.Thelast snapshotatt = 1.5s
showsa complexstructurewith featuresovera wide rangeof
W.Schmidtetal.:AlocalisedsubgridscalemodelforfluiddynamicalsimulationsII 9
Fig.5.ThesameplotasinFig.4att=0.45s.
scales.Thereappeartobefiveorsixmajormodeswhicheven- Note that integrating each PDF over the decade logarithm of
tuallyprevail.Theresultinglayeringofnuclearspeciesisillus- the speed yields unity. Also shown are the PDFs of q and
sgs
tratedbythecontourplotsofthecorrespondingmassdensities v (∆ ).Therelationbetweentheturbulentflamespeeds and
RT eff t
inFig.7.Nickelisconcentratedinthecentralregion,whereas theSGSturbulencevelocityq isformulatedinequation(2).
sgs
both magnesium and unburned carbon and oxygen are found During the first tenth of a second, s is basically givenby the
t
further outside. The outermost layers and the narrow down- laminarflamespeed.Thentheflamepropagationbecomesin-
draftsbetween the convectivefingers of nuclear ash are com- creasingly affected by SGS turbulence. From about 0.3s on-
posedalmostexclusivelyofcarbonandoxygen.Thestratifica- wards, s is dominatedbyq . Atlater times,onecansee the
t sgs
tion of the nuclear species in the explosionejecta is reflected asymptotic relation s 2q /√3. The Rayleigh-Taylor ve-
t sgs
≃
inthecorrespondingmassdensityfunctionsdM/dvr,wherevr locityscalevRT(∆eff)isinitiallymuchsmallerthanthelaminar
is the radial velocity component. In particular, Fig. 8 shows burning velocity. As the flame propagates outwards, both the
that little carbon and oxygen is found for velocities less than gravity and the density contrast at the flame surface become
3000kms−1inthelatephaseofalmosthomologousexpansion. largerandvRT(∆eff)increases.Eventually,thePDFofvRT(∆eff)
tendstowardarathernarrowpeakaround107cms 1.ThePDF
−
Tounderstandtheflamedynamics,itisinstructivetocon-
of q , on the other hand, extends over a substantially wider
sgs
sidertheprobabilitydensityfunction(PDF)ofthelogarithmof
range.Forthisreason,thelocalisedSGSmodelgeneratesmore
the flame propagation speed s over the surface of the flame.
t variationinthepropagationoftheflamefrontincomparisonto
The PDFs for several instants of time are plotted in Fig. 9.
10 W.Schmidtetal.:AlocalisedsubgridscalemodelforfluiddynamicalsimulationsII
Fig.6.ThesameplotasinFig.4att=0.75s.
theSharp-Wheelermodel.ThisisexpectedbecausetheSharp- in the interior ash regions.Consequently,there is pronounced
Wheeler relation is ignorant of the interaction between sub- anisotropy at the flame surface which can be tackled by the
gridandresolvedscalesandtheeffectsofnon-localtransport. localisedSGSmodelonly.TheArchimedianforcecontributes
Nevertheless, v (∆ ) is seemingly a velocity scale which is noticeably to the turbulent flame speed, particularly once the
RT eff
representative for the magnitude of q at the flame surface flame surfacehas grownsubstantially.However,the dominat-
sgs
duringmostoftheburningprocess. ingeffectistheenergytransferthroughtheturbulencecascade.
Inthelate stageoftheexplosion,sustainedturbulenceenergy
4. Conclusion comesfromtherescaling,whilethemajordynamicalcontribu-
tionisSGSdissipation.Furthermore,itappearsthatnumerical
WeappliedtheSGSturbulenceenergymodeltothelargeeddy gridswith morethan N = 2563 cells in one octantare neces-
simulation of turbulent deflagration in thermonuclear super- sary in order to sufficiently resolve the turbulent dynamicsin
nova explosions. The novel features of this model are a lo-
theburningregionsandtoobtainconvergedresults.
calised closure for the rate of energy transfer, an additional
Archimedian force term which accounts for buoyancy effects An investigation of probability density functions over the
on unresolvedscales and the rescaling of the SGS turbulence flame fronts (see Fig. 9) reveals that the Rayleigh-Taylor ve-
energyduetotheshiftofthecutofflengthinsimulationswith locity scale v (∆ ) given by the Sharp-Wheelerrelation (1)
RT eff
a co-expanding grid. We found that the production of turbu- is not negligible compared to the SGS turbulence velocity
lenceislargelyconfinedtotheregionsneartheflamefrontsand q , once the regime of fully turbulent burning has been en-
sgs