Con en s lis s a ailable a ScienceDi ec
Compu e s and Fluids
jou nal homepage: www.else ie .com/loca e/comp luid
Ene gy-consis en disc e iza ion o iscous dissipa ion wi h applica ion o
na u al con ec ion low
B. Sande se a,∗, F.X. T ias b
aCen um Wiskunde & In o ma ica, Science Pa k 123, Ams e dam, The Ne he lands
bHea and Mass T ans e Technological Cen e , Technical Uni e si y o Ca alonia, ESEIAAT, c/ Colom 11, 08222 Te assa (Ba celona), Spain
ARTICLE INFO
Da ase link: h ps://gi hub.com/bsande se/IN
S2D,h ps://gi hub.com/agdes ein/Incomp ess
ibleNa ie S okes.jl
Keywo ds:
Viscous dissipa ion
Ene gy conse a ion
S agge ed g id
Na u al con ec ion
Rayleigh–Béna d
Gebha numbe
ABSTRACT
A new ene gy-consis en disc e iza ion o he iscous dissipa ion unc ion in incomp essible lows is p oposed.
I is implied by choosing a disc e iza ion o he di usi e e ms and a disc e iza ion o he local kine ic
ene gy equa ion and by equi ing ha con inuous iden i ies like he p oduc ule a e mimicked disc e ely. The
p oposed iscous dissipa ion unc ion has a quad a ic, s ic ly dissipa i e o m, o bo h simpli ied (cons an
iscosi y) s ess enso s and gene al s ess enso s. The p oposed exp ession is no only use ul in e alua ing
ene gy budge s in u bulen lows, bu also in na u al con ec ion lows, whe e i appea s in he in e nal ene gy
equa ion and is esponsible o iscous hea ing. The iscous dissipa ion unc ion is such ha a consis en o al
ene gy balance is ob ained: he ‘implied’ p esence as sink in he kine ic ene gy equa ion is exac ly balanced by
explici ly adding i as sou ce e m in he in e nal ene gy equa ion.
Nume ical expe imen s o Rayleigh–Béna d con ec ion (RBC) and Rayleigh–Taylo ins abili ies con i m ha
wi h he p oposed dissipa ion unc ion, he ene gy exchange be ween kine ic and in e nal ene gy is exac ly
p ese ed. The expe imen s show u he mo e ha iscous dissipa ion does no a ec he c i ical Rayleigh
numbe a which ins abili ies o m, bu i does signi ican ly impac he de elopmen o ins abili ies once hey
occu . Consequen ly, he alue o he Nussel numbe on he cold pla e becomes la ge han on he ho pla e,
wi h he di e ence inc easing wi h inc easing Gebha numbe . Finally, 3D simula ions o u bulen RBC show
ha ene gy balances a e exac ly sa is ied e en o e y coa se g ids. The e o e, he p oposed disc e iza ion
also o ms an excellen s a ing poin o es ing sub-g id scale models and is a use ul ool o assess ene gy
budge s in any u bulence simula ion, wi h o wi hou he p esence o na u al con ec ion.
1. In oduc ion and p oblem desc ip ion
In his a icle we s udy he iscous dissipa ion unc ion and i s ole
in na u al con ec ion lows desc ibed by he incomp essible Na ie –
S okes equa ions, wi h buoyancy e ec s modeled by he Boussinesq
app oxima ion [1]. These ‘Boussinesq‘ o ‘Obe beck-Boussinesq’ equa-
ions ha e a ac ed much scien i ic in e es o e se e al decades [2],
no only because o hei physical ele ance, bu also o hei in igu-
ing ma hema ical p ope ies. An impo an es case s udied wi h he
Boussinesq sys em is ha o Rayleigh–Béna d con ec ion [3], in which
a box o luid is hea ed om he bo om and cooled om he op,
gi ing ise o con ec ion cells. The Boussinesq equa ions also desc ibe
a (miscible) o m o Rayleigh–Taylo ins abili y, which occu s when a
hea y (cold) luid is posi ioned abo e a ligh (wa m) luid.
A common assump ion in many incomp essible na u al con ec ion
s udies is ha he e ec o iscous dissipa ion on he in e nal ene gy
(e ec i ely on he empe a u e) is neglec ed. This assump ion is no
∗Co esponding au ho .
E-mail add ess: [email p o ec ed] (B. Sande se).
always alid, o example when conside ing na u al con ec ion in he
Ea h man le, when conside ing highly iscous liquids, when la ge
leng h scales a e in ol ed, o in de ices ope a ing a high o a ional
speed [4–11]. O cou se, when conside ing comp essible lows, e.g.
high-speed lows, including hea ing by iscous dissipa ion is known o
be impo an , and se e al benchma king s udies ha e been pe o med
ela ed o modeling na u al con ec ion in he Ea h man le [12,13].
These s udies ypically assume in ini e P and l numbe s, and igno e
he uns eady and con ec i e e ms in he momen um equa ions. In
his pape , we will es ic ou sel es o he incomp essible si ua ion
wi h he Boussinesq app oxima ion. Ne e heless, we an icipa e ha
ou idea o disc e izing he iscous dissipa ion e m in an ene gy-
consis en manne has a b oade scope o applicabili y since i is also
applicable o non-Obe beck–Boussinesq [14] and comp essible lows
(see e.g. [15,16]).
In he incomp essible case, Os ach [11], Gebha [10] and Tu co e
e al. [9] should be explici ly men ioned, being among he i s o
h ps://doi.o g/10.1016/j.comp luid.2024.106473
Recei ed 20 June 2024; Recei ed in e ised o m 24 Sep embe 2024; Accep ed 1 No embe 2024
Compu e s and Fluids 286 (2025) 106473
A ailable online 13 No embe 2024
0045-7930/© 2024 The Au ho (s). Published by Else ie L d. This is an open access a icle unde he CC BY license (
h p://c ea i ecommons.o g/licenses/by/4.0/ ).
B. Sande se and F.X. T ias
add ess he ole o iscous dissipa ion and o in oduce nex o he well-
known Rayleigh and P and l numbe s ano he dimensionless quan i y,
which is known as he dissipa ion numbe o he Gebha numbe . In
addi ion, we men ion he wo k o Ba le a and co-au ho s [17–20],
who conside ed he ole o iscous dissipa ion in na u al con ec ion in
se e al pape s, s udying he co ec ma hema ical o mula ion o he
p oblem and linea s abili y analysis o di e en geome ies. Tu co e
e al. [9] we e p obably one o he i s o pe o m nume ical expe i-
men s o incomp essible na u al con ec ion lows ha include iscous
dissipa ion. They pe o med simula ions on coa se g ids (10 ×10)
and low Rayleigh numbe s (Ra = 104,105) o di e en alues o he
dissipa ion numbe and concluded ha Rayleigh–Béna d con ec ion
was signi ican ly a ec ed when he dissipa ion numbe was o o de
uni y. The main quan i y o in e es is he Nussel numbe , which is a
measu e o he hea ans e a he walls.
F om an ene gy pe spec i e, he iscous dissipa ion sou ce e m in
he in e nal ene gy equa ion occu s as a sink in he kine ic ene gy
equa ion, which cancel each o he when conside ing he o al ene gy
equa ion. Howe e , mos ene gy analyses, especially o incomp essible
low, ocus on he ole o he po en ial ene gy e m and i s spli in o
a ailable and backg ound po en ial ene gy [21–23], o on he kine ic
ene gy budge [24]. To he au ho ’s knowledge, he ole o iscous
dissipa ion in he kine ic ene gy equa ion and i s nume ical ea men
o he in e nal ene gy equa ion ha e no been explo ed in de ail. In
addi ion, e en o he cases whe e iscous dissipa ion does no ha e a
s ong e ec on he empe a u e, ha ing an accu a e and consis en way
o e alua ing he iscous and he mal dissipa ion budge s (as p esen
in o example G ossmann–Lohse heo y [3]) is ano he bene i o ou
p oposed disc e iza ion scheme.
In his pape , he main no el y is ha we p opose a disc e iza ion
o he iscous dissipa ion unc ion and apply i o he con ex o na u al
con ec ion low, whe e i appea s as a sou ce e m in he in e nal
ene gy equa ion. Ou disc e iza ion is such ha we ge a co ec global
ene gy balance, on con inuous, semi-disc e e, and ully disc e e le el.
Fi s , on he con inuous le el, a non-dimensionaliza ion is p oposed
ha makes he in e nal and kine ic ene gy scaling consis en . Second,
on he semi-disc e e le el, we p opose a disc e e dissipa ion ope a o ,
and show ha i canno be chosen eely bu is implied by he dis-
c e iza ion o he iscous e ms in he momen um equa ions and by
he de ini ion o he kine ic ene gy. Thi d, on he ully disc e e le el,
we p opose a ime in eg a ion me hod ha p ese es he o al ene gy
balance upon ime ma ching.
Impo an ly, he disc e e dissipa ion ope a o ha we p opose he e
is no es ic ed o he con ex o na u al con ec ion lows. Fo exam-
ple, when es ima ing he dissipa ion o kine ic ene gy in DNS o LES
simula ions o u bulen lows, a consis en exp ession o he ene gy
dissipa ion is c ucial in e alua ing ene gy budge s and de eloping sub-
g id scale models. Ano he example o which ou dissipa ion ope a o
is impo an is he case o incomp essible Taylo –Coue e low ( he low
be ween wo o a ing cylinde s): he powe supplied o he cylinde s is
con e ed in o hea ing o he luid, which can be e y signi ican (1
K/min o he se -up wi h wa e epo ed in [25]). In o de o p edic
he co ec empe a u e inc ease, he iscous dissipa ion ope a o ha
we p opose in his wo k is needed. In exis ing simula ions o Taylo –
Coue e low his e ec is igno ed (p obably because in expe imen al
se -ups ac i e cooling is used o keep he empe a u e unde con ol).
Wi h ou dissipa ion ope a o , one can pe o m mo e ealis ic s udies, in
which in e nal hea ing h ough iscous dissipa ion and cooling h ough
he bounda ies can bo h be included.
The pape is s uc u ed as ollows. Sec ion 2in oduces he go -
e ning equa ions, ene gy balances, and new non-dimensionaliza ion.
Sec ions 3and 4desc ibe he ene gy-consis en spa ial and empo al
disc e iza ion. Sec ion 5desc ibes s eady-s a e esul s o Rayleigh–
Béna d con ec ion including iscous dissipa ion, and Sec ion 6de-
sc ibes ene gy-conse ing simula ions o Rayleigh–Taylo ins abili ies
including iscous dissipa ion. Sec ion 7shows he e ec o iscous
dissipa ion in 3D DNS o Rayleigh–Béna d con ec ion.
Fig. 1. P oblem se -up o Rayleigh–Béna d con ec ion.
2. Ene gy-conse ing o mula ion
2.1. Go e ning equa ions
The Boussinesq app oxima ion s a es ha densi y a ia ions a e
small and can be igno ed in all e ms o he Na ie –S okes (NS) equa-
ions, excep in he one pe aining o he g a i y e m. The NS equa ions
desc ibing conse a ion o mass and momen um hen ead
∇⋅𝒖= 0,(1)
𝜌0(𝜕𝒖
𝜕 𝑡+ ∇⋅(𝒖⊗𝒖))= −∇𝑝+𝜇∇2𝒖+𝜌𝒈,(2)
whe e 𝒖(𝒙, 𝑡)is he eloci y ield, 𝑝(𝒙, 𝑡) he p essu e, 𝜇 he dynamic is-
cosi y, 𝜌(𝒙, 𝑡) he densi y and 𝜌0a e e ence densi y. Wi hou loss o gen-
e ali y, we conside a wo-dimensional (ins ead o h ee-dimensional)
domain 𝛺, wi h he g a i y ec o poin ing in he nega i e 𝑦-di ec ion
so ha 𝒈= −𝑔𝒆𝑦. An example o he domain as used in he Rayleigh–
Béna d p oblem, including he bounda y condi ions, is gi en in Fig. 1.
In he esul s sec ion we will also conside he Rayleigh–Taylo p ob-
lem, which has adiaba ic bounda ies on op and bo om, ins ead o
iso he mal as in case o Rayleigh–Béna d.
The densi y 𝜌is assumed o a y only wi h empe a u e 𝑇(𝒙, 𝑡),
acco ding o 𝜌(𝑇) =𝜌0−𝛽 𝜌0(𝑇−𝑇0), whe e 𝛽is he isoba ic coe icien o
he mal expansion (𝛽= −1
𝜌(𝜕 𝜌
𝜕 𝑇)𝑝). The NS equa ions a e hen w i en
as
𝜌0(𝜕𝒖
𝜕 𝑡+ ∇⋅(𝒖⊗𝒖))= −∇𝑝′+𝜇∇2𝒖−𝛽 𝜌0(𝑇−𝑇0)𝒈,(3)
whe e 𝑝=𝑝′−𝜌0𝑔 𝑦and ∇𝑝= ∇𝑝′−𝜌0𝑔𝒆𝑦.
The equa ion o he in e nal ene gy 𝑒𝑖desc ibes he empe a u e
e olu ion acco ding o
𝜕
𝜕 𝑡(𝜌0𝑐 𝑇
⏟⏟⏟
𝑒𝑖
) + ∇⋅(𝒖(𝜌0𝑐 𝑇)) =𝛷+𝜆∇2𝑇 ,(4)
whe e 𝜆is he he mal conduc i i y and 𝑐equals 𝑐𝑣in case o an ideal
gas ( he speci ic hea a cons an olume), and equals 𝑐𝑝−𝑝𝛽
𝜌 o a eal
gas [26]). The con ibu ion o p essu e wo k o he change in in e nal
ene gy, 𝑝∇⋅𝒖, has been disca ded in Eq. (3) because o Eq. (1).
The iscous dissipa ion unc ion
𝛷∶=
𝝉∶ ∇𝒖=𝜇[2(𝜕 𝑢
𝜕 𝑥)2
+ 2(𝜕 𝑣
𝜕 𝑦)2
+(𝜕 𝑢
𝜕 𝑦+𝜕 𝑣
𝜕 𝑥)2]≥0.(5)
is he key quan i y in his wo k, whe e he s ess enso is gi en by
𝝉=𝜇(∇𝒖+ (∇𝒖)𝑇). Exp ession (5) holds in 2D and is easily gene alized
Compu e s and Fluids 286 (2025) 106473
2
B. Sande se and F.X. T ias
o 3D. Since he luid is incomp essible and 𝜇is assumed cons an , we
ha e ∇⋅
𝝉=𝜇∇⋅(∇𝒖+ (∇𝒖)𝑇) =𝜇∇2𝒖+𝜇∇(∇ ⋅𝒖) =𝜇∇2𝒖, which is he
o m o he di usi e e ms used in Eq. (3). The simpli ied o m could
be in e p e ed as ∇⋅𝝉=𝜇∇2𝒖, wi h 𝝉=𝜇∇𝒖, al hough 𝝉is no a p ope
s ess enso (i is no symme ic). Rema kably, he simpli ied o m o
he di usi e e ms implies adi e en dissipa ion unc ion, namely
𝛷∶= 𝜇‖∇𝒖‖2≥0.(6)
whe e ‖∇𝒖‖2= ∇𝒖∶ ∇𝒖( he F obenius inne p oduc ). The de ails
ega ding he di e ence be ween 𝛷and
𝛷a e gi en in Appendix A. In
2D and Ca esian coo dina es he iscous dissipa ion can be w i en as
𝛷=𝜇[(𝜕 𝑢
𝜕 𝑥)2
+(𝜕 𝑢
𝜕 𝑦)2
+(𝜕 𝑣
𝜕 𝑥)2
+(𝜕 𝑣
𝜕 𝑦)2].(7)
In his wo k we ocus mainly on he disc e iza ion o exp ession (7),
bu we will also explain he disc e iza ion o he mo e gene al o m (5),
see Appendix B, Eq. (B.21).
2.2. To al ene gy conse a ion
Conse a ion o kine ic ene gy ollows by aking he do p oduc
o Eq. (3) wi h 𝒖:
𝜕
𝜕 𝑡(1
2𝜌0|𝒖|2
⏟⏞⏟⏞⏟
𝑒𝑘
) + ∇⋅(1
2𝜌0|𝒖|2𝒖) = −𝒖⋅∇𝑝′+𝜇∇⋅(𝒖⋅∇𝒖) −𝜇‖∇𝒖‖2+𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣,
(8)
whe e 𝒈⋅𝒖= −𝑔 𝑣and we ha e used he iden i y
𝒖⋅∇2𝒖= −‖∇𝒖‖2+ ∇⋅(𝒖⋅∇𝒖).(9)
Upon adding he kine ic and in e nal ene gy Eqs. (4) and (8), he
iscous dissipa ion e m cancels and we a i e a he equa ion o he
o al ene gy 𝑒=𝑒𝑘+𝑒𝑖:
𝜕
𝜕 𝑡(𝑒𝑘+𝑒𝑖) + ∇⋅((𝑒𝑘+𝑒𝑖)𝒖) = −∇⋅(𝑝′𝒖) +𝜇∇⋅(𝒖⋅∇𝒖) +𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣+𝜆∇2𝑇 .
(10)
All e ms a e in conse a i e (di e gence) o m, excep he po en ial
ene gy e m. Upon in eg a ing o e he domain 𝛺and assuming no-slip
condi ions 𝒖=𝟎on all bounda ies, we ob ain he global balances
d𝐸𝑘
d𝑡= −∫𝛺
𝛷d𝛺+∫𝛺
𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣d𝛺 ,(11)
d𝐸𝑖
d𝑡=∫𝛺
𝛷d𝛺+∫𝜕 𝛺
𝜆∇𝑇⋅𝒏d𝑆 ,(12)
d𝐸
d𝑡=d𝐸𝑘
d𝑡+d𝐸𝑖
d𝑡=∫𝛺
𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣d𝛺+∫𝜕 𝛺
𝜆∇𝑇⋅𝒏d𝑆 ,(13)
whe e 𝐸=∫𝛺𝑒d𝛺=𝐸𝑘+𝐸𝑖. In case he bounda y condi ions a e
adiaba ic (∇𝑇⋅𝒏= 0), he las e m in (13) anishes and he o al
ene gy equa ion exp esses ha he sum o in e nal and kine ic ene gy
changes due o he buoyancy lux ∫𝛺𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣d𝛺— his case will
be deal wi h in he Rayleigh–Taylo se -up in Sec ion 6. No e ha in
comp essible lows, he buoyancy lux can be w i en in e ms o he
ime de i a i e o he po en ial ene gy — see Appendix D. In such a
case, he Boussinesq sys em wi h he iscous dissipa ion unc ion in
he in e nal ene gy equa ion and wi h adiaba ic bounda ies can be
conside ed o be uly ene gy-conse ing, wi h he o al ene gy being
he sum o kine ic, in e nal and po en ial ene gy. In he incomp essible
case, his does no hold, so we use he e m ‘ene gy-consis en ’ o
indica e ha we ha e included he iscous dissipa ion unc ion in he
in e nal ene gy equa ion in such a way ha he o al ene gy equa ion
is no a ec ed by i .
Table 1
Di e en non-dimensional o ms esul ing om di e en choices o 𝑢 e .
𝑢 e 𝛼1=𝜈
𝑢 e 𝐻𝛼2=𝛽 𝑔 𝛥𝑇 𝐻
𝑢2
e
𝛼3=𝜈 𝑢 e
𝑐 𝛥𝑇 𝐻𝛼4=𝜅
𝑢 e 𝐻𝛾=𝛼1
𝛼3
I√𝛽 𝑔 𝛥𝑇 𝐻√P
Ra 1Ge√P
Ra
1
√P Ra
1
Ge
II 𝜅
𝐻P P Ra Ge
Ra 1P Ra
Ge
III √𝑐 𝛥𝑇 √P Ge
Ra Ge √P Ge
Ra √Ge
P Ra 1
In mos s udies o Rayleigh–Béna d con ec ion he dissipa ion unc-
ion 𝛷is le ou om he in e nal ene gy Eq. (4), while i s co espond-
ing coun e pa in he momen um equa ion (𝜇∇2𝒖) is s ill included. As
a consequence, he ene gy los in he kine ic ene gy equa ion is no
balanced by he hea gene a ed in he in e nal ene gy equa ion, so ha
he o al ene gy equa ion ea u es a dissipa ion e m, which des oys
he global ene gy balance.
2.3. Non-dimensionaliza ion
We non-dimensionalize Eqs. (1),(2) and (4) by aking a e e ence
leng h 𝐻(ca i y heigh ), a e e ence empe a u e di e ence 𝛥𝑇 (di -
e ence be ween he cold and ho pla es), and a e e ence eloci y 𝑢 e
ye o be speci ied. These choices de e mine he ime scale 𝐻∕𝑢 e and
he p essu e scale 𝜌0𝑢2
e . An impo an ques ion, which we will add ess
he e, is how he choice o non-dimensionaliza ion changes he o al
ene gy equa ion. The non-dimensional equa ions a e w i en as ( o
de ails, see Appendix C):
∇⋅
𝒖= 0,(14)
𝜕
𝒖
𝜕
𝑡+
∇⋅(
𝒖⊗
𝒖)= −
∇𝑝′+𝛼1
∇2
𝒖+𝛼2
𝑇𝒆𝑦,(15)
𝜕
𝑇
𝜕
𝑡+
∇⋅(
𝒖
𝑇) =𝛼3
𝛷+𝛼4
∇2
𝑇 ,(16)
whe e he pa ame e s 𝛼𝑖,𝑖= 1 … 4a e a unc ion o he Rayleigh num-
be Ra =𝛽 𝑔 𝛥𝑇 𝐻3
𝜈 𝜅, he P and l numbe P =𝜈
𝜅and he Gebha numbe
Ge =𝛽 𝑔 𝐻
𝑐(also known as he dissipa ion numbe [6]). In Table 1we
p esen h ee di e en op ions o 𝑢 e wi h he co esponding alues
o 𝛼. Choices I and II a e common in li e a u e, see o example [27]
o choice I and [1,18,28] o choice II; hey co espond o a ee-
all eloci y scale and he he mal di usi i y scale, espec i ely. O he
choices a e also possible, e.g. 𝑢 e =𝛽 𝑔 𝛥𝑇 𝐻2∕𝜈[5], bu his choice does
no lead o a ‘clean’ exp ession in e ms o he dimensionless numbe s
de ined abo e. To ou bes knowledge, choice III is new and inspi ed
by he o m o he o al ene gy equa ion. Physically, his choice can
be in e p e ed as he eloci y ha is ob ained when in e nal ene gy is
ans o med in o kine ic ene gy.
The non-dimensional o m o he o al ene gy equa ion ollows by
aking he do p oduc o (15) wi h
𝒖and add he in e nal ene gy
Eq. (16). The global ene gy balances in non-dimensional o m ead
d
𝐸𝑘
d
𝑡= −𝛼1
𝛬∫
𝛺
𝛷d
𝛺+𝛼2
𝛬∫
𝛺
𝑇 𝑣 d
𝛺 ,(17)
d
𝐸𝑖
d
𝑡=𝛼3
𝛬∫
𝛺
𝛷d
𝛺+𝛼4
𝛬∫𝜕
𝛺
∇
𝑇⋅𝒏d
𝑆 ,(18)
d
𝐸
d
𝑡=d
𝐸𝑘
d𝑡+𝛾d
𝐸𝑖
d𝑡=𝛼2
𝛬∫
𝛺
𝑇 𝑣 d
𝛺+𝛾 𝛼4
𝛬∫𝜕
𝛺
∇
𝑇⋅𝒏d
𝑆 ,(19)
whe e
𝐸=1
𝛬∫
𝛺𝑒 d
𝛺,𝛬=𝐿∕𝐻is he aspec a io o he box, and
𝛾=𝛼1
𝛼3is a weigh ing ac o , which is epo ed in Table 1 o di e en
choices o 𝑢 e . Fo de ini ions o 𝑒𝑘,𝑒𝑖and 𝑒, see Appendix C. The
p oposed choice III is he only choice ha ea u es 𝛾= 1, meaning ha
he dimensionless kine ic and in e nal ene gy equa ion a e consis en
wi h each o he and do no equi e a weigh ing ac o in o de o he
iscous dissipa ion e m o cancel.
The choice o a pa icula e e ence eloci y ypically depends on
he p oblem a hand. Choices I and II ha e he ad an age ha in case
Compu e s and Fluids 286 (2025) 106473
3
B. Sande se and F.X. T ias
o Ge = 0(mos commonly in es iga ed in li e a u e), one ob ains
𝛼3= 0and he dissipa ion e ms simply d ops om he in e nal ene gy
equa ion. Howe e , when Ge is small bu nonze o, he weigh ac o
𝛾becomes e y la ge o choices I and II. Choice III does no su e
om his issue, because 𝛾= 1independen o Ge, so kine ic ene gy and
in e nal ene gy can be summed independen o Ge. Howe e , choice III
has he disad an age ha i does no wo k in he case Ge = 0, since
i leads o 𝛼𝑖= 0 o all 𝑖. In summa y: o Ge = 0, choices I and II
a e p e e ed; o small bu nonze o Ge, choice III is p e e ed; in o he
cases, all choices a e ine.
The discussion in he nex sec ions will be agnos ic o he choice
o 𝑢 e , and exp essed in e ms o he gene al pa ame e s 𝛼𝑖. No e ha
in he simula ions in Sec ions 5–7, we will employ choice I. Choices II
and III gi e equi alen esul s apa om scaling ac o s.
2.4. E ec o iscous dissipa ion on Nussel numbe and he mal dissipa ion
A main quan i y o in e es in na u al con ec ion lows is he Nussel
numbe Nu and we will in es iga e how i changes upon including
iscous dissipa ion in he in e nal ene gy equa ion. Fi s , de ine he
a e age o he sum o con ec i e and conduc i e luxes h ough a
ho izon al plane 𝑦=𝑦′by
𝐹(𝑦′) ∶= 1
𝐿∫𝐿
0(𝜌0𝑐 𝑇 𝑣−𝜆𝜕 𝑇
𝜕 𝑦)(𝑥,𝑦′)
d𝑥. (20)
Then, he Nussel numbe based on 𝐹 ollows as [3]:
Nu( 𝑦′) ∶= 𝐹(𝑦′)
𝜆𝛥𝑇 ∕𝐻=1
𝛬∫𝛬
0(1
𝛼4
𝑇 𝑣 −𝜕
𝑇
𝜕 𝑦 )(𝑥, 𝑦′)
d𝑥. (21)
Fo s eady s a e o s a is ically s eady s a e (using a sui able a e age),
and in he absence o iscous dissipa ion, i is s aigh o wa d o show
om he in e nal ene gy equa ion ha Nu( 𝑦) = Nu( 𝑦 = 0) = Nu, which is
a cons an , independen o 𝑦′[1,28]. Howe e , upon including iscous
dissipa ion, his ela ion no longe holds ue and ins ead he s eady
in e nal ene gy equa ion yields
𝛼4(Nu( 𝑦′) − Nu(0)) =𝛼3𝜖𝑈(𝑦′),(22)
whe e he in eg a ed dissipa ion unc ion is gi en by
𝜖𝑈(𝑦′) ∶= 1
𝛬∫𝑦′
0∫𝛬
0
𝛷d𝑥d𝑦. (23)
Eq. (22) is an impo an ela ion which shows ha ( aking 𝑦′= 1)
𝛼4(Nu(1) − Nu(0)) =𝛼3𝜖𝑈(1),(24)
so he Nussel numbe o he uppe pla e is always la ge han o equal o
he Nussel numbe o he lowe pla e.
A second ela ion be ween Nussel numbe and iscous dissipa ion
can be ob ained om he global kine ic ene gy balance, Eq. (17). The
second e m in he igh -hand side o Eq. (17) can be ew i en wi h
Eq. (22), ollowing he analysis in [28]:
𝛼2
𝛬∫
𝛺
𝑇 𝑣 d
𝛺=𝛼2
𝛬∫1
0∫𝛬
0
𝑇 𝑣 d𝑥d𝑦 =𝛼2𝛼4∫1
0
Nu( 𝑦) d𝑦
+𝛼2𝛼4
𝛬∫𝛬
0∫1
0
𝜕
𝑇
𝜕 𝑦 d𝑦d𝑥
=𝛼2𝛼4Nu(0) +𝛼2𝛼3∫1
0
𝜖𝑈(𝑦) d𝑦
+𝛼2𝛼4
𝛬∫𝛬
0
(
𝑇(𝑥, 𝑦 = 1) −
𝑇(𝑥, 𝑦 = 0))d 𝑥
=𝛼2𝛼4(Nu(0) − 1) +𝛼2𝛼3∫1
0
𝜖𝑈(𝑦) d𝑦.
(25)
Fo (s a is ically) s eady low, his e m equals he i s e m in he
igh -hand side o Eq. (17), yielding he second ela ion be ween he
Nussel numbe and he iscous dissipa ion 𝜖𝑈
𝛼2𝛼4(Nu(0) − 1) =𝛼1𝜖𝑈(1) −𝛼2𝛼3∫1
0
𝜖𝑈(𝑦) d𝑦. (26)
We ecognize he well-known equa ion 𝛼2𝛼4(Nu(0) − 1) =𝛼1𝜖𝑈(1), see
e.g. [1], bu wi h he addi ional nega i e e m −𝛼2𝛼3∫1
0𝜖𝑈(𝑦) d𝑦.
Las ly, we link he he mal dissipa ion 𝜖𝑇 o he Nussel numbe and
he iscous dissipa ion unc ion. The non-dimensional in e nal ene gy
equa ion, Eq. (16), is mul iplied by
𝑇, and a e in eg a ing by pa s,
using he skew-symme y o he con ec i e ope a o , and employing
he bounda y condi ion
𝑇(𝑦 = 1) = 0, one ob ains
1
𝛬
d
d𝑡∫
𝛺
1
2
𝑇2d
𝛺=𝛼3
𝛬∫
𝛺
𝑇
𝛷d
𝛺−𝛼4
𝛬∫𝛬
0(
𝑇𝜕
𝑇
𝜕 𝑦 )𝑦=0
d𝑥
−𝛼4
𝛬∫
𝛺‖
∇
𝑇‖2d
𝛺 .(27)
Wi h he bounda y condi ion
𝑇(𝑦 = 0) = 1, and he assump ion o
(s a is ically) s eady low, his ela ion is u he simpli ied o
𝛼4Nu(0) =𝛼4𝜖𝑇−𝛼3
𝛬∫
𝛺
𝑇
𝛷d
𝛺 ,(28)
whe e
𝜖𝑇∶= 1
𝛬∫
𝛺‖
∇
𝑇‖2d
𝛺 .(29)
Since
𝑇≥0,
𝛷≥0, we conclude ha iscous dissipa ion lowe s he
Nussel numbe o he lowe pla e. In absence o iscous dissipa ion in
he in e nal ene gy equa ion, one ob ains he amilia ela ion Nu =𝜖𝑇.
In combina ion wi h Eq. (24), we ob ain o he Nussel numbe o he
uppe pla e:
𝛼4Nu(1) =𝛼4𝜖𝑇+𝛼3
𝛬∫(1 −
𝑇)
𝛷d
𝛺 .(30)
Assuming ha he empe a u e sa is ies 0≤
𝑇≤1, we ind ha iscous
dissipa ion inc eases he Nussel numbe o he uppe pla e. In o he wo ds,
he he mal dissipa ion lies in be ween he wo Nussel numbe s:
Nu(0) ≤𝜖𝑇≤Nu(1).(31)
The h ee ela ions (24),(26) and (28) a e summa ized in Table 2and
will be con i med in he nume ical expe imen s in Sec ion 5.
3. Ene gy-consis en spa ial disc e iza ion
3.1. Mass, momen um and kine ic ene gy equa ion
To disc e ize he non-dimensional mass and momen um Eqs. (14)
and (15), we use he s agge ed-g id ene gy-conse ing ini e olume
me hod desc ibed in [29], ex ended by including he buoyancy e m
in he momen um equa ions. This leads o he ollowing semi-disc e e
equa ions:
𝑀 𝑉ℎ(𝑡) = 0,(32)
𝛺𝑉
d𝑉ℎ(𝑡)
d𝑡= −𝐶𝑉(𝑉ℎ(𝑡)) −𝐺 𝑝ℎ(𝑡) +𝛼1𝐷𝑉𝑉ℎ(𝑡) +𝛼2(𝐴𝑇ℎ(𝑡) +𝑦𝑇).(33)
He e, 𝑉ℎ∈R𝑁𝑉a e he eloci y unknowns, 𝑝ℎ∈R𝑁𝑝 he p essu e
unknowns, and 𝑇ℎ∈R𝑁𝑝 he empe a u e unknowns; see Fig. 2 o
hei posi ioning. 𝑀∈R𝑁𝑝×𝑁𝑉is he disc e ized di e gence ope a o ,
𝐺= −𝑀𝑇∈R𝑁𝑉×𝑁𝑝 he disc e ized g adien ope a o , 𝛺𝑉∈R𝑁𝑉×𝑁𝑉
a ma ix wi h he ‘ eloci y’ ini e olume sizes on i s diagonal, and 𝐶𝑉
and 𝐷𝑉cons i u e cen al di e ence app oxima ions o he con ec i e
and di usi e e ms. 𝐴is a ma ix ha a e ages he empe a u e om
he cen e o he ‘ empe a u e’ ini e olumes o cen e o he ‘ eloci y’
ini e olumes, and he ec o 𝑦𝑇inco po a es he nonze o bounda y
condi ion o he empe a u e a he lowe pla e.
The ene gy-conse ing na u e o ou ini e olume me hod is c ucial
in de i ing an ene gy-consis en disc e iza ion o iscous dissipa ion.
The ene gy-conse ing p ope y means ha , in absence o bounda y
con ibu ions, he disc e ized con ec i e and p essu e g adien ope a-
o s do no con ibu e o he kine ic ene gy balance: 𝑉𝑇
ℎ𝐶𝑉(𝑉ℎ) = 0and
𝑉𝑇
ℎ𝐺 𝑝ℎ= 0, jus like in he con inuous case. This is achie ed by using a
Compu e s and Fluids 286 (2025) 106473
4
B. Sande se and F.X. T ias
Table 2
S eady-s a e Nussel numbe ela ions, wi h and wi hou iscous dissipa ion.
O igin Wi hou iscous dissipa ion Wi h iscous dissipa ion
In e nal Nu(1) = Nu(0) 𝛼4(Nu(1) − Nu(0)) =𝛼3𝜖𝑈(1)
Kine ic 𝛼2𝛼4(Nu(0) − 1) =𝛼1𝜖𝑈(1) 𝛼2𝛼4(Nu(0) − 1) =𝛼1𝜖𝑈(1) −𝛼2𝛼3∫1
0𝜖𝑈(𝑦) d𝑦
In e nal ene gy ×𝑇Nu(0) =𝜖𝑇𝛼4Nu(0) =𝛼4𝜖𝑇−𝛼3
𝛬∫
𝛺
𝑇
𝛷d
𝛺
Fig. 2. S agge ed g id wi h posi ioning o unknowns a ound a p essu e olume.
skew-symme ic con ec ion ope a o and he compa ibili y be ween 𝑀
and 𝐺 ia 𝐺= −𝑀𝑇. The disc e e kine ic ene gy balance hen eads:
d𝐸𝑘,ℎ
d𝑡= −𝛼1𝜖𝑈 ,ℎ +𝛼2𝑉𝑇
ℎ(𝐴𝑇ℎ+𝑦𝑇),(34)
whe e 𝐸𝑘,ℎ =1
2𝑉𝑇
ℎ𝛺𝑉𝑉ℎ. The global iscous dissipa ion (i.e. summed
o e he en i e domain) is gi en by 𝜖𝑈 ,ℎ =‖𝑄𝑉ℎ‖2
2>0, whe e 𝑄s ems
om decomposing he symme ic nega i e-de ini e di usi e ope a o
as 𝐷𝑉= −𝑄𝑇𝑄. Eq. (34) is he semi-disc e e coun e pa o Eq. (17).
3.2. P oposed iscous dissipa ion unc ion
Gi en a disc e iza ion ha sa is ies a disc e e kine ic ene gy bal-
ance, he key s ep is o design a disc e iza ion scheme o he in e nal
ene gy Eq. (16) which is such ha disc e e e sions o he global bal-
ances (12) and (13) a e ob ained. In pa icula , he iscous dissipa ion
in he in e nal ene gy equa ion should cancel he iscous dissipa ion
e m in he kine ic ene gy equa ion, whe e he la e is ully de e -
mined by he choice o he di usion ope a o and he exp ession o he
local kine ic ene gy. The choice o he di usion ope a o (second-o de
cen al di e encing) is s aigh o wa d. The choice o he exp ession
o he local kine ic ene gy on a s agge ed g id is howe e no ob ious.
We p opose he ollowing de ini ion:
𝑘𝑖,𝑗 ∶= 1
4𝑢2
𝑖+1∕2,𝑗 +1
4𝑢2
𝑖−1∕2,𝑗 +1
4𝑣2
𝑖,𝑗+1∕2 +1
4𝑣2
𝑖,𝑗−1∕2.(35)
This choice gi es a local kine ic ene gy equa ion ha is consis en wi h
he con inuous equa ions, as is de ailed in Appendix B, and consis en
wi h he global ene gy de ini ion.
The exp ession o 𝛷ℎ hen ollows om di e en ia ing he ex-
p ession o 𝑘𝑖𝑗 in ime, subs i u ing he momen um equa ions, and
ew i ing he e ms in ol ing he di usi e ope a o (see Appendix B).
The implied dissipa ion hen ollows by cons uc ing a disc e e e -
sion o (9). As example, we cons uc he disc e e e sion o 𝑢𝜕2𝑢
𝜕 𝑥2=
−(𝜕 𝑢
𝜕 𝑥)2
+𝜕
𝜕 𝑥(𝑢𝜕 𝑢
𝜕 𝑥), being
𝑢𝑖+1∕2,𝑗
𝛥𝑥 (𝑢𝑖+3∕2,𝑗 −𝑢𝑖+1∕2,𝑗
𝛥𝑥 −𝑢𝑖+1∕2,𝑗 −𝑢𝑖−1∕2,𝑗
𝛥𝑥 )
= −1
2(𝑢𝑖+3∕2,𝑗 −𝑢𝑖+1∕2,𝑗
𝛥𝑥 )2
−1
2(𝑢𝑖+1∕2,𝑗 −𝑢𝑖−1∕2,𝑗
𝛥𝑥 )2
+1
𝛥𝑥 (1
2(𝑢𝑖+3∕2,𝑗 +𝑢𝑖+1∕2,𝑗 )𝑢𝑖+3∕2,𝑗 −𝑢𝑖+1∕2,𝑗
𝛥𝑥
−1
2(𝑢𝑖+1∕2,𝑗 +𝑢𝑖−1∕2,𝑗 )𝑢𝑖+1∕2,𝑗 −𝑢𝑖−1∕2,𝑗
𝛥𝑥 ).(36)
The i s wo e ms on he igh -hand side con ibu e o he iscous
dissipa ion unc ion. Repea ing his p ocess o he o he componen s
(𝑢𝜕2𝑢
𝜕 𝑦2,𝑣𝜕2𝑣
𝜕 𝑥2,𝑣𝜕2𝑣
𝜕 𝑦2), as ou lined in Appendix B.2, yields he ollowing
no el exp ession o he local dissipa ion unc ion:
𝛷𝑖,𝑗 =1
2𝛷𝑢
𝑖+1∕2,𝑗 +1
2𝛷𝑢
𝑖−1∕2,𝑗 +1
2𝛷𝑣
𝑖,𝑗+1∕2 +1
2𝛷𝑣
𝑖,𝑗−1∕2 ,(37)
whe e
𝛷𝑢
𝑖+1∕2,𝑗 = −1
2(𝑢𝑖+3∕2,𝑗 −𝑢𝑖+1∕2,𝑗
𝛥𝑥 )2
−1
2(𝑢𝑖+1∕2,𝑗 −𝑢𝑖−1∕2,𝑗
𝛥𝑥 )2
−1
2(𝑢𝑖+1∕2,𝑗+1 −𝑢𝑖+1∕2,𝑗
𝛥𝑦 )2
−1
2(𝑢𝑖+1∕2,𝑗 −𝑢𝑖+1∕2,𝑗−1
𝛥𝑦 )2
,(38)
𝛷𝑣
𝑖,𝑗+1∕2 = −1
2(𝑣𝑖+1,𝑗+1∕2 −𝑣𝑖,𝑗+1∕2
𝛥𝑥 )2
−1
2(𝑣𝑖+1,𝑗−1∕2 −𝑣𝑖,𝑗−1∕2
𝛥𝑥 )2
−1
2(𝑣𝑖,𝑗+3∕2 −𝑣𝑖,𝑗+1∕2
𝛥𝑦 )2
−1
2(𝑣𝑖,𝑗+1∕2 −𝑣𝑖,𝑗−1∕2
𝛥𝑦 )2
.(39)
A bounda ies, an adap a ion o 𝛷ℎis equi ed in o de o ha e a
disc e e equi alen o Eq. (9). This is de ailed in Eq. (B.15).
No e ha 𝛷ℎis de i ed based on local ene gy conside a ion which
upon summa ion equals he global dissipa ion, jus like Eq. (23):
1𝑇𝛺𝑝𝛷ℎ=𝜖𝑈 ,ℎ.(40)
3.3. In e nal ene gy equa ion
Ha ing p oposed a consis en exp ession o 𝛷ℎ, he spa ial dis-
c e iza ion o he in e nal ene gy Eq. (16) eads:
𝛺𝑝
d𝑇ℎ
d𝑡= −𝐶𝑇(𝑉ℎ, 𝑇ℎ) +𝛼3𝛺𝑝𝛷ℎ(𝑉ℎ) +𝛼4(𝐷𝑇𝑇ℎ+𝑦𝑇),(41)
whe e
[𝐶𝑇(𝑉ℎ, 𝑇ℎ)]𝑖,𝑗 =𝛥𝑦 (𝑢𝑖+1∕2,𝑗
1
2(𝑇𝑖+1,𝑗 +𝑇𝑖,𝑗 ) −𝑢𝑖−1∕2,𝑗
1
2(𝑇𝑖,𝑗 +𝑇𝑖−1,𝑗 ))+
𝛥𝑥 (𝑣𝑖,𝑗+1∕2
1
2(𝑇𝑖,𝑗+1 +𝑇𝑖,𝑗 ) −𝑣𝑖,𝑗−1∕2
1
2(𝑇𝑖,𝑗 +𝑇𝑖,𝑗−1))(42)
is he con ec ion ope a o . The con ec ion ope a o has a disc e e
skew-symme y p ope y which will be used in he de i a ion o he
he mal dissipa ion balance in he nex subsec ion. 𝐷𝑇 he s anda d
second-o de di e ence s encil wi h bounda y condi ions encoded in
𝑦𝑇.
The o al in e nal ene gy is gi en by 𝐸𝑖,ℎ = 1𝑇𝛺𝑝𝑇ℎ(simply sum-
ming o e all ini e olumes). Due o he no-slip bounda y condi ions
on he eloci y ield, he con ec i e ope a o sa is ies 1𝑇𝐶𝑇(𝑉ℎ, 𝑇ℎ) = 0.
The summa ion o e he di usi e ope a o can be w i en in e ms o
he Nussel numbe s (de ailed in he nex sec ion). The o al in e nal
ene gy equa ion hus eads
d𝐸𝑖,ℎ
d𝑡=𝛼31𝑇𝛺𝑝𝛷ℎ+𝛼41𝑇(𝐷𝑇𝑇ℎ+𝑦𝑇),
=𝛼31𝑇𝛺𝑝𝛷ℎ+𝛼4(Nu𝐻− Nu𝐶),
(43)
Compu e s and Fluids 286 (2025) 106473
5
B. Sande se and F.X. T ias
Fig. 3. S eady-s a e empe a u e ield o Ra = 105on a 128 ×128 g id, o di e en Ge.
whe e in he second line he Nussel numbe s a e ins an aneous Nussel
numbe s. Upon adding he o al kine ic ene gy Eq. (34), and using
p ope y (40), he global ene gy balance esul s:
d𝐸ℎ
d𝑡=d𝐸𝑘,ℎ
d𝑡+𝛾d𝐸𝑖,ℎ
d𝑡=𝛼2𝑉𝑇
ℎ(𝐴𝑇ℎ+𝑦𝑇) +𝛾 𝛼41𝑇(𝐷𝑇𝑇ℎ+𝑦𝑇),
=𝛼2𝑉𝑇
ℎ(𝐴𝑇ℎ+𝑦𝑇) +𝛾 𝛼4(Nu𝐻− Nu𝐶),
(44)
which is he semi-disc e e coun e pa o Eq. (19). In o he wo ds,
we ha e p oposed a disc e e iscous dissipa ion unc ion ha leads o
a co ec exp ession o he o al ene gy equa ion, namely such ha
he iscous dissipa ion om he kine ic and in e nal ene gy equa ions
exac ly balances, independen o he mesh size. No e ha in he case o
homogeneous Neumann bounda y condi ions o he empe a u e on all
bounda ies, he las e m disappea s.
3.4. Disc e e global balances and Nussel numbe ela ions
We now de i e disc e e e sions o he Nussel ela ions ha inco -
po a e he iscous dissipa ion unc ion, i.e. ela ions (24) and (28). Ou
symme y-p ese ing spa ial disc e iza ion is such ha exac disc e e
ela ions can be de i ed. I is impo an o ealize ha he disc e e
app oxima ion o he Nussel numbe canno be chosen independen ly
(when he goal is o ha e exac disc e e global balances) bu is implic-
i ly de ined once he disc e iza ion o he di usi e ope a o is chosen.
Conside he disc e ized global in e nal ene gy equa ion o s eady
condi ions,
𝛼31𝑇𝛺𝑝𝛷ℎ(𝑉ℎ) +𝛼41𝑇(𝐷𝑇𝑇ℎ+𝑦𝑇) = 0.(45)
The second e m can be simpli ied as
1𝑇(𝐷𝑇𝑇ℎ+𝑦𝑇) = −
𝑁𝑥
∑
𝑖=1
𝑇𝑖,1−𝑇𝐻
1
2𝛥𝑦
𝛥𝑥+
𝑁𝑥
∑
𝑖=1
𝑇𝐶−𝑇𝑖,𝑁𝑦
1
2𝛥𝑦
𝛥𝑥 = Nu𝐻− Nu𝐶,(46)
whe e he Nussel numbe s on he lowe (ho ) and uppe (cold) pla e
a e de ined as
Nu𝐻∶= −
𝑁𝑥
∑
𝑖=1
𝑇𝑖,1−𝑇𝐻
1
2𝛥𝑦
𝛥𝑥, (47)
Nu𝐶∶= −
𝑁𝑥
∑
𝑖=1
𝑇𝐶−𝑇𝑖,𝑁𝑦
1
2𝛥𝑦
𝛥𝑥. (48)
This leads o he disc e e e sion o (24):
𝛼4(Nu𝐶− Nu𝐻) =𝛼31𝑇𝛺𝑝𝛷ℎ(𝑉ℎ).(49)
The disc e e e sion o (28) ollows by conside ing he inne p oduc
o Eq. (41) wi h 𝑇𝑇
ℎins ead o 1𝑇. An impo an p ope y o he
con ec i e disc e iza ion (42) is ha
𝑇𝑇
ℎ𝐶𝑇(𝑉ℎ, 𝑇ℎ) = 0,∀𝑇ℎ,i 𝑀 𝑉ℎ= 0.(50)
This p ope y is mos easily de i ed by ecognizing ha 𝐶ℎ(𝑉ℎ, 𝑇ℎ)can
be w i en in e ms o a ma ix– ec o p oduc
𝐶𝑇(𝑉ℎ)𝑇ℎ, whe e
𝐶𝑇(𝑉ℎ)
is skew-symme ic i 𝑀 𝑉ℎ= 0. In addi ion, he inne p oduc o 𝑇ℎwi h
he di usi e e ms can be w i en as
𝑇𝑇
ℎ(𝐷𝑇𝑇ℎ+𝑦𝑇) =
𝑁𝑥
∑
𝑖=1 ⎛⎜⎜⎝
−𝑇𝐻
𝑇𝑖,1−𝑇𝐻
1
2𝛥𝑦
+𝑇𝐶
𝑇𝐶−𝑇𝑖,𝑁𝑦
1
2𝛥𝑦 ⎞⎟⎟⎠
𝛥𝑥 −𝜖𝑇 ,ℎ,(51)
whe e
𝜖𝑇 ,ℎ ∶=
𝑁𝑥
∑
𝑖=1 ⎛⎜⎜⎜⎝
1
2⎛⎜⎜⎝
𝑇𝑖,1−𝑇𝐻
1
2𝛥𝑦 ⎞⎟⎟⎠
2
+
𝑁𝑦
∑
𝑗=2 (𝑇𝑖,𝑗 −𝑇𝑖,𝑗−1
𝛥𝑦 )2
+1
2⎛⎜⎜⎝
𝑇𝐶−𝑇𝑖,𝑁𝑦
1
2𝛥𝑦 ⎞⎟⎟⎠
2⎞⎟⎟⎟⎠
×𝛥𝑥𝛥𝑦 +
𝑁𝑦
∑
𝑗=1
𝑁𝑥
∑
𝑖=2 (𝑇𝑖,𝑗 −𝑇𝑖−1,𝑗
𝛥𝑥 )2
𝛥𝑥𝛥𝑦
(52)
is he disc e e analogue o (27) and Eq. (51) is he disc e e e sion o
∫𝑇d2𝑇
d𝑦2= [𝑇d𝑇
d𝑦] −∫(d𝑇
d𝑦)2. Wi h he bounda y condi ion 𝑇𝐻= 1,𝑇𝐶= 0,
we ge he balance
𝛼4Nu𝐻=𝛼4𝜖𝑇 ,ℎ −𝛼3𝑇𝑇
ℎ𝛺𝑝𝛷ℎ(𝑉ℎ),(53)
which is he disc e e e sion o Eq. (28).
4. Ene gy-consis en empo al disc e iza ion
The sys em o Eqs. (32),(33) and (41) needs o be in eg a ed in
ime wi h a sui able me hod in o de o p ese e a ime-disc e e e sion
o he global ene gy balance (44). A common choice is o use an
explici me hod (e.g. Adams–Bash o h) o he nonlinea con ec i e
e ms and an implici me hod (e.g. C ank–Nicolson) o he (s i ) linea
di usion e ms [23,28,30], o an explici me hod o bo h con ec ion
and di usion [31,32]. In such an app oach, he empe a u e equa ion is
ypically sol ed i s (gi en eloci y ields a p e ious ime ins ances),
and hen he mass and momen um equa ions a e sol ed wi h a p essu e-
co ec ion app oach. Howe e , hese me hods do no p ese e he
global ene gy balance as hey iola e he ene gy-conse ing na u e o
he nonlinea e ms when ma ching in ime [33].
Ins ead, we show he e ha he implici midpoin me hod can be em-
ployed o achie e ene gy-consis en ime in eg a ion. The ully disc e e
sys em eads:
𝑀 𝑉𝑛+1∕2
ℎ= 0,(54)
𝛺𝑉
𝑉𝑛+1
ℎ−𝑉𝑛
ℎ
𝛥𝑡 = −𝐶𝑉(𝑉𝑛+1∕2
ℎ) −𝐺 𝑝𝑛+1∕2
ℎ+𝛼1𝐷𝑉𝑉𝑛+1∕2
ℎ
+𝛼2(𝐴𝑇 𝑛+1∕2
ℎ+𝑦𝑇),(55)
Compu e s and Fluids 286 (2025) 106473
6
B. Sande se and F.X. T ias
𝛺𝑝
𝑇𝑛+1
ℎ−𝑇𝑛
ℎ
𝛥𝑡 = −𝐶𝑇(𝑉𝑛+1∕2
ℎ, 𝑇𝑛+1∕2
ℎ) +𝛼3𝛺𝑝𝛷(𝑉𝑛+1∕2
ℎ)
+𝛼4(𝐷𝑇𝑇𝑛+1∕2
ℎ+𝑦𝑇).(56)
He e 𝑉𝑛+1∕2
ℎ=1
2(𝑉𝑛
ℎ+𝑉𝑛+1
ℎ)and 𝑇𝑛+1∕2
ℎ=1
2(𝑇𝑛
ℎ+𝑇𝑛+1
ℎ). Upon
mul iplying (55) by (𝑉𝑛+1∕2
ℎ)𝑇and (56) by 1𝑇, and adding he wo
esul ing equa ions, we ge he disc e e ene gy balance,
𝐸𝑛+1
ℎ−𝐸𝑛
ℎ
𝛥𝑡 =
𝐸𝑛+1
𝑘,ℎ −𝐸𝑛
𝑘,ℎ
𝛥𝑡 +𝛾
𝐸𝑛+1
𝑖,ℎ −𝐸𝑛
𝑖,ℎ
𝛥𝑡
=𝛼2(𝑉𝑛+1∕2
ℎ)𝑇(𝐴𝑇 𝑛+1∕2
ℎ+𝑦𝑇) +𝛾 𝛼41𝑇(𝐷𝑇𝑇𝑛+1∕2
ℎ+𝑦𝑇),(57)
which is he ully-disc e e coun e pa o Eqs. (19) and (44). The de i a-
ions hinges again on skew-symme y o he con ec ion ope a o 𝐶𝑉,
he compa ibili y be ween 𝑀and 𝐺(𝐺= −𝑀𝑇), and he consis ency
equi emen on he iscous dissipa ion unc ion, Eq. (40). Again, we
s ess ha he disc e e ene gy equa ion esul s om exac ly balanc-
ing he iscous dissipa ion be ween he kine ic and in e nal ene gy
equa ions, independen o he mesh size and he ime s ep.
The sys em o Eqs. (54)–(56) leads o a la ge sys em o nonlinea
equa ions which has a saddle poin s uc u e due o he di e gence-
ee cons ain . We sol e he sys em in a seg ega ed ashion and i e a e
a each ime s ep wi h a s anda d p essu e-co ec ion me hod un il he
esidual o he en i e sys em is below a p esc ibed ole ance.
We will compa e his ene gy-conse ing ime in eg a ion app oach
o an explici one-leg me hod commonly used o di ec nume ical
simula ions [31,32] in Sec ion 6. This one-leg scheme eads
𝑀 𝑉𝑛+1
ℎ= 0,(58)
𝛺𝑉
(𝑏+1
2)𝑉𝑛+1
ℎ− 2𝑏𝑉 𝑛
ℎ+ (𝑏−1
2)𝑉𝑛−1
ℎ
𝛥𝑡 = −𝐶𝑉(𝑉∗
ℎ) −𝐺 𝑝𝑛+1
ℎ
+𝛼1𝐷𝑉𝑉∗
ℎ+𝛼2(𝐴𝑇 ∗
ℎ+𝑦𝑇),
(59)
𝛺𝑝
(𝑏+1
2)𝑇𝑛+1
ℎ− 2𝑏𝑇 𝑛
ℎ+ (𝑏−1
2)𝑇𝑛−1
ℎ
𝛥𝑡 = −𝐶𝑇(𝑉∗
ℎ, 𝑇∗
ℎ)
+𝛼3𝛺𝑝𝛷(𝑉∗
ℎ) +𝛼4(𝐷𝑇𝑇∗
ℎ+𝑦𝑇),
(60)
whe e 𝑉∗
ℎ= (1 +𝑏)𝑉𝑛
ℎ−𝑏𝑉 𝑛−1
ℎand 𝑇∗
ℎ= (1 +𝑏)𝑇𝑛
ℎ−𝑏𝑇 𝑛−1
ℎ, and we will
ake 𝑏=1
2. The ene gy balance o his scheme ollows again om
mul iplying he momen um equa ions by (𝑉𝑛+1∕2
ℎ)𝑇and he in e nal
ene gy equa ion by 1𝑇, and adding he wo, leading o
𝐸𝑛+1
ℎ−𝐸𝑛
ℎ
𝛥𝑡 = −(𝑉𝑛+1∕2
ℎ)𝑇𝐶𝑉(𝑉∗
ℎ) +𝛼1(𝑉𝑛+1∕2
ℎ𝐷𝑉𝑉∗
ℎ+ 1𝑇𝛺𝑝𝛷(𝑉∗
ℎ))
+𝛼2(𝑉𝑛+1∕2
ℎ)𝑇(𝐴𝑇 ∗
ℎ+𝑦𝑇) +𝛾 𝛼41𝑇(𝐷𝑇𝑇𝑛+1∕2
ℎ+𝑦𝑇).(61)
We obse e ha he explici na u e o he one-leg scheme in oduces
wo e o s in he ene gy equa ion. Fi s , he con ec i e e ms do no
cancel om he ene gy equa ion, because (𝑉𝑛+1∕2
ℎ)𝑇𝐶𝑉(𝑉∗
ℎ)is no equal
o ze o. Second, a con ibu ion om he iscous dissipa ion unc ion
appea s as (𝑉𝑛+1∕2
ℎ)𝑇𝐷𝑉𝑉∗
ℎdoes no exac ly cancel 1𝑇𝛺𝑝𝛷(𝑉∗
ℎ).
5. S eady s a e esul s (Rayleigh–Béna d)
The concep o ene gy consis ency is bes demons a ed h ough
ime-dependen simula ions. Howe e , we s a wi h s eady-s a e e-
sul s in o de o alida e he spa ial disc e iza ion me hod and o ge
in ui ion o he e ec o he Gebha numbe on he Nussel numbe .
Fo he esul s epo ed he e we employ a di ec sol e ha sol es he
en i e coupled non-linea sys em o equa ions ha a ises om spa ial
disc e iza ion. As ini ial guess we ake he ollowing di e gence- ee
eloci y ield:
𝑢(𝑥, 𝑦) = −64𝑥2(𝑥− 1)2𝑦(𝑦− 1)(2𝑦− 1),(62)
𝑣(𝑥, 𝑦) = 64𝑥(𝑥− 1)(2𝑥− 1)𝑦2(𝑦− 1)2,(63)
Table 3
Con e gence o Nussel numbe (47) wi h g id e inemen o di e en Rayleigh
numbe s and Ge = 0.
G id Ra = 103Ra = 104Ra = 105
3221.000 2.170 3.933
6421.000 2.161 3.916
12821.000 2.159 3.912
25621.000 2.158 3.911
Cai e al. [35] (2562) 1.000 2.158 3.911
Table 4
Con e gence o Nussel numbe s (47) and (48) wi h g id e inemen o di e en
Rayleigh and di e en Gebha numbe s.
(a) Ra = 104(b) Ra = 105
G id Ge = 0.1 Ge = 1G id Ge = 0.1 Ge = 1
Nu𝐻Nu𝐶Nu𝐻Nu𝐶Nu𝐻Nu𝐶Nu𝐻Nu𝐶
3222.111 2.228 1.582 2.729 3223.786 4.080 2.448 5.319
6422.103 2.219 1.578 2.716 6423.770 4.062 2.441 5.299
12822.101 2.217 1.576 2.713 12823.766 4.057 2.439 5.293
25622.100 2.216 1.576 2.712 25623.765 4.056 2.439 5.292
which is inspi ed by he egula ized d i en ca i y p oblem [34]. Fo
he empe a u e we ake a andom ield (be ween 0 and 1). The idea
behind his choice o ini ial condi ion is o a oid he non-linea sol e
o be s uck in he i ial solu ion (𝒖= 0). No e ha in all simula ions in
his a icle, we will se P = 0.71 (ai ), and use non-dimensionaliza ion
choice I. Choices II and III gi e equi alen esul s apa om scaling
ac o s.
5.1. G id con e gence s udy o no-dissipa ion case (Ge =0)
Fig. 3(a) shows he empe a u e ield when iscous dissipa ion is
no included (Ge = 0). The esul ing Nussel numbe s as a unc ion
o g id e inemen a e displayed in Table 3and indica e excellen
ag eemen wi h li e a u e [35]. We no e ha he Nussel numbe s as
de ined by (47) and (48) a e i s -o de app oxima ions. Mo e accu a e
app oxima ions can be cons uc ed by including mo e in e io poin s.
We a e no using such high-o de accu a e app oxima ions as hey
would no sa is y he disc e e global balance (49). No e also ha we
only epo Nu𝐻since Nu𝐶= Nu𝐻up o machine p ecision.
5.2. G id con e gence s udy o iscous dissipa ion case (Ge >0)
When including iscous dissipa ion (Ge >0) in he in e nal ene gy
equa ion, he low ield changes quali a i ely and loses i s symme ic
na u e, as can be obse ed in Figs. 3(b)–3(c). The Nussel numbe s a
he ho and cold pla e s a o de ia e om each o he , hei di e ence
being equal o he dissipa ion unc ion, acco ding o Eq. (49) (o
(24)). This is epo ed in Table 4and Fig. 4(a). The c i ical Rayleigh
numbe ha we ind om he bi u ca ion diag am is Ra𝑐≈ 2585,
which is in excellen ag eemen wi h he alue o 2585.02 epo ed in
li e a u e [36,37]. I is independen o he alue o he P and l numbe ,
as shown in [36], and also independen o he alue o he Gebha
numbe . This la e ac ollows by ex ending he linea s abili y anal-
ysis in [36] and ealizing ha he e m ∇𝒖∶ ∇𝒖wi h 𝒖=𝒖0+𝜀𝒖′
and backg ound s a e 𝒖0= 0leads o he e m 𝜀2∇𝒖′∶ ∇𝒖′, which
disappea s when ga he ing e ms o (𝜀). The esul s in Fig. 4(a) show
indeed ha he bi u ca ion poin is he same o di e en alues o Ge.
Fig. 4(b) shows a di e en in e p e a ion o he Nussel numbe ,
indica ing he ela ion wi h he he mal dissipa ion and iscous dis-
sipa ion acco ding o Eq. (53) (o (28)). The esul s con i m ha he
he mal dissipa ion lies in be ween he Nussel numbe o he ho and
cold pla e.
Compu e s and Fluids 286 (2025) 106473
7
B. Sande se and F.X. T ias
Fig. 4. Bi u ca ion diag am o Rayleigh–Béna d p oblem including iscous dissipa ion.
6. Time-dependen , ene gy-conse ing simula ion (Rayleigh–
Taylo )
The p e ious sec ion con i med he (disc e e) s eady-s a e Nussel
numbe balances. In his sec ion we conside he co e idea o his
a icle: achie ing exac ene gy conse a ion in a ime-dependen simu-
la ion. Exac ene gy conse a ion equi es ha all con ibu ions om
bounda y e ms disappea , which we achie e by p esc ibing no-slip
condi ions 𝒖= 0and adiaba ic condi ions 𝜕 𝑇
𝜕 𝑛= 0on all bounda ies ( he
p essu e does no equi e bounda y condi ions). The ene gy balance
hen ep esen s a pu e exchange o kine ic, in e nal and po en ial
ene gy acco ding o
𝐸𝑛+1
ℎ−𝐸𝑛
ℎ
𝛥𝑡 =
𝐸𝑛+1
𝑘,ℎ −𝐸𝑛
𝑘,ℎ
𝛥𝑡 +𝛾
𝐸𝑛+1
𝑖,ℎ −𝐸𝑛
𝑖,ℎ
𝛥𝑡 =𝛼2(𝑉𝑛+1∕2
ℎ)𝑇(𝐴𝑇 𝑛+1∕2
ℎ+𝑦𝑇).
(64)
Howe e , wi h adiaba ic bounda y condi ions we canno simula e he
classic Rayleigh–Béna d p oblem. Ins ead, we u n o he well-known
Rayleigh–Taylo p oblem, ea u ing a cold (hea y) luid on op o a
wa m (ligh ) luid. A ske ch o he se -up is shown in Fig. 5. The ene gy-
conse ing implici midpoin (‘IM’) me hod de ailed in Sec ion 4will be
compa ed o he explici one-leg (‘OL’) me hod commonly used in DNS
s udies [31,32] (whe e we ake 𝑏=1
2and a ixed ime s ep).
The domain size is 1×2, he g id is 64 ×128, he ime s ep
𝛥𝑡 = 5⋅10−3 and he end ime 𝑇= 50. We conside he case P = 0.71,
Ra = 106and Ge = {0.1,1}. The ins abili y does na u ally a ise due
o g ow h o ound-o e o s, bu his akes a he long, so ins ead
a pe u ba ion is added o he ini ial in e ace: 𝑦= 1 + 0.05 sin(2𝜋 𝑥).
The ins abili y quickly de elops and an asymme y in he solu ion
appea s, igge ing a sequence o well-known ‘mush oom’ ype plumes:
ho plumes ising upwa d and cold plumes sinking downwa d (Fig. 6).
The de elopmen o he ins abili y is essen ially he same o IM and
OL — see also Fig. 7(a) o a mo e quan i a i e compa ison. No e ha i
no pe u ba ion is added, he onse o s abili y is sensi i e o he choice
o ime in eg a ion me hod, due o di e ences in he accumula ion o
ound-o e o s. Fig. 6also shows ha he ini ial de elopmen o he
ins abili y is insensi i e o he alue o Ge, jus like he bi u ca ion poin
in he s eady s a e Rayleigh–Béna d simula ion was insensi i e o he
alue o Ge.
Since he e is no d i ing o ce and all bounda y condi ions a e
homogeneous, iscosi y damps he eloci y ield back o a homoge-
neous s eady s a e, while a he same ime inc easing he empe a u e
h ough dissipa ion. This inc ease in empe a u e is clea om Fig. 7(a),
whe e he a e age empe a u e is displayed. Compa ed o he ini ial
empe a u e di e ence 𝛥𝑇 = 1, he ela i e empe a u e inc ease is
Fig. 5. P oblem se -up wi h ini ial condi ion o Rayleigh Taylo p oblem.
abou 2% o Ge = 0.1and mo e han 20% o Ge = 1. No e ha
many exis ing na u al con ec ion models, which igno e he iscous
dissipa ion e m, would no p edic any empe a u e inc ease. Wi h ou
p oposed ene gy-consis en iscous dissipa ion unc ion, he empe a-
u e inc ease exac ly ma ches he kine ic ene gy loss h ough iscous
dissipa ion. This is con i med in Fig. 7(b), which shows he ene gy e o
𝜀𝐸∶= ||||||
𝐸𝑛+1
𝑘,ℎ −𝐸𝑛
𝑘,ℎ
𝛥𝑡 +𝛾
𝐸𝑛+1
𝑖,ℎ −𝐸𝑛
𝑖,ℎ
𝛥𝑡 −𝛼2(𝑉𝑛+1∕2
ℎ)𝑇(𝐴𝑇 𝑛+1∕2
ℎ+𝑦𝑇)||||||
.(65)
Fo IM he e o emains a he ole ance wi h which we sol e he
sys em o nonlinea equa ions (10−12). Fo OL, he e o is a ound
(10−6)when he ins abili y is mos p onounced (a ound 𝑡= 5, see
Fig. 6), and dec eases when he low se les back o a s eady s a e. This
small ene gy e o o he OL scheme seems accep able gi en ha OL
is oughly 4–5×less expensi e han IM, because IM equi es oughly
4–5 i e a ions (Poisson sol es) pe ime s ep, ins ead o only 1 o OL.
Consequen ly, OL will be employed o he 3D simula ions in he nex
sec ion. No e ha his balance o accu acy e sus compu a ional cos s
depends on he de ails o he low p oblem and migh di e in o he
es cases.
Compu e s and Fluids 286 (2025) 106473
8
B. Sande se and F.X. T ias
Fig. 6. Rayleigh–Taylo empe a u e ields a 𝑡= 5 o di e en Ge and di e en ime in eg a ion me hods (IM =Implici Midpoin , OL =One-Leg scheme).
Fig. 7. Rayleigh–Taylo esul s, IM =Implici Midpoin , OL =One-Leg scheme.
7. Ene gy-conse ing simula ion o a u bulen low
As a inal es -case, we conside he nume ical simula ion o an
ai - illed (P = 0.71) Rayleigh–Béna d low a wo di e en Rayleigh
numbe s, Ra = 108and 1010. Di ec nume ical simula ions (DNS) we e
ca ied ou and analyzed in p e ious s udies [38,39] wi hou aking
in o accoun he iscous dissipa ion e ec s (Ge = 0). He e, he esul s
a e ex ended o Ge = 0.1and Ge = 1keeping he same domain size
(𝜋× 1 × 1) and mesh esolu ion (400 × 208 × 208 o Ra = 108, and
1024 ×768 ×768 o Ra = 1010). G ids a e cons uc ed wi h a uni o m
g id spacing in he pe iodic 𝑥-di ec ion whe eas wall-no mal poin s
(𝑦and 𝑧di ec ions) a e dis ibu ed ollowing a hype bolic- angen
unc ion as ollows (iden ical o he 𝑧-di ec ion)
𝑦𝑖=1
2(1 + anh (𝛾𝑦(2(𝑖− 1)∕𝑁𝑦− 1))
anh 𝛾𝑦), 𝑖= 1,…, 𝑁𝑦+ 1,(66)
whe e 𝑁𝑦and 𝛾𝑦a e he numbe o con ol olumes and he concen-
a ion ac o in he 𝑦-di ec ion, espec i ely. In ou case, 𝛾𝑦=𝛾𝑧= 1.4
o Ra = 108and 𝛾𝑦=𝛾𝑧= 1.6 o Ra = 1010. Fo u he de ails, he
eade is e e ed o ou p e ious wo ks [38,39].
Ins an aneous empe a u e ields co esponding o he s a is ically
s eady s a e a e displayed in Fig. 8. As expec ed, he mal dissipa ion
e ec s a Ge = 1lead o a signi ican inc ease in he a e age ca i y
empe a u e which is clea ly isible o bo h Rayleigh numbe s. As in
2D, he low symme y (in a e age sense) wi h espec o he mid-
heigh plane is los o Ge >0leading o highe (lowe ) Nussel numbe
o he op (bo om) wall. Subsequen ly, he op (bo om) he mal
bounda y laye becomes hinne ( hicke ) wi h espec o he case a
Ge = 0. This implies ha mesh esolu ion equi emen s in he nea -
wall egion a e also asymme ical; howe e , in his wo k, o he sake
o simplici y, he g id spacing a he wo walls is he same ega dless
o he Gebha numbe .
All simula ions ha e been ca ied ou o 500 ime-uni s s a ing
om a ze o eloci y ield and uni o mly dis ibu ed andom empe -
a u es be ween 𝑇𝐶and 𝑇𝐻. As he luid se s in mo ion, ini ially he
disc e e kine ic ene gy o he sys em inc eases. Then, a e a su icien ly
long pe iod o ime (a ound 50 ime-uni s) a s a is ically s eady s a e is
eached. This is clea ly obse ed in Fig. 9whe e he ime-e olu ion
o a ious a e-o -changes o ene gy a e shown. Resul s co espond o
Ra = 108and Ge = 1using a e y ine (400 × 208 × 208 ≈ 17.3M) and
a e y coa se mesh. Simila esul s a e ob ained o he o he es ed
con igu a ions. As expec ed, once a s a is ically s eady s a e is eached,
he kine ic ene gy luc ua es a ound i s mean alue and he e o e i s
a e-o -change d𝐸𝑘,ℎ∕d𝑡(in ed) luc ua es a ound ze o. Only wo e ms
con ibu e o he global kine ic ene gy o he sys em (see Eq. (34)):
he global iscous dissipa ion, 𝜖𝑢,ℎ (in yellow), and he con ibu ion o
he buoyancy o ces gi en by 𝛼2𝑉𝑇
ℎ(𝐴𝑇ℎ(𝑡) +𝑦𝑇)(in blue). These wo
con ibu ions cancel each o he on a e age when a s a is ically s eady
Compu e s and Fluids 286 (2025) 106473
9
B. Sande se and F.X. T ias
[8] Hewi JM, Mckenzie DP, Weiss NO. Dissipa i e hea ing in con ec-
i e lows. J Fluid Mech 1975;68(4):721–38. h p://dx.doi.o g/10.1017/
S002211207500119X.
[9] Tu co e DL, Hsui AT, To ance KE, Schube G. In luence o iscous dissipa ion
on Béna d con ec ion. J Fluid Mech 1974;64(2):369–74. h p://dx.doi.o g/10.
1017/S0022112074002448.
[10] Gebha B. E ec s o iscous dissipa ion in na u al con ec ion. J Fluid Mech
1962;14(2):225–32. h p://dx.doi.o g/10.1017/S0022112062001196.
[11] Os ach S. Lamina na u al-con ec ion low and hea ans e o luids wi h and
wi hou hea sou ces in channels wi h cons an wall empe a u es. Tech. ep.
NACA-TN-2863, Lewis Fligh P opulsion Lab., NACA; 1952.
[12] Blankenbach B, Busse F, Ch is ensen U, Cse epes L, Gunkel D, Hansen U,
Ha de H, Ja is G, Koch M, Ma qua G, Moo e D, Olson P, Schmeling H,
Schnaubel T. A benchma k compa ison o man le con ec ion codes. Geophys J
In 1989;98(1):23–38. h p://dx.doi.o g/10.1111/j.1365-246X.1989. b05511.x.
[13] King SD, Lee C, an Keken PE, Leng W, Zhong S, Tan E, Tosi N, Kameyama MC.
A communi y benchma k o 2-D ca esian comp essible con ec ion in he Ea h’s
man le. Geophys J In 2010;180(1):73–87. h p://dx.doi.o g/10.1111/j.1365-
246X.2009.04413.x.
[14] Sugiyama K, Calza a ini E, G ossmann S, Lohse D. Flow o ganiza ion in wo-
dimensional non-Obe beck-Boussinesq Rayleigh-Béna d con ec ion in wa e . J
Fluid Mech 2009;637:105–35.
[15] Najm HN, Wycko PS, Knio OM. A semi-implici nume ical scheme o eac ing
low: I. S i chemis y. J Compu Phys 1998;143(2):381–402. h p://dx.doi.o g/
10.1006/jcph.1997.5856.
[16] Nema i H. Di ec nume ical simula ion o u bulen hea ans e o luids a
supe c i ical p essu es (Ph.D. hesis), Del Uni e si y o Technology; 2016.
[17] Ba le a A. Commen s on a pa adox o iscous dissipa ion and i s ela ion o he
Obe beck–Boussinesq app oach. In J Hea Mass T ans e 2008;51(25–26):6312–
6. h p://dx.doi.o g/10.1016/j.ijhea mass ans e .2007.10.044.
[18] Ba le a A, Nield D. E ec o p essu e wo k and iscous dissipa ion in he analysis
o he Rayleigh–Béna d p oblem. In J Hea Mass T ans e 2009;52(13–14):3279–
89. h p://dx.doi.o g/10.1016/j.ijhea mass ans e .2009.02.005.
[19] Ba le a A, Celli M, Nield DA. On he onse o dissipa ion he mal ins abili y
o he Poiseuille low o a highly iscous luid in a ho izon al channel. J Fluid
Mech 2011;681:499–514. h p://dx.doi.o g/10.1017/j m.2011.213.
[20] Ba le a A, Celli M, B andão PV. On mixed con ec ion in a ho izon al channel,
iscous dissipa ion and low duali y. Fluids 2022;7(5):170. h p://dx.doi.o g/10.
3390/ luids7050170.
[21] Win e s KB, Lomba d PN, Riley JJ, D’Asa o EA. A ailable po en ial ene gy
and mixing in Densi y-S a i ied luids. J Fluid Mech 1995;289:115–28. h p:
//dx.doi.o g/10.1017/S002211209500125X.
[22] Hughes GO, Gayen B, G i i hs RW. A ailable po en ial ene gy in Rayleigh–
béna d con ec ion. J Fluid Mech 2013;729:R3. h p://dx.doi.o g/10.1017/j m.
2013.353.
[23] Gayen B, Hughes GO, G i i hs RW. Comple ing he mechanical ene gy pa hways
in u bulen Rayleigh-Béna d con ec ion. Phys Re Le 2013;111(12):124301.
h p://dx.doi.o g/10.1103/PhysRe Le .111.124301.
[24] Pe schel K, S ellmach S, Wilczek M, Lül J, Hansen U. Kine ic ene gy anspo
in Rayleigh–Béna d con ec ion. J Fluid Mech 2015;773:395–417. h p://dx.doi.
o g/10.1017/j m.2015.216.
[25] an Gils DPM, B ugge G-W, La h op DP, Sun C, Lohse D. The wen e u bulen
Taylo –Coue e (T3C) acili y: S ongly u bulen (mul iphase) low be ween
wo independen ly o a ing cylinde s. Re Sci Ins um 2011;82(2):025105. h p:
//dx.doi.o g/10.1063/1.3548924.
[26] Ba le a A. Local ene gy balance, speci ic hea s and he Obe beck–Boussinesq
app oxima ion. In J Hea Mass T ans e 2009;52(21–22):5266–70. h p://dx.
doi.o g/10.1016/j.ijhea mass ans e .2009.06.006.
[27] an de Poel EP, S e ens RJAM, Lohse D. Compa ison be ween wo- and h ee-
dimensional Rayleigh–Béna d con ec ion. J Fluid Mech 2013;736:177–94. h p:
//dx.doi.o g/10.1017/j m.2013.488.
[28] Hepwo h BJ. Nonlinea wo-dimensional Rayleigh-béna d con ec ion (Ph.D.
hesis), Uni e si y o Leeds; 2014.
[29] Sande se B. Ene gy-conse ing disc e iza ion me hods o he incomp essible
Na ie -S okes equa ions:applica ion o he simula ion o wind- u bine wakes
(Ph.D. hesis), Technische Uni e si ei Eindho en; 2013.
[30] Sugiyama K, Calza a ini E, G ossmann S, Lohse D. Flow o ganiza ion in wo-
dimensional non-Obe beck–Boussinesq Rayleigh–Béna d con ec ion in wa e . J
Fluid Mech 2009;637:105–35. h p://dx.doi.o g/10.1017/S0022112009008027.
[31] Ve s appen R, Veldman A. Symme y-p ese ing disc e iza ion o u bulen low.
J Compu Phys 2003;187(1):343–68. h p://dx.doi.o g/10.1016/S0021-9991(03)
00126-8.
[32] T ias FX, Lehmkuhl O. A sel -adap i e s a egy o he ime in eg a ion o Na ie -
S okes equa ions. Nume Hea T ans e B 2011;60(2):116–34. h p://dx.doi.o g/
10.1080/10407790.2011.594398.
[33] Sande se B. Ene gy-conse ing Runge–Ku a me hods o he incomp essible
Na ie –S okes equa ions. J Compu Phys 2013;233:100–31. h p://dx.doi.o g/
10.1016/j.jcp.2012.07.039.
[34] Shih TM, Tan CH, Hwang BC. E ec s o g id s agge ing on nume ical schemes.
In e na J Nume Me hods Fluids 1989;9(2):193–212. h p://dx.doi.o g/10.1002/
ld.1650090206.
[35] Cai W, Ma H, Wang Y, Chen J, Zheng X, Zhang H. De elopmen o POD educed-
o de model and i s closu e scheme o 2D Rayleigh–Béna d con ec ion. Appl
Ma h Model 2019;66:562–75. h p://dx.doi.o g/10.1016/j.apm.2018.09.031.
[36] Gel ga AY. Di e en modes o Rayleigh–Béna d Ins abili y in wo- and h ee-
dimensional ec angula enclosu es. J Compu Phys 1999;156(2):300–24. h p:
//dx.doi.o g/10.1006/jcph.1999.6363.
[37] Ven u i D, Wan X, Ka niadakis GE. S ochas ic bi u ca ion analysis o Rayleigh–
Béna d con ec ion. J Fluid Mech 2010;650:391–413. h p://dx.doi.o g/10.1017/
S0022112009993685.
[38] Dabbagh F, T ias FX, Go obe s A, Oli a A. On he e olu ion o low opology in
u bulen Rayleigh-Béna d con ec ion. Phys Fluids 2016;28:115105.
[39] Dabbagh F, T ias FX, Go obe s A, Oli a A. Flow opology dynamics in a h ee-
dimensional phase space o u bulen Rayleigh-Béna d con ec ion. Phys Re
Fluids 2020;5:024603.
[40] T ias FX, Ve s appen RWCP, Go obe s A, So ia M, Oli a A. Pa ame e - ee
symme y-p ese ing egula iza ion modeling o a u bulen di e en ially hea ed
ca i y. Compu & Fluids 2010;39:1815–31.
[41] Smi h W. All hings low. O egon S a e Uni e si y; 2019.
[42] Tailleux R. On he ene ge ics o s a i ied u bulen mixing, i e e sible he mo-
dynamics, Boussinesq models and he ocean hea engine con o e sy. J Fluid
Mech 2009;638:339–82. h p://dx.doi.o g/10.1017/S002211200999111X.
Compu e s and Fluids 286 (2025) 106473
16
