# Energy-consistent discretization of viscous dissipation with application to
natural convection flow

B. Sanderse  F.X. Trias

###### Abstract

A new energy-consistent discretization of the viscous dissipation function in
incompressible flows is proposed. It is implied by choosing a discretization
of the diffusive terms and a discretization of the local kinetic energy
equation and by requiring that continuous identities like the product rule are
mimicked discretely. The proposed viscous dissipation function has a
quadratic, strictly dissipative form, for both simplified (constant viscosity)
stress tensors and general stress tensors. The proposed expression is not only
useful in evaluating energy budgets in turbulent flows, but also in natural
convection flows, where it appears in the internal energy equation and is
responsible for viscous heating. The viscous dissipation function is such that
a consistent total energy balance is obtained: the ‘implied’ presence as sink
in the kinetic energy equation is exactly balanced by explicitly adding it as
source term in the internal energy equation.

Numerical experiments of Rayleigh-Bénard convection (RBC) and Rayleigh-Taylor
instabilities confirm that with the proposed dissipation function, the energy
exchange between kinetic and internal energy is exactly preserved. The
experiments show furthermore that viscous dissipation does not affect the
critical Rayleigh number at which instabilities form, but it does
significantly impact the development of instabilities once they occur.
Consequently, the value of the Nusselt number on the cold plate becomes larger
than on the hot plate, with the difference increasing with increasing Gebhart
number. Finally, 3D simulations of turbulent RBC show that energy balances are
exactly satisfied even for very coarse grids; therefore, we consider that the
proposed discretization forms an excellent starting point for testing sub-grid
scale models.

###### keywords:

viscous dissipation, energy conservation, staggered grid, natural convection,
Rayleigh-Bénard, Gebhart number

$\displaystyle{\linenomathAMS}{}{}\pretocmd\begin{gathered}{\linenomathAMS}{}{}\apptocmd\end{gathered}{\endlinenomath}{}{}$
${\linenomathAMS}{}{}\pretocmd\begin{multline*}{\linenomathAMS}{}{}\apptocmd\end{multline*}{\linenomathAMS}{}{}\apptocmd{\endlinenomath}{}{}\apptocmd$
(4) $\displaystyle{\linenomathAMS}{}{}\pretocmd\begin{aligned}
{\linenomathAMS}{}{}\apptocmd\end{aligned}{\endlinenomath}{}{}\apptocmd$ (5)
$\displaystyle{}{}\pretocmd\begin{aligned}
{}{}\apptocmd\end{aligned}{\endlinenomath}{}{}\apptocmd$ (6)
$\displaystyle{\linenomathAMS}{}{}\pretocmd\begin{aligned}
{\linenomathAMS}{}{}\apptocmd\end{aligned}{\endlinenomath}{}{}\apptocmd$ (7)

[1]organization=Centrum Wiskunde & Informatica,addressline=Science Park 123,
city=Amsterdam, country=The Netherlands [2]organization=Heat and Mass Transfer
Technological Center, Technical University of Catalonia, ESEIAAT,
addressline=c/ Colom 11, 08222 Terrassa (Barcelona), country=Spain

##  1 Introduction and problem description

In this article we study the viscous dissipation function and its role in
natural convection flows described by the incompressible Navier-Stokes
equations, with buoyancy effects modelled by the Boussinesq approximation [1].
These ‘Boussinesq‘ or ‘Oberbeck-Boussinesq’ equations have attracted much
scientific interest over several decades [2], not only because of their
physical relevance, but also of their intriguing mathematical properties. An
important test case studied with the Boussinesq system is that of Rayleigh-
Bénard convection [3], in which a box of fluid is heated from the bottom and
cooled from the top, giving rise to convection cells. The Boussinesq equations
also describe a (miscible) form of Rayleigh-Taylor instability, which occurs
when a heavy (cold) fluid is positioned above a light (warm) fluid.

A common assumption in many incompressible natural convection studies is that
the effect of viscous dissipation on the internal energy (effectively on the
temperature) is neglected. This assumption is not always valid, for example
when considering natural convection in the Earth mantle, when considering
highly viscous liquids, when large length scales are involved, or in devices
operating at high rotational speed [4, 5, 6, 7, 8, 9, 10, 11]. Of course, when
considering compressible flows, e.g. high-speed flows, including heating by
viscous dissipation is known to be important, and several benchmarking studies
have been performed related to modelling natural convection in the Earth
mantle [12, 13]. These studies typically assume infinite Prandtl numbers, and
ignore the unsteady and convective terms in the momentum equations. In this
paper we will restrict ourselves to the incompressible situation, for which
this effect is less studied. In this incompressible case, Ostrach [11],
Gebhart [10] and Turcotte et al. [9] should be explicitly mentioned, being
among the first to address the role of viscous dissipation and to introduce
next to the well-known Rayleigh and Prandtl numbers another dimensionless
quantity, which is known as the dissipation number or the Gebhart number. In
addition, we mention the work of Barletta and co-authors [14, 15, 16, 17], who
considered the role of viscous dissipation in natural convection in several
papers, studying the correct mathematical formulation of the problem and
linear stability analysis for different geometries. Turcotte et al. [9] were
probably one of the first to perform numerical experiments of incompressible
natural convection flows that include viscous dissipation. They performed
simulations on coarse grids ($10\times 10$) and low Rayleigh numbers
($\mathrm{Ra}=10^{4},10^{5}$) for different values of the dissipation number
and concluded that Rayleigh-Bénard convection was significantly affected when
the dissipation number was of order unity.

From an energy perspective, the viscous dissipation source term in the
internal energy equation occurs as a sink in the kinetic energy equation,
which cancel each other when considering the total energy equation. However,
most energy analyses, especially for incompressible flow, focus on the role of
the potential energy term and its split into available and background
potential energy [18, 19, 20], or on the kinetic energy budget [21]. To the
author’s knowledge, the role of viscous dissipation in the kinetic energy
equation and its numerical treatment for the internal energy equation have not
been explored in detail.

In this paper, the main novelty is that we propose a discretization of the
viscous dissipation function and apply it to the context of natural convection
flow, where it appears as a source term in the internal energy equation. Our
discretization is such that we get a correct global energy balance, on
continuous, semi-discrete, and fully discrete level. First, on the continuous
level, a non-dimensionalization is proposed that makes the internal and
kinetic energy scaling consistent. Second, on the semi-discrete level, we
propose a discrete dissipation operator, and show that it cannot be chosen
freely but is implied by the discretization of the viscous terms in the
momentum equations and by the definition of the kinetic energy. This discrete
dissipation operator is not only of use in the internal energy equation, but
also useful beyond the context of natural convection flows, e.g. when
estimating the dissipation of kinetic energy in turbulent flows in a numerical
simulation. Third, on the fully discrete level, we propose a time integration
method that preserves the total energy balance upon time marching.

The paper is structured as follows. Section 2 introduces the governing
equations, energy balances, and new non-dimensionalization. Sections 3 and 4
describe the energy-consistent spatial and temporal discretization. Section 5
describes steady-state results of Rayleigh-Bénard convection including viscous
dissipation, and section 6 describes energy-conserving simulations of
Rayleigh-Taylor instabilities including viscous dissipation. Section 7 shows
the effect of viscous dissipation in 3D DNS of Rayleigh-Bénard convection.

##  2 Energy-conserving formulation

###  2.1 Governing equations

The Boussinesq approximation states that density variations are small and can
be ignored in all terms of the Navier-Stokes (NS) equations, except in the one
pertaining to the gravity term. The NS equations describing conservation of
mass and momentum then read

$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0,$ (8)
$\displaystyle\rho_{0}\left(\frac{\partial\bm{u}}{\partial
t}+\nabla\cdot\left(\bm{u}\otimes\bm{u}\right)\right)$ $\displaystyle=-\nabla
p+\mu\nabla^{2}\bm{u}+\rho\bm{g},$ (9)

where $\bm{u}(\bm{x},t)$ is the velocity field, $p(\bm{x},t)$ the pressure,
$\mu$ the dynamic viscosity, $\rho(\bm{x},t)$ the density and $\rho_{0}$ a
reference density. Without loss of generality, we consider a two-dimensional
domain $\Omega$, with the gravity vector pointing in the negative
$y$-direction so that $\bm{g}=-g\bm{e}_{y}$. An example of the domain as used
in the Rayleigh-Bénard problem, including the boundary conditions, is given in
figure 1. In the results section we will also consider the Rayleigh-Taylor
problem, which has adiabatic boundaries on top and bottom, instead of
isothermal as in case of Rayleigh-Bénard.

$L$
---
$H$
---
$\frac{\partial T}{\partial x}=0$
---
$\bm{u}=0$
---
$\frac{\partial T}{\partial x}=0$
---
$\bm{u}=0,T=1$
---
$\bm{u}=0,T=0$
---
$y$
---
$x$
---
$\bm{u}=0$
---
Figure 1: Problem set-up for Rayleigh-Bénard convection.

The density $\rho$ is assumed to vary only with temperature $T(\bm{x},t)$,
according to $\rho(T)=\rho_{0}-\beta\rho_{0}(T-T_{0})$, where $\beta$ is the
isobaric coefficient of thermal expansion
($\beta=-\frac{1}{\rho}\left(\frac{\partial\rho}{\partial T}\right)_{p}$). The
NS equations are then written as

$\rho_{0}\left(\frac{\partial\bm{u}}{\partial
t}+\nabla\cdot\left(\bm{u}\otimes\bm{u}\right)\right)=-\nabla
p^{\prime}+\mu\nabla^{2}\bm{u}-\beta\rho_{0}(T-T_{0})\bm{g},$ (20)

where $p=p^{\prime}-\rho_{0}gy$ and $\nabla p=\nabla
p^{\prime}-\rho_{0}g\bm{e}_{y}$.

The equation for the internal energy $e_{i}$ describes the temperature
evolution according to

$\frac{\partial}{\partial
t}(\underbrace{\rho_{0}cT}_{e_{i}})+\nabla\cdot(\bm{u}(\rho_{0}cT))=\Phi+\lambda\nabla^{2}T,$
(21)

where $\lambda$ is the thermal conductivity and $c$ equals $c_{v}$ in case of
an ideal gas (the specific heat at constant volume), and equals
$c_{p}-\frac{p\beta}{\rho}$ for a real gas [22]). The contribution of pressure
work to the change in internal energy, $p\nabla\cdot\bm{u}$, has been
discarded in equation (20) because of equation (8). The viscous dissipation
function

$\Phi:=\mu\|\nabla\bm{u}\|^{2}\geq 0,$ (22)

is the key quantity in this work, where
$\|\nabla\bm{u}\|^{2}=\nabla\bm{u}:\nabla\bm{u}$ (the Frobenius inner
product). In 2D and Cartesian coordinates it can be written as

$\Phi=\mu\left[\left(\frac{\partial u}{\partial
x}\right)^{2}+\left(\frac{\partial u}{\partial
y}\right)^{2}+\left(\frac{\partial v}{\partial
x}\right)^{2}+\left(\frac{\partial v}{\partial y}\right)^{2}\right].$ (23)

######  Remark 1.

The viscous dissipation expression (22) is only valid if the diffusive terms
are written as $\nabla\cdot\bm{\tau}=\mu\nabla^{2}\bm{u}$, with
$\bm{\tau}=\nabla\bm{u}$. This simplified form follows from the more general
stress tensor expression

$\nabla\cdot\hat{\bm{\tau}},\qquad\hat{\bm{\tau}}=\mu(\nabla\bm{u}+(\nabla\bm{u})^{T}),$
(24)

by assuming constant $\mu$ and incompressibility, so that the identity
$\nabla\cdot(\mu\nabla\bm{u})^{T}=\mu\nabla(\nabla\cdot\bm{u})=0$ holds. For
the more general stress tensor $\hat{\bm{\tau}}$ we have instead

$\hat{\Phi}=\hat{\bm{\tau}}:\nabla\bm{u}=\mu\left[2\left(\frac{\partial
u}{\partial x}\right)^{2}+2\left(\frac{\partial v}{\partial
y}\right)^{2}+\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial
x}\right)^{2}\right],$ (25)

which is notably different from expression (23). The dissipation function is
therefore linked to the form of the stress tensor used in the momentum
equations. In this work we mainly use expression (23), but we will also
explain the extension to the more general form (25), see B, equation (236).
For an alternative approach, see also [23].

###  2.2 Total energy conservation

Conservation of kinetic energy follows by taking the dot product of equation
(20) with $\bm{u}$:

$\frac{\partial}{\partial
t}(\underbrace{\frac{1}{2}\rho_{0}|\bm{u}|^{2}}_{e_{k}})+\nabla\cdot(\frac{1}{2}\rho_{0}|\bm{u}|^{2}\bm{u})=-\bm{u}\cdot\nabla
p^{\prime}+\mu\nabla\cdot(\bm{u}\cdot\nabla\bm{u})-\mu\|\nabla\bm{u}\|^{2}+\beta
g\rho_{0}(T-T_{0})v,$ (26)

where $\bm{g}\cdot\bm{u}=-gv$ and we have used the identity

$\bm{u}\cdot\nabla^{2}\bm{u}=-\|\nabla\bm{u}\|^{2}+\nabla\cdot(\bm{u}\cdot\nabla\bm{u}).$
(27)

Upon adding the kinetic and internal energy equations (21) and (26), the
viscous dissipation term cancels and we arrive at the equation for the total
energy $e=e_{k}+e_{i}$:

$\frac{\partial}{\partial
t}(e_{k}+e_{i})+\nabla\cdot((e_{k}+e_{i})\bm{u})=-\nabla\cdot(p^{\prime}\bm{u})+\mu\nabla\cdot(\bm{u}\cdot\nabla\bm{u})+\beta
g\rho_{0}(T-T_{0})v+\lambda\nabla^{2}T.$ (28)

All terms are in conservative (divergence) form, except the potential energy
term. Upon integrating over the domain $\Omega$ and assuming no-slip
conditions $\bm{u}=\bm{0}$ on all boundaries, we obtain the global balances

$\displaystyle\frac{\mathrm{d}E_{k}}{\mathrm{d}t}$
$\displaystyle=-\int_{\Omega}\Phi\mathrm{d}\Omega+\int_{\Omega}\beta
g\rho_{0}(T-T_{0})v\mathrm{d}\Omega,$ (29)
$\displaystyle\frac{\mathrm{d}E_{i}}{\mathrm{d}t}$
$\displaystyle=\int_{\Omega}\Phi\mathrm{d}\Omega+\int_{\partial\Omega}\lambda\nabla
T\cdot\bm{n}\,\mathrm{d}S,$ (30)
$\displaystyle\frac{\mathrm{d}E}{\mathrm{d}t}$
$\displaystyle=\frac{\mathrm{d}E_{k}}{\mathrm{d}t}+\frac{\mathrm{d}E_{i}}{\mathrm{d}t}=\int_{\Omega}\beta
g\rho_{0}(T-T_{0})v\mathrm{d}\Omega+\int_{\partial\Omega}\lambda\nabla
T\cdot\bm{n}\,\mathrm{d}S,$ (31)

where $E=\int_{\Omega}e\,\mathrm{d}\Omega=E_{k}+E_{i}$. In case the boundary
conditions are adiabatic ($\nabla T\cdot\bm{n}=0$), the last term in (31)
vanishes and the total energy equation expresses that the sum of internal and
kinetic energy changes due to the buoyancy flux $\int_{\Omega}\beta
g\rho_{0}(T-T_{0})v\,\mathrm{d}\Omega$ – this case will be dealt with in the
Rayleigh-Taylor set-up in section 6.

In most studies of Rayleigh-Bénard convection the dissipation function $\Phi$
is left out from the internal energy equation (21), while its corresponding
counterpart in the momentum equation ($\mu\nabla^{2}\bm{u}$) is still
included. As a consequence, the energy lost in the kinetic energy equation is
not balanced by the heat generated in the internal energy equation, so that
the total energy equation features a dissipation term, which destroys the
global energy balance.

######  Remark 2.

Equations (29) and (31) feature the buoyancy flux $\int\beta
g\rho_{0}(T-T_{0})v\mathrm{d}\Omega$ (stemming from the term
$\int\rho\bm{g}\cdot\bm{u}\mathrm{d}\Omega$). In general compressible fluids,
i.e. those that satisfy

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\bm{u})=0,$ (32)

one can show that the buoyancy flux can be written as the time derivative of
the potential energy $E_{p}=\int\rho gy\mathrm{d}\Omega$ (see [24], section
6.4.2; [7], section 3.8). In that case, one could define
$\hat{E}=E_{k}+E_{i}+E_{p}$ and have a total energy conservation statement of
the form [25]:

$\frac{\mathrm{d}\hat{E}}{\mathrm{d}t}=\int_{\partial\Omega}\lambda\nabla
T\cdot\bm{n}\,\mathrm{d}S.$ (33)

However, in Boussinesq fluids, equation (32) is not satisfied; instead we have

$\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\bm{u})=\frac{\partial\rho}{\partial
t}+(\bm{u}\cdot\nabla)\rho+\underbrace{\rho\nabla\cdot\bm{u}}_{=0}=-\rho_{0}\beta\left[\frac{\partial
T}{\partial t}+(\bm{u}\cdot\nabla)T\right].$ (34)

The right-hand side can be written in terms of the sum of thermal diffusion
and viscous dissipation, see equation (21), and is generally nonzero. As a
consequence, the time derivative of the potential energy includes not only the
buoyancy flux, but also additional terms [18, 19]. Therefore, for Boussinesq
fluids additional terms appear in the right-hand side of equation (33)
(independent of whether viscous dissipation is included in the internal energy
equation). In this paper, the meaning ‘energy-consistent’ thus refers to the
exchange between internal and kinetic energy, and not to the total energy
(kinetic + potential + internal), which is not conserved under the Boussinesq
approximation.

###  2.3 Non-dimensionalization

The study of the Rayleigh-Bénard convection problem is simplified by
introducing dimensionless quantities. As explained in [4], p. 46, three
dimensionless groups (or ‘similarity parameters’) are needed to fully describe
the problem. An important question that we address here is how the choice of
non-dimensionalization changes the total energy equation.

We non-dimensionalize equations (8), (9) and (21) by taking a reference length
$H$ (cavity height), a reference temperature difference $\Delta T$ (difference
between the cold and hot plates), and a reference velocity yet to be
specified. From these choices we find the time scale $H/u_{\mathrm{ref}}$ and
the pressure scale $\rho_{0}u_{\mathrm{ref}}^{2}$. The non-dimensional
quantities are thus

$\displaystyle\bm{\tilde{x}}$ $\displaystyle=\frac{\bm{x}}{H},$
$\displaystyle\tilde{t}$ $\displaystyle=\frac{tu_{\mathrm{ref}}}{H},$
$\displaystyle\bm{\tilde{u}}$ $\displaystyle=\frac{\bm{u}}{u_{\mathrm{ref}}},$
$\displaystyle\tilde{T}$ $\displaystyle=\frac{T-T_{0}}{\Delta T},$
$\displaystyle\tilde{p}^{\prime}$
$\displaystyle=\frac{p^{\prime}}{\rho_{0}u_{\mathrm{ref}}^{2}},$ (35)

and the non-dimensional equations read

$\displaystyle\tilde{\nabla}\cdot\bm{\tilde{u}}$ $\displaystyle=0,$ (36)
$\displaystyle\frac{\partial\bm{\tilde{u}}}{\partial\tilde{t}}+\tilde{\nabla}\cdot\left(\bm{\tilde{u}}\otimes\bm{\tilde{u}}\right)$
$\displaystyle=-\tilde{\nabla}\tilde{p}^{\prime}+\frac{\mu}{\rho_{0}u_{\mathrm{ref}}H}\tilde{\nabla}^{2}\bm{\tilde{u}}+\frac{\beta
g\Delta TH}{u_{\mathrm{ref}}^{2}}\tilde{T}\bm{e}_{y},$ (37)
$\displaystyle\frac{\partial\tilde{T}}{\partial\tilde{t}}+\tilde{\nabla}\cdot(\bm{\tilde{u}}\tilde{T})$
$\displaystyle=\frac{\nu u_{\mathrm{ref}}}{cH\Delta
T}\tilde{\Phi}+\frac{\kappa}{u_{\mathrm{ref}}H}\tilde{\nabla}^{2}\tilde{T},$
(38)

where $\nu=\mu/\rho_{0}$ and $\kappa=\lambda/(\rho_{0}c)$. The two latter
equations are re-written by introducing the parameters $\alpha_{i}$,
$i=1\ldots 4$, as

$\displaystyle\frac{\partial\bm{\tilde{u}}}{\partial\tilde{t}}+\tilde{\nabla}\cdot\left(\bm{\tilde{u}}\otimes\bm{\tilde{u}}\right)$
$\displaystyle=-\tilde{\nabla}\tilde{p}^{\prime}+\alpha_{1}\tilde{\nabla}^{2}\bm{\tilde{u}}+\alpha_{2}\tilde{T}\bm{e}_{y},$
(39)
$\displaystyle\frac{\partial\tilde{T}}{\partial\tilde{t}}+\tilde{\nabla}\cdot(\bm{\tilde{u}}\tilde{T})$
$\displaystyle=\alpha_{3}\tilde{\Phi}+\alpha_{4}\tilde{\nabla}^{2}\tilde{T}.$
(40)

The $\alpha_{i}$’s can be expressed in terms of three dimensionless numbers,
being the Rayleigh number $\mathrm{Ra}$, the Prandtl number $\Pr$ and the
Gebhart number $\mathrm{Ge}$ (also known as the dissipation number [6]):

$\displaystyle\mathrm{Ra}$ $\displaystyle=\frac{\beta g\Delta
TH^{3}}{\nu\kappa},$ (41) $\displaystyle\mathrm{Pr}$
$\displaystyle=\frac{\nu}{\kappa},$ (42) $\displaystyle\mathrm{Ge}$
$\displaystyle=\frac{\beta gH}{c}.$ (43)

Alternatively, one can employ the Grashof number $\text{Gr}=\mathrm{Ra}/\Pr$
[4]. In table 1 we present three different options for $u_{\mathrm{ref}}$ with
the corresponding values of $\alpha$. Choices I and II are common in
literature, see for example [26] for choice I and [15, 1, 27] for choice II.
Other choices are also possible, e.g. $u_{\mathrm{ref}}=\beta g\Delta
TH^{2}/\nu$ [5], but this choice does not lead to a ‘clean’ expression in
terms of the dimensionless numbers defined above. To our best knowledge,
choice III is new and inspired by the form of the total energy equation, as we
will explain below.

| $u_{\mathrm{ref}}$ | $\alpha_{1}=\frac{\nu}{u_{\mathrm{ref}}H}$ | $\alpha_{2}=\frac{\beta g\Delta TH}{u_{\mathrm{ref}}^{2}}$ | $\alpha_{3}=\frac{\nu u_{\mathrm{ref}}}{c\Delta TH}$ | $\alpha_{4}=\frac{\kappa}{u_{\mathrm{ref}}H}$ | $\gamma=\frac{\alpha_{1}}{\alpha_{3}}$
---|---|---|---|---|---|---
I | $\sqrt{\beta g\Delta TH}$ | $\sqrt{\frac{\mathrm{Pr}}{\mathrm{Ra}}}$ | 1 | $\mathrm{Ge}\sqrt{\frac{\mathrm{Pr}}{\mathrm{Ra}}}$ | $\frac{1}{\sqrt{\mathrm{Pr}\mathrm{Ra}}}$ | $\frac{1}{\mathrm{Ge}}$
II | $\frac{\kappa}{H}$ | $\mathrm{Pr}$ | $\mathrm{Pr}\mathrm{Ra}$ | $\frac{\mathrm{Ge}}{\mathrm{Ra}}$ | 1 | $\frac{\mathrm{Pr}\mathrm{Ra}}{\mathrm{Ge}}$
III | $\sqrt{c\Delta T}$ | $\sqrt{\frac{\mathrm{Pr}\mathrm{Ge}}{\mathrm{Ra}}}$ | $\mathrm{Ge}$ | $\sqrt{\frac{\mathrm{Pr}\mathrm{Ge}}{\mathrm{Ra}}}$ | $\sqrt{\frac{\mathrm{Ge}}{\mathrm{Pr}\mathrm{Ra}}}$ | 1
Table 1: Different non-dimensional forms resulting from different choices of
$u_{\mathrm{ref}}$.

It is important to realize that the time scales and the velocity fields
corresponding to numerical simulations with choices I, II and III are
different. The time scales are related as
$\frac{\tilde{t}_{I}}{u_{\text{ref},{I}}}=\frac{\tilde{t}_{II}}{u_{\text{ref},{II}}}=\frac{\tilde{t}_{III}}{u_{\text{ref},{III}}}$,
so $\tilde{t}_{III}=\tilde{t}_{I}/\sqrt{\mathrm{Ge}}$ and
$\tilde{t}_{III}=\tilde{t}_{II}\sqrt{\mathrm{Ra}\Pr}/\sqrt{\mathrm{Ge}}$. The
velocity fields are related as
$\bm{\tilde{u}}_{I}u_{\text{ref},{I}}=\bm{\tilde{u}}_{II}u_{\text{ref},{II}}=\bm{\tilde{u}}_{III}u_{\text{ref},{III}}$,
so that $\bm{\tilde{u}}_{III}=\bm{\tilde{u}}_{I}\sqrt{\mathrm{Ge}}$, and
$\bm{\tilde{u}}_{III}=\bm{\tilde{u}}_{II}\sqrt{\mathrm{Ge}}/\sqrt{\mathrm{Ra}\Pr}$.
On the other hand, the temperature fields corresponding to each choice are
equivalent, and consequently the Nusselt numbers are the same.

To obtain the non-dimensional form of the total energy equation we take the
dot product of (39) with $\bm{\tilde{u}}$ and add the internal energy equation
(40). In order for the dissipation function of the kinetic energy equation to
cancel with the internal energy equation, we require $\alpha_{1}=\alpha_{3}$.
This requirement is satisfied by $u_{\mathrm{ref}}=\sqrt{c\Delta T}$, i.e. our
proposed choice III in table 1. For the other choices (I and II), a weighting
of the kinetic and internal energy equations is needed in order to cancel the
dissipation function in the non-dimensional total energy equation. The
weighting factor depends on the definition of the non-dimensional total
energy. First define the dimensionless kinetic and internal energy as

$\displaystyle\tilde{e}_{k}$
$\displaystyle:=\frac{e_{k}}{\rho_{0}u_{\mathrm{ref}}^{2}}=\frac{\frac{1}{2}\rho_{0}|\bm{u}|^{2}}{\rho_{0}u_{\mathrm{ref}}^{2}}=\frac{\frac{1}{2}\rho_{0}u_{\mathrm{ref}}^{2}|\bm{\tilde{u}}|^{2}}{\rho_{0}u_{\mathrm{ref}}^{2}}=\frac{1}{2}|\bm{\tilde{u}}|^{2},$
(48) $\displaystyle\tilde{e}_{i}$ $\displaystyle:=\frac{e_{i}}{\rho_{0}c\Delta
T}=\frac{\rho_{0}cT}{\rho_{0}c\Delta T}=\frac{\rho_{0}c\Delta
T(\tilde{T}+T_{0}/\Delta T)}{\rho_{0}c\Delta T}=(\tilde{T}+T_{0}/\Delta T),$
(49)

so that

$e=e_{k}+e_{i}=\rho_{0}u_{\mathrm{ref}}^{2}\tilde{e}_{k}+\rho_{0}c\Delta
T\tilde{e}_{i}=\rho_{0}u_{\mathrm{ref}}^{2}\left(\tilde{e}_{k}+\frac{c\Delta
T}{u_{\mathrm{ref}}^{2}}\tilde{e}_{i}\right).$ (50)

By choosing the non-dimensional total energy as
$\tilde{e}=e/\rho_{0}u_{\mathrm{ref}}^{2}$, we obtain

$\tilde{e}=\tilde{e}_{k}+\frac{c\Delta
T}{u_{\mathrm{ref}}^{2}}\tilde{e}_{i}=\tilde{e}_{k}+\frac{\alpha_{1}}{\alpha_{3}}\tilde{e}_{i}=\tilde{e}_{k}+\gamma\tilde{e}_{i}.$
(51)

Here $\gamma=\frac{\alpha_{1}}{\alpha_{3}}$ is the weighting factor, which is
reported in table 1 for different choices of $u_{\mathrm{ref}}$. The global
energy balances in non-dimensional form read

$\displaystyle\frac{\mathrm{d}\tilde{E}_{k}}{\mathrm{d}\tilde{t}}$
$\displaystyle=-\frac{\alpha_{1}}{\Lambda}\int_{\tilde{\Omega}}\tilde{\Phi}\,\mathrm{d}\tilde{\Omega}+\frac{\alpha_{2}}{\Lambda}\int_{\tilde{\Omega}}\tilde{T}\tilde{v}\,\mathrm{d}\tilde{\Omega},$
(52) $\displaystyle\frac{\mathrm{d}\tilde{E}_{i}}{\mathrm{d}\tilde{t}}$
$\displaystyle=\frac{\alpha_{3}}{\Lambda}\int_{\tilde{\Omega}}\tilde{\Phi}\,\mathrm{d}\tilde{\Omega}+\frac{\alpha_{4}}{\Lambda}\int_{\partial\tilde{\Omega}}\tilde{\nabla}\tilde{T}\cdot\bm{n}\,\mathrm{d}\tilde{S},$
(53) $\displaystyle\frac{\mathrm{d}\tilde{E}}{\mathrm{d}\tilde{t}}$
$\displaystyle=\frac{\mathrm{d}\tilde{E}_{k}}{\mathrm{d}t}+\gamma\frac{\mathrm{d}\tilde{E}_{i}}{\mathrm{d}t}=\frac{\alpha_{2}}{\Lambda}\int_{\tilde{\Omega}}\tilde{T}\tilde{v}\,\mathrm{d}\tilde{\Omega}+\frac{\gamma\alpha_{4}}{\Lambda}\int_{\partial\tilde{\Omega}}\tilde{\nabla}\tilde{T}\cdot\bm{n}\,\mathrm{d}\tilde{S},$
(54)

where we define
$\tilde{E}=\frac{1}{\Lambda}\int_{\tilde{\Omega}}\tilde{e}\,\mathrm{d}\tilde{\Omega}$,
and $\Lambda=L/H$ is the aspect ratio of the box.

The choice for a particular reference velocity typically depends on the
problem at hand. Choices I and II have the advantage that in case of
$\mathrm{Ge}=0$ (most commonly investigated in literature), one obtains
$\alpha_{3}=0$ and the dissipation terms simply drops from the internal energy
equation. However, when $\mathrm{Ge}$ is small but nonzero, the weight factor
$\gamma$ becomes very large for choices I and II. Choice III does not suffer
from this issue, because $\gamma=1$ independent of $\mathrm{Ge}$. However,
choice III has the disadvantage that it does not work in the case
$\mathrm{Ge}=0$, since it leads to $\alpha_{i}=0$ for all $i$. In summary: for
$\mathrm{Ge}=0$, choices I and II are preferred; for small but nonzero
$\mathrm{Ge}$, choice III is preferred; in other cases, all choices are fine.

###  2.4 Effect of viscous dissipation on Nusselt number and thermal
dissipation

A main quantity of interest in natural convection flows is the Nusselt number
$\mathrm{Nu}$ and we will investigate how it changes upon including viscous
dissipation in the internal energy equation. First, define the average of the
sum of convective and conductive fluxes through a horizontal plane
$y=y^{\prime}$ by

$\overline{F}(y^{\prime}):=\frac{1}{L}\int_{0}^{L}\left(\rho_{0}cTv-\lambda\frac{\partial
T}{\partial y}\right)_{(x,y^{\prime})}\,\mathrm{d}x.$ (55)

Then, the Nusselt number based on $\overline{F}$ follows as [3]:

$\mathrm{Nu}(\tilde{y}^{\prime}):=\frac{\overline{F}(y^{\prime})}{\lambda\Delta
T/H}=\frac{1}{\Lambda}\int_{0}^{\Lambda}\left(\frac{1}{\alpha_{4}}\tilde{T}\tilde{v}-\frac{\partial\tilde{T}}{\partial\tilde{y}}\right)_{(\tilde{x},\tilde{y}^{\prime})}\,\mathrm{d}\tilde{x}.$
(56)

For steady state or statistically steady state (using a suitable average), and
in the absence of viscous dissipation, it is straightforward to show from the
internal energy equation that
$\mathrm{Nu}(\tilde{y})=\mathrm{Nu}(\tilde{y}=0)=\mathrm{Nu}$, which is a
constant, independent of $\tilde{y}^{\prime}$ [1, 27]. However, upon including
viscous dissipation, this relation no longer holds true and instead the steady
internal energy equation yields

$\alpha_{4}(\mathrm{Nu}(\tilde{y}^{\prime})-\mathrm{Nu}(0))=\alpha_{3}\epsilon_{U}(\tilde{y}^{\prime}),$
(57)

where the integrated dissipation function is given by

$\epsilon_{U}(\tilde{y}^{\prime}):=\frac{1}{\Lambda}\int_{0}^{\tilde{y}^{\prime}}\int_{0}^{\Lambda}\tilde{\Phi}\,\mathrm{d}\tilde{x}\mathrm{d}\tilde{y}.$
(58)

Equation (57) is an important relation which shows that (taking
$\tilde{y}^{\prime}=1$)

$\alpha_{4}(\mathrm{Nu}(1)-\mathrm{Nu}(0))=\alpha_{3}\epsilon_{U}(1),$ (59)

so the Nusselt number of the upper plate is always larger than or equal to the
Nusselt number of the lower plate.

A second relation between Nusselt number and viscous dissipation can be
obtained from the global kinetic energy balance, equation (52). The second
term in the right-hand side of equation (52) can be rewritten with equation
(57), following the analysis in [27]:

$\begin{split}\frac{\alpha_{2}}{\Lambda}\int_{\tilde{\Omega}}\tilde{T}\tilde{v}\,\mathrm{d}\tilde{\Omega}&=\frac{\alpha_{2}}{\Lambda}\int_{0}^{1}\int_{0}^{\Lambda}\tilde{T}\tilde{v}\,\mathrm{d}\tilde{x}\mathrm{d}\tilde{y}=\alpha_{2}\alpha_{4}\int_{0}^{1}\mathrm{Nu}(\tilde{y})\,\mathrm{d}\tilde{y}+\frac{\alpha_{2}\alpha_{4}}{\Lambda}\int_{0}^{\Lambda}\int_{0}^{1}\frac{\partial\tilde{T}}{\partial\tilde{y}}\,\mathrm{d}\tilde{y}\mathrm{d}\tilde{x}\\\
&=\alpha_{2}\alpha_{4}\mathrm{Nu}(0)+\alpha_{2}\alpha_{3}\int_{0}^{1}\epsilon_{U}(\tilde{y})\,\mathrm{d}\tilde{y}+\frac{\alpha_{2}\alpha_{4}}{\Lambda}\int_{0}^{\Lambda}(\tilde{T}(\tilde{x},\tilde{y}=1)-\tilde{T}(\tilde{x},\tilde{y}=0))\mathrm{d}\tilde{x}\\\
&=\alpha_{2}\alpha_{4}(\mathrm{Nu}(0)-1)+\alpha_{2}\alpha_{3}\int_{0}^{1}\epsilon_{U}(\tilde{y})\,\mathrm{d}\tilde{y}.\end{split}$
(60)

For (statistically) steady flow, this term equals the first term in the right-
hand side of equation (52), yielding the second relation between the Nusselt
number and the viscous dissipation $\epsilon_{U}$

$\alpha_{2}\alpha_{4}(\mathrm{Nu}(0)-1)=\alpha_{1}\epsilon_{U}(1)-\alpha_{2}\alpha_{3}\int_{0}^{1}\epsilon_{U}(\tilde{y})\,\mathrm{d}\tilde{y}.$
(61)

We recognize the well-known equation
$\alpha_{2}\alpha_{4}(\mathrm{Nu}(0)-1)=\alpha_{1}\epsilon_{U}(1)$, see e.g.
[1], but with the additional negative term
$-\alpha_{2}\alpha_{3}\int_{0}^{1}\epsilon_{U}(\tilde{y})\,\mathrm{d}\tilde{y}$.

Lastly, we link the thermal dissipation $\epsilon_{T}$ to the Nusselt number
and the viscous dissipation function. The non-dimensional internal energy
equation, equation (40), is multiplied by $\tilde{T}$, and after integrating
by parts, using the skew-symmetry of the convective operator, and employing
the boundary condition $\tilde{T}(\tilde{y}=1)=0$, one obtains

$\frac{1}{\Lambda}\frac{\mathrm{d}}{\mathrm{d}t}\int_{\tilde{\Omega}}\frac{1}{2}\tilde{T}^{2}\,\mathrm{d}\tilde{\Omega}=\frac{\alpha_{3}}{\Lambda}\int_{\tilde{\Omega}}\tilde{T}\tilde{\Phi}\,\mathrm{d}\tilde{\Omega}-\frac{\alpha_{4}}{\Lambda}\int_{0}^{\Lambda}\left(\tilde{T}\frac{\partial\tilde{T}}{\partial\tilde{y}}\right)_{\tilde{y}=0}\,\mathrm{d}\tilde{x}-\frac{\alpha_{4}}{\Lambda}\int_{\tilde{\Omega}}\|\tilde{\nabla}\tilde{T}\|^{2}\,\mathrm{d}\tilde{\Omega}.$
(62)

With the boundary condition $\tilde{T}(\tilde{y}=0)=1$, and the assumption of
(statistically) steady flow, this relation is further simplified to

$\alpha_{4}\mathrm{Nu}(0)=\alpha_{4}\epsilon_{T}-\frac{\alpha_{3}}{\Lambda}\int_{\tilde{\Omega}}\tilde{T}\tilde{\Phi}\,\mathrm{d}\tilde{\Omega},$
(63)

where

$\epsilon_{T}:=\frac{1}{\Lambda}\int_{\tilde{\Omega}}\|\tilde{\nabla}\tilde{T}\|^{2}\,\mathrm{d}\tilde{\Omega}.$
(64)

Since $\tilde{T}\geq 0$, $\tilde{\Phi}\geq 0$, we conclude that viscous
dissipation lowers the Nusselt number of the lower plate. In absence of
viscous dissipation in the internal energy equation, one obtains the familiar
relation $\mathrm{Nu}=\epsilon_{T}$. In combination with equation (59), we
obtain for the Nusselt number of the upper plate:

$\alpha_{4}\mathrm{Nu}(1)=\alpha_{4}\epsilon_{T}+\frac{\alpha_{3}}{\Lambda}\int(1-\tilde{T})\tilde{\Phi}\,\mathrm{d}\tilde{\Omega}.$
(65)

Assuming that the temperature satisfies $0\leq\tilde{T}\leq 1$, we find that
viscous dissipation increases the Nusselt number of the upper plate. In other
words, the thermal dissipation lies in between the two Nusselt numbers:

$\mathrm{Nu}(0)\leq\epsilon_{T}\leq\mathrm{Nu}(1).$ (66)

The three relations (59), (61) and (63) are summarized in table 2 and will be
confirmed in the numerical experiments in section 5.

origin | without viscous dissipation | with viscous dissipation
---|---|---
internal | $\mathrm{Nu}(1)=\mathrm{Nu}(0)$ | $\alpha_{4}(\mathrm{Nu}(1)-\mathrm{Nu}(0))=\alpha_{3}\epsilon_{U}(1)$
kinetic | $\alpha_{2}\alpha_{4}(\mathrm{Nu}(0)-1)=\alpha_{1}\epsilon_{U}(1)$ | $\alpha_{2}\alpha_{4}(\mathrm{Nu}(0)-1)=\alpha_{1}\epsilon_{U}(1)-\alpha_{2}\alpha_{3}\int_{0}^{1}\epsilon_{U}(\tilde{y})\,\mathrm{d}\tilde{y}$
internal energy $\times T$ | $\mathrm{Nu}(0)=\epsilon_{T}$ | $\alpha_{4}\mathrm{Nu}(0)=\alpha_{4}\epsilon_{T}-\frac{\alpha_{3}}{\Lambda}\int_{\tilde{\Omega}}\tilde{T}\tilde{\Phi}\,\mathrm{d}\tilde{\Omega}$
Table 2: Steady-state Nusselt number relations, with and without viscous
dissipation.

##  3 Energy-consistent spatial discretization

###  3.1 Mass, momentum and kinetic energy equation

To discretize the non-dimensional mass and momentum equations (36) and (39),
we use the staggered-grid energy-conserving finite volume method described in
[28], extended by including the buoyancy term in the momentum equations. This
leads to the following semi-discrete equations:

$\displaystyle MV_{h}(t)$ $\displaystyle=0,$ (71)
$\displaystyle\Omega_{V}\frac{\mathrm{d}V_{h}(t)}{\mathrm{d}t}$
$\displaystyle=-C_{V}(V_{h}(t))-Gp_{h}(t)+\alpha_{1}D_{V}V_{h}(t)+\alpha_{2}(AT_{h}(t)+y_{T}).$
(72)

Here, $V_{h}\in\mathbb{R}^{N_{V}}$ are the velocity unknowns,
$p_{h}\in\mathbb{R}^{N_{p}}$ the pressure unknowns, and
$T_{h}\in\mathbb{R}^{N_{p}}$ the temperature unknowns; see figure 2 for their
positioning. $M\in\mathbb{R}^{N_{p}\times N_{V}}$ is the discretized
divergence operator, $G=-M^{T}\in\mathbb{R}^{N_{V}\times N_{p}}$ the
discretized gradient operator, $\Omega_{V}\in\mathbb{R}^{N_{V}\times N_{V}}$ a
matrix with the ‘velocity’ finite volume sizes on its diagonal, and $C_{V}$
and $D_{V}$ constitute central difference approximations of the convective and
diffusive terms. $A$ is a matrix that averages the temperature from the center
of the ‘temperature’ finite volumes to center of the ‘velocity’ finite
volumes, and the vector $y_{T}$ incorporates the nonzero boundary condition
for the temperature at the lower plate.

The energy-conserving nature of our finite volume method is crucial in
deriving an energy-consistent discretization of viscous dissipation. The
energy-conserving property means that, in absence of boundary contributions,
the discretized convective and pressure gradient operators do not contribute
to the kinetic energy balance: $V_{h}^{T}C_{V}(V_{h})=0$ and
$V_{h}^{T}Gp_{h}=0$, just like in the continuous case. This is achieved by
using a skew-symmetric convection operator and the compatibility between $M$
and $G$ via $G=-M^{T}$. The discrete kinetic energy balance then reads:

$\frac{\mathrm{d}E_{k,h}}{\mathrm{d}t}=-\alpha_{1}\epsilon_{U,h}+\alpha_{2}V_{h}^{T}(AT_{h}+y_{T}),$
(73)

where $E_{k,h}=\frac{1}{2}V_{h}^{T}\Omega_{V}V_{h}$. The global viscous
dissipation (i.e. summed over the entire domain) is given by
$\epsilon_{U,h}=\|QV_{h}\|_{2}^{2}>0$, where $Q$ stems from decomposing the
symmetric negative-definite diffusive operator as $D_{V}=-Q^{T}Q$. Equation
(73) is the semi-discrete counterpart of equation (52).

$v_{i,j+1/2}$
---
$p_{i,j}$
---
$u_{i+1/2,j}$
---
$v_{i,j-1/2}$
---
$\Delta x$
---
$\Delta y$
---
$u_{i-1/2,j}$
---
$T_{i,j}$
---
Figure 2: Staggered grid with positioning of unknowns around a pressure
volume.

###  3.2 Proposed viscous dissipation function

Given a discretization that satisfies a discrete kinetic energy balance, the
key step is to design a discretization scheme of the internal energy equation
(40) which is such that discrete versions of the global balances (30) and (31)
are obtained. In particular, the viscous dissipation in the internal energy
equation should cancel the viscous dissipation term in the kinetic energy
equation, where the latter is fully determined by the choice of the diffusion
operator and the expression for the local kinetic energy. The choice for the
diffusion operator (second-order central differencing) is straightforward. The
choice for the expression of the local kinetic energy on a staggered grid is
however not obvious. We propose the following definition:

$k_{i,j}:=\frac{1}{4}u_{i+1/2,j}^{2}+\frac{1}{4}u_{i-1/2,j}^{2}+\frac{1}{4}v_{i,j+1/2}^{2}+\frac{1}{4}v_{i,j-1/2}^{2}.$
(82)

This choice gives a local kinetic energy equation that is consistent with the
continuous equations, as is detailed in B, and consistent with the global
energy definition.

The expression for $\Phi_{h}$ then follows from differentiating the expression
for $k_{ij}$ in time, substituting the momentum equations, and rewriting the
terms involving the diffusive operator (see B). The implied dissipation then
follows by constructing a discrete version of (27). As example, we construct
the discrete version of $u\frac{\partial^{2}u}{\partial
x^{2}}=-\left(\frac{\partial u}{\partial
x}\right)^{2}+\frac{\partial}{\partial x}\left(u\frac{\partial u}{\partial
x}\right)$, being

$\frac{u_{i+1/2,j}}{\Delta x}\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}-\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)=-\frac{1}{2}\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}\right)^{2}-\frac{1}{2}\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)^{2}\\\ +\frac{1}{\Delta
x}\left(\frac{1}{2}(u_{i+3/2,j}+u_{i+1/2,j})\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}-\frac{1}{2}(u_{i+1/2,j}+u_{i-1/2,j})\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right).$ (85)

The first two terms on the right-hand side contribute to the viscous
dissipation function. Repeating this process for the other components
($u\frac{\partial^{2}u}{\partial y^{2}}$, $v\frac{\partial^{2}v}{\partial
x^{2}}$, $v\frac{\partial^{2}v}{\partial y^{2}}$), as outlined in B.2, yields
the following novel expression for the local dissipation function:

$\boxed{\Phi_{i,j}=\frac{1}{2}\Phi^{u}_{i+1/2,j}+\frac{1}{2}\Phi^{u}_{i-1/2,j}+\frac{1}{2}\Phi^{v}_{i,j+1/2}+\frac{1}{2}\Phi^{v}_{i,j-1/2}},$
(86)

where

$\Phi^{u}_{i+1/2,j}=-\frac{1}{2}\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}\right)^{2}-\frac{1}{2}\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)^{2}-\frac{1}{2}\left(\frac{u_{i+1/2,j+1}-u_{i+1/2,j}}{\Delta
y}\right)^{2}-\frac{1}{2}\left(\frac{u_{i+1/2,j}-u_{i+1/2,j-1}}{\Delta
y}\right)^{2},$ (88)
$\Phi^{v}_{i,j+1/2}=-\frac{1}{2}\left(\frac{v_{i+1,j+1/2}-v_{i,j+1/2}}{\Delta
x}\right)^{2}-\frac{1}{2}\left(\frac{v_{i+1,j-1/2}-v_{i,j-1/2}}{\Delta
x}\right)^{2}-\frac{1}{2}\left(\frac{v_{i,j+3/2}-v_{i,j+1/2}}{\Delta
y}\right)^{2}-\frac{1}{2}\left(\frac{v_{i,j+1/2}-v_{i,j-1/2}}{\Delta
y}\right)^{2}.$ (90)

At boundaries, an adaptation of $\Phi_{h}$ is required in order to have a
discrete equivalent of equation (27). This is detailed in equation (215).

Note that $\Phi_{h}$ is derived based on local energy consideration which upon
summation equals the global dissipation, just like equation (58):

$1^{T}\Omega_{p}\Phi_{h}=\epsilon_{U,h}.$ (91)

###  3.3 Internal energy equation

Having proposed a consistent expression for $\Phi_{h}$, the spatial
discretization of the internal energy equation (40) reads:

$\Omega_{p}\frac{\mathrm{d}T_{h}}{\mathrm{d}t}=-C_{T}(V_{h},T_{h})+\alpha_{3}\Omega_{p}\Phi_{h}(V_{h})+\alpha_{4}(D_{T}T_{h}+\hat{y}_{T}),$
(92)

where

$\begin{split}[C_{T}(V_{h},T_{h})]_{i,j}=&\Delta
y\left(u_{i+1/2,j}\frac{1}{2}(T_{i+1,j}+T_{i,j})-u_{i-1/2,j}\frac{1}{2}(T_{i,j}+T_{i-1,j})\right)+\\\
&\Delta
x\left(v_{i,j+1/2}\frac{1}{2}(T_{i,j+1}+T_{i,j})-v_{i,j-1/2}\frac{1}{2}(T_{i,j}+T_{i,j-1})\right)\end{split}$
(93)

is the convection operator. The convection operator has a discrete skew-
symmetry property which will be used in the derivation of the thermal
dissipation balance in the next subsection. $D_{T}$ the standard second-order
difference stencil with boundary conditions encoded in $\hat{y}_{T}$.

The total internal energy is given by $E_{i,h}=1^{T}\Omega_{p}T_{h}$ (simply
summing over all finite volumes). Due to the no-slip boundary conditions on
the velocity field, the convective operator satisfies
$1^{T}C_{T}(V_{h},T_{h})=0$. The summation over the diffusive operator can be
written in terms of the Nusselt numbers (detailed in the next section). The
total internal energy equation thus reads

$\begin{split}\frac{\mathrm{d}E_{i,h}}{\mathrm{d}t}&=\alpha_{3}1^{T}\Omega_{p}\Phi_{h}+\alpha_{4}1^{T}(D_{T}T_{h}+\hat{y}_{T}),\\\
&=\alpha_{3}1^{T}\Omega_{p}\Phi_{h}+\alpha_{4}(\mathrm{Nu}_{H}-\mathrm{Nu}_{C}),\end{split}$
(94)

where in the second line the Nusselt numbers are instantaneous Nusselt
numbers. Upon adding the total kinetic energy equation (73), and using
property (91), the global energy balance results:

$\begin{split}\frac{\mathrm{d}E_{h}}{\mathrm{d}t}=\frac{\mathrm{d}E_{k,h}}{\mathrm{d}t}+\gamma\frac{\mathrm{d}E_{i,h}}{\mathrm{d}t}&=\alpha_{2}V_{h}^{T}(AT_{h}+y_{T})+\gamma\alpha_{4}1^{T}(D_{T}T_{h}+\hat{y}_{T}),\\\
&=\alpha_{2}V_{h}^{T}(AT_{h}+y_{T})+\gamma\alpha_{4}(\mathrm{Nu}_{H}-\mathrm{Nu}_{C}),\end{split}$
(95)

which is the semi-discrete counterpart of equation (54). In other words, we
have proposed a discrete viscous dissipation function that leads to a correct
expression for the total energy equation, namely such that the viscous
dissipation from the kinetic and internal energy equations exactly balances.
Note that in the case of homogeneous Neumann boundary conditions for the
temperature on all boundaries, the last term disappears.

###  3.4 Discrete global balances and Nusselt number relations

We now derive discrete versions of the Nusselt relations that incorporate the
viscous dissipation function, i.e. relations (59) and (63). Our symmetry-
preserving spatial discretization is such that exact discrete relations can be
derived. It is important to realize that the discrete approximation for the
Nusselt number cannot be chosen independently (when the goal is to have exact
discrete global balances) but is implicitly defined once the discretization of
the diffusive operator is chosen. Consider the discretized global internal
energy equation for steady conditions,

$\alpha_{3}1^{T}\Omega_{p}\Phi_{h}(V_{h})+\alpha_{4}1^{T}(D_{T}T_{h}+\hat{y}_{T})=0.$
(96)

The second term can be simplified as

$1^{T}(D_{T}T_{h}+\hat{y}_{T})=-\sum_{i=1}^{N_{x}}\frac{T_{i,1}-T_{H}}{\frac{1}{2}\Delta
y}\Delta x+\sum_{i=1}^{N_{x}}\frac{T_{C}-T_{i,N_{y}}}{\frac{1}{2}\Delta
y}\Delta x=\mathrm{Nu}_{H}-\mathrm{Nu}_{C},$ (97)

where the Nusselt numbers on the lower (hot) and upper (cold) plate are
defined as

$\displaystyle\mathrm{Nu}_{H}$
$\displaystyle:=-\sum_{i=1}^{N_{x}}\frac{T_{i,1}-T_{H}}{\frac{1}{2}\Delta
y}\Delta x,$ (98) $\displaystyle\mathrm{Nu}_{C}$
$\displaystyle:=-\sum_{i=1}^{N_{x}}\frac{T_{C}-T_{i,N_{y}}}{\frac{1}{2}\Delta
y}\Delta x.$ (99)

This leads to the discrete version of (59):

$\alpha_{4}(\mathrm{Nu}_{C}-\mathrm{Nu}_{H})=\alpha_{3}1^{T}\Omega_{p}\Phi_{h}(V_{h}).$
(100)

The discrete version of (63) follows by considering the inner product of
equation (92) with $T_{h}^{T}$ instead of $1^{T}$. An important property of
the convective discretization (93) is that

$T_{h}^{T}C_{T}(V_{h},T_{h})=0,\qquad\forall\,T_{h},\quad\mathrm{if}\quad
MV_{h}=0.$ (101)

This property is most easily derived by recognizing that $C_{h}(V_{h},T_{h})$
can be written in terms of a matrix-vector product
$\tilde{C}_{T}(V_{h})T_{h}$, where $\tilde{C}_{T}(V_{h})$ is skew-symmetric if
$MV_{h}=0$. In addition, the inner product of $T_{h}$ with the diffusive terms
can be written as

$T_{h}^{T}(D_{T}T_{h}+\hat{y}_{T})=\sum_{i=1}^{N_{x}}\left(-T_{H}\frac{T_{i,1}-T_{H}}{\frac{1}{2}\Delta
y}+T_{C}\frac{T_{C}-T_{i,N_{y}}}{\frac{1}{2}\Delta y}\right)\Delta
x-\epsilon_{T,h},$ (102)

where

$\begin{split}\epsilon_{T,h}:=&\sum_{i=1}^{N_{x}}\left(\frac{1}{2}\left(\frac{T_{i,1}-T_{H}}{\frac{1}{2}\Delta
y}\right)^{2}+\sum_{j=2}^{N_{y}}\left(\frac{T_{i,j}-T_{i,j-1}}{\Delta
y}\right)^{2}+\frac{1}{2}\left(\frac{T_{C}-T_{i,N_{y}}}{\frac{1}{2}\Delta
y}\right)^{2}\right)\Delta x\Delta
y+\sum_{j=1}^{N_{y}}\sum_{i=2}^{N_{x}}\left(\frac{T_{i,j}-T_{i-1,j}}{\Delta
x}\right)^{2}\Delta x\Delta y\end{split}$ (103)

is the discrete analogue of (62) and equation (102) is the discrete version of
$\int
T\frac{\mathrm{d}^{2}T}{\mathrm{d}y^{2}}=[T\frac{\mathrm{d}T}{\mathrm{d}y}]-\int(\frac{\mathrm{d}T}{\mathrm{d}y})^{2}$.
With the boundary condition $T_{H}=1$, $T_{C}=0$, we get the balance

$\alpha_{4}\mathrm{Nu}_{H}=\alpha_{4}\epsilon_{T,h}-\alpha_{3}T_{h}^{T}\Omega_{p}\Phi_{h}(V_{h}),$
(104)

which is the discrete version of equation (63).

##  4 Energy-consistent temporal discretization

The system of equations (71), (72) and (92) needs to be integrated in time
with a suitable method in order to preserve a time-discrete version of the
global energy balance (95). A common choice is to use an explicit method (e.g.
Adams-Bashforth) for the nonlinear convective terms and an implicit method
(e.g. Crank-Nicolson) for the (stiff) linear diffusion terms [20, 27, 29], or
an explicit method for both convection and diffusion [30, 31]. In such an
approach, the temperature equation is typically solved first (given velocity
fields at previous time instances), and then the mass and momentum equations
are solved with a pressure-correction approach. However, these methods do not
preserve the global energy balance as they violate the energy-conserving
nature of the nonlinear terms when marching in time [32].

Instead, we show here that the implicit midpoint method can be employed to
achieve energy-consistent time integration. The fully discrete system reads:

$\displaystyle MV_{h}^{n+1/2}$ $\displaystyle=0,$ (105)
$\displaystyle\Omega_{V}\frac{V^{n+1}_{h}-V^{n}_{h}}{\Delta t}$
$\displaystyle=-C_{V}(V_{h}^{n+1/2})-Gp_{h}^{n+1/2}+\alpha_{1}D_{V}V_{h}^{n+1/2}+\alpha_{2}(AT_{h}^{n+1/2}+y_{T}),$
(106) $\displaystyle\Omega_{p}\frac{T^{n+1}_{h}-T^{n}_{h}}{\Delta t}$
$\displaystyle=-C_{T}(V_{h}^{n+1/2},T_{h}^{n+1/2})+\alpha_{3}\Omega_{p}\Phi(V_{h}^{n+1/2})+\alpha_{4}(D_{T}T_{h}^{n+1/2}+\hat{y}_{T}).$
(107)

Here $V_{h}^{n+1/2}=\frac{1}{2}(V_{h}^{n}+V_{h}^{n+1})$ and
$T_{h}^{n+1/2}=\frac{1}{2}(T_{h}^{n}+T_{h}^{n+1})$. Upon multiplying (106) by
$(V_{h}^{n+1/2})^{T}$ and (107) by $1^{T}$, and adding the two resulting
equations, we get the discrete energy balance,

$\frac{E_{h}^{n+1}-E_{h}^{n}}{\Delta
t}=\frac{E_{k,h}^{n+1}-E_{k,h}^{n}}{\Delta
t}+\gamma\frac{E_{i,h}^{n+1}-E_{i,h}^{n}}{\Delta
t}=\alpha_{2}(V_{h}^{n+1/2})^{T}(AT_{h}^{n+1/2}+y_{T})+\gamma\alpha_{4}1^{T}(D_{T}T_{h}^{n+1/2}+\hat{y}_{T}),$
(108)

which is the fully-discrete counterpart of equations (54) and (95). The
derivations hinges again on skew-symmetry of the convection operator $C_{V}$,
the compatibility between $M$ and $G$ ($G=-M^{T}$), and the consistency
requirement on the viscous dissipation function, equation (91).

The system of equations (105) - (107) leads to a large system of nonlinear
equations which has a saddle point structure due to the divergence-free
constraint. We solve the system in a segregated fashion and iterate at each
time step with a standard pressure-correction method until the residual of the
entire system is below a prescribed tolerance. We will compare this energy-
conserving time integration approach to an explicit one-leg method [30, 31] in
section 6.

##  5 Steady state results (Rayleigh-Bénard)

The concept of energy consistency is best demonstrated through time-dependent
simulations. However, we start with steady-state results in order to validate
the spatial discretization method and to get intuition for the effect of the
Gebhart number on the Nusselt number. For the results reported here we employ
a direct solver that solves the entire coupled non-linear system of equations
that arises from spatial discretization. As initial guess we take the
following divergence-free velocity field:

$\displaystyle u(x,y)$ $\displaystyle=-64x^{2}(x-1)^{2}y(y-1)(2y-1),$ (109)
$\displaystyle v(x,y)$ $\displaystyle=64x(x-1)(2x-1)y^{2}(y-1)^{2},$ (110)

which is inspired by the regularized driven cavity problem [33]. For the
temperature we take a random field (between 0 and 1). The idea behind this
choice of initial condition is to avoid the non-linear solver to be stuck in
the trivial solution ($\bm{u}=0$). Note that in all simulations in this
article, we will set $\Pr=0.71$ (air), and use non-dimensionalization choice
I. Choices II and III give equivalent results apart from scaling factors.

###  5.1 Grid convergence study for no-dissipation case (Ge = 0)

Figure 3(a) shows the temperature field when viscous dissipation is not
included ($\mathrm{Ge}=0$). The resulting Nusselt numbers as a function of
grid refinement are displayed in Table 3 and indicate excellent agreement with
literature [34]. We note that the Nusselt numbers as defined by (98) and (99)
are first-order approximations. More accurate approximations can be
constructed by including more interior points. We are not using such high-
order accurate approximations as they would not satisfy the discrete global
balance (100). Note also that we only report $\mathrm{Nu}_{H}$ since
$\mathrm{Nu}_{C}=\mathrm{Nu}_{H}$ up to machine precision.

(a) $\mathrm{Ge}=0$.

(b) $\mathrm{Ge}=0.1$.

(c) $\mathrm{Ge}=1$.

Figure 3: Steady-state temperature field for $\mathrm{Ra}=10^{5}$ on a 128 $\times$ 128 grid, for different $\mathrm{Ge}$.  grid | $\mathrm{Ra}=10^{3}$ | $\mathrm{Ra}=10^{4}$ | $\mathrm{Ra}=10^{5}$
---|---|---|---
$32^{2}$ | 1.000 | 2.170 | 3.933
$64^{2}$ | 1.000 | 2.161 | 3.916
$128^{2}$ | 1.000 | 2.159 | 3.912
$256^{2}$ | 1.000 | 2.158 | 3.911
Cai et al. [34] ($256^{2}$) | 1.000 | 2.158 | 3.911
Table 3: Convergence of Nusselt number (98) with grid refinement for different
Rayleigh numbers and $\mathrm{Ge}=0$.

###  5.2 Grid convergence study for viscous dissipation case (Ge¿ 0)

When including viscous dissipation ($\mathrm{Ge}>0$) in the internal energy
equation, the flow field changes qualitatively and loses its symmetric nature,
as can be observed in figures 3(b)-3(c). The Nusselt numbers at the hot and
cold plate start to deviate from each other, their difference being equal to
the dissipation function, according to equation (100) (or (59)). This is
reported in table 4 and figure 4(a). The critical Rayleigh number that we find
from the bifurcation diagram is $\mathrm{Ra}_{c}\approx 2585$, which is in
excellent agreement with the value of 2585.02 reported in literature [35, 36].
It is independent of the value of the Prandtl number, as shown in [35], and
also independent of the value of the Gebhart number. This latter fact follows
by extending the linear stability analysis in [35] and realizing that the term
$\nabla\bm{u}:\nabla\bm{u}$ with
$\bm{u}=\bm{u}_{0}+\varepsilon\bm{u}^{\prime}$ and background state
$\bm{u}_{0}=0$ leads to the term
$\varepsilon^{2}\nabla\bm{u}^{\prime}:\nabla\bm{u}^{\prime}$, which disappears
when gathering terms of $\mathcal{O}(\varepsilon)$. The results in figure 4(a)
show indeed that the bifurcation point is the same for different values of
$\mathrm{Ge}$.

Figure 4(b) shows a different interpretation of the Nusselt number, indicating
the relation with the thermal dissipation and viscous dissipation according to
equation (104) (or (63)). The results confirm that the thermal dissipation
lies in between the Nusselt number of the hot and cold plate.

grid | $\mathrm{Ge}=0.1$ | $\mathrm{Ge}=1$
---|---|---
| $\mathrm{Nu}_{H}$ | $\mathrm{Nu}_{C}$ | $\mathrm{Nu}_{H}$ | $\mathrm{Nu}_{C}$
$32^{2}$ | 2.111 | 2.228 | 1.582 | 2.729
$64^{2}$ | 2.103 | 2.219 | 1.578 | 2.716
$128^{2}$ | 2.101 | 2.217 | 1.576 | 2.713
$256^{2}$ | 2.100 | 2.216 | 1.576 | 2.712
(a) $\mathrm{Ra}=10^{4}$.

grid | $\mathrm{Ge}=0.1$ | $\mathrm{Ge}=1$
---|---|---
| $\mathrm{Nu}_{H}$ | $\mathrm{Nu}_{C}$ | $\mathrm{Nu}_{H}$ | $\mathrm{Nu}_{C}$
$32^{2}$ | 3.786 | 4.080 | 2.448 | 5.319
$64^{2}$ | 3.770 | 4.062 | 2.441 | 5.299
$128^{2}$ | 3.766 | 4.057 | 2.439 | 5.293
$256^{2}$ | 3.765 | 4.056 | 2.439 | 5.292
(b) $\mathrm{Ra}=10^{5}$.

Table 4: Convergence of Nusselt numbers (98) and (99) with grid refinement for
different Rayleigh and different Gebhart numbers.

(a) Viscous dissipation does not change the critical Rayleigh number, but
leads to a difference between the Nusselt number of the upper and lower plate.

(b) Relation between thermal dissipation and viscous dissipation for
$\mathrm{Ge}=1$.

Figure 4: Bifurcation diagram for Rayleigh-Bénard problem including viscous
dissipation.

##  6 Time-dependent, energy-conserving simulation (Rayleigh-Taylor)

The previous section confirmed the (discrete) steady-state Nusselt number
balances. In this section we consider the core idea of this article: achieving
exact energy conservation in a time-dependent simulation. Exact energy
conservation requires that all contributions from boundary terms disappear,
which we achieve by prescribing no-slip conditions $\bm{u}=0$ and adiabatic
conditions $\frac{\partial T}{\partial n}=0$ on all boundaries (the pressure
does not require boundary conditions). The energy balance then represents a
pure exchange of kinetic, internal and potential energy according to

$\frac{E_{h}^{n+1}-E_{h}^{n}}{\Delta
t}=\frac{E_{k,h}^{n+1}-E_{k,h}^{n}}{\Delta
t}+\gamma\frac{E_{i,h}^{n+1}-E_{i,h}^{n}}{\Delta
t}=\alpha_{2}(V_{h}^{n+1/2})^{T}(AT_{h}^{n+1/2}+y_{T}).$ (129)

However, with adiabatic boundary conditions we cannot simulate the classic
Rayleigh-Bénard problem. Instead, we turn to the well-known Rayleigh-Taylor
problem, featuring a cold (heavy) fluid on top of a warm (light) fluid. A
sketch of the set-up is shown in figure 5. The energy-conserving implicit
midpoint (‘IM’) method detailed in section 4 will be compared to the explicit
one-leg (‘OL’) method commonly used in DNS studies [30, 31] (where we take
$\kappa=\frac{1}{2}$ and a fixed time step).

$\frac{\partial T}{\partial x}=0$
---
$\bm{u}=0$
---
$\frac{\partial T}{\partial x}=0$
---
$\bm{u}=0,\frac{\partial T}{\partial y}=0$
---
$\bm{u}=0,\frac{\partial T}{\partial y}=0$
---
$y$
---
$x$
---
$\bm{u}=0$
---
$T=1$
---
$T=0$
---
Figure 5: Problem set-up with initial condition for Rayleigh Taylor problem.

The domain size is $1\times 2$, the grid is $64\times 128$, the time step
$\Delta t=5\cdot 10^{-3}$ and the end time $T=100$. We consider the case
$\mathrm{Ra}=10^{6}$ and $\mathrm{Ge}=\\{0.1,1\\}$ . The instability naturally
arises due to growth of round-off errors (no perturbation is added in the
initial condition). After the initial instability has developed, an asymmetry
in the solution appears, triggering a sequence of well-known ‘mushroom’ type
plumes: hot plumes rising upward and cold plumes sinking downward (figure 6).
Noteworthy are the differences in the development of the instability due to
different time integration methods: IM predicts an earlier onset (around
$t=23$) than OL (around $t=33$), and the evolution stays symmetric for a much
longer period of time in case of IM. For both methods, the time of onset of
instability is insensitive to the value of $\mathrm{Ge}$, just like the
bifurcation point in the steady state Rayleigh-Bénard simulation was
insensitive to the value of $\mathrm{Ge}$. The differences between the methods
might be attributed to the absence of artificial dissipation in the IM scheme
compared to OL, as well as its more symmetric nature. However, one should note
that the problem is chaotic, and similar differences in the solution can also
be obtained by adding minute perturbations to the initial condition.

(a) $t=30$, $\mathrm{Ge}=0.1$, IM

(b) $t=40$, $\mathrm{Ge}=0.1$, OL

(c) $t=30$, $\mathrm{Ge}=1$, IM

(d) $t=40$, $\mathrm{Ge}=1$, OL

Figure 6: Rayleigh-Taylor temperature fields.

Since there is no driving force and all boundary conditions are homogeneous,
viscosity damps the velocity field back to a homogeneous steady state, while
at the same time increasing the temperature through dissipation. This increase
in temperature is clear from figure 7(a), where the average temperature is
displayed. Compared to the initial temperature difference $\Delta T=1$, the
relative temperature increase is about $2\%$ for $\mathrm{Ge}=0.1$ and more
than $20\%$ for $\mathrm{Ge}=1$. Note that many existing natural convection
models, which ignore the viscous dissipation term, would not predict any
temperature increase. With our proposed energy-consistent viscous dissipation
function, the temperature increase exactly matches the kinetic energy loss
through viscous dissipation. This is confirmed in figure 7(b), which shows the
energy error

$\varepsilon_{E}:=\left|\frac{E_{k,h}^{n+1}-E_{k,h}^{n}}{\Delta
t}+\gamma\frac{E_{i,h}^{n+1}-E_{i,h}^{n}}{\Delta
t}-\alpha_{2}(V_{h}^{n+1/2})^{T}(AT_{h}^{n+1/2}+y_{T})\right|.$ (140)

For IM the error remains at the tolerance with which we solve the system of
nonlinear equations ($10^{-12}$). For OL, the error is initially at a similar
level but increases to $\mathcal{O}(10^{-6})$ when the instability develops
(for $t>30$). However, one could argue that this advantage is offset by the
fact that IM is roughly 4-5$\times$ as expensive because it requires roughly
4-5 iterations (Poisson solves) per time step, making OL much faster to run.
Consequently, OL will be employed for the 3D simulations in the next section.
Note that this balance of accuracy versus computational costs depends on the
details of the flow problem and might differ in other test cases.

(a) Average temperature and velocity.

(b) Energy error.

Figure 7: Rayleigh-Taylor results, IM = Implicit Midpoint, OL = One-Leg
scheme.

(a) $\mathrm{Ra}=10^{8}$ and $\mathrm{Ge}=0$.

(b) $\mathrm{Ra}=10^{8}$ and $\mathrm{Ge}=1$.

(c) $\mathrm{Ra}=10^{10}$ and $\mathrm{Ge}=0$.

(d) $\mathrm{Ra}=10^{10}$ and $\mathrm{Ge}=1$.

Figure 8: Instantaneous temperature fields for 3D RBC at different Rayleigh
and Gebhart numbers. For a visualization of the 3D time-dependent simulation
results, we refer to the supplementary material.

##  7 Energy-conserving simulation of a turbulent flow

As a final test-case, we consider the numerical simulation of an air-filled
($\Pr=0.71$) Rayleigh–Bénard flow at two different Rayleigh numbers,
$\mathrm{Ra}=10^{8}$ and $10^{10}$. Direct numerical simulations (DNS) were
carried out and analyzed in previous studies [37, 38] without taking into
account the viscous dissipation effects ($\mathrm{Ge}=0$). Here, the results
are extended to $\mathrm{Ge}=0.1$ and $\mathrm{Ge}=1$ keeping the same domain
size ($\pi\times 1\times 1$) and mesh resolution ($400\times 208\times 208$
for $\mathrm{Ra}=10^{8}$, and $1024\times 768\times 768$ for
$\mathrm{Ra}=10^{10}$). Grids are constructed with a uniform grid spacing in
the periodic $x$-direction whereas wall-normal points ($y$ and $z$ directions)
are distributed following a hyperbolic-tangent function as follows (identical
for the $z$-direction)

$y_{i}=\frac{1}{2}\left(1+\frac{\tanh\left(\gamma_{y}(2(i-1)/N_{y}-1)\right)}{\tanh\gamma_{y}}\right),\hskip
17.07164pti=1,\dots,N_{y}+1,$ (141)

where $N_{y}$ and $\gamma_{y}$ are the number of control volumes and the
concentration factor in the $y$-direction, respectively. In our case,
$\gamma_{y}=\gamma_{z}=1.4$ for $\mathrm{Ra}=10^{8}$ and
$\gamma_{y}=\gamma_{z}=1.6$ for $\mathrm{Ra}=10^{10}$. For further details,
the reader is referred to our previous works [37, 38].

Instantaneous temperature fields corresponding to the statistically steady
state are displayed in Figure 8. As expected, thermal dissipation effects at
$\mathrm{Ge}=1$ lead to a significant increase in the average cavity
temperature which is clearly visible for both Rayleigh numbers. As in 2D, the
flow symmetry (in average sense) with respect to the mid-height plane is lost
for $\mathrm{Ge}>0$ leading to higher (lower) Nusselt number for the top
(bottom) wall. Subsequently, the top (bottom) thermal boundary layer becomes
thinner (thicker) with respect to the case at $\mathrm{Ge}=0$. This implies
that mesh resolution requirements in the near-wall region are also
asymmetrical; however, in this work, for the sake of simplicity, the grid
spacing at the two walls is the same regardless of the Gebhart number.

Figure 9: Time-evolution of the most relevant energy contributions for
$\mathrm{Ge}=1$; (left) Finest grid: $400\times 208\times 208\approx
17.3\text{M}$; (right) Coarsest grid: $50\times 26\times 26\approx
0.034\text{M}$.

All simulations have been carried out for $500$ time-units starting from a
zero velocity field and uniformly distributed random temperatures between
$T_{C}$ and $T_{H}$. As the fluid sets in motion, initially the discrete
kinetic energy of the system increases. Then, after a sufficiently long period
of time (around $50$ time-units) a statistically steady state is reached. This
is clearly observed in Figure 9 where the time-evolution of various rate-of-
changes of energy are shown. Results correspond to $\mathrm{Ra}=10^{8}$ and
$\mathrm{Ge}=1$ using a very fine ($400\times 208\times 208\approx
17.3\text{M}$) and a very coarse mesh. Similar results are obtained for the
other tested configurations. As expected, once a statistically steady state is
reached, the kinetic energy fluctuates around its mean value and therefore its
rate-of-change $\mathrm{d}E_{k,h}/\mathrm{d}t$ (in red) fluctuates around
zero. Only two terms contribute to the global kinetic energy of the system
(see equation (73)): the global viscous dissipation, $\epsilon_{u,h}$ (in
yellow), and the contribution of the buoyancy forces given by
$\alpha_{2}V_{h}^{T}(AT_{h}(t)+y_{T})$ (in blue). These two contributions
cancel each other on average when a statistically steady state is reached. The
former is transferred into internal energy, $E_{i,h}$, whereas the latter can
be viewed as a transfer from potential to kinetic energy. In addition, the
total energy of the system is exactly in balance with the buoyancy term and
the heat conduction through the top and bottom boundaries (green line), as
given by (95), repeated here for convenience:

$\frac{\mathrm{d}E_{k,h}}{\mathrm{d}t}+\gamma\frac{\mathrm{d}E_{i,h}}{\mathrm{d}t}-\alpha_{2}V_{h}^{T}(AT_{h}+y_{T})-\gamma\alpha_{4}(\mathrm{Nu}_{H}-\mathrm{Nu}_{C})=0.$
(142)

This proofs that the viscous dissipation function has indeed been discretized
correctly, since an imbalance between the viscous dissipation implied by the
kinetic energy equation and the explicitly added viscous dissipation in the
internal energy equation would otherwise show up. These energy balances are
exactly satisfied for any grid, so even for very coarse grids (see Figure 9,
right). Notice again that it is important that the Nusselt numbers are
evaluated consistently with the discretization of the diffusive terms in the
internal energy equation, as explained in Section 3.4.

In addition to these instantaneous balances, we show in Figure 10 that the
time averages of the exact relations given in equations (59), (61) and (63)
are preserved at the discrete level, similar to what was shown in steady-state
in 2D (see Figure 4(a)). However, here we display time-averaged Nusselt
numbers and consider a wide range of meshes. The finest meshes correspond to
the DNS simulations shown in Figure 8 whereas coarser and coarser meshes have
been generated by reducing the number of grid points in each spatial direction
by factors of approximately $\sqrt{2}$. Hence, after six successive mesh
coarsenings, the total number of grid points is reduced by approximately
$((\sqrt{2})^{6})^{3}=2^{9}=512$. This under-resolution causes a pile-up of
(kinetic) energy close to the smallest resolved scales, that leads to higher
values of $\epsilon_{U}$ and, therefore, an increase of both $\mathrm{Nu}_{H}$
(see equation (61)) and $\mathrm{Nu}_{C}-\mathrm{Nu}_{H}$ (see equation (59)).
Although the solution is surely less accurate at coarse grids, the fact that
an energy balance is still satisfied, makes our approach an excellent starting
point for developing or testing sub-grid scale models, as the additional
dissipation that is introduced by the sub-grid scale model can be exactly
quantified.

(a) $\mathrm{Ra}=10^{8}$.

(b) $\mathrm{Ra}=10^{10}$.

Figure 10: Time-averaged Nusselt numbers at lower and upper plate for a set of
meshes at $\mathrm{Ge}=0$, $\mathrm{Ge}=0.1$ and $\mathrm{Ge}=1$.

##  8 Conclusions

In this paper we have proposed a new energy-consistent discretization of the
viscous dissipation function. The viscous dissipation function is an important
quantity, for example in turbulent flow computations, where it is critical to
assess the global energy balances, or in natural convection flows, where it
leads to internal heating. This latter case has been the focus of this
article, and we have shown that including the viscous dissipation function in
the internal energy equation leads to a consistent total energy balance:
viscous dissipation acts as a sink in the kinetic energy equation and as a
source in the internal energy equation, such that the sum of internal and
kinetic energy only changes due to buoyancy and thermal diffusion.

Our key result is a new discretization of the local viscous dissipation
function that abides by the total energy balance. We have shown that it is
determined by two choices, namely the discretization of the diffusive terms in
the momentum equations and the expression for the local kinetic energy. The
discretization of the diffusive terms is detailed for both general (non-
constant viscosity) and simplified (constant viscosity) stress tensor
expressions. The proposed expression for the local kinetic energy is such that
a discrete local kinetic energy equation is satisfied, and leads to a
quadratic, strictly dissipative form of the viscous dissipation function, also
for general stress tensors. Near boundaries we have proposed corrections to
the viscous dissipation function to keep the dissipative property.

The numerical experiments in 2D and 3D show that viscous dissipation does not
affect the critical Rayleigh number at which instabilities form, but it does
significantly impact the development of instabilities once they occur, leading
to a significant difference between the Nusselt numbers on the cold and hot
plates. Moreover, simulations of turbulent Rayleigh-Bénard convection have
clearly shown that the proposed discretization is stable even for very coarse
grids. Namely, the numerical discretization does not interfere with the energy
balances and, therefore, we consider that the proposed method is an excellent
starting point for testing sub-grid scale models.

The analysis in this paper has been performed for the classic finite-volume
staggered grid method. Extensions to other discretization methods, such as
finite differences or finite elements, are in principle possible provided that
a discrete local kinetic energy balance mimicking the continuous balance can
be identified. Another limitation of this work is the assumption of
incompressible flow, which might seem restrictive given the fact that viscous
dissipation typically becomes important for compressible flows. However, the
idea of discretizing the viscous dissipation term in an energy-consistent
manner is also applicable to compressible flows, see e.g. [39, 40], and we
expect our work can therefore be extended in that direction.

As mentioned, an important avenue for future work lies in the assessment of
subgrid-scale models for turbulent flows, including those driven by buoyancy.
For example, in large-eddy simulation, the kinetic energy equation of the
resolved scales and of the subgrid-scales features viscous dissipation terms,
and the current work provides a basis for proper discrete representations of
these terms.

## Data availability statement

The incompressible Navier-Stokes code is available at
https://github.com/bsanderse/INS2D (Matlab version). A Julia version is
available from https://github.com/agdestein/IncompressibleNavierStokes.jl. The
data generated in this work is available upon request.

## CRediT

Benjamin Sanderse: Conceptualization, Methodology, Writing - Original Draft,
Software; Xavier Trias: Writing - Original Draft, Writing - Review & Editing,
Software

## Acknowledgements

This publication is part of the project ”Discretize first, reduce next” (with
project number VI.Vidi.193.105) of the research programme NWO Talent Programme
Vidi which is (partly) financed by the Dutch Research Council (NWO). F.X.T. is
supported by the Ministerio de Economía y Competitividad, Spain, RETOtwin
project (PDC2021-120970-I00). Turbulent Rayleigh-Bénard simulations were
carried out on MareNostrum 4 supercomputer at BSC. The authors thankfully
acknowledge these institutions.

##  Appendix A Forms of the dissipation function

In this appendix we expalin why the dissipation function changes depending on
which form of the stress tensor is used. The stress tensor for an
incompressible fluid with non-constant viscosity is given by

$\hat{\bm{\tau}}=\mu(\nabla\bm{u}+(\nabla\bm{u})^{T}).$ (143)

In the case of constant viscosity $\mu$, the divergence of the stress tensor
can be simplified:

$\nabla\cdot\hat{\bm{\tau}}=\mu\nabla\cdot(\nabla\bm{u}+(\nabla\bm{u})^{T})=\mu\nabla\cdot\nabla\bm{u}+\mu\nabla\cdot(\nabla\bm{u})^{T}=\mu\nabla^{2}\bm{u}+\mu\nabla(\nabla\cdot\bm{u})=\mu\nabla^{2}\bm{u}=:\nabla\cdot\bm{\tau},$
(144)

where

$\bm{\tau}=\mu\nabla\bm{u}=\hat{\bm{\tau}}-\mu(\nabla\bm{u})^{T}.$ (145)

Note that $\bm{\tau}$ is not a proper stress tensor, since it is not
symmetric. We stress that $\nabla\cdot\bm{\tau}=\nabla\cdot\hat{\bm{\tau}}$,
even though $\bm{\tau}\neq\hat{\bm{\tau}}$.

In the kinetic energy equation the divergence of the stress tensor is
multiplied by $\bm{u}$: $\bm{u}\cdot(\nabla\cdot\bm{\tau})$. Since
$\nabla\cdot\bm{\tau}=\nabla\cdot\hat{\bm{\tau}}$, we also have

$\bm{u}\cdot(\nabla\cdot\bm{\tau})=\bm{u}\cdot(\nabla\cdot\hat{\bm{\tau}}).$
(146)

Expanding both the left-hand and right-hand side with a vector identity (note:
also valid for non-symmetric $\bm{\tau}$) gives:

$\nabla\cdot(\bm{\tau}\cdot\bm{u})-\bm{\tau}:\nabla\bm{u}=\nabla\cdot(\hat{\bm{\tau}}\cdot\bm{u})-\hat{\bm{\tau}}:\nabla\bm{u}.$
(147)

or

$\nabla\cdot(\bm{\tau}\cdot\bm{u})-\Phi=\nabla\cdot(\hat{\bm{\tau}}\cdot\bm{u})-\hat{\Phi}.$
(148)

The crucial point is that, even though equation (148) holds, the individual
terms are not equal, i.e. $\Phi\neq\hat{\Phi}$ and
$\nabla\cdot(\bm{\tau}\cdot\bm{u})\neq\nabla\cdot(\hat{\bm{\tau}}\cdot\bm{u})$.

This insight can be further clarified by evaluating these expressions in 2D
Cartesian coordinates:

$\displaystyle\nabla\cdot(\hat{\bm{\tau}}\cdot\bm{u})$
$\displaystyle=\mu\left[\frac{\partial}{\partial x}\left(2u\frac{\partial
u}{\partial x}\right)+\frac{\partial}{\partial y}\left(u\left(\frac{\partial
u}{\partial y}+\frac{\partial v}{\partial
x}\right)\right)+\frac{\partial}{\partial x}\left(v\left(\frac{\partial
u}{\partial y}+\frac{\partial v}{\partial
x}\right)\right)+\frac{\partial}{\partial y}\left(2v\frac{\partial v}{\partial
y}\right)\right],$ (149) $\displaystyle\hat{\Phi}$
$\displaystyle=\mu\left[2\left(\frac{\partial u}{\partial
x}\right)^{2}+\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial
x}\right)^{2}+2\left(\frac{\partial v}{\partial y}\right)^{2}\right],$ (150)
$\displaystyle\nabla\cdot(\bm{\tau}\cdot\bm{u})$
$\displaystyle=\mu\left[\frac{\partial}{\partial x}\left(u\frac{\partial
u}{\partial x}\right)+\frac{\partial}{\partial y}\left(u\frac{\partial
u}{\partial y}\right)+\frac{\partial}{\partial x}\left(v\frac{\partial
v}{\partial x}\right)+\frac{\partial}{\partial y}\left(v\frac{\partial
v}{\partial y}\right)\right],$ (151) $\displaystyle\Phi$
$\displaystyle=\mu\left[\left(\frac{\partial u}{\partial
x}\right)^{2}+\left(\frac{\partial u}{\partial
y}\right)^{2}+\left(\frac{\partial v}{\partial
x}\right)^{2}+\left(\frac{\partial v}{\partial y}\right)^{2}\right].$ (152)

Note that in a closed domain ($\bm{u}=0$ on all boundaries), we have the
relation

$\int_{\Omega}\Phi\,\mathrm{d}\Omega=\int_{\Omega}\hat{\Phi}\,\mathrm{d}\Omega.$
(153)

##  Appendix B Discrete dissipation operator from local kinetic energy
equation

###  B.1 Momentum equations and choice of local kinetic energy

The energy-conserving discretization presented in equation (72) can be written
component-wise as:

$\displaystyle\frac{\mathrm{d}u_{i+1/2,j}}{\mathrm{d}t}$
$\displaystyle=-\text{conv}^{u}_{i+1/2,j}-\frac{p_{i+1,j}-p_{i,j}}{\Delta
x}+\alpha_{1}\text{diff}^{u}_{i+1/2,j},$ (154)
$\displaystyle\frac{\mathrm{d}v_{i,j+1/2}}{\mathrm{d}t}$
$\displaystyle=-\text{conv}^{v}_{i+1/2,j}-\frac{p_{i,j+1}-p_{i,j}}{\Delta
y}+\alpha_{1}\text{diff}^{v}_{i,j+1/2}+\alpha_{2}\frac{1}{2}(T_{i,j}+T_{i,j+1}).$
(155)

The convective terms are discretized starting from the divergence form, and
due to discrete mass conservation they can be written in skew-symmetric form,
which is energy-conserving. These terms are not the main focus of this work
and we refer to [30, 32] for details.

The aim here is to find a local kinetic energy equation expression and the
exact form of the dissipation terms. The local kinetic energy should be such
that it results in the well-known global kinetic energy balance [30] upon
integration over the entire domain. This global kinetic energy equation is
obtained by taking the inner product of all momentum equations with the full
velocity vector $V_{h}$ (containing $u_{i+1/2,j}$ and $v_{i,j+1/2}$ at all
locations). This resulting global kinetic energy definition
$\frac{1}{2}V_{h}^{T}\Omega_{h}V_{h}$ still leaves room for the definition of
the local kinetic energy.

Our proposal is to choose for the local kinetic energy the definition

$k_{i,j}:=\frac{1}{4}u_{i+1/2,j}^{2}+\frac{1}{4}u_{i-1/2,j}^{2}+\frac{1}{4}v_{i,j+1/2}^{2}+\frac{1}{4}v_{i,j-1/2}^{2}.$
(156)

Upon differentiating,

$\frac{\mathrm{d}k_{i,j}}{\mathrm{d}t}=\frac{1}{2}u_{i+1/2,j}\frac{\mathrm{d}u_{i+1/2,j}}{\mathrm{d}t}+\frac{1}{2}u_{i-1/2,j}\frac{\mathrm{d}u_{i-1/2,j}}{\mathrm{d}t}+\frac{1}{2}v_{i,j+1/2}\frac{\mathrm{d}v_{i+1/2,j}}{\mathrm{d}t}+\frac{1}{2}v_{i,j-1/2}\frac{\mathrm{d}v_{i,j-1/2}}{\mathrm{d}t},$
(157)

and substituting the momentum equations, our proposed definition gives a local
kinetic energy equation that is consistent with the continuous equations. The
stencil of points required to evaluate (157) is shown in figure 11.

$v_{i,j+1/2}$
---
$p_{i,j}$
---
$u_{i+1/2,j}$
---
$v_{i,j-1/2}$
---
$u_{i-1/2,j}$
---
$u_{i+3/2,j}$
---
$u_{i-3/2,j}$
---
$u_{i+1/2,j+1}$
---
$u_{i-1/2,j+1}$
---
$u_{i+1/2,j-1}$
---
$u_{i-1/2,j-1}$
---
$v_{i+1,j+1/2}$
---
$v_{i+1,j-1/2}$
---
$v_{i-1,j+1/2}$
---
$v_{i-1,j-1/2}$
---
$p_{i+1,j}$
---
$p_{i-1,j}$
---
$p_{i,j+1}$
---
$p_{i,j-1}$
---
$v_{i,j+3/2}$
---
$v_{i,j-3/2}$
---
Figure 11: Stencil of velocity and pressure points involved in the local
kinetic energy equation.

The choice (156) is inspired by the fact that it naturally allows a discrete
equivalent of $\bm{u}\cdot\nabla p=\nabla\cdot(p\bm{u})-p\nabla\cdot\bm{u}$:

$\frac{1}{2}u_{i+1/2,j}\frac{p_{i+1,j}-p_{i,j}}{\Delta
x}+\frac{1}{2}u_{i-1/2,j}\frac{p_{i,j}-p_{i-1,j}}{\Delta
x}+\frac{1}{2}v_{i,j+1/2}\frac{p_{i,j+1}-p_{i,j}}{\Delta
y}+\frac{1}{2}v_{i,j-1/2}\frac{p_{i,j}-p_{i,j-1}}{\Delta y}=\\\
\frac{u_{i+1/2,j}\frac{1}{2}(p_{i+1,j}+p_{i,j})-\frac{1}{2}u_{i-1/2,j}(p_{i,j}+p_{i-1,j})}{\Delta
x}+\frac{v_{i,j+1/2}\frac{1}{2}(p_{i,j+1}+p_{i,j})-v_{i,j-1/2}\frac{1}{2}(p_{i,j}+p_{i,j-1})}{\Delta
y}\\\ -p_{i,j}\underbrace{\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}+\frac{v_{i,j+1/2}-v_{i,j-1/2}}{\Delta y}\right)}_{\text{div}(u)_{i,j}}.$
(182)

Furthermore, choice (156) for the local kinetic energy leads to a consistent
quadratic dissipation form in the case of a general stress tensor, as will be
shown in B.4.

###  B.2 Diffusion and dissipation

We continue to investigate the dissipation implied by the diffusive term in
the momentum equation (154) and the kinetic energy choice (156). Restricting
ourselves momentarily to the term $\frac{\partial^{2}u}{\partial x^{2}}$, we
are looking for a discrete equivalent of the relation

$u\frac{\partial^{2}u}{\partial x^{2}}=-\left(\frac{\partial u}{\partial
x}\right)^{2}+\frac{\partial}{\partial x}\left(u\frac{\partial u}{\partial
x}\right).$ (183)

This is given by $u_{i+1/2,j}\cdot\text{diff}^{u}_{i+1/2,j}$:

$\frac{u_{i+1/2,j}}{\Delta x}\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}-\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)=-\frac{1}{2}\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}\right)^{2}-\frac{1}{2}\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)^{2}\\\ +\frac{1}{\Delta
x}\left(\frac{1}{2}(u_{i+3/2,j}+u_{i+1/2,j})\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}-\frac{1}{2}(u_{i+1/2,j}+u_{i-1/2,j})\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right).$ (186)

Equation (186) is important because the discrete local dissipation expression
is explicitly needed in the internal energy equation.

The analysis for the term $\frac{\partial^{2}u}{\partial y^{2}}$ is completely
analogous, and hence we can define the following discrete function that
describes the dissipation implied by the discretized diffusion term of the
momentum equation for $u_{i+1/2,j}$:

$\Phi^{u}_{i+1/2,j}=-\frac{1}{2}\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}\right)^{2}-\frac{1}{2}\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)^{2}\\\ -\frac{1}{2}\left(\frac{u_{i+1/2,j+1}-u_{i+1/2,j}}{\Delta
y}\right)^{2}-\frac{1}{2}\left(\frac{u_{i+1/2,j}-u_{i+1/2,j-1}}{\Delta
y}\right)^{2}.$ (189)

Similarly, the dissipation implied by the discretized diffusion term of the
momentum equation for $v_{i,j+1/2}$ is:

$\Phi^{v}_{i,j+1/2}=-\frac{1}{2}\left(\frac{v_{i+1,j+1/2}-v_{i,j+1/2}}{\Delta
x}\right)^{2}-\frac{1}{2}\left(\frac{v_{i+1,j-1/2}-v_{i,j-1/2}}{\Delta
x}\right)^{2}\\\ -\frac{1}{2}\left(\frac{v_{i,j+3/2}-v_{i,j+1/2}}{\Delta
y}\right)^{2}-\frac{1}{2}\left(\frac{v_{i,j+1/2}-v_{i,j-1/2}}{\Delta
y}\right)^{2}.$ (192)

The entire dissipation term appearing in the kinetic energy equation for
$\frac{\mathrm{d}k_{i,j}}{\mathrm{d}t}$ is then

$\boxed{\Phi_{i,j}=\frac{1}{2}\Phi^{u}_{i+1/2,j}+\frac{1}{2}\Phi^{u}_{i-1/2,j}+\frac{1}{2}\Phi^{v}_{i,j+1/2}+\frac{1}{2}\Phi^{v}_{i,j-1/2}.}$
(193)

###  B.3 Boundary conditions

The analysis in the previous section ignored the effect of boundary
conditions. Upon integrating (183) over the domain, we get

$\int u\frac{\partial^{2}u}{\partial
x^{2}}\mathrm{d}x=-\int\left(\frac{\partial u}{\partial
x}\right)^{2}\mathrm{d}x+\underbrace{\left[u\frac{\partial u}{\partial
x}\right]}_{\text{boundary term}},$ (194)

and the boundary term vanishes in case of homogeneous Dirichlet, homogeneous
Neumann or periodic conditions. The discrete version should mimic this
behavior.

$v_{1,j+1/2}$
---
$p_{1,j}$
---
$u_{3/2,j}$
---
$u_{b,j}$
---
$v_{1,j-1/2}$
---
$\Delta x$
---
$\Delta y$
---
$v_{i,3/2}$
---
$p_{i,1}$
---
$u_{i+1/2,1}$
---
$u_{i+1/2,b}$
---
$v_{i,b}$
---
$v_{i+1,3/2}$
---
$v_{i+1,b}$
---
$u_{i+1/2,2}$
---
Figure 12: Staggered grid near vertical (left) and horizontal (right)
boundary.

Consider the case where the solution on the boundary is given by
$u_{1/2,j}=u_{b,j}$ (figure 12, left). Then the first unknown for which the
momentum equation (154) is solved is $u_{3/2,j}$, and equation (186) becomes

$\frac{u_{3/2,j}}{\Delta x}\left(\frac{u_{5/2,j}-u_{3/2,j}}{\Delta
x}-\frac{u_{3/2,j}-u_{b,j}}{\Delta x}\right)=\frac{1}{\Delta
x}\left(\frac{1}{2}(u_{5/2,j}+u_{3/2,j})\frac{u_{5/2,j}-u_{3/2,j}}{\Delta
x}-\frac{1}{2}(u_{3/2,j}+u_{b,j})\frac{u_{3/2,j}-u_{b,j}}{\Delta x}\right)\\\
-\frac{1}{2}\left(\frac{u_{5/2,j}-u_{3/2,j}}{\Delta
x}\right)^{2}-\frac{1}{2}\left(\frac{u_{3/2,j}-u_{b,j}}{\Delta x}\right)^{2}.$
(212)

In case where $u_{b,j}=0$, we want the boundary term to vanish, like the term
$u\frac{\partial u}{\partial x}$ in the continuous case. However, when setting
$u_{b,j}=0$, the term that corresponds to $u\frac{\partial u}{\partial x}$
reads:

$-\frac{1}{2}(u_{3/2,j}+u_{b,j})\frac{u_{3/2,j}-u_{b,j}}{\Delta
x}=-\frac{1}{2}\frac{u_{3/2,j}^{2}}{\Delta x},$ (213)

and the discrete boundary contribution does not vanish for $u_{b,j}=0$. This
issue is caused by the fact that the finite volumes do not cover the entire
domain, because there is no momentum equation to be solved for $u_{b,j}$ (as
it is given by the boundary data). We resolve this issue by splitting instead
as

$-\frac{u_{3/2,j}}{\Delta x}\left(\frac{u_{3/2,j}-u_{b,j}}{\Delta
x}\right)=-\underbrace{\frac{u_{b,j}}{\Delta
x}\left(\frac{u_{3/2,j}-u_{b,j}}{\Delta x}\right)}_{\textrm{boundary
contribution}}-\underbrace{\frac{u_{3/2,j}-u_{b,j}}{\Delta
x}\left(\frac{u_{3/2,j}-u_{b,j}}{\Delta x}\right)}_{\textrm{dissipation
contribution}},$ (214)

so that the contribution to the dissipation function is

$-\left(\frac{u_{3/2,j}-u_{b},j}{\Delta x}\right)^{2},$ (215)

instead of $-\frac{1}{2}\left(\frac{u_{3/2,j}-u_{b,j}}{\Delta x}\right)^{2}$.

For the discretization of $\frac{\partial^{2}u}{\partial y^{2}}$, we have a
different situation, because the solution points are not aligned with the
boundary. The first unknown is $u_{i+1/2,1}$, which is situated at a distance
$\frac{1}{2}\Delta y$ above the lower boundary. In this case we can write

$\frac{u_{i+1/2,j}}{\Delta y}\left(\frac{u_{i+1/2,j+1}-u_{i+1/2,j}}{\Delta
y}-\frac{u_{i+1/2,j}-u_{i+1/2,j-1}}{\Delta
y}\right)\overset{j=1}{=}\frac{u_{i+1/2,1}}{\Delta
y}\left(\frac{u_{i+1/2,2}-u_{i+1/2,1}}{\Delta
y}-\frac{u_{i+1/2,1}-u_{i+1/2,b}}{\frac{1}{2}\Delta y}\right)\\\
=\frac{1}{\Delta
y}\left(\frac{1}{2}(u_{i+1/2,2}+u_{i+1/2,1})\frac{u_{i+1/2,2}-u_{i+1/2,1}}{\Delta
y}-u_{i+1/2,b}\frac{u_{i+1/2,2}-u_{i+1/2,1}}{\frac{1}{2}\Delta y}\right)\\\
-\frac{1}{2}\left(\frac{u_{i+1/2,2}-u_{i+1/2,1}}{\Delta
y}\right)^{2}-\frac{1}{2}\left(\frac{u_{i+1/2,1}-u_{i+1/2,b}}{\frac{1}{2}\Delta
y}\right)^{2},$ (219)

and we have a correct discrete equivalent of the continuous expression, and no
correction to $\Phi$ is needed.

The analysis for the $v$-component follows in a similar fashion. A correction
is needed in the expression for $\Phi$ associated to
$\frac{\partial^{2}v}{\partial x^{2}}$, but not for
$\frac{\partial^{2}v}{\partial y^{2}}$.

###  B.4 Extension to non-constant viscosity: general stress tensor

For the case of non-constant $\mu$, the discretization of the diffusion terms
in the momentum equation changes to

$\begin{split}\text{diff}^{u}_{i+1/2,j}=&\frac{1}{\Delta
x}\left[\left(2\mu_{i+1,j}\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}\right)-\left(2\mu_{i,j}\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)\right]+\\\ &\frac{1}{\Delta
y}\left[\left(\mu_{i+1/2,j+1/2}\frac{u_{i+1/2,j+1}-u_{i+1/2,j}}{\Delta
y}\right)-\left(\mu_{i+1/2,j-1/2}\frac{u_{i+1/2,j}-u_{i+1/2,j-1}}{\Delta
y}\right)\right]+\\\ &\frac{1}{\Delta
y}\left[\left(\mu_{i+1/2,j+1/2}\frac{v_{i+1,j+1/2}-v_{i,j+1/2}}{\Delta
x}\right)-\left(\mu_{i+1/2,j-1/2}\frac{v_{i+1,j-1/2}-v_{i,j-1/2}}{\Delta
x}\right)\right].\end{split}$ (220)

Importantly, we first show that this form reduces to the expression in
equation (154) for constant $\mu$. In the continuous equations, this happens
because

$\frac{\partial}{\partial x}\left(2\frac{\partial u}{\partial
x}\right)+\frac{\partial}{\partial y}\left(\frac{\partial u}{\partial
y}+\frac{\partial v}{\partial x}\right)=\frac{\partial}{\partial
x}\left(\frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial
y}\left(\frac{\partial u}{\partial y}\right)+\frac{\partial}{\partial
x}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial
y}\right)=\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial
y^{2}}.$ (221)

The derivation hinges on the divergence-freeness of $\bm{u}$ and interchanging
of differentiation in $x-$ and $y$-directions. Discretely, the same derivation
holds, which can be shown by rewriting as follows:

$\begin{split}&\frac{1}{\Delta
x}\left[\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}\right)-\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)\right]+\frac{1}{\Delta
y}\left[\left(\frac{u_{i+1/2,j+1}-u_{i+1/2,j}}{\Delta
y}\right)-\left(\frac{u_{i+1/2,j}-u_{i+1/2,j-1}}{\Delta y}\right)\right]+\\\
&\frac{1}{\Delta x}\left[\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}\right)-\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)+\left(\frac{v_{i+1,j+1/2}-v_{i+1,j-1/2}}{\Delta
y}\right)-\left(\frac{v_{i,j+1/2}-v_{i,j-1/2}}{\Delta
y}\right)\right],\end{split}$ (222)

and the second line evaluates to zero, as it contains the difference of the
divergence associated to volumes $(i+1,j)$ and $(i,j)$.

We continue to derive the dissipation function. As explain in Remark 1 and in
A, the dissipation function changes when the generic stress tensor for non-
constant $\mu$ is considered. Multiplying (220) with $u_{i+1/2,j}$ and
rewriting leads to

$\hat{\Phi}^{u}_{i+1/2,j}=u_{i+1/2,j}\cdot\text{diff}^{u}_{i+1/2,j}=\\\
\frac{1}{\Delta
x}\left(\mu_{i+1,j}(u_{i+3/2,j}+u_{i+1/2,j})\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}-\mu_{i,j}(u_{i+1/2,j}+u_{i-1/2,j})\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)\\\ -\mu_{i+1,j}\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}\right)^{2}-\mu_{i,j}\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)^{2}+\\\ \frac{1}{\Delta
y}\Bigg{(}\mu_{i+1/2,j+1/2}\frac{u_{i+1/2,j+1}+u_{i+1/2,j}}{2}\left[\frac{u_{i+1/2,j+1}-u_{i+1/2,j}}{\Delta
y}+\frac{v_{i+1,j+1/2}-v_{i,j+1/2}}{\Delta x}\right]\\\
-\mu_{i+1/2,j-1/2}\frac{u_{i+1/2,j}+u_{i+1/2,j-1}}{2}\left[\frac{u_{i+1/2,j}-u_{i+1/2,j-1}}{\Delta
y}+\frac{v_{i+1,j-1/2}-v_{i,j-1/2}}{\Delta x}\right]\Bigg{)}\\\
-\frac{\mu_{i+1/2,j+1/2}}{2}\left(\frac{u_{i+1/2,j+1}-u_{i+1/2,j}}{\Delta
y}\right)\left(\frac{u_{i+1/2,j+1}-u_{i+1/2,j}}{\Delta
y}+\frac{v_{i+1,j+1/2}-v_{i,j+1/2}}{\Delta x}\right)\\\
-\frac{\mu_{i+1/2,j-1/2}}{2}\left(\frac{u_{i+1/2,j}-u_{i+1/2,j-1}}{\Delta
y}\right)\left(\frac{u_{i+1/2,j}-u_{i+1/2,j-1}}{\Delta
y}+\frac{v_{i+1,j-1/2}-v_{i,j-1/2}}{\Delta x}\right).$ (230)

The last two terms are not in quadratic form yet. The quadratic form results
upon considering the full kinetic energy expression (156), i.e. adding
$\hat{\Phi}^{u}_{i-1/2,j}=u_{i-1/2,j}\cdot\text{diff}^{u}_{i-1/2,j}$,
$\hat{\Phi}^{v}_{i,j+1/2}=v_{i,j+1/2}\cdot\text{diff}^{v}_{i,j+1/2}$ and
$\hat{\Phi}^{v}_{i,j-1/2}=v_{i,j-1/2}\cdot\text{diff}^{v}_{i,j-1/2}$. The full
dissipation function then reads

$\hat{\Phi}_{i,j}=\frac{1}{2}\hat{\Phi}^{u}_{i+1/2,j}+\frac{1}{2}\hat{\Phi}^{u}_{i-1/2,j}+\frac{1}{2}\hat{\Phi}^{v}_{i,j+1/2}+\frac{1}{2}\hat{\Phi}^{v}_{i,j-1/2}=\\\
-\frac{\mu_{i+1,j}}{2}\left(\frac{u_{i+3/2,j}-u_{i+1/2,j}}{\Delta
x}\right)^{2}-\mu_{i,j}\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta
x}\right)^{2}-\frac{\mu_{i-1,j}}{2}\left(\frac{u_{i-1/2,j}-u_{i-3/2,j}}{\Delta
x}\right)^{2}\\\
-\frac{\mu_{i+1/2,j+1/2}}{4}\left(\frac{u_{i+1/2,j+1}-u_{i+1/2,j}}{\Delta
y}+\frac{v_{i+1,j+1/2}-v_{i,j+1/2}}{\Delta
x}\right)^{2}-\frac{\mu_{i+1/2,j-1/2}}{4}\left(\frac{u_{i+1/2,j}-u_{i+1/2,j-1}}{\Delta
y}+\frac{v_{i+1,j-1/2}-v_{i,j-1/2}}{\Delta x}\right)^{2}\\\
-\frac{\mu_{i-1/2,j+1/2}}{4}\left(\frac{u_{i-1/2,j+1}-u_{i-1/2,j}}{\Delta
y}+\frac{v_{i,j+1/2}-v_{i-1,j+1/2}}{\Delta
x}\right)^{2}-\frac{\mu_{i-1/2,j-1/2}}{4}\left(\frac{u_{i-1/2,j}-u_{i-1/2,j-1}}{\Delta
y}+\frac{v_{i,j-1/2}-v_{i-1,j-1/2}}{\Delta x}\right)^{2}\\\
-\frac{\mu_{i,j+1}}{2}\left(\frac{v_{i,j+3/2}-v_{i,j+1/2}}{\Delta
y}\right)^{2}-\mu_{i,j}\left(\frac{v_{i,j+1/2}-v_{i,j-1/2}}{\Delta
y}\right)^{2}-\frac{\mu_{i,j-1}}{2}\left(\frac{v_{i,j-1/2}-v_{i,j-3/2}}{\Delta
y}\right)^{2}.$ (236)

## References

  * [1] E. D. Siggia, High Rayleigh Number Convection, Annual Review of Fluid Mechanics 26 (1994) 137–168.  doi:10.1146/annurev.fl.26.010194.001033.
  * [2] G. K. Batchelor, Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures, Quart. Appl. Math. 12 (3) (1954) 209–233.  doi:10.1090/qam/64563.
  * [3] S. Grossmann, D. Lohse, Scaling in thermal convection: A unifying theory, J. Fluid Mech. 407 (2000) 27–56.  doi:10.1017/S0022112099007545.
  * [4] G. I. Barenblatt, Scaling, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 2003.
  * [5] D. P. Mckenzie, J. M. Roberts, N. O. Weiss, Convection in the earth’s mantle: Towards a numerical simulation, J. Fluid Mech. 62 (03) (1974) 465.  doi:10.1017/S0022112074000784.
  * [6] G. Schubert, D. L. Turcotte, P. Olson, Mantle Convection in the Earth and Planets, Cambridge University Press, 2001.
  * [7] R. J. Kee, M. E. Coltrin, P. Glarborg, Chemically Reacting Flow: Theory and Practice, Wiley-Interscience, Hoboken, N.J, 2003.
  * [8] J. M. Hewitt, D. P. Mckenzie, N. O. Weiss, Dissipative heating in convective flows, J. Fluid Mech. 68 (4) (1975) 721–738.  doi:10.1017/S002211207500119X.
  * [9] D. L. Turcotte, A. T. Hsui, K. E. Torrance, G. Schubert, Influence of viscous dissipation on Bénard convection, J. Fluid Mech. 64 (2) (1974) 369–374.  doi:10.1017/S0022112074002448.
  * [10] B. Gebhart, Effects of viscous dissipation in natural convection, J. Fluid Mech. 14 (2) (1962) 225–232.  doi:10.1017/S0022112062001196.
  * [11] S. Ostrach, Laminar natural-convection flow and heat transfer of fluids with and without heat sources in channels with constant wall temperatures, Tech. Rep. NACA-TN-2863, Lewis Flight Propulsion Lab., NACA (Dec. 1952).
URL https://ntrs.nasa.gov/api/citations/19930083632/downloads/19930083632.pdf

  * [12] B. Blankenbach, F. Busse, U. Christensen, L. Cserepes, D. Gunkel, U. Hansen, H. Harder, G. Jarvis, M. Koch, G. Marquart, D. Moore, P. Olson, H. Schmeling, T. Schnaubelt, A benchmark comparison for mantle convection codes, Geophysical Journal International 98 (1) (1989) 23–38.  doi:10.1111/j.1365-246X.1989.tb05511.x.
  * [13] S. D. King, C. Lee, P. E. van Keken, W. Leng, S. Zhong, E. Tan, N. Tosi, M. C. Kameyama, A community benchmark for 2-D Cartesian compressible convection in the Earth’s mantle, Geophysical Journal International 180 (1) (2010) 73–87.  doi:10.1111/j.1365-246X.2009.04413.x.
  * [14] A. Barletta, Comments on a paradox of viscous dissipation and its relation to the Oberbeck–Boussinesq approach, International Journal of Heat and Mass Transfer 51 (25-26) (2008) 6312–6316.  doi:10.1016/j.ijheatmasstransfer.2007.10.044.
  * [15] A. Barletta, D. Nield, Effect of pressure work and viscous dissipation in the analysis of the Rayleigh–Bénard problem, International Journal of Heat and Mass Transfer 52 (13-14) (2009) 3279–3289.  doi:10.1016/j.ijheatmasstransfer.2009.02.005.
  * [16] A. Barletta, M. Celli, D. A. Nield, On the onset of dissipation thermal instability for the Poiseuille flow of a highly viscous fluid in a horizontal channel, J. Fluid Mech. 681 (2011) 499–514.  doi:10.1017/jfm.2011.213.
  * [17] A. Barletta, M. Celli, P. V. Brandão, On Mixed Convection in a Horizontal Channel, Viscous Dissipation and Flow Duality, Fluids 7 (5) (2022) 170.  doi:10.3390/fluids7050170.
  * [18] K. B. Winters, P. N. Lombard, J. J. Riley, E. A. D’Asaro, Available potential energy and mixing in density-stratified fluids, J. Fluid Mech. 289 (1995) 115–128.  doi:10.1017/S002211209500125X.
  * [19] G. O. Hughes, B. Gayen, R. W. Griffiths, Available potential energy in Rayleigh–Bénard convection, J. Fluid Mech. 729 (2013) R3.  doi:10.1017/jfm.2013.353.
  * [20] B. Gayen, G. O. Hughes, R. W. Griffiths, Completing the Mechanical Energy Pathways in Turbulent Rayleigh-Bénard Convection, Phys. Rev. Lett. 111 (12) (2013) 124301.  doi:10.1103/PhysRevLett.111.124301.
  * [21] K. Petschel, S. Stellmach, M. Wilczek, J. Lülff, U. Hansen, Kinetic energy transport in Rayleigh–Bénard convection, J. Fluid Mech. 773 (2015) 395–417.  doi:10.1017/jfm.2015.216.
  * [22] A. Barletta, Local energy balance, specific heats and the Oberbeck–Boussinesq approximation, International Journal of Heat and Mass Transfer 52 (21-22) (2009) 5266–5270.  doi:10.1016/j.ijheatmasstransfer.2009.06.006.
  * [23] F. Trias, A. Gorobets, A. Oliva, A simple approach to discretize the viscous term with spatially varying (eddy-)viscosity, Journal of Computational Physics 253 (2013) 405–417.  doi:10.1016/j.jcp.2013.07.021.
  * [24] W. Smith, All Things Flow, Oregon State University, 2019.
  * [25] R. Tailleux, On the energetics of stratified turbulent mixing, irreversible thermodynamics, Boussinesq models and the ocean heat engine controversy, J. Fluid Mech. 638 (2009) 339–382.  doi:10.1017/S002211200999111X.
  * [26] E. P. van der Poel, R. J. A. M. Stevens, D. Lohse, Comparison between two- and three-dimensional Rayleigh–Bénard convection, J. Fluid Mech. 736 (2013) 177–194.  doi:10.1017/jfm.2013.488.
  * [27] B. J. Hepworth, Nonlinear two-dimensional Rayleigh-Bénard convection, Ph.D. thesis, University of Leeds (2014).
URL https://etheses.whiterose.ac.uk/6822/1/thesis_bjh.pdf

  * [28] B. Sanderse, Energy-conserving discretization methods for the incompressible Navier-Stokes equations:application to the simulation of wind-turbine wakes, Ph.D. thesis, Technische Universiteit Eindhoven (2013).
URL https://research.tue.nl/en/publications/energyconserving-discretization-
methods-for-the-incompressible-navierstokes-equations--application-to-the-
simulation-of-windturbine-wakes(d188919d-72bb-4bbd-a15f-20360ec22e95).html

  * [29] K. Sugiyama, E. Calzavarini, S. Grossmann, D. Lohse, Flow organization in two-dimensional non-Oberbeck–Boussinesq Rayleigh–Bénard convection in water, J. Fluid Mech. 637 (2009) 105–135.  doi:10.1017/S0022112009008027.
  * [30] R. Verstappen, A. Veldman, Symmetry-preserving discretization of turbulent flow, Journal of Computational Physics 187 (1) (2003) 343–368.  doi:10.1016/S0021-9991(03)00126-8.
  * [31] F. X. Trias, O. Lehmkuhl, A Self-Adaptive Strategy for the Time Integration of Navier-Stokes Equations, Numerical Heat Transfer, Part B: Fundamentals 60 (2) (2011) 116–134.  doi:10.1080/10407790.2011.594398.
  * [32] B. Sanderse, Energy-conserving Runge–Kutta methods for the incompressible Navier–Stokes equations, Journal of Computational Physics 233 (2013) 100–131.  doi:10.1016/j.jcp.2012.07.039.
  * [33] T. M. Shih, C. H. Tan, B. C. Hwang, Effects of grid staggering on numerical schemes, Int. J. Numer. Meth. Fluids 9 (2) (1989) 193–212.  doi:10.1002/fld.1650090206.
  * [34] W. Cai, H. Ma, Y. Wang, J. Chen, X. Zheng, H. Zhang, Development of POD reduced-order model and its closure scheme for 2D Rayleigh–Bénard convection, Applied Mathematical Modelling 66 (2019) 562–575.  doi:10.1016/j.apm.2018.09.031.
  * [35] A. Y. Gelfgat, Different Modes of Rayleigh–Bénard Instability in Two- and Three-Dimensional Rectangular Enclosures, Journal of Computational Physics 156 (2) (1999) 300–324.  doi:10.1006/jcph.1999.6363.
  * [36] D. Venturi, X. Wan, G. E. Karniadakis, Stochastic bifurcation analysis of Rayleigh–Bénard convection, J. Fluid Mech. 650 (2010) 391–413.  doi:10.1017/S0022112009993685.
  * [37] F. Dabbagh, F. X. Trias, A. Gorobets, A. Oliva, On the evolution of flow topology in turbulent Rayleigh-Bénard convection, Physics of Fluids 28 (2016) 115105.
  * [38] F. Dabbagh, F. X. Trias, A. Gorobets, A. Oliva, Flow topology dynamics in a three-dimensional phase space for turbulent Rayleigh-Bénard convection, Physical Review Fluids 5 (2020) 024603.
  * [39] H. N. Najm, P. S. Wyckoff, O. M. Knio, A Semi-implicit Numerical Scheme for Reacting Flow: I. Stiff Chemistry, Journal of Computational Physics 143 (2) (1998) 381–402.  doi:10.1006/jcph.1997.5856.
  * [40] H. Nemati, Direct numerical simulation of turbulent heat transfer to fluids at supercritical pressures, Ph.D. thesis, Delft University of Technology (2016).
URL
https://repository.tudelft.nl/islandora/object/uuid%3Ab3fcb3a0-28b1-4335-a7a1-3fa1629e03d9