Anelastic Versus Fully Compressible Turbulent Rayleigh-B\'enard Convection

1.9 MB·English

by  Jan Verhoeven

ANELASTIC VERSUS FULLY COMPRESSIBLE TURBULENT
´
RAYLEIGH-BENARD CONVECTION
Jan Verhoeven, Thomas Wieseh¨ofer and Stephan Stellmach
5 Institut fu¨r Geophysik, Westf¨alische Wilhelms Universit¨at Mu¨nster, Germany
1 [email protected]
0
2

r
a
M
ABSTRACT

3 Numerical simulations of turbulent Rayleigh-B´enard convection in an ideal gas, using either
1 the anelastic approximation or the fully compressible equations, are compared. Theoretically,

the anelastic approximation is expected to hold in weakly superadiabatic systems with (cid:15) =
]
R ∆T/T (cid:28) 1, where ∆T denotes the superadiabatic temperature drop over the convective layer
r
S and Tr the bottom temperature. Using direct numerical simulations, a systematic comparison
. of anelastic and fully compressible convection is carried out. With decreasing superadiabaticity
h
(cid:15), the fully compressible results are found to converge linearly to the anelastic solution with
p
- larger density contrasts generally improving the match. We conclude that in many solar and
o
planetary applications, where the superadiabaticity is expected to be vanishingly small, results
r
t obtained with the anelastic approximation are in fact more accurate than fully compressible
s
a computations, which typically fail to reach small (cid:15) for numerical reasons. On the other hand, if
[ the astrophysical system studied contains (cid:15) ∼ O(1) regions, such as the solar photosphere, fully

compressible simulations have the advantage of capturing the full physics. Interestingly, even in
2
v weaklysuperadiabaticregions,likethebulkofthesolarconvectionzone,theerrorsintroducedby
7 usingartificiallylargevaluesfor(cid:15)forefficiencyreasonsremainmoderate. Ifquantitativeerrorsof
3 the order of 10% are acceptable in such low (cid:15) regions, our work suggests that fully compressible
2
simulations can indeed be computationally more efficient than their anelastic counterparts.
1
0 Subject headings: convection — Earth — planets and satellites: gaseous planets — Sun: interior —
.
1 Turbulence
0
5 1. INTRODUCTION The convective regions in stellar and planetary
1 objects typically feature a non-negligible density
: Thermalconvectionisofprimaryimportancein stratification, and the flows are often subsonic. In
v
astrophysicalobjects. Itcarriestheheatflowover
i this paper, we compare two approaches that are
X largeregionsinstellarandplanetaryinteriors,and
commonly used for modeling convection in these
r is one of the major sources of mechanical mixing systems numerically—the fully compressible ap-
a
in these objects. Furthermore, some of the most
proach and the so-called anelastic approximation.
striking, large-scale features of stars and planets
Ourgoalistoquantifytheaccuracyandefficiency
are powered by convective motions, such as in-
ofbothmethodsinagivensituation,guidingmod-
trinsic dynamo-generated magnetic fields, plate
elersinmakingtherightchoicefortheirparticular
tectonics on Earth and possibly also the zonal
problem at hand.
winds on Jupiter and other giant planets (e.g.
The fully compressible equations are the most
Brun et al. 2004; Brandenburg and Subramanian
fundamental equations governing thermal convec-
2005; Trompert and Hansen 1998; Tackley 2000;
tion. They can directly be derived from first prin-
Bercovici 2003; Heimpel et al. 2005; Verhoeven
ciples of physics, such as mass, energy, and mo-
and Stellmach 2014).
mentum conservation, equipped with constitutive
1
relationsthatcharacterizethefluid. Theresulting to study regions where the Mach number is not
equationsarethusverygeneralandencompassthe small.
full range of physical behavior, from the temporal Among the sound-proof approaches described
evolution of the convective motions to the prop- above,theanelasticapproximationistheonemost
agation of sound waves. On the one hand, this commonlydeployedformodelingstellarandplan-
allows to study regions such as the outermost lay- etaryinteriors(e.g. GlatzmaierandRoberts1996;
ers of the Sun, where the Mach number, i.e. the Miesch et al. 2008; Brun et al. 2011; Jones et al.
ratio of convective velocity to the sound speed, 2011). The anelastic equations are theoretically
becomes O(1). On the other hand, problems arise predicted to hold for low Mach number systems
in low Mach number regions where the flow ve- inwhichonlyslightthermodynamicperturbations
locities are much slower than the sound speed, from a hydrostatic background state occur (e.g.
which is typically the case in the bulk of deep Gough 1969). In convective systems, the back-
stellar and planetary interiors. Even though the ground state is typically assumed to be adiabatic.
convective motions in such regions occur on time The above conditions are believed to be satisfied
scales which are many orders of magnitude larger in the deep interiors of giant planets and in the
than the acoustic time scale, standard numerical bulk of the solar convection zone, but break down
schemeshavetoexplicitlyresolvethesoundwaves in their outermost parts which feature relatively
for stability reasons. This forces modelers to as- small sound speeds (Ulrich 1970; Bahcall and Ul-
sume artificially large Mach numbers, which re- rich1988;Christensen-Dalsgaardetal.1996;Guil-
duces the differences between the convective and lotetal.2004). Thedynamicsoftheseouterlayers
acoustic time scales to numerically tractable val- thus cannot be accounted for within the anelas-
ues (e.g. Tobias et al. 1998; Brummell et al. 2002; ticframework, andmodelersareforcedtoexclude
K¨apyl¨a et al. 2010). Errors introduced by this themfromthesimulationdomain. Thedynamical
procedure occur as an unavoidable side-effect in consequences of neglecting these regions remain
the fully compressible framework. Still, most of unclear.
the numerical resources typically go into captur-
In summary, both approaches have advantages
ing acoustic wave propagation phenomena, which
and drawbacks. While the fully compressible ap-
are generally believed to be irrelevant for the in-
proach is the method of choice for modeling O(1)
vestigated convection dynamics (but see Bogdan
Machnumberflowsinnear-surfaceregionsofstel-
et al. 1993; Meakin and Arnett 2006).
lar objects, the anelastic approximation seems to
To circumvent the problems arising from the be beneficial in the deep interiors where the Mach
numerical stiffness of the fully compressible equa- numbers are usually very small and where the
tions, different ”sound-proof” models, such as the thermodynamic state is close to the adiabat. Un-
low Mach number approach (e.g. Majda and fortunately, in many astrophysical applications, it
Sethian 1985; Bell et al. 2004; Almgren et al. remains unclear which approach performs best,
2006), the pseudo-incompressible approximation with anelastic and fully compressible models be-
(Durran 1989) or the anelastic approximation ing used side by side. The main goal of this study
(Batchelor 1953; Ogura and Phillips 1962; Gough is thus twofold: First, we aim to quantify and
1969; Gilman and Glatzmaier 1981; Lantz and comparetheerrorsinherentinmodelingturbulent
Fan 1999) have been developed. Instead of pre- convection in both approaches. Secondly, we seek
scribing artificially large Mach numbers, all these tocomparetheircomputationalefficiency,thereby
approachestaketheoppositeroutebyconsidering guiding modelers in minimizing the tradeoff be-
thesmallMachnumberlimitofthefullycompress- tween accuracy and efficiency for any given situa-
ible equations. The same time scale disparities tion.
which make solving the fully compressible equa-
Perhaps somewhat surprisingly, comparing re-
tions numerically challenging are thus exploited
sults from the anelastic models currently used in
to substantially simplify the equations. As a re-
astro- and geophysics to standard fully compress-
sult, the pressure field adapts instantaneously,
ible simulations is non-trivial. This is because
whicheffectivelyfiltersoutthesoundwaves. This
theanelasticmodelsusuallyparameterizethetur-
comes, however, at the price of loosing the ability
bulent, subgrid-scale entropy flux, while similar
2
turbulence models are typically not used in fully
compressible models. The popularity of turbu-
lence modeling in the anelastic framework stems
from the fact that it allows further simplifications
of the governing equations, which eases the nu-
g
merical implementation considerably. Typically,
molecular heat conduction is neglected and re-
zˆ
placed by an artificial eddy diffusion model that
yˆ
represents turbulent mixing of entropy (Gilman
and Glatzmaier 1981; Glatzmaier 1984; Bragin-
sky and Roberts 1995; Lantz and Fan 1999). This xˆ
turbulententropydiffusionmodel,however,isnot
Fig. 1.— Compressible convection is modeled in
mandatoryfortheactualanelasticapproximation,
Rayleigh-B´enardgeometry,i.e. inaCartesianbox
andanelasticequationshavebeenformulatedthat
that is cooled from above and heated from below.
do not rely on parameterizations of the subgrid-
Gravity g is pointing downward, antiparallel to
scale transport (Gough 1969). These equations
the z-axis.
have not found widespread use so far. In order to
providedirectcomparabilitybetweentheanelastic
andthefullycompressibleapproach, inthisstudy comparability of the anelastic and fully compress-
wewillrestrictourselvestomolecularthermalheat ible influences.
conduction in both cases.
In this paper, we present the first systematic
While direct comparisons of anelastic and fully one-to-one comparisons between fully compress-
compressible gravity wave dynamics in stably ibleandanelasticnumericalsimulationsofconvec-
stratified set-ups have been performed in several tion in the fully nonlinear, turbulent regime. As
studies (e.g. Davies et al. 2003; Klein et al. 2010; a starting point, we neglect important ingredients
Brown et al. 2012), the unstable thermal convec- of stellar convection, such as spherical geometry,
tion case considered in this paper has received rotation, compositional inhomogeneities, nuclear
less attention so far. The work of Berkoff et al. reactions, magnetic fields, penetration and over-
(2010) focussed on linear magnetoconvection and shooting in stably stratified layers, and the corre-
found good agreement between both approaches sponding wave-emission. This allows us to quan-
for the weakly superadiabatic case. Subsequently, tify the respective errors, as well as the compu-
Lecoanet et al. (2014) studied differences between tationalefficiencyencounteredinbothapproaches
temperature and entropy diffusion, while Calkins in the simplest setup possible. The influences of
et al. (2014, 2015) focussed on the influence of the above physical processes will be investigated
rotation on the onset of anelastic and fully com- in future studies.
pressible convection. Their linear study identified
The paper is organized as follows: In section 2,
shortcomingsoftheanelasticequationsforrapidly
we start with defining our idealized model, which
rotating, low Prandtl number fluids, where fast
is followed by discussing the fully compressible
density oscillations were found to become non-
equations along with the anelastic approximation
negligible. Calkins et al. (2014) conclude that
in section 3. A brief overview of the applied nu-
fully non-linear studies tracing the validity range
merical methods is given in section 4, while a di-
of the anelastic approximation are crucial in both
rectcomparisonofanelasticandfullycompressible
rotating and non-rotating systems, especially in
results and the computational efficiencies of both
the turbulent regime characterized by a broad-
approachesarediscussedinsection5. Finally,gen-
band frequency spectrum. A first step in this
eral conclusions are drawn in section 6.
direction has been taken by Meakin and Arnett
(2007), who compared non-linear anelastic and
2. MODEL
fully compressible simulations of stellar oxygen
burning. Thedifferingphysicalprocessesincluded Fully compressible and anelastic convection in
in each model, however, precluded a one-to-one an ideal gas are modeled in a plane fluid layer of
3
depth d confined between rigid, horizontal plates, way
as displayed in figure 1. Gravity g = −gzˆ is con-
stant and pointing downward, antiparallel to the cvρ[∂tT +(v·∇)T]+p(∇·v)
unit vector zˆ. The fluid is cooled from above and =c ρ[∂ T +(v·∇)T]−[∂ p+(v·∇)p] (5)
p t t
heatedfrombelowbymaintainingaconstant,pre-
scribed temperature difference across the layer. byusingequations(1)and(4)(detailsaregivenin
The remaining boundary conditions are periodic appendix A). As usual, the thermodynamic quan-
inthehorizontaldirectionsandnoslipatthebot- tities are decomposed into a time-independent,
tom and the top boundary. The ideal gas is char- verticallyvarying,hydrostaticandadiabaticback-
acterized by constant dynamic viscosity µ = νρ, groundstate(indexA)1 andasuperadiabaticpart
heat conductivity k = c ρκ and specific heat ca- (index S), which is allowed to vary in time and
p
pacities at fixed volume and pressure, c and c . space,
v p
The quantities ν and κ are the kinematic viscos-
T(t,x)=T (z)+T (t,x), (6)
ityandthethermaldiffusivity,respectively,which A S
consistently vary across the fluid layer inversely ρ(t,x)=ρA(z)+ρS(t,x), (7)
proportional to the density. p(t,x)=p (z)+p (t,x). (8)
A S
The governing equations for fully compressible
convection describing the temporal evolution of While for subadiabatic or stably stratified fluids a
density ρ, temperature T, pressure p and veloc- conductive background state is the better choice,
ity v are the assumption of approximate adiabaticity (i.e.
isentropy) is justified in superadiabatic regions,
∂ ρ+∇·(ρv)=0, (1) where convection turbulently mixes the fluid and
t
thus homogenizes entropy. The background pro-
filecanbederivedfromhydrostaticity∇p=−ρgzˆ
ρ[∂ v+(v·∇)v]=−∇p−ρgzˆ
t (i.e. equation (2) with v = 0) and the thermody-
(cid:20) (cid:21)
+µ ∇2v+ 1∇(∇·v) , (2) namic relation
3
ρTds=c ρdT −δ dp, (9)
p p
c ρ[∂ T +(v·∇)T]+p(∇·v)=
v t with s denoting specific entropy and the dimen-
(cid:20) 1 (cid:21)2 sionless thermal expansion coefficient being de-
k∇2T +2µ e − (∇·v)δ , (3)
ij 3 ij fined as δp = −(∂lnρ/∂lnT). Note that for an
ideal gas, δ = 1, see (4). By further assuming
p
adiabaticity (i.e. ds=0), the background state is
p=(cp−cv)ρT, (4) found to be
(cid:18) (cid:19)
g
winigththtedsetnraoitninrgatteimteenasnodr.eEijq=uat12io(∂njsv(i1+-3∂)ievxjp)rbeses- TA(z)=Tr 1− cpTrz , (10)
the conservation of mass, momentum and energy, (cid:18) g (cid:19)cv/(cp−cv)
respectively,whileequation(4)istheidealgaslaw. ρA(z)=ρr 1− c T z , (11)
p r
3. FULLY COMPRESSIBLE AND AN- (cid:18) g (cid:19)cp/(cp−cv)
p (z)=(c −c )ρ T 1− z ,
ELASTIC EQUATIONS A p v r r cpTr
(12)
In the following, the anelastic and fully com-
pressible equations are discussed in detail. where the index r in T and ρ refers to reference
r r
values, here defined as the adiabatic values at the
3.1. ReformulationandNon-dimensionali-
zation 1We use the term ”adiabatic” here for constant entropy
states. More accurately, the background state may be
We begin with reformulating the left-hand-side
called”isentropic”,whichhoweverappearstobelesscom-
of equation (3) in the more ”anelastic-friendly” monintheliterature.
4
bottom boundary. Applying the decomposition of peradiabatic pressure p extracts kinetic energy
S
the thermodynamic variables (6-8) to equations from the vertical flows todrive the horizontal mo-
(1-4) results in tions(seee.g. Gough1969). Thenon-dimensional
thermodynamic quantities2 then read
∂ (ρ +ρ )+∇·[(ρ +ρ )v]=0, (13)
t A S A S
T(t,x)=T (z)+(cid:15)T (t,x), (17)
A S
ρ(t,x)=ρ (z)+(cid:15)ρ (t,x), (18)
(ρ +ρ )[∂ v+(v·∇)v]=−∇(p +p ) A S
A S t A S
(cid:20) 1 (cid:21) p(t,x)=pA(z)+(cid:15)pS(t,x), (19)
−(ρ +ρ )gzˆ+µ ∇2v+ ∇(∇·v) , (14)
A S 3
with (cid:15) = ∆T/T and the adiabatic background
r
state being
c (ρ +ρ )[∂ (T +T )+(v·∇)(T +T )]
p A S t A S A S
T (z)=(1−Dz), (20)
A
−[∂ (p +p )+(v·∇)(p +p )]=
t A S A S
ρ (z)=(1−Dz)1/(γ−1), (21)
(cid:20) 1 (cid:21)2 A
k∇2(TA+TS)+2µ eij − 3(∇·v)δij , (15) pA(z)=(1−Dz)γ/(γ−1). (22)
Upon dividing equations (13-16) by ρ v /d, ρ g,
r f r
(pA+pS)=(cp−cv)(ρA+ρS)(TA+TS). (16) cpρrvfTr/d,and(cp−cv)ρrTr,respectively,weob-
tain
Within the anelastic approximation, insignifi-
cant terms in the above equations are neglected. (cid:15)∂tρS +∇·[(ρA+(cid:15)ρS)v]=0, (23)
To judge which terms are negligible, the magni-
tude of each single term needs to be estimated,
(cid:15)(ρ +(cid:15)ρ )[∂ v+(v·∇)v]=
A S t
which is best done after a proper rescaling. If not
(cid:32)1− 1 (cid:33)
stated otherwise, from now on non-dimensional −∇ γp +(cid:15)p −(ρ +(cid:15)ρ )zˆ
variables will be used. All spatial variables are D A S A S
scaled with the domain depths d and velocity is (cid:114) (cid:20) (cid:21)
Pr 1
non-dimensionalizedwithaconvectivefree-fallve- +(cid:15) ∇2v+ ∇(∇·v) , (24)
locity v = (cid:112)∆ρgd/ρ . Correspondingly, time Ra 3
f r
is scaled with the free-fall time t = d/v =
f f
(cid:112)
ρrd/(∆ρg). In choosing these units, we implic- (ρA+(cid:15)ρS)[(cid:15)∂tTS +(v·∇)(TA+(cid:15)TS)]
itly assume that shorter timescales, as for exam- (cid:26) (cid:20)(cid:18) (cid:19) (cid:21)(cid:27)
1
plecausedbysoundwaves,playaminorrole. The − (cid:15)D∂ p +(v·∇) 1− p +(cid:15)Dp
t s γ A s
scalefortemperatureT andadiabatictemperature
1
TA isTr,i.e. thetemperatureatthebottomofthe = √ ∇2(T +(cid:15)T )
A S
domain,whilethesuperadiabatictemperaturedif- RaPr
ference ∆T, as dictated by the thermal boundary (cid:114)Pr (cid:20) 1 (cid:21)2
conditions, scales the superadiabatic temperature +2(cid:15)D Ra eij − 3(∇·v)δij , (25)
T . Since temperature and density perturbations
S
are usually assumed to be closely correlated (see
e.g. Clayton 1968), the superadiabatic density D
p +(cid:15) p =(ρ +(cid:15)ρ )(T +(cid:15)T ). (26)
ρ is scaled with the approximate superadiabatic A 1− 1 S A S A S
S γ
density jump ∆ρ = ρ ∆T/T . Consistently, den-
r r
sity ρ and adiabatic background density ρ are Due to the non-dimensionalization with charac-
A
scaled with ρ , which is the adiabatic density at teristicscales,allvariablesandoperatorsareO(1)
r
the bottom of the fluid layer. While pressure p
andadiabaticpressurep arenon-dimensionalized 2Notethatasthetemperatureatthebottomofthedomain,
A
with (c −c )ρ T as suggested by the ideal gas which is dictated by the boundary conditions, is used to
p v r r
scalethetemperature,itfollowsthatT(z=0)=1. There-
law, theappropriatesuperadiabaticpressurescale
fore,T isgenerallynegativeforasuperadiabaticallystrat-
S
∆ρgd can be inferred from the fact that the su- ifiedsystemasconsideredhere.
5
and the magnitude of each term in equations (17- can be interpreted in several different ways. Its
26) can be estimated by its prefactor. The non- name originates from the fact that it constraints
dimensional control parameters (cid:15), Ra, Pr, γ, and how much internal energy can be generated by
D occurring in these coefficients are discussed in viscous dissipation, i.e. D is a measure for the
the following section. significance of viscous heating with 0 ≤ D ≤ 1.
This becomes evident from (27) and (30), which
3.2. Control Parameters and Magnitude allow to rearrange the dissipation number to D =
of the Terms 1 gd∆ρ = 1 Epot . This alternative formula-
γρrcv∆T γ∆Eint
tionrevealsthatthedissipationnumberispropor-
Seven control parameters determine the fate of
tionaltotheratiooftheavailablepotentialenergy
the convection system governed by (20-22) and
E =gd∆ρ, whichdrivesconvection, tothetyp-
(23-26). The superadiabaticity pot
ical internal energy variations ∆E = ρ c ∆T.
int r v
∆T ∆ρ As viscous heating results from the dissipation of
(cid:15)= = (27)
convective kinetic energy (for which E defines
T ρ pot
r r
the upper limit), viscous heating can only signif-
compares the superadiabatic temperature differ- icantly contribute to internal energy variations if
ence as dictated by the boundary conditions to a Epot is of the same order of magnitude as ∆Eint.
typicalreferencetemperaturethatis—asallother D can also be interpreted to be the inverse adi-
reference values—evaluated at the bottom. We abatic temperature scale height evaluated at the
will show later that (cid:15) constrains the typical Mach bottomboundary. Finally,thedissipationnumber
number M, which is defined as the ratio of a typ- is directly linked to the density contrast χ that
ical convective free-fall velocity and the speed of may serve as an alternative parameter. It is de-
sound. The Rayleigh number fined as the ratio of the adiabatic density at the
bottom and at the top,
gd3∆T gd3(cid:15)
Ra= = (28)
νrκrTr νrκr χ= ρbAot = ρA(z =0) =(1−D)−1/(γ−1). (32)
ρtAop ρA(z =1)
controls the vigor of convection with large g, d,
and (cid:15) enhancing and large diffusivities ν and κ The total mass of the fluid, as determined by
reducingtheconvectivevigor. MoreformallyRais the initial conditions, and the aspect ratio of the
the ratio of the product of the diffusive timescales periodicboxformthelasttwocontrolparameters.
d2/κ d2/ν tothesquareofthefree-falltimescale
r r The scaled equations (23-26), which still repre-
t2. The Prandtl number
f sent the full compressible dynamics, can be fur-
ther simplified by noting that the (cid:15)0 terms in
ν
Pr = , (29) equations (24-26) balance exactly. In the mo-
κ
mentum equation (24), the (cid:15)0 terms simply rep-
which for the setup chosen here is constant with resent the hydrostatic balance of the reference
depth, denotes the ratio of momentum diffusivity state, i.e. −(1−1/γ)/D∇pA −ρAzˆ = 0. Like-
to the thermal diffusivity and therefore is a ma- wise, the first two (cid:15)0 terms in the energy equa-
terial property. It effectively controls the impor- tion (25) ρAvz∂zTA −(1−1/γ)vz∂zpA = 0 bal-
tanceofinertiainthesystem,withPr (cid:28)1leading ancebec√auseof(20-22). Notethattheconduction
to strong and Pr (cid:29)1 leading to weak inertial ef- term 1/ RaPr∇2TA drops out here because the
fects. The ratio of the heat capacities defines the adiabatic temperature gradient is constant in our
parameter simple model configuration3. Finally, in equation
c
γ = p, (30)
cv 3Forgeneraldepthdependentheatconductivitieskandadi-
abatic temperature gradients ∇T , this term must be re-
while the Dissipation number A
tained. It then effectively acts as a heat source or sink
anddrivesorhindersconvectionwitht√hemagnitudebeing
gd
estimated by the term’s prefactor 1/ RaPr. For astro-
D = (31)
c T physical systems that exhibit large Rayleigh numbers this
p r
6
(26), the (cid:15)0 terms p = ρ T just represent the
A A A
ideal gas law for the reference state.
By subtracting the (cid:15)0 terms from (24-26) and ∇·(ρAv)=0, (37)
dividing by (cid:15), we arrive at
ρ [∂ v+(v·∇)v]=−∇p
A t S
(cid:15)∂ ρ +∇·[(ρ +(cid:15)ρ )v]=0, (33)
t S A S (cid:114) (cid:20) (cid:21)
Pr 1
−ρ zˆ+ ∇2v+ ∇(∇·v) , (38)
S Ra 3
(ρ +(cid:15)ρ )[∂ v+(v·∇)v]=−∇p
A S t S
(cid:114) (cid:20) (cid:21)
Pr 1
−ρSzˆ+ Ra ∇2v+ 3∇(∇·v) , (34) ρA[∂tTS +(v·∇)TS]−Dρsvz
1
−D[∂ p +(v·∇)p ]= √ ∇2T
t s s S
RaPr
(ρA+(cid:15)ρS)[∂tTS +(v·∇)TS]−Dρsvz (cid:114)Pr (cid:20) 1 (cid:21)2
−D[∂ p +(v·∇)p ]= √ 1 ∇2T +2D Ra eij − 3(∇·v)δij , (39)
t s s S
RaPr
(cid:114)Pr (cid:20) 1 (cid:21)2
+2D Ra eij − 3(∇·v)δij , (35) D pS = TS + ρS. (40)
1− 1 p T ρ
γ A A A
D p T ρ ρ T Note that for the setup chosen here, the supera-
S = S + S +(cid:15) S S (36)
1− 1 p T ρ ρ T diabaticity parameter (cid:15) drops out of the non-
γ A A A A A dimensional anelastic equations4. Furthermore,
the well-known Boussinesq equations describing
which describe fully compressible convection as
shallow convection follow in the limit D →0.
perturbations from the adiabatic, hydrostatic
background state. In a very simple manner the above equations
illustrate the neglected physical processes in the
3.3. Anelastic Approximation anelastic and the Boussinesq approximation: The
continuity equation (33) reveals that by letting
Theenergyconservinganelasticequations,that
(cid:15) → 0, sound waves are effectively filtered out as
can also be derived more formally by a rigorous
the time derivative term becomes negligible. Fur-
amplitudeexpansion(e.g. Gough1969;Lantzand
thermore, unpleasant nonlinearities disappear in
Fan 1999), follow from (33-36) in the limit (cid:15)→0,
(34-36). In the Boussinesq limit D → 0, the en-
resulting in
ergyequation(35)uncoversthatpressurelosesits
role in the energy budget, while viscous heating
magnitude is typically very small and may be comparable
can be neglected as the available potential energy
orevensmallerthanthemagnitudeofthe(cid:15)1 termsrepre-
is much smaller than internal energy variations
senting the usual convective perturbation. For numerical
simulations that do not reach realistic parameter values, (c.f. section 3.2). Equation (36) further shows
the diffusion of adiabatic background temperature, how- thatthesuperadiabaticdensityisdirectlypropor-
ever,maybeofsignificance.
tional to the superadiabatic temperature in the
Boussinesq limit. Finally, the Mach number
(cid:112) (cid:115)
v ∆ρgd/ρ (cid:15)D
M = f = r = , (41)
(cid:112)
vs cp(cp−cv)Tr/cv γ−1
based on the free-fall velocity v and the speed
f
of sound v at the bottom of the domain, can be
s
4For the general case that contains the diffusion of back-
ground temperature, the (cid:15)-parameter controls the signifi-
canceofthisprocessandisthusretained.
7
estimated from the input parameters. Obviously,
it is considered to be small in both, the anelas-
(cid:15) Ra χ Resolution t Re
run
tic and the Boussinesq approximation. Note that 0 104 2.72 1442×129 519 25.0
the Mach number can serve as an alternative con- 0.01 104 2.72 1283 146 25.0
trol parameter that replaces (cid:15). In solar and giant 0.05 104 2.72 1283 326 25.6
planets’interiors,whereD andγ−1cantypically 0.1 104 2.72 1283 463 26.3
be assumed to be O(1), the square of the Mach 0.15 104 2.72 1283 148 27.0
number is crudely approximated by the superadi- 0.2 104 2.72 1283 151 27.7
abaticity, 0.25 104 2.72 1283 77.8 28.4
0.3 104 2.72 1283 185 29.1
M2 ≈(cid:15), (42)
0.35 104 2.72 1283 92.0 30.0
0.4 104 2.72 1283 248 30.9
which suggests that the anelastic approximation
0 105 2.72 1442×129 4180 99.7
holds for M (cid:28)1.
0.1 105 2.72 1283 3519 102
0.2 105 2.72 1283 2990 104
4. NUMERICAL REALIZATION
0.3 105 2.72 1283 2937 107
Theequationsgoverningfullycompressiblecon- 0.4 105 2.72 1283 4260 111
vection (33-36) are solved on a collocated grid us- 0 106 2.72 1922×193 3003 316
ingsecondorderfinitedifferencesandathirdorder 0.1 106 2.72 1923 1364 322
upwind method for the advection terms. A semi- 0.2 106 2.72 1923 1063 330
implicit time stepping scheme based on a third 0.3 106 2.72 1923 2384 339
orderAdams-Bashforth/backward-differencefor- 0.4 106 2.72 1923 2041 350
mula (AB3/BDF3) is applied (e.g. Boyd 2001; 0 107 2.72 2882×257 1913 954
Peyret 2002). All terms except for the vertical 0.1 107 2.72 2563 1373 973
diffusion terms are treated explicitly. 0.2 107 2.72 2563 1262 991
Theanelasticsimulationsthatwillbepresented 0.3 107 2.72 2563 1878 1016
in this paper are performed with an anelastic 0.4 107 2.72 2563 1066 1050
code, which is a modified version of the Boussi- 0 106 4.48 1922×193 2744 300
nesq code by Stellmach and Hansen (2008). It 0.1 106 4.48 1923 1205 313
uses a mixed pseudo-spectral fourth order finite- 0 106 7.39 1922×193 2535 294
difference discretization of the spatial deriva- 0.1 106 7.39 1923 1422 299
tives and an AB3/BDF3 time integration scheme, 0 106 12.18 1922×193 2811 280
which treats all linear terms implicitly. Instead 0.1 106 12.18 1923 1347 286
of using (33-36) directly, for numerical reasons it 0 106 20.1 1922×193 2945 269
turns out to be beneficial to use an equivalent 0.1 106 20.1 1923 1464 271
formulationbasedonentropyratherthantemper-
Table 1: Overview of the simulations carried out
ature. Therelevantequations(C5-C7)arederived
for this study, with Pr = 0.7 and γ = 5/3 apply-
in detail in appendix C.
ing to all simulations. The horizontal dimensions
of the simulation domain are l =l =2d, result-
5. RESULTS x y
ing in an aspect ratio of two. While the spatial
In this section results from a suite of anelas- resolution is given in the number of x, y, and z
tic and fully compressible direct numerical simu- grid points, trun denotes the run time measured
(cid:112)
lations (DNS) are presented in order to test the in free-fall times, and Re = vrms/ Pr/Ra is
accuracy and efficiency of both approaches in the theapproximatedReynoldsnumber,withthenon-
fully nonlinear regime of convection. dimensionalroot-mean-squarevelocityvrms being
definedinequation(47). WhileallRa=104 cases
Equations (33-36) are solved for various su-
result in stationary solutions, the remaining sim-
peradiabaticities (0 ≤ (cid:15) ≤ 0.4), density contrasts
ulations stay time-dependent.
(2.72 ≈ exp(1) ≤ χ ≤ exp(3) ≈ 20.1 or analo-
gously 0.49 ≤ D ≤ 0.86) and Rayleigh numbers
8
a) b)
Fig. 2.—TypicalvolumerenderingsofthesuperadiabatictemperatureT (a)andverticalvelocityv (b)for
S z
an anelastic simulation run that reached statistical equilibrium. Red colors denote warm, buoyant material
and positive v , blue signifies cold fluid and negative v , and yellow structures refer to intermediate values
z z
ofT andv . Thecorrespondingparametersare(cid:15)=0, Ra=107, Pr =0.7, γ =5/3andχ=exp(1)≈2.72.
S z
Correspondingsnapshotstakenfromnumericalsimulationsoffullycompressibleconvectionlookqualitatively
the same and cannot be visually distinguished from the displayed example. A stereoscopic 3-d version of
these volume renderings, which reflects the full 3-d structures when wearing red-cyan filter glasses, is shown
in figure 9 in the appendix.
(104 ≤Ra≤107). See Table 1 for an overview of isthepolytropicindex,determinesthetotalmass.
all simulations. The simulation runs with (cid:15) = 0 Note that the polytropic index n is often used as
are carried out with the anelastic code, while the an alternative parameter to the superadiabaticity
remainingsimulationsareexecutedwithourinde- (cid:15) (e.g. Brummell et al. 1996, 1998; Berkoff et al.
pendent code for fully compressible convection as 2010). For χ ≈ 2.72 and 0.1 ≤ (cid:15) ≤ 0.4, which are
described in section 4. The remaining four con- typical parameters for this study, the polytropic
trol parameters are kept constant for all simula- index varies within the range 1.07≥n≥0.37.
tions, with the Prandtl number set to Pr = 0.7 To give the reader a feeling for the level of tur-
and the ratio of specific heats chosen to repre- bulence reached in our simulations, figure 2 shows
sent a monoatomic ideal gas, γ = 5. The hor- a typical snapshot of an anelastic simulation run
3
izontal dimensions of the simulation domain are that reached statistical equilibrium. Correspond-
lx = ly = 2d, resulting in an aspect ratio of two. ing snapshots taken from numerical simulations
The total mass of the fluid is determined by the of fully compressible convection look qualitatively
initial state, for which we choose the hydrostatic, similar.
conductive solution with
5.1. Comparison of fully compressible and
T(t=0,z)=TA+(cid:15)TS(t=0) anelastic results
=[1−(D+(cid:15))z]. (43)
Global diagnostic quantities can provide a first
impressionastowhatextenttheanelasticapprox-
Theintegraloverthecorrespondinginitialdensity
imation holds. Initially, we vary (cid:15) and Ra, while
distribution
keeping χ = exp(1) ≈ 2.72 constant. Covering
several orders of magnitude in Rayleigh number
ρ(t=0,z)=ρ +(cid:15)ρ (t=0)
A S
Ra, figure 3 shows three different global diagnos-
=[1−(D+(cid:15))z]n, (44)
tics plotted against the superadiabaticity param-
eter (cid:15). From left to right, the graphs in the top
where
γ D row display the heat flux in terms of the Nusselt
n= −1 (45)
γ−1D+(cid:15)
9
20 0.32 0.05
n
16 ms Eki 0.04
uxNu 12 cityvr 0.28 dens. 0.03
fl o y
Heat 8 msvel 0.24 energ 0.02
4 R Kin.
0.20 0.01
1.5 Ra=104
1.3 Ra=105 2.2
Ra=106
1.4 Ra=107
ǫ=0u 1.3 ǫ=0vrms 1.2 ǫ=0Ekin 1.8
N / /
Nu/ 1.2 vrms Ekin
1.1 1.4
1.1
1.0 1.0 1.0
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4
Superadiabaticityǫ Superadiabaticityǫ Superadiabaticityǫ
Fig. 3.— Global diagnostic quantities are plotted against the superadiabaticity parameter (cid:15). From left to
right,thegraphsinthetoprowdisplaytheheatfluxintermsofaNusseltnumberNu,therootmeansquare
velocityv ,andthekineticenergydensityE . Thebottomrowshowsthesamequantitiesnormalizedto
rms kin
thecorrespondinganelasticvalues((cid:15)=0). ForallRayleighnumbers, thefullycompressibleresultsconverge
to the anelastic values for (cid:15)→0. For large Rayleigh numbers Ra and superadiabaticities (cid:15) smaller than 0.3,
theoutputsfromcompressibleconvectiondeviatebynomorethan30%fromtheassociatedanelasticvalues.
The error bars given for the Nusselt numbers are estimates based on the difference between temporally
averaged Nusselt numbers computed at the top and bottom boundary. In all cases, χ=2.72.
10