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

Contents lists available at ScienceDirect Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid Energy-consistent discretization of viscous dissipation with application to natural convection flow B. Sanderse a,∗, F.X. Trias b aCentrum Wiskunde & Informatica, Science Park 123, Amsterdam, The Netherlands bHeat and Mass Transfer Technological Center, Technical University of Catalonia, ESEIAAT, c/ Colom 11, 08222 Terrassa (Barcelona), Spain ARTICLE INFO Dataset link: https://github.com/bsanderse/IN S2D,https://github.com/agdestein/Incompress ibleNavierStokes.jl Keywords: Viscous dissipation Energy conservation Staggered grid Natural convection Rayleigh–Bénard Gebhart number 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, the proposed discretization also forms an excellent starting point for testing sub-grid scale models and is a useful tool to assess energy budgets in any turbulence simulation, with or without the presence of natural convection. 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 modeled 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 ∗Corresponding author. E-mail address: [email protected] (B. Sanderse). 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–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 modeling 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 with the Boussinesq approximation. Nevertheless, we anticipate that our idea of discretizing the viscous dissipation term in an energyconsistent manner has a broader scope of applicability since it is also applicable to non-Oberbeck–Boussinesq [14] and compressible flows (see e.g. [15,16]). In the incompressible case, Ostrach [11], Gebhart [10] and Turcotte et al. [9] should be explicitly mentioned, being among the first to https://doi.org/10.1016/j.compfluid.2024.106473 Received 20 June 2024; Received in revised form 24 September 2024; Accepted 1 November 2024 Computers and Fluids 286 (2025) 106473 Available online 13 November 2024 0045-7930/© 2024 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). B. Sanderse and F.X. Trias address the role of viscous dissipation and to introduce next to the wellknown 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 [17–20], 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 ×10) and low Rayleigh numbers (Ra = 104,105) 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. The main quantity of interest is the Nusselt number, which is a measure for the heat transfer at the walls. 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 [21–23], or on the kinetic energy budget [24]. 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 addition, even for the cases where viscous dissipation does not have a strong effect on the temperature, having an accurate and consistent way of evaluating the viscous and thermal dissipation budgets (as present in for example Grossmann–Lohse theory [3]) is another benefit of our proposed discretization scheme. 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. Third, on the fully discrete level, we propose a time integration method that preserves the total energy balance upon time marching. Importantly, the discrete dissipation operator that we propose here is not restricted to the context of natural convection flows. For example, when estimating the dissipation of kinetic energy in DNS or LES simulations of turbulent flows, a consistent expression for the energy dissipation is crucial in evaluating energy budgets and developing subgrid scale models. Another example for which our dissipation operator is important is the case of incompressible Taylor–Couette flow (the flow between two rotating cylinders): the power supplied to the cylinders is converted into heating of the fluid, which can be very significant (1 K/min for the set-up with water reported in [25]). In order to predict the correct temperature increase, the viscous dissipation operator that we propose in this work is needed. In existing simulations of Taylor– Couette flow this effect is ignored (probably because in experimental set-ups active cooling is used to keep the temperature under control). With our dissipation operator, one can perform more realistic studies, in which internal heating through viscous dissipation and cooling through the boundaries can both be included. The paper is structured as follows. Section 2introduces the governing equations, energy balances, and new non-dimensionalization. Sections 3and 4describe the energy-consistent spatial and temporal discretization. Section 5describes steady-state results of Rayleigh– Bénard convection including viscous dissipation, and Section 6describes energy-conserving simulations of Rayleigh–Taylor instabilities including viscous dissipation. Section 7shows the effect of viscous dissipation in 3D DNS of Rayleigh–Bénard convection. Fig. 1. Problem set-up for 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 ∇⋅𝒖= 0,(1) 𝜌0(𝜕𝒖 𝜕 𝑡+ ∇⋅(𝒖⊗𝒖))= −∇𝑝+𝜇∇2𝒖+𝜌𝒈,(2) where 𝒖(𝒙, 𝑡)is the velocity field, 𝑝(𝒙, 𝑡)the pressure, 𝜇the dynamic viscosity, 𝜌(𝒙, 𝑡)the density and 𝜌0a reference density. Without loss of generality, we consider a two-dimensional (instead of three-dimensional) domain 𝛺, with the gravity vector pointing in the negative 𝑦-direction so that 𝒈= −𝑔𝒆𝑦. An example of the domain as used in the Rayleigh– Bénard problem, including the boundary conditions, is given in Fig. 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. The density 𝜌is assumed to vary only with temperature 𝑇(𝒙, 𝑡), according to 𝜌(𝑇) =𝜌0−𝛽 𝜌0(𝑇−𝑇0), where 𝛽is the isobaric coefficient of thermal expansion (𝛽= −1 𝜌(𝜕 𝜌 𝜕 𝑇)𝑝). The NS equations are then written as 𝜌0(𝜕𝒖 𝜕 𝑡+ ∇⋅(𝒖⊗𝒖))= −∇𝑝′+𝜇∇2𝒖−𝛽 𝜌0(𝑇−𝑇0)𝒈,(3) where 𝑝=𝑝′−𝜌0𝑔 𝑦and ∇𝑝= ∇𝑝′−𝜌0𝑔𝒆𝑦. The equation for the internal energy 𝑒𝑖describes the temperature evolution according to 𝜕 𝜕 𝑡(𝜌0𝑐 𝑇 ⏟⏟⏟ 𝑒𝑖 ) + ∇⋅(𝒖(𝜌0𝑐 𝑇)) =𝛷+𝜆∇2𝑇 ,(4) where 𝜆is the thermal conductivity and 𝑐equals 𝑐𝑣in case of an ideal gas (the specific heat at constant volume), and equals 𝑐𝑝−𝑝𝛽 𝜌for a real gas [26]). The contribution of pressure work to the change in internal energy, 𝑝∇⋅𝒖, has been discarded in Eq. (3) because of Eq. (1). The viscous dissipation function  𝛷∶=  𝝉∶ ∇𝒖=𝜇[2(𝜕 𝑢 𝜕 𝑥)2 + 2(𝜕 𝑣 𝜕 𝑦)2 +(𝜕 𝑢 𝜕 𝑦+𝜕 𝑣 𝜕 𝑥)2]≥0.(5) is the key quantity in this work, where the stress tensor is given by  𝝉=𝜇(∇𝒖+ (∇𝒖)𝑇). Expression (5) holds in 2D and is easily generalized Computers and Fluids 286 (2025) 106473 2 B. Sanderse and F.X. Trias to 3D. Since the fluid is incompressible and 𝜇is assumed constant, we have ∇⋅ 𝝉=𝜇∇⋅(∇𝒖+ (∇𝒖)𝑇) =𝜇∇2𝒖+𝜇∇(∇ ⋅𝒖) =𝜇∇2𝒖, which is the form of the diffusive terms used in Eq. (3). The simplified form could be interpreted as ∇⋅𝝉=𝜇∇2𝒖, with 𝝉=𝜇∇𝒖, although 𝝉is not a proper stress tensor (it is not symmetric). Remarkably, the simplified form of the diffusive terms implies adifferent dissipation function, namely 𝛷∶= 𝜇‖∇𝒖‖2≥0.(6) where ‖∇𝒖‖2= ∇𝒖∶ ∇𝒖(the Frobenius inner product). The details regarding the difference between 𝛷and  𝛷are given in Appendix A. In 2D and Cartesian coordinates the viscous dissipation can be written as 𝛷=𝜇[(𝜕 𝑢 𝜕 𝑥)2 +(𝜕 𝑢 𝜕 𝑦)2 +(𝜕 𝑣 𝜕 𝑥)2 +(𝜕 𝑣 𝜕 𝑦)2].(7) In this work we focus mainly on the discretization for expression (7), but we will also explain the discretization of the more general form (5), see Appendix B, Eq. (B.21). 2.2. Total energy conservation Conservation of kinetic energy follows by taking the dot product of Eq. (3) with 𝒖: 𝜕 𝜕 𝑡(1 2𝜌0|𝒖|2 ⏟⏞⏟⏞⏟ 𝑒𝑘 ) + ∇⋅(1 2𝜌0|𝒖|2𝒖) = −𝒖⋅∇𝑝′+𝜇∇⋅(𝒖⋅∇𝒖) −𝜇‖∇𝒖‖2+𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣, (8) where 𝒈⋅𝒖= −𝑔 𝑣and we have used the identity 𝒖⋅∇2𝒖= −‖∇𝒖‖2+ ∇⋅(𝒖⋅∇𝒖).(9) Upon adding the kinetic and internal energy Eqs. (4) and (8), the viscous dissipation term cancels and we arrive at the equation for the total energy 𝑒=𝑒𝑘+𝑒𝑖: 𝜕 𝜕 𝑡(𝑒𝑘+𝑒𝑖) + ∇⋅((𝑒𝑘+𝑒𝑖)𝒖) = −∇⋅(𝑝′𝒖) +𝜇∇⋅(𝒖⋅∇𝒖) +𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣+𝜆∇2𝑇 . (10) All terms are in conservative (divergence) form, except the potential energy term. Upon integrating over the domain 𝛺and assuming no-slip conditions 𝒖=𝟎on all boundaries, we obtain the global balances d𝐸𝑘 d𝑡= −∫𝛺 𝛷d𝛺+∫𝛺 𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣d𝛺 ,(11) d𝐸𝑖 d𝑡=∫𝛺 𝛷d𝛺+∫𝜕 𝛺 𝜆∇𝑇⋅𝒏d𝑆 ,(12) d𝐸 d𝑡=d𝐸𝑘 d𝑡+d𝐸𝑖 d𝑡=∫𝛺 𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣d𝛺+∫𝜕 𝛺 𝜆∇𝑇⋅𝒏d𝑆 ,(13) where 𝐸=∫𝛺𝑒d𝛺=𝐸𝑘+𝐸𝑖. In case the boundary conditions are adiabatic (∇𝑇⋅𝒏= 0), the last term in (13) vanishes and the total energy equation expresses that the sum of internal and kinetic energy changes due to the buoyancy flux ∫𝛺𝛽 𝑔 𝜌0(𝑇−𝑇0)𝑣d𝛺— this case will be dealt with in the Rayleigh–Taylor set-up in Section 6. Note that in compressible flows, the buoyancy flux can be written in terms of the time derivative of the potential energy — see Appendix D. In such a case, the Boussinesq system with the viscous dissipation function in the internal energy equation and with adiabatic boundaries can be considered to be truly energy-conserving, with the total energy being the sum of kinetic, internal and potential energy. In the incompressible case, this does not hold, so we use the term ‘energy-consistent’ to indicate that we have included the viscous dissipation function in the internal energy equation in such a way that the total energy equation is not affected by it. Table 1 Different non-dimensional forms resulting from different choices of 𝑢r ef . 𝑢r ef 𝛼1=𝜈 𝑢r ef 𝐻𝛼2=𝛽 𝑔 𝛥𝑇 𝐻 𝑢2 r ef 𝛼3=𝜈 𝑢r ef 𝑐 𝛥𝑇 𝐻𝛼4=𝜅 𝑢r ef 𝐻𝛾=𝛼1 𝛼3 I√𝛽 𝑔 𝛥𝑇 𝐻√Pr Ra 1Ge√Pr Ra 1 √Pr Ra 1 Ge II 𝜅 𝐻Pr Pr Ra Ge Ra 1Pr Ra Ge III √𝑐 𝛥𝑇 √Pr Ge Ra Ge √Pr Ge Ra √Ge Pr Ra 1 In most studies of Rayleigh–Bénard convection the dissipation function 𝛷is left out from the internal energy Eq. (4), while its corresponding counterpart in the momentum equation (𝜇∇2𝒖) 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. 2.3. Non-dimensionalization We non-dimensionalize Eqs. (1),(2) and (4) by taking a reference length 𝐻(cavity height), a reference temperature difference 𝛥𝑇 (difference between the cold and hot plates), and a reference velocity 𝑢r ef yet to be specified. These choices determine the time scale 𝐻∕𝑢r ef and the pressure scale 𝜌0𝑢2 r ef . An important question, which we will address here, is how the choice of non-dimensionalization changes the total energy equation. The non-dimensional equations are written as (for details, see Appendix C):  ∇⋅ 𝒖= 0,(14) 𝜕 𝒖 𝜕 𝑡+ ∇⋅( 𝒖⊗ 𝒖)= − ∇𝑝′+𝛼1 ∇2 𝒖+𝛼2 𝑇𝒆𝑦,(15) 𝜕 𝑇 𝜕 𝑡+ ∇⋅( 𝒖 𝑇) =𝛼3 𝛷+𝛼4 ∇2 𝑇 ,(16) where the parameters 𝛼𝑖,𝑖= 1 … 4are a function of the Rayleigh number Ra =𝛽 𝑔 𝛥𝑇 𝐻3 𝜈 𝜅, the Prandtl number Pr =𝜈 𝜅and the Gebhart number Ge =𝛽 𝑔 𝐻 𝑐(also known as the dissipation number [6]). In Table 1we present three different options for 𝑢r ef with the corresponding values of 𝛼. Choices I and II are common in literature, see for example [27] for choice I and [1,18,28] for choice II; they correspond to a freefall velocity scale and the thermal diffusivity scale, respectively. Other choices are also possible, e.g. 𝑢r ef =𝛽 𝑔 𝛥𝑇 𝐻2∕𝜈[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. Physically, this choice can be interpreted as the velocity that is obtained when internal energy is transformed into kinetic energy. The non-dimensional form of the total energy equation follows by taking the dot product of (15) with  𝒖and add the internal energy Eq. (16). The global energy balances in non-dimensional form read d 𝐸𝑘 d 𝑡= −𝛼1 𝛬∫ 𝛺  𝛷d 𝛺+𝛼2 𝛬∫ 𝛺  𝑇 𝑣 d 𝛺 ,(17) d 𝐸𝑖 d 𝑡=𝛼3 𝛬∫ 𝛺  𝛷d 𝛺+𝛼4 𝛬∫𝜕 𝛺  ∇ 𝑇⋅𝒏d 𝑆 ,(18) d 𝐸 d 𝑡=d 𝐸𝑘 d𝑡+𝛾d 𝐸𝑖 d𝑡=𝛼2 𝛬∫ 𝛺  𝑇 𝑣 d 𝛺+𝛾 𝛼4 𝛬∫𝜕 𝛺  ∇ 𝑇⋅𝒏d 𝑆 ,(19) where  𝐸=1 𝛬∫ 𝛺𝑒 d 𝛺,𝛬=𝐿∕𝐻is the aspect ratio of the box, and 𝛾=𝛼1 𝛼3is a weighting factor, which is reported in Table 1for different choices of 𝑢r ef . For definitions of 𝑒𝑘,𝑒𝑖and 𝑒, see Appendix C. The proposed choice III is the only choice that features 𝛾= 1, meaning that the dimensionless kinetic and internal energy equation are consistent with each other and do not require a weighting factor in order for the viscous dissipation term to cancel. The choice for a particular reference velocity typically depends on the problem at hand. Choices I and II have the advantage that in case Computers and Fluids 286 (2025) 106473 3 B. Sanderse and F.X. Trias of Ge = 0(most commonly investigated in literature), one obtains 𝛼3= 0and the dissipation terms simply drops from the internal energy equation. However, when Ge is small but nonzero, the weight factor 𝛾becomes very large for choices I and II. Choice III does not suffer from this issue, because 𝛾= 1independent of Ge, so kinetic energy and internal energy can be summed independent of Ge. However, choice III has the disadvantage that it does not work in the case Ge = 0, since it leads to 𝛼𝑖= 0for all 𝑖. In summary: for Ge = 0, choices I and II are preferred; for small but nonzero Ge, choice III is preferred; in other cases, all choices are fine. The discussion in the next sections will be agnostic for the choice of 𝑢ref, and expressed in terms of the general parameters 𝛼𝑖. Note that in the simulations in Sections 5–7, we will employ choice I. Choices II and III give equivalent results apart from scaling factors. 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 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 𝑦=𝑦′by 𝐹(𝑦′) ∶= 1 𝐿∫𝐿 0(𝜌0𝑐 𝑇 𝑣−𝜆𝜕 𝑇 𝜕 𝑦)(𝑥,𝑦′) d𝑥. (20) Then, the Nusselt number based on 𝐹follows as [3]: Nu( 𝑦′) ∶= 𝐹(𝑦′) 𝜆𝛥𝑇 ∕𝐻=1 𝛬∫𝛬 0(1 𝛼4  𝑇 𝑣 −𝜕 𝑇 𝜕 𝑦 )(𝑥, 𝑦′) d𝑥. (21) 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 Nu( 𝑦) = Nu( 𝑦 = 0) = Nu, which is a constant, independent of 𝑦′[1,28]. However, upon including viscous dissipation, this relation no longer holds true and instead the steady internal energy equation yields 𝛼4(Nu( 𝑦′) − Nu(0)) =𝛼3𝜖𝑈(𝑦′),(22) where the integrated dissipation function is given by 𝜖𝑈(𝑦′) ∶= 1 𝛬∫𝑦′ 0∫𝛬 0  𝛷d𝑥d𝑦. (23) Eq. (22) is an important relation which shows that (taking 𝑦′= 1) 𝛼4(Nu(1) − Nu(0)) =𝛼3𝜖𝑈(1),(24) 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, Eq. (17). The second term in the right-hand side of Eq. (17) can be rewritten with Eq. (22), following the analysis in [28]: 𝛼2 𝛬∫ 𝛺  𝑇 𝑣 d 𝛺=𝛼2 𝛬∫1 0∫𝛬 0  𝑇 𝑣 d𝑥d𝑦 =𝛼2𝛼4∫1 0 Nu( 𝑦) d𝑦 +𝛼2𝛼4 𝛬∫𝛬 0∫1 0 𝜕 𝑇 𝜕 𝑦 d𝑦d𝑥 =𝛼2𝛼4Nu(0) +𝛼2𝛼3∫1 0 𝜖𝑈(𝑦) d𝑦 +𝛼2𝛼4 𝛬∫𝛬 0 ( 𝑇(𝑥, 𝑦 = 1) − 𝑇(𝑥, 𝑦 = 0))d 𝑥 =𝛼2𝛼4(Nu(0) − 1) +𝛼2𝛼3∫1 0 𝜖𝑈(𝑦) d𝑦. (25) For (statistically) steady flow, this term equals the first term in the right-hand side of Eq. (17), yielding the second relation between the Nusselt number and the viscous dissipation 𝜖𝑈 𝛼2𝛼4(Nu(0) − 1) =𝛼1𝜖𝑈(1) −𝛼2𝛼3∫1 0 𝜖𝑈(𝑦) d𝑦. (26) We recognize the well-known equation 𝛼2𝛼4(Nu(0) − 1) =𝛼1𝜖𝑈(1), see e.g. [1], but with the additional negative term −𝛼2𝛼3∫1 0𝜖𝑈(𝑦) d𝑦. Lastly, we link the thermal dissipation 𝜖𝑇to the Nusselt number and the viscous dissipation function. The non-dimensional internal energy equation, Eq. (16), is multiplied by  𝑇, and after integrating by parts, using the skew-symmetry of the convective operator, and employing the boundary condition  𝑇(𝑦 = 1) = 0, one obtains 1 𝛬 d d𝑡∫ 𝛺 1 2 𝑇2d 𝛺=𝛼3 𝛬∫ 𝛺  𝑇 𝛷d 𝛺−𝛼4 𝛬∫𝛬 0( 𝑇𝜕 𝑇 𝜕 𝑦 )𝑦=0 d𝑥 −𝛼4 𝛬∫ 𝛺‖ ∇ 𝑇‖2d 𝛺 .(27) With the boundary condition  𝑇(𝑦 = 0) = 1, and the assumption of (statistically) steady flow, this relation is further simplified to 𝛼4Nu(0) =𝛼4𝜖𝑇−𝛼3 𝛬∫ 𝛺  𝑇 𝛷d 𝛺 ,(28) where 𝜖𝑇∶= 1 𝛬∫ 𝛺‖ ∇ 𝑇‖2d 𝛺 .(29) Since  𝑇≥0, 𝛷≥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 Nu =𝜖𝑇. In combination with Eq. (24), we obtain for the Nusselt number of the upper plate: 𝛼4Nu(1) =𝛼4𝜖𝑇+𝛼3 𝛬∫(1 − 𝑇) 𝛷d 𝛺 .(30) Assuming that the temperature satisfies 0≤ 𝑇≤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: Nu(0) ≤𝜖𝑇≤Nu(1).(31) The three relations (24),(26) and (28) are summarized in Table 2and will be confirmed in the numerical experiments in Section 5. 3. Energy-consistent spatial discretization 3.1. Mass, momentum and kinetic energy equation To discretize the non-dimensional mass and momentum Eqs. (14) and (15), we use the staggered-grid energy-conserving finite volume method described in [29], extended by including the buoyancy term in the momentum equations. This leads to the following semi-discrete equations: 𝑀 𝑉ℎ(𝑡) = 0,(32) 𝛺𝑉 d𝑉ℎ(𝑡) d𝑡= −𝐶𝑉(𝑉ℎ(𝑡)) −𝐺 𝑝ℎ(𝑡) +𝛼1𝐷𝑉𝑉ℎ(𝑡) +𝛼2(𝐴𝑇ℎ(𝑡) +𝑦𝑇).(33) Here, 𝑉ℎ∈R𝑁𝑉are the velocity unknowns, 𝑝ℎ∈R𝑁𝑝the pressure unknowns, and 𝑇ℎ∈R𝑁𝑝the temperature unknowns; see Fig. 2for their positioning. 𝑀∈R𝑁𝑝×𝑁𝑉is the discretized divergence operator, 𝐺= −𝑀𝑇∈R𝑁𝑉×𝑁𝑝the discretized gradient operator, 𝛺𝑉∈R𝑁𝑉×𝑁𝑉 a matrix with the ‘velocity’ finite volume sizes on its diagonal, and 𝐶𝑉 and 𝐷𝑉constitute central difference approximations of the convective and diffusive terms. 𝐴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 𝑦𝑇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: 𝑉𝑇 ℎ𝐶𝑉(𝑉ℎ) = 0and 𝑉𝑇 ℎ𝐺 𝑝ℎ= 0, just like in the continuous case. This is achieved by using a Computers and Fluids 286 (2025) 106473 4 B. Sanderse and F.X. Trias Table 2 Steady-state Nusselt number relations, with and without viscous dissipation. Origin Without viscous dissipation With viscous dissipation Internal Nu(1) = Nu(0) 𝛼4(Nu(1) − Nu(0)) =𝛼3𝜖𝑈(1) Kinetic 𝛼2𝛼4(Nu(0) − 1) =𝛼1𝜖𝑈(1) 𝛼2𝛼4(Nu(0) − 1) =𝛼1𝜖𝑈(1) −𝛼2𝛼3∫1 0𝜖𝑈(𝑦) d𝑦 Internal energy ×𝑇Nu(0) =𝜖𝑇𝛼4Nu(0) =𝛼4𝜖𝑇−𝛼3 𝛬∫ 𝛺 𝑇 𝛷d 𝛺 Fig. 2. Staggered grid with positioning of unknowns around a pressure volume. skew-symmetric convection operator and the compatibility between 𝑀 and 𝐺via 𝐺= −𝑀𝑇. The discrete kinetic energy balance then reads: d𝐸𝑘,ℎ d𝑡= −𝛼1𝜖𝑈 ,ℎ +𝛼2𝑉𝑇 ℎ(𝐴𝑇ℎ+𝑦𝑇),(34) where 𝐸𝑘,ℎ =1 2𝑉𝑇 ℎ𝛺𝑉𝑉ℎ. The global viscous dissipation (i.e. summed over the entire domain) is given by 𝜖𝑈 ,ℎ =‖𝑄𝑉ℎ‖2 2>0, where 𝑄stems from decomposing the symmetric negative-definite diffusive operator as 𝐷𝑉= −𝑄𝑇𝑄. Eq. (34) is the semi-discrete counterpart of Eq. (17). 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 Eq. (16) which is such that discrete versions of the global balances (12) and (13) 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: 𝑘𝑖,𝑗 ∶= 1 4𝑢2 𝑖+1∕2,𝑗 +1 4𝑢2 𝑖−1∕2,𝑗 +1 4𝑣2 𝑖,𝑗+1∕2 +1 4𝑣2 𝑖,𝑗−1∕2.(35) This choice gives a local kinetic energy equation that is consistent with the continuous equations, as is detailed in Appendix B, and consistent with the global energy definition. The expression for 𝛷ℎthen follows from differentiating the expression for 𝑘𝑖𝑗 in time, substituting the momentum equations, and rewriting the terms involving the diffusive operator (see Appendix B). The implied dissipation then follows by constructing a discrete version of (9). As example, we construct the discrete version of 𝑢𝜕2𝑢 𝜕 𝑥2= −(𝜕 𝑢 𝜕 𝑥)2 +𝜕 𝜕 𝑥(𝑢𝜕 𝑢 𝜕 𝑥), being 𝑢𝑖+1∕2,𝑗 𝛥𝑥 (𝑢𝑖+3∕2,𝑗 −𝑢𝑖+1∕2,𝑗 𝛥𝑥 −𝑢𝑖+1∕2,𝑗 −𝑢𝑖−1∕2,𝑗 𝛥𝑥 ) = −1 2(𝑢𝑖+3∕2,𝑗 −𝑢𝑖+1∕2,𝑗 𝛥𝑥 )2 −1 2(𝑢𝑖+1∕2,𝑗 −𝑢𝑖−1∕2,𝑗 𝛥𝑥 )2 +1 𝛥𝑥 (1 2(𝑢𝑖+3∕2,𝑗 +𝑢𝑖+1∕2,𝑗 )𝑢𝑖+3∕2,𝑗 −𝑢𝑖+1∕2,𝑗 𝛥𝑥 −1 2(𝑢𝑖+1∕2,𝑗 +𝑢𝑖−1∕2,𝑗 )𝑢𝑖+1∕2,𝑗 −𝑢𝑖−1∕2,𝑗 𝛥𝑥 ).(36) The first two terms on the right-hand side contribute to the viscous dissipation function. Repeating this process for the other components (𝑢𝜕2𝑢 𝜕 𝑦2,𝑣𝜕2𝑣 𝜕 𝑥2,𝑣𝜕2𝑣 𝜕 𝑦2), as outlined in Appendix B.2, yields the following novel expression for the local dissipation function: 𝛷𝑖,𝑗 =1 2𝛷𝑢 𝑖+1∕2,𝑗 +1 2𝛷𝑢 𝑖−1∕2,𝑗 +1 2𝛷𝑣 𝑖,𝑗+1∕2 +1 2𝛷𝑣 𝑖,𝑗−1∕2 ,(37) where 𝛷𝑢 𝑖+1∕2,𝑗 = −1 2(𝑢𝑖+3∕2,𝑗 −𝑢𝑖+1∕2,𝑗 𝛥𝑥 )2 −1 2(𝑢𝑖+1∕2,𝑗 −𝑢𝑖−1∕2,𝑗 𝛥𝑥 )2 −1 2(𝑢𝑖+1∕2,𝑗+1 −𝑢𝑖+1∕2,𝑗 𝛥𝑦 )2 −1 2(𝑢𝑖+1∕2,𝑗 −𝑢𝑖+1∕2,𝑗−1 𝛥𝑦 )2 ,(38) 𝛷𝑣 𝑖,𝑗+1∕2 = −1 2(𝑣𝑖+1,𝑗+1∕2 −𝑣𝑖,𝑗+1∕2 𝛥𝑥 )2 −1 2(𝑣𝑖+1,𝑗−1∕2 −𝑣𝑖,𝑗−1∕2 𝛥𝑥 )2 −1 2(𝑣𝑖,𝑗+3∕2 −𝑣𝑖,𝑗+1∕2 𝛥𝑦 )2 −1 2(𝑣𝑖,𝑗+1∕2 −𝑣𝑖,𝑗−1∕2 𝛥𝑦 )2 .(39) At boundaries, an adaptation of 𝛷ℎis required in order to have a discrete equivalent of Eq. (9). This is detailed in Eq. (B.15). Note that 𝛷ℎis derived based on local energy consideration which upon summation equals the global dissipation, just like Eq. (23): 1𝑇𝛺𝑝𝛷ℎ=𝜖𝑈 ,ℎ.(40) 3.3. Internal energy equation Having proposed a consistent expression for 𝛷ℎ, the spatial discretization of the internal energy Eq. (16) reads: 𝛺𝑝 d𝑇ℎ d𝑡= −𝐶𝑇(𝑉ℎ, 𝑇ℎ) +𝛼3𝛺𝑝𝛷ℎ(𝑉ℎ) +𝛼4(𝐷𝑇𝑇ℎ+𝑦𝑇),(41) where [𝐶𝑇(𝑉ℎ, 𝑇ℎ)]𝑖,𝑗 =𝛥𝑦 (𝑢𝑖+1∕2,𝑗 1 2(𝑇𝑖+1,𝑗 +𝑇𝑖,𝑗 ) −𝑢𝑖−1∕2,𝑗 1 2(𝑇𝑖,𝑗 +𝑇𝑖−1,𝑗 ))+ 𝛥𝑥 (𝑣𝑖,𝑗+1∕2 1 2(𝑇𝑖,𝑗+1 +𝑇𝑖,𝑗 ) −𝑣𝑖,𝑗−1∕2 1 2(𝑇𝑖,𝑗 +𝑇𝑖,𝑗−1))(42) 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. 𝐷𝑇the standard second-order difference stencil with boundary conditions encoded in 𝑦𝑇. The total internal energy is given by 𝐸𝑖,ℎ = 1𝑇𝛺𝑝𝑇ℎ(simply summing over all finite volumes). Due to the no-slip boundary conditions on the velocity field, the convective operator satisfies 1𝑇𝐶𝑇(𝑉ℎ, 𝑇ℎ) = 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 d𝐸𝑖,ℎ d𝑡=𝛼31𝑇𝛺𝑝𝛷ℎ+𝛼41𝑇(𝐷𝑇𝑇ℎ+𝑦𝑇), =𝛼31𝑇𝛺𝑝𝛷ℎ+𝛼4(Nu𝐻− Nu𝐶), (43) Computers and Fluids 286 (2025) 106473 5 B. Sanderse and F.X. Trias Fig. 3. Steady-state temperature field for Ra = 105on a 128 ×128 grid, for different Ge. where in the second line the Nusselt numbers are instantaneous Nusselt numbers. Upon adding the total kinetic energy Eq. (34), and using property (40), the global energy balance results: d𝐸ℎ d𝑡=d𝐸𝑘,ℎ d𝑡+𝛾d𝐸𝑖,ℎ d𝑡=𝛼2𝑉𝑇 ℎ(𝐴𝑇ℎ+𝑦𝑇) +𝛾 𝛼41𝑇(𝐷𝑇𝑇ℎ+𝑦𝑇), =𝛼2𝑉𝑇 ℎ(𝐴𝑇ℎ+𝑦𝑇) +𝛾 𝛼4(Nu𝐻− Nu𝐶), (44) which is the semi-discrete counterpart of Eq. (19). 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, independent of the mesh size. 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 (24) and (28). 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, 𝛼31𝑇𝛺𝑝𝛷ℎ(𝑉ℎ) +𝛼41𝑇(𝐷𝑇𝑇ℎ+𝑦𝑇) = 0.(45) The second term can be simplified as 1𝑇(𝐷𝑇𝑇ℎ+𝑦𝑇) = − 𝑁𝑥 ∑ 𝑖=1 𝑇𝑖,1−𝑇𝐻 1 2𝛥𝑦 𝛥𝑥+ 𝑁𝑥 ∑ 𝑖=1 𝑇𝐶−𝑇𝑖,𝑁𝑦 1 2𝛥𝑦 𝛥𝑥 = Nu𝐻− Nu𝐶,(46) where the Nusselt numbers on the lower (hot) and upper (cold) plate are defined as Nu𝐻∶= − 𝑁𝑥 ∑ 𝑖=1 𝑇𝑖,1−𝑇𝐻 1 2𝛥𝑦 𝛥𝑥, (47) Nu𝐶∶= − 𝑁𝑥 ∑ 𝑖=1 𝑇𝐶−𝑇𝑖,𝑁𝑦 1 2𝛥𝑦 𝛥𝑥. (48) This leads to the discrete version of (24): 𝛼4(Nu𝐶− Nu𝐻) =𝛼31𝑇𝛺𝑝𝛷ℎ(𝑉ℎ).(49) The discrete version of (28) follows by considering the inner product of Eq. (41) with 𝑇𝑇 ℎinstead of 1𝑇. An important property of the convective discretization (42) is that 𝑇𝑇 ℎ𝐶𝑇(𝑉ℎ, 𝑇ℎ) = 0,∀𝑇ℎ,if 𝑀 𝑉ℎ= 0.(50) This property is most easily derived by recognizing that 𝐶ℎ(𝑉ℎ, 𝑇ℎ)can be written in terms of a matrix–vector product  𝐶𝑇(𝑉ℎ)𝑇ℎ, where  𝐶𝑇(𝑉ℎ) is skew-symmetric if 𝑀 𝑉ℎ= 0. In addition, the inner product of 𝑇ℎwith the diffusive terms can be written as 𝑇𝑇 ℎ(𝐷𝑇𝑇ℎ+𝑦𝑇) = 𝑁𝑥 ∑ 𝑖=1 ⎛⎜⎜⎝ −𝑇𝐻 𝑇𝑖,1−𝑇𝐻 1 2𝛥𝑦 +𝑇𝐶 𝑇𝐶−𝑇𝑖,𝑁𝑦 1 2𝛥𝑦 ⎞⎟⎟⎠ 𝛥𝑥 −𝜖𝑇 ,ℎ,(51) where 𝜖𝑇 ,ℎ ∶= 𝑁𝑥 ∑ 𝑖=1 ⎛⎜⎜⎜⎝ 1 2⎛⎜⎜⎝ 𝑇𝑖,1−𝑇𝐻 1 2𝛥𝑦 ⎞⎟⎟⎠ 2 + 𝑁𝑦 ∑ 𝑗=2 (𝑇𝑖,𝑗 −𝑇𝑖,𝑗−1 𝛥𝑦 )2 +1 2⎛⎜⎜⎝ 𝑇𝐶−𝑇𝑖,𝑁𝑦 1 2𝛥𝑦 ⎞⎟⎟⎠ 2⎞⎟⎟⎟⎠ ×𝛥𝑥𝛥𝑦 + 𝑁𝑦 ∑ 𝑗=1 𝑁𝑥 ∑ 𝑖=2 (𝑇𝑖,𝑗 −𝑇𝑖−1,𝑗 𝛥𝑥 )2 𝛥𝑥𝛥𝑦 (52) is the discrete analogue of (27) and Eq. (51) is the discrete version of ∫𝑇d2𝑇 d𝑦2= [𝑇d𝑇 d𝑦] −∫(d𝑇 d𝑦)2. With the boundary condition 𝑇𝐻= 1,𝑇𝐶= 0, we get the balance 𝛼4Nu𝐻=𝛼4𝜖𝑇 ,ℎ −𝛼3𝑇𝑇 ℎ𝛺𝑝𝛷ℎ(𝑉ℎ),(53) which is the discrete version of Eq. (28). 4. Energy-consistent temporal discretization The system of Eqs. (32),(33) and (41) needs to be integrated in time with a suitable method in order to preserve a time-discrete version of the global energy balance (44). 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 [23,28,30], or an explicit method for both convection and diffusion [31,32]. 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 pressurecorrection 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 [33]. Instead, we show here that the implicit midpoint method can be employed to achieve energy-consistent time integration. The fully discrete system reads: 𝑀 𝑉𝑛+1∕2 ℎ= 0,(54) 𝛺𝑉 𝑉𝑛+1 ℎ−𝑉𝑛 ℎ 𝛥𝑡 = −𝐶𝑉(𝑉𝑛+1∕2 ℎ) −𝐺 𝑝𝑛+1∕2 ℎ+𝛼1𝐷𝑉𝑉𝑛+1∕2 ℎ +𝛼2(𝐴𝑇 𝑛+1∕2 ℎ+𝑦𝑇),(55) Computers and Fluids 286 (2025) 106473 6 B. Sanderse and F.X. Trias 𝛺𝑝 𝑇𝑛+1 ℎ−𝑇𝑛 ℎ 𝛥𝑡 = −𝐶𝑇(𝑉𝑛+1∕2 ℎ, 𝑇𝑛+1∕2 ℎ) +𝛼3𝛺𝑝𝛷(𝑉𝑛+1∕2 ℎ) +𝛼4(𝐷𝑇𝑇𝑛+1∕2 ℎ+𝑦𝑇).(56) Here 𝑉𝑛+1∕2 ℎ=1 2(𝑉𝑛 ℎ+𝑉𝑛+1 ℎ)and 𝑇𝑛+1∕2 ℎ=1 2(𝑇𝑛 ℎ+𝑇𝑛+1 ℎ). Upon multiplying (55) by (𝑉𝑛+1∕2 ℎ)𝑇and (56) by 1𝑇, and adding the two resulting equations, we get the discrete energy balance, 𝐸𝑛+1 ℎ−𝐸𝑛 ℎ 𝛥𝑡 = 𝐸𝑛+1 𝑘,ℎ −𝐸𝑛 𝑘,ℎ 𝛥𝑡 +𝛾 𝐸𝑛+1 𝑖,ℎ −𝐸𝑛 𝑖,ℎ 𝛥𝑡 =𝛼2(𝑉𝑛+1∕2 ℎ)𝑇(𝐴𝑇 𝑛+1∕2 ℎ+𝑦𝑇) +𝛾 𝛼41𝑇(𝐷𝑇𝑇𝑛+1∕2 ℎ+𝑦𝑇),(57) which is the fully-discrete counterpart of Eqs. (19) and (44). The derivations hinges again on skew-symmetry of the convection operator 𝐶𝑉, the compatibility between 𝑀and 𝐺(𝐺= −𝑀𝑇), and the consistency requirement on the viscous dissipation function, Eq. (40). Again, we stress that the discrete energy equation results from exactly balancing the viscous dissipation between the kinetic and internal energy equations, independent of the mesh size and the time step. The system of Eqs. (54)–(56) leads to a large system of nonlinear equations which has a saddle point structure due to the divergencefree 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 commonly used for direct numerical simulations [31,32] in Section 6. This one-leg scheme reads 𝑀 𝑉𝑛+1 ℎ= 0,(58) 𝛺𝑉 (𝑏+1 2)𝑉𝑛+1 ℎ− 2𝑏𝑉 𝑛 ℎ+ (𝑏−1 2)𝑉𝑛−1 ℎ 𝛥𝑡 = −𝐶𝑉(𝑉∗ ℎ) −𝐺 𝑝𝑛+1 ℎ +𝛼1𝐷𝑉𝑉∗ ℎ+𝛼2(𝐴𝑇 ∗ ℎ+𝑦𝑇), (59) 𝛺𝑝 (𝑏+1 2)𝑇𝑛+1 ℎ− 2𝑏𝑇 𝑛 ℎ+ (𝑏−1 2)𝑇𝑛−1 ℎ 𝛥𝑡 = −𝐶𝑇(𝑉∗ ℎ, 𝑇∗ ℎ) +𝛼3𝛺𝑝𝛷(𝑉∗ ℎ) +𝛼4(𝐷𝑇𝑇∗ ℎ+𝑦𝑇), (60) where 𝑉∗ ℎ= (1 +𝑏)𝑉𝑛 ℎ−𝑏𝑉 𝑛−1 ℎand 𝑇∗ ℎ= (1 +𝑏)𝑇𝑛 ℎ−𝑏𝑇 𝑛−1 ℎ, and we will take 𝑏=1 2. The energy balance for this scheme follows again from multiplying the momentum equations by (𝑉𝑛+1∕2 ℎ)𝑇and the internal energy equation by 1𝑇, and adding the two, leading to 𝐸𝑛+1 ℎ−𝐸𝑛 ℎ 𝛥𝑡 = −(𝑉𝑛+1∕2 ℎ)𝑇𝐶𝑉(𝑉∗ ℎ) +𝛼1(𝑉𝑛+1∕2 ℎ𝐷𝑉𝑉∗ ℎ+ 1𝑇𝛺𝑝𝛷(𝑉∗ ℎ)) +𝛼2(𝑉𝑛+1∕2 ℎ)𝑇(𝐴𝑇 ∗ ℎ+𝑦𝑇) +𝛾 𝛼41𝑇(𝐷𝑇𝑇𝑛+1∕2 ℎ+𝑦𝑇).(61) We observe that the explicit nature of the one-leg scheme introduces two errors in the energy equation. First, the convective terms do not cancel from the energy equation, because (𝑉𝑛+1∕2 ℎ)𝑇𝐶𝑉(𝑉∗ ℎ)is not equal to zero. Second, a contribution from the viscous dissipation function appears as (𝑉𝑛+1∕2 ℎ)𝑇𝐷𝑉𝑉∗ ℎdoes not exactly cancel 1𝑇𝛺𝑝𝛷(𝑉∗ ℎ). 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: 𝑢(𝑥, 𝑦) = −64𝑥2(𝑥− 1)2𝑦(𝑦− 1)(2𝑦− 1),(62) 𝑣(𝑥, 𝑦) = 64𝑥(𝑥− 1)(2𝑥− 1)𝑦2(𝑦− 1)2,(63) Table 3 Convergence of Nusselt number (47) with grid refinement for different Rayleigh numbers and Ge = 0. Grid Ra = 103Ra = 104Ra = 105 3221.000 2.170 3.933 6421.000 2.161 3.916 12821.000 2.159 3.912 25621.000 2.158 3.911 Cai et al. [35] (2562) 1.000 2.158 3.911 Table 4 Convergence of Nusselt numbers (47) and (48) with grid refinement for different Rayleigh and different Gebhart numbers. (a) Ra = 104(b) Ra = 105 Grid Ge = 0.1 Ge = 1Grid Ge = 0.1 Ge = 1 Nu𝐻Nu𝐶Nu𝐻Nu𝐶Nu𝐻Nu𝐶Nu𝐻Nu𝐶 3222.111 2.228 1.582 2.729 3223.786 4.080 2.448 5.319 6422.103 2.219 1.578 2.716 6423.770 4.062 2.441 5.299 12822.101 2.217 1.576 2.713 12823.766 4.057 2.439 5.293 25622.100 2.216 1.576 2.712 25623.765 4.056 2.439 5.292 which is inspired by the regularized driven cavity problem [34]. 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 (𝒖= 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) Fig. 3(a) shows the temperature field when viscous dissipation is not included (Ge = 0). The resulting Nusselt numbers as a function of grid refinement are displayed in Table 3and indicate excellent agreement with literature [35]. We note that the Nusselt numbers as defined by (47) and (48) 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 (49). Note also that we only report Nu𝐻since Nu𝐶= Nu𝐻up to machine precision. 5.2. Grid convergence study for viscous dissipation case (Ge >0) When including viscous dissipation (Ge >0) in the internal energy equation, the flow field changes qualitatively and loses its symmetric nature, as can be observed in Figs. 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 Eq. (49) (or (24)). This is reported in Table 4and Fig. 4(a). The critical Rayleigh number that we find from the bifurcation diagram is Ra𝑐≈ 2585, which is in excellent agreement with the value of 2585.02 reported in literature [36,37]. It is independent of the value of the Prandtl number, as shown in [36], and also independent of the value of the Gebhart number. This latter fact follows by extending the linear stability analysis in [36] and realizing that the term ∇𝒖∶ ∇𝒖with 𝒖=𝒖0+𝜀𝒖′ and background state 𝒖0= 0leads to the term 𝜀2∇𝒖′∶ ∇𝒖′, which disappears when gathering terms of (𝜀). The results in Fig. 4(a) show indeed that the bifurcation point is the same for different values of Ge. Fig. 4(b) shows a different interpretation of the Nusselt number, indicating the relation with the thermal dissipation and viscous dissipation according to Eq. (53) (or (28)). The results confirm that the thermal dissipation lies in between the Nusselt number of the hot and cold plate. Computers and Fluids 286 (2025) 106473 7 B. Sanderse and F.X. Trias Fig. 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 𝒖= 0and adiabatic conditions 𝜕 𝑇 𝜕 𝑛= 0on 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 𝐸𝑛+1 ℎ−𝐸𝑛 ℎ 𝛥𝑡 = 𝐸𝑛+1 𝑘,ℎ −𝐸𝑛 𝑘,ℎ 𝛥𝑡 +𝛾 𝐸𝑛+1 𝑖,ℎ −𝐸𝑛 𝑖,ℎ 𝛥𝑡 =𝛼2(𝑉𝑛+1∕2 ℎ)𝑇(𝐴𝑇 𝑛+1∕2 ℎ+𝑦𝑇). (64) 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 Fig. 5. The energyconserving implicit midpoint (‘IM’) method detailed in Section 4will be compared to the explicit one-leg (‘OL’) method commonly used in DNS studies [31,32] (where we take 𝑏=1 2and a fixed time step). The domain size is 1×2, the grid is 64 ×128, the time step 𝛥𝑡 = 5⋅10−3 and the end time 𝑇= 50. We consider the case Pr = 0.71, Ra = 106and Ge = {0.1,1}. The instability does naturally arise due to growth of round-off errors, but this takes rather long, so instead a perturbation is added to the initial interface: 𝑦= 1 + 0.05 sin(2𝜋 𝑥). The instability quickly develops and an asymmetry in the solution appears, triggering a sequence of well-known ‘mushroom’ type plumes: hot plumes rising upward and cold plumes sinking downward (Fig. 6). The development of the instability is essentially the same for IM and OL — see also Fig. 7(a) for a more quantitative comparison. Note that if no perturbation is added, the onset of stability is sensitive to the choice of time integration method, due to differences in the accumulation of round-off errors. Fig. 6also shows that the initial development of the instability is insensitive to the value of Ge, just like the bifurcation point in the steady state Rayleigh–Bénard simulation was insensitive to the value of Ge. 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 Fig. 7(a), where the average temperature is displayed. Compared to the initial temperature difference 𝛥𝑇 = 1, the relative temperature increase is Fig. 5. Problem set-up with initial condition for Rayleigh Taylor problem. about 2% for Ge = 0.1and more than 20% for 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 Fig. 7(b), which shows the energy error 𝜀𝐸∶= |||||| 𝐸𝑛+1 𝑘,ℎ −𝐸𝑛 𝑘,ℎ 𝛥𝑡 +𝛾 𝐸𝑛+1 𝑖,ℎ −𝐸𝑛 𝑖,ℎ 𝛥𝑡 −𝛼2(𝑉𝑛+1∕2 ℎ)𝑇(𝐴𝑇 𝑛+1∕2 ℎ+𝑦𝑇)|||||| .(65) For IM the error remains at the tolerance with which we solve the system of nonlinear equations (10−12). For OL, the error is around (10−6)when the instability is most pronounced (around 𝑡= 5, see Fig. 6), and decreases when the flow settles back to a steady state. This small energy error of the OL scheme seems acceptable given that OL is roughly 4–5×less expensive than IM, because IM requires roughly 4–5 iterations (Poisson solves) per time step, instead of only 1 for OL. 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. Computers and Fluids 286 (2025) 106473 8 B. Sanderse and F.X. Trias Fig. 6. Rayleigh–Taylor temperature fields at 𝑡= 5for different Ge and different time integration methods (IM =Implicit Midpoint, OL =One-Leg scheme). Fig. 7. Rayleigh–Taylor results, IM =Implicit Midpoint, OL =One-Leg scheme. 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, Ra = 108and 1010. Direct numerical simulations (DNS) were carried out and analyzed in previous studies [38,39] without taking into account the viscous dissipation effects (Ge = 0). Here, the results are extended to Ge = 0.1and Ge = 1keeping the same domain size (𝜋× 1 × 1) and mesh resolution (400 × 208 × 208 for Ra = 108, and 1024 ×768 ×768 for Ra = 1010). Grids are constructed with a uniform grid spacing in the periodic 𝑥-direction whereas wall-normal points (𝑦and 𝑧directions) are distributed following a hyperbolic-tangent function as follows (identical for the 𝑧-direction) 𝑦𝑖=1 2(1 +t anh (𝛾𝑦(2(𝑖− 1)∕𝑁𝑦− 1)) t anh 𝛾𝑦), 𝑖= 1,…, 𝑁𝑦+ 1,(66) where 𝑁𝑦and 𝛾𝑦are the number of control volumes and the concentration factor in the 𝑦-direction, respectively. In our case, 𝛾𝑦=𝛾𝑧= 1.4 for Ra = 108and 𝛾𝑦=𝛾𝑧= 1.6for Ra = 1010. For further details, the reader is referred to our previous works [38,39]. Instantaneous temperature fields corresponding to the statistically steady state are displayed in Fig. 8. As expected, thermal dissipation effects at Ge = 1lead 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 midheight plane is lost for Ge >0leading 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 Ge = 0. This implies that mesh resolution requirements in the nearwall 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. All simulations have been carried out for 500 time-units starting from a zero velocity field and uniformly distributed random temperatures between 𝑇𝐶and 𝑇𝐻. 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 Fig. 9where the time-evolution of various rate-of-changes of energy are shown. Results correspond to Ra = 108and Ge = 1using a very fine (400 × 208 × 208 ≈ 17.3M) 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 d𝐸𝑘,ℎ∕d𝑡(in red) fluctuates around zero. Only two terms contribute to the global kinetic energy of the system (see Eq. (34)): the global viscous dissipation, 𝜖𝑢,ℎ (in yellow), and the contribution of the buoyancy forces given by 𝛼2𝑉𝑇 ℎ(𝐴𝑇ℎ(𝑡) +𝑦𝑇)(in blue). These two contributions cancel each other on average when a statistically steady Computers and Fluids 286 (2025) 106473 9 B. Sanderse and F.X. Trias [8] Hewitt JM, Mckenzie DP, Weiss NO. Dissipative heating in convective flows. J Fluid Mech 1975;68(4):721–38. http://dx.doi.org/10.1017/ S002211207500119X. [9] Turcotte DL, Hsui AT, Torrance KE, Schubert G. Influence of viscous dissipation on Bénard convection. J Fluid Mech 1974;64(2):369–74. http://dx.doi.org/10. 1017/S0022112074002448. [10] Gebhart B. Effects of viscous dissipation in natural convection. J Fluid Mech 1962;14(2):225–32. http://dx.doi.org/10.1017/S0022112062001196. [11] Ostrach S. 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; 1952. [12] Blankenbach B, Busse F, Christensen U, Cserepes L, Gunkel D, Hansen U, Harder H, Jarvis G, Koch M, Marquart G, Moore D, Olson P, Schmeling H, Schnaubelt T. A benchmark comparison for mantle convection codes. Geophys J Int 1989;98(1):23–38. http://dx.doi.org/10.1111/j.1365-246X.1989.tb05511.x. [13] King SD, Lee C, van Keken PE, Leng W, Zhong S, Tan E, Tosi N, Kameyama MC. A community benchmark for 2-D cartesian compressible convection in the Earth’s mantle. Geophys J Int 2010;180(1):73–87. http://dx.doi.org/10.1111/j.1365246X.2009.04413.x. [14] Sugiyama K, Calzavarini E, Grossmann S, Lohse D. Flow organization in twodimensional non-Oberbeck-Boussinesq Rayleigh-Bénard convection in water. J Fluid Mech 2009;637:105–35. [15] Najm HN, Wyckoff PS, Knio OM. A semi-implicit numerical scheme for reacting flow: I. Stiff chemistry. J Comput Phys 1998;143(2):381–402. http://dx.doi.org/ 10.1006/jcph.1997.5856. [16] Nemati H. Direct numerical simulation of turbulent heat transfer to fluids at supercritical pressures (Ph.D. thesis), Delft University of Technology; 2016. [17] Barletta A. Comments on a paradox of viscous dissipation and its relation to the Oberbeck–Boussinesq approach. Int J Heat Mass Transfer 2008;51(25–26):6312– 6. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2007.10.044. [18] Barletta A, Nield D. Effect of pressure work and viscous dissipation in the analysis of the Rayleigh–Bénard problem. Int J Heat Mass Transfer 2009;52(13–14):3279– 89. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2009.02.005. [19] Barletta A, Celli M, Nield DA. On the onset of dissipation thermal instability for the Poiseuille flow of a highly viscous fluid in a horizontal channel. J Fluid Mech 2011;681:499–514. http://dx.doi.org/10.1017/jfm.2011.213. [20] Barletta A, Celli M, Brandão PV. On mixed convection in a horizontal channel, viscous dissipation and flow duality. Fluids 2022;7(5):170. http://dx.doi.org/10. 3390/fluids7050170. [21] Winters KB, Lombard PN, Riley JJ, D’Asaro EA. Available potential energy and mixing in Density-Stratified fluids. J Fluid Mech 1995;289:115–28. http: //dx.doi.org/10.1017/S002211209500125X. [22] Hughes GO, Gayen B, Griffiths RW. Available potential energy in Rayleigh– bénard convection. J Fluid Mech 2013;729:R3. http://dx.doi.org/10.1017/jfm. 2013.353. [23] Gayen B, Hughes GO, Griffiths RW. Completing the mechanical energy pathways in turbulent Rayleigh-Bénard convection. Phys Rev Lett 2013;111(12):124301. http://dx.doi.org/10.1103/PhysRevLett.111.124301. [24] Petschel K, Stellmach S, Wilczek M, Lülff J, Hansen U. Kinetic energy transport in Rayleigh–Bénard convection. J Fluid Mech 2015;773:395–417. http://dx.doi. org/10.1017/jfm.2015.216. [25] van Gils DPM, Bruggert G-W, Lathrop DP, Sun C, Lohse D. The twente turbulent Taylor–Couette (T3C) facility: Strongly turbulent (multiphase) flow between two independently rotating cylinders. Rev Sci Instrum 2011;82(2):025105. http: //dx.doi.org/10.1063/1.3548924. [26] Barletta A. Local energy balance, specific heats and the Oberbeck–Boussinesq approximation. Int J Heat Mass Transfer 2009;52(21–22):5266–70. http://dx. doi.org/10.1016/j.ijheatmasstransfer.2009.06.006. [27] van der Poel EP, Stevens RJAM, Lohse D. Comparison between twoand threedimensional Rayleigh–Bénard convection. J Fluid Mech 2013;736:177–94. http: //dx.doi.org/10.1017/jfm.2013.488. [28] Hepworth BJ. Nonlinear two-dimensional Rayleigh-bénard convection (Ph.D. thesis), University of Leeds; 2014. [29] Sanderse B. 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. [30] Sugiyama K, Calzavarini E, Grossmann S, Lohse D. Flow organization in twodimensional non-Oberbeck–Boussinesq Rayleigh–Bénard convection in water. J Fluid Mech 2009;637:105–35. http://dx.doi.org/10.1017/S0022112009008027. [31] Verstappen R, Veldman A. Symmetry-preserving discretization of turbulent flow. J Comput Phys 2003;187(1):343–68. http://dx.doi.org/10.1016/S0021-9991(03) 00126-8. [32] Trias FX, Lehmkuhl O. A self-adaptive strategy for the time integration of NavierStokes equations. Numer Heat Transfer B 2011;60(2):116–34. http://dx.doi.org/ 10.1080/10407790.2011.594398. [33] Sanderse B. Energy-conserving Runge–Kutta methods for the incompressible Navier–Stokes equations. J Comput Phys 2013;233:100–31. http://dx.doi.org/ 10.1016/j.jcp.2012.07.039. [34] Shih TM, Tan CH, Hwang BC. Effects of grid staggering on numerical schemes. Internat J Numer Methods Fluids 1989;9(2):193–212. http://dx.doi.org/10.1002/ fld.1650090206. [35] Cai W, Ma H, Wang Y, Chen J, Zheng X, Zhang H. Development of POD reducedorder model and its closure scheme for 2D Rayleigh–Bénard convection. Appl Math Model 2019;66:562–75. http://dx.doi.org/10.1016/j.apm.2018.09.031. [36] Gelfgat AY. Different modes of Rayleigh–Bénard Instability in twoand threedimensional rectangular enclosures. J Comput Phys 1999;156(2):300–24. http: //dx.doi.org/10.1006/jcph.1999.6363. [37] Venturi D, Wan X, Karniadakis GE. Stochastic bifurcation analysis of Rayleigh– Bénard convection. J Fluid Mech 2010;650:391–413. http://dx.doi.org/10.1017/ S0022112009993685. [38] Dabbagh F, Trias FX, Gorobets A, Oliva A. On the evolution of flow topology in turbulent Rayleigh-Bénard convection. Phys Fluids 2016;28:115105. [39] Dabbagh F, Trias FX, Gorobets A, Oliva A. Flow topology dynamics in a threedimensional phase space for turbulent Rayleigh-Bénard convection. Phys Rev Fluids 2020;5:024603. [40] Trias FX, Verstappen RWCP, Gorobets A, Soria M, Oliva A. Parameter-free symmetry-preserving regularization modeling of a turbulent differentially heated cavity. Comput & Fluids 2010;39:1815–31. [41] Smith W. All things flow. Oregon State University; 2019. [42] Tailleux R. On the energetics of stratified turbulent mixing, irreversible thermodynamics, Boussinesq models and the ocean heat engine controversy. J Fluid Mech 2009;638:339–82. http://dx.doi.org/10.1017/S002211200999111X. Computers and Fluids 286 (2025) 106473 16