Determination of Numerical Errors in Under-Resolved DNS of Turbulent Non-isothermal Flows

Abstract

The methodology to quantify the numerical dissipation in Under-resolved DNS (UDNS) based on the balance of the kinetic energy equation by Schranner et al. (Comput Fluids 114:84–97, 2015), has been examined in this study. Furthermore, this methodology has been extended for considering active scalars, based on both balance of the kinetic energy and thermal variance equations for assessing the quality of UDNS. As a first step towards the analysis of reactive flows, where temperature is actively coupled to the momentum equation, the turbulent Rayleigh–Bénard convection in a cubical enclosure is considered here. The simulations have been carried out for different Prandtl numbers ($Pr=0.1, 1.0, 10$) to analyse different levels of momentum-temperature coupling for a representative value of nominal Rayleigh number ($Ra=10^7$) in order to investigate the effects of grid resolution and discretisation schemes on the numerical dissipation. It has been found that the numerical dissipation for both kinetic energy ($R_{k}$) and thermal variance ($R_{t}$) can range between positive or negative values depending on the combined effects of temporal and spatial discretisation schemes and grid spacing. The positive values of $R_{k}$ and $R_{t}$ indicate a loss of the kinetic energy and thermal variance due to numerical errors, whereas the negative values of $R_{k}$ and $R_{t}$ result in a moderate attenuation of these quantities. Accordingly, the numerical dissipation results obtained here (i.e. $R_{k}$ and $R_{t}$) have been utilised for evaluating the effective values of the Rayleigh and Prandtl number (i.e. $Ra_{eff}$, $Pr_{eff}$) for UDNS. Using the scaling of the mean Nusselt number ($Nu \sim {Ra^a} {Pr^b}$), it has been shown that the error relative to exact mean Nusselt number can be estimated based on $R_{k}$, $R_{t}$ and the effective mean Nusselt number of the UDNS ($Nu_{UDNS}$). Overall, this methodology has been found a promising tool for the quality assessment of UDNS for the applications of non-isothermal flows (e.g. Rayleigh–Bénard convection).

1 Introduction

Turbulent non-isothermal flow is a phenomenon where there is a two-way coupling between fluid flow and heat transfer due to density variations arising from temperature differences. This phenomenon is prevalent in heat exchangers, chemical reactors, combustion and atmospheric flows. Rayleigh–Bénard convection (i.e. fluid between horizontal walls heated from below and cooled from above) is a canonical type of non-isothermal flow which is significantly important for numerous applications in environment and technology. Industrial processes such as chemical and food processing, heating and ventilation of indoor spaces, cooling of electronics, solar and nuclear power systems are few important examples (Bodenschatz et al. 2000).

Direct Numerical Simulation (DNS) of Rayleigh–Bénard convection is not feasible especially at high Rayleigh (i.e. ratio of thermal transport due to buoyancy and viscous forces) and Prandtl (i.e. ratio of momentum and thermal diffusivities) numbers due to high resolution requirements. Therefore, Large-Eddy Simulation (LES) can be used as a reliable tool for investigating turbulent Rayleigh–Bénard convection at high Rayleigh and Prandtl numbers for near future applications. In this work we consider UDNS which can be considered as implicit LES i.e. an LES without explicit subgrid scale model. LES is also increasingly used for the design of combustion systems (Janicka and Sadiki 2005) which represents a different kind of non-isothermal flow. However, LES quality assessment is immensely important in order to obtain predictive LES simulations. Vervisch et al. (2008) state that energy conservation for velocity and scalar fields on unstructured meshes is a major issue for turbulent combustion LES. Due to interaction of numerical and modelling errors, it is much more complex than the verification of solutions from the Reynolds-averaged Navier–Stokes equations. Another problem in the context of implicit filtering is, that the governing equations change, whenever the grid is changed (Klein 2005). An overview of different methodologies to assess the quality of LES simulations can be found in Celik et al. (2009) and Roache (1998). Recently, Schranner et al. (2015) proposed a method on computing the residual of the kinetic-energy evolution in the physical space which can serve as a post-processing tool for the quantification of the numerical errors in Computational Fluid Dynamics (CFD) simulations. According to the taxonomy for obtaining error estimates, as suggested by Roache (1998), this method belongs to the category of “Auxiliary algebraic evaluations on the same grid”, more precisely “Non-conservation of higher order moments”, and in their particular case: turbulent kinetic energy. Komen et al. (2017) proposed a complementary method to Schranner et al. (2015) based on the computation of the residuals in the transport equation for the mean turbulent kinetic energy.

Following the work of Schranner et al. (2015), their methodology has been tested for quantifying the effective numerical dissipation rate in UDNS for the flow around an air foil (Castiglioni and Domaradzki 2015) and a sphere (Cadieux et al. 2017). It has been reported that the methodology can be used as promising utility for estimation of numerical dissipation in turbulent flows (Castiglioni and Domaradzki 2015; Cadieux et al. 2017). However, this methodology has not been tested so far for non-isothermal flows such as turbulent Rayleigh–Bénard convection in enclosures. Therefore, this gap has been addressed here by extending the aforementioned approach for scalar transport (i.e. thermal variance) in addition to kinetic energy since momentum and thermal transport are coupled due to buoyancy forces in the Rayleigh–Bénard convection. Accordingly the main objectives of this study are:

• (i) To analyse the effects of grid resolution and discretisation schemes on the numerical dissipation in UDNS of the turbulent Rayleigh–Bénard convection in a cubical enclosure.

• (ii) To assess the quality of simulations based on the numerical dissipation evaluation.

2 Methodology

The balances of both kinetic energy ($e_{kin}=u_i^2/2$) and thermal variance ($e_t=T^2/2$) have been taken into account as shown in Eqs. 1 and 2 for the buoyancy driven flows of incompressible Newtonian fluids:

$$\begin{aligned} \varepsilon _n^k&= - \iiint _V \Bigg [ {\partial e_{kin} \over \partial t} + \underbrace{\partial (e_{kin}u_i) \over \partial x_i}_{F_{e_{kin}}} + \underbrace{{1 \over \rho }{\partial (pu_i) \over \partial x_i}}_{F_{p}} \nonumber \\&\quad -\underbrace{{\partial \over \partial x_i}\left( \nu {\partial e_{kin} \over \partial x_i}\right) }_{F_{\nu }} +\underbrace{\nu {\partial u_j \over \partial x_i}{\partial u_j \over \partial x_i}}_{\varepsilon _\nu } - \underbrace{u_ig_i(1-\beta (T-T_{ref}))}_{F_{s}} \Bigg ] dV, \end{aligned}$$

(1)

$$\begin{aligned} \varepsilon _n^{t}&= - \iiint _V \Bigg [ {\partial e_t \over \partial t} + \underbrace{\partial (e_t u_i) \over \partial x_i}_{F_{e_{t}}} - \underbrace{{\partial \over \partial x_i} \left( \alpha {\partial e_t \over {\partial x_i}}\right) }_{F_{\alpha }} +\underbrace{\alpha {\partial T \over \partial x_i}{\partial T \over \partial x_i}}_{\varepsilon _\alpha } \Bigg ] dV. \end{aligned}$$

(2)

Conservative transport of momentum (temperature) does not necessarily imply conservation of kinetic energy (thermal variance) (Morinishi et al. 1998). Hence, the residuals of Eqs. 1 and 2, represent according to Cadieux et al. (2017) the numerical dissipation of kinetic energy (i.e. $\varepsilon _n^k$), numerical dissipation of thermal variance (i.e. $\varepsilon _n^{t}$) and $\varepsilon _\nu$, $\varepsilon _\alpha$ denote the physical viscous dissipation and physical thermal dissipation, respectively. These equations are passively evaluated on the fly, using an extended version of the source code provided in Cadieux et al. (2017), based on resolved velocity and temperature fields for different grid resolutions (see Sect. 4 for the numerical details of attaining resolved velocity and temperature fields). All the terms of Eqs. 1 and 2 were computed on the fly (see Fig. 1 a), during the run time of the solver for a cubic domain from sequential snapshots of the velocity and temperature fields. Accordingly, the rate of the change in the kinetic energy and thermal variance is computed using a second-order central difference in time time step index n as follows:

$$\begin{aligned} {\partial e_{\phi } \over \partial t} = \frac{{e_{\phi }}^{n+1} - {e_{\phi }}^{n-1}}{2 \Delta t} \, \end{aligned}$$

(3)

where $e_\phi =e_{kin}$ or $e_\phi =e_t$. Following Cadieux et al. (2017) the difference between the rate of change evaluated based on transported velocity and temperature and the remaining terms (involving all spatial derivatives and $F_{s}$), that are passively evaluated according to Eqs. 1 and 2, without explicitly solving a transport equation for kinetic energy and thermal variance, can be considered as a numerical error and following the original work will henceforth be referred to as numerical dissipation.

Figure 1b shows (note that all terms are shown without leading sign) that the buoyancy term has the biggest contribution to the energy balance and this is countered by viscous dissipation. These are followed in magnitude by the viscous flux and the pressure work. The same qualitative trends have been observed for all simulations and mesh resolutions considered in this work.

3 Non-dimensional Numbers

In the Rayleigh–Bénard convection, the mean Nusselt number ($Nu=hH/k$) characterizes the heat transfer rates at the walls which depends on the non-dimensional parameters such as Rayleigh number ($Ra=g\beta (T_{H}-T_{C})H^3/\nu \alpha$), Prandtl number ($Pr=\nu /\alpha$), and the aspect ratios of the enclosure. In the Rayleigh, Prandtl and Nusselt number definitions, $g, \beta , T_{H}, T_{C}, H, \nu , \alpha , k$ represent the acceleration due to gravity, thermal expansion coefficient, hot wall temperature, cold wall temperature, height of the enclosure, kinematic viscosity, thermal diffusivity and thermal conductivity respectively.

The mean Nusselt number is defined as a dimensionless heat flux averaged over horizontal walls and over time as follows:

$$\begin{aligned} h = |{-k \,(\partial T / \partial y)_{wf} }\,/\,{(T_{wall} - T_{ref})}| \, , \quad {Nu}_{wall}=({\langle {h}\rangle _{(A,t)} \,H})/\,{k} \, , \end{aligned}$$

(4)

where subscript ’wf’ refers to the condition of the fluid in contact with the wall, $T_{wall}$ is the wall temperature and $T_{ref}$ is the appropriate reference temperature, which can be taken to be the temperature of the hot or cold wall respectively. In Eq. 4, the mean value of the Nusselt number on the walls (i.e. $Nu_{wall}$) is obtained by averaging convective heat transfer coefficients ($\langle {h}\rangle _{(A,t)}$) over horizontal active walls and over time here. In the current analysis, $\varepsilon _n^k$ , $\varepsilon _n^{t},$ as well as all other quantities, were evaluated after the simulations reached a quasi-steady state (Yigit et al. 2020), i.e. when the mean values of ${Nu}_{wall}$ for top and bottom walls converge to each other and become time independent. This is not explicitly shown here but the quasi-steady state can be observed in Fig. 1b for the terms of balance of kinetic energy. Cadieux et al. (2017) normalised numerical dissipation with physical dissipation for kinetic energy and the same is done for thermal variance,

$$\begin{aligned}&R_{k}=\frac{{\langle \varepsilon _n^k\rangle }_{(V,t)}}{{\langle \varepsilon _\nu \rangle }_{(V,t)}} \, ,\quad R_{t}=\frac{{\langle \varepsilon _n^{t}\rangle }_{(V,t)}}{{\langle \varepsilon _\alpha \rangle }_{(V,t)}} \, \end{aligned}$$

(5)

in order to evaluate the overall accuracy of the simulations where the average value (i.e. $\langle ...\rangle _{(V,t)}$) is over the volume of the cubical enclosure and time.

Similarly, numerical dissipation results in this work are also presented in terms of $R_{k}$ and $R_{t}$. It is admitted that the evaluation of $R_k$ and $R_t$ is unavoidably affected by interpolation and discretization errors. It is also worth noting that the ratio of resolved to numerical dissipation of kinetic energy does not provide the error itself. However, it has been suggested by several other authors that it is indicative of LES errors and this follows exactly the ideas presented in (Schranner et al. 2015; Castiglioni and Domaradzki 2015; Cadieux et al. 2017). The resolved dissipation rate (of turbulent kinetic energy) is also used in Komen et al. (2017) for normalisation of the numerical dissipation. Also, Geurts and Fröhlich (2002) define a subgrid activity parameter based on turbulent and molecular dissipation. Finally, Celik et al. (2005) follow a similar philosophy.

Moreover, Calzavarini et al. (2015) suggested two alternative definitions of the mean Nusselt number from volume and time averaged viscous (i.e. $\varepsilon _\nu = \nu {(\partial _{i}u_{j})^2}$) and thermal (i.e. $\varepsilon _\alpha = \alpha {(\partial _{i} T)^2}$) dissipation rates in the following manner:

$$\begin{aligned} {Nu}_{\varepsilon _{\nu }}=\frac{{\langle \varepsilon _\nu \rangle _{(V,t)}} \,H^4\, Pr^2}{Ra \, {\nu }^3}+1 \, , \quad {Nu}_{\varepsilon _{\alpha }}=\frac{{\langle \varepsilon _\alpha \rangle _{(V,t)} \,H^2}}{{\alpha \, ({\Delta T})^2}} \, . \end{aligned}$$

(6)

Accordingly, one can numerically compute the mean Nusselt number in three different ways based on Eqs. 4 and 6. For Nusselt numbers sufficiently larger than unity, Eq. 6 suggests that numerical dissipation errors are expected to directly result in mean Nusselt number errors in the following manner:

$$\begin{aligned} {R_{k}}=\frac{\langle \varepsilon _n^k\rangle _{(V,t)}}{\langle \varepsilon _\nu \rangle _{(V,t)}} \approx {\frac{{Nu}_{\varepsilon _n^k}}{{Nu}_{\varepsilon _{\nu }}}} \, , \quad {R_{t}}= \frac{{\langle \varepsilon _n^t\rangle }_{(V,t)}}{{\langle \varepsilon _\alpha \rangle }_{(V,t)}} \approx {\frac{{Nu}_{\varepsilon _n^{t}}}{{Nu}_{\varepsilon _{\alpha }}}} \, , \end{aligned}$$

(7)

where ${Nu}_{\varepsilon _n^k}$ and ${Nu}_{\varepsilon _n^t}$ characterize changes in the mean Nusselt number due to numerical dissipation in the velocity and temperature fields, respectively. Thus, the variation of $R_{k}$ and $R_{t}$ versus ${Nu}_{\varepsilon _{\nu }}$ and ${Nu}_{\varepsilon _{\alpha }}$ will be analysed in this study with the motivation to asses the quality of UDNS.

4 Numerical Implementation

DNS have been performed for Rayleigh–Bénard convection in order to validate UDNS results. The simulations have been carried out by using the finite volume code OpenFOAM for three different Pr values (i.e. 0.1, 1, 10) at a representative value of Ra (i.e. $Ra=10^7$) in a three-dimensional cubical enclosure. Figure 1a demonstrates that the differentially heated horizontal walls of the cubic enclosure ($H=L=W$) are subjected to constant wall temperature boundary conditions. The bottom wall is taken to be at higher temperature than the top wall ($T_{bottom}>T_{top}$) and all other walls are considered to be adiabatic. No-slip and impermeability conditions are specified for all walls. The pressure-velocity coupling has been addressed by the use of the Merged PISO-SIMPLE (PIMPLE) algorithm, which has been optimised for unsteady problems. Convective and diffusive fluxes are evaluated by the second-order upwind and central differencing schemes, respectively. Temporal advancement has been achieved by the second-order Crank-Nicolson scheme with adaptive time-stepping. It has been ensured that the CFL number is always sufficiently below 0.5 so that the underlying physics is captured with sufficient temporal resolution. The convergence criteria were set to $10^{-7}$ for all the relative (scaled) residuals. Based on the grid resolution criteria by Grötzbach (1983), uniform Cartesian grids of up to $250^3$ have been used for the analysis of different values of Pr at $Ra=10^7$. It has been checked a-posteriori that the chosen grid resolution in terms of non-dimensional wall units $y^{+}$ (where $y^{+}=u_{\tau } y / \nu$ with $u_{\tau }=\sqrt{\tau _{w}/\rho }$, $\tau _{w}$ and y being the friction velocity, wall shear stress and wall normal distance of the wall adjacent grid point, respectively) remains 1.12 for $Pr=0.1$ and smaller than unity (i.e. 0.8) for $Pr=1$ at the highest resolution (i.e. $250^3$). Correspondingly, $y^{+}\sqrt{Pr}$ also remains smaller than unity (i.e. 0.8) for $Pr=10$ at the highest resolution.

Fig. 1

a Schematic diagram of the three-dimensional simulation setup, b temporal evolution ($t^*=t\alpha \sqrt{RaPr}/H^2$) of volume averaged terms of Eq. 1 at $Ra=10^7$ and $Pr=1$ for a $40^3$ uniform grid

The standard OpenFOAM solver (i.e. buoyantBoussinesqPimpleFoam) has been modified for this study since it has been widely used in industrial LES. However, it is worth mentioning that the solver has a collocated arrangement of the flow variables, thus the standard PIMPLE corrector loop includes a pressure Poisson equation which is formulated with the fluxes calculated with the predicted velocity. This pressure projection method requires Rhie-Chow type interpolation which can dissipate kinetic energy of the flow (Shashank et al. 2010). Also, in the pressure correction loop the mass flux on time step n is of the form ${{\phi _{PIMPLE}}^n}=\phi _{1}+\phi _{2}$. The first term is the standard mass flux ${\phi _{1}}=u^n \cdot n_{c}\,dS$ where $n_{c}$, dS are the face normal and the face area of the cell. The second term $\phi _{2}$ is a correction term which is specific to OpenFOAM. This correction term is not implemented in the current solver since it is extremely diffusive, as reported by Vourinen et al. (2014).

5 Results and Discussion

One, two and three-dimensional analysis will be reported here for better understanding of numerical dissipation in Rayleigh–Bénard convection. The one-dimensional analysis will provide a basic understanding of the methodology for a generic test case. The two dimensional Rayleigh–Bénard convection analysis serves the purpose of demonstrating the convergence of numerical dissipation towards zero for very fine grids, which is not possible in the three dimensional case due to excessive computational cost. Finally, three-dimensional numerical dissipation results will be presented in detail in the last part of the results section.

5.1 One-Dimensional Analysis

In order to better understand the methodology of evaluating the numerical dissipation on the fly by using a balance equation for kinetic energy, the analysis is started with a simple 1D linear convection-diffusion equation in a periodic domain of length L. Figure 2a shows the initial pseudo turbulent profile used for the convection diffusion equation together with the quasi-exact solution obtained on a very fine grid. Six different mesh resolutions (i.e. $L/ \Delta x=50,100,400,800,1600,3200$) have been used for the analysis. Despite the simplicity of this test case several important features are apparent from the results of this one-dimensional problem.

Fig. 2

1D case: a Initial profile for the linear convection-diffusion equation together with the quasi-exact solution after one revolution. b Relative error in percent between numerical and exact solution versus $R_{k}$ for different spatial discretisation schemes combined with the Euler explicit time integration. c Relative error in percent between numerical and exact solution versus $R_{k}$ for a central differencing scheme combined with different time integration methods

Figure 2b shows the relative error in percent between numerical and exact solution versus $R_{k}$ for different spatial discretisation schemes (first order upwind differencing (UDS), central differencing (CDS), a total variation diminishing (TVD) scheme with CHARM limiter, a third order weighted essentially non oscillatory (WENO) scheme and the quadratic upwind interpolation for convective kinematics (QUICK) scheme) combined with the Euler explicit time integration and a CFL number of 0.1. It becomes obvious for all schemes that the numerical dissipation is well correlated with the relative error according to the maximum norm (i.e. $||S_{qe}- S_{num}||_{\infty } \,/\, ||S_{qe}||_{\infty }$ where $S_{qe}$ and $S_{num}$ represent quasi-exact and numerical solutions, respectively). For small mesh sizes all schemes except UDS converge towards the CDS solution and this convergence is the fastest for the QUICK scheme followed by the TVD and then the WENO scheme. It is also important to note that the numerical dissipation changes its sign (and it is non monotonic) for the WENO and TVD schemes starting from positive values and then switching to negative values for finer meshes. The negative dissipation is induced by the Euler explicit scheme. Although obvious it is worth mentioning the possibility of negative numerical dissipation explicitly. Negative values of $R_{k}$ are also reported in Cadieux et al. (2017), but were not discussed further, and they can locally be seen in the DNS results reported in Komen et al. (2017).

Figure 2c shows the influence of the time discretization scheme (Euler explicit, Euler implicit, Crank Nicolson) for a central differencing scheme used for convective transport. A second order Runge-Kutta method (Heun) has been tested as well, but the results are for small time steps very similar to the Crank Nicolson scheme. The results indicate that the value and sign of $R_{k}$ is a result of the combined effects of the temporal and spatial discretization schemes and grid spacing.

It is important to mention here that the term numerical dissipation, which was not coined by the authors, is to some extent misleading as it suggests positive $R_k$ values. While numerical dissipation typically refers to the convective term (implicit time discretization can also have numerical dissipation), here it is related to the imbalance of a whole equation.

Cadieux et al. (2017) suggest a value of $R_{k}$ smaller than 1% in order to obtain good simulation results. Figure 2 demonstrates that this is only possible for very fine grid resolution, at least with second-order accurate spatial and temporal discretization.

5.2 Two-Dimensional Analysis

After the simple introductory one-dimensional example, two-dimensional Rayleigh Bénard convection will be examined in this section. This is, because two-dimensional Rayleigh Bénard convection features similar physics as three-dimensional Rayleigh Bénard convection but allows for finer meshes.

The variations of $R_{k}$ and $R_{t}$ for the two-dimensional Rayleigh Bénard convection are shown in Fig. 3a, b for different discretisation schemes and grid resolutions at $Pr=1$. Figure 3 shows that $R_{k}$ and $R_{t}$ values remain positive and decrease with decreasing non-dimensional grid size for the first order upwind discretisation scheme. Very high values of numerical dissipation are obtained on the coarsest grid for the first order upwind scheme (${R_{k}=105}$ % and ${R_{t}=27.5}$ %). They approach moderately small values with decreasing non-dimensional grid size reaching ${R_{k}=6.1}$ % and ${R_{t}=4.82}$ % on the very fine grid with ${\Delta x /H}=0.001$ corresponding to a $1000^2$ grid.

Fig. 3

2D case: the variations of a $R_{k}$, b $R_{t}$ values versus non-dimensional grid size (i.e. $\Delta x /H$) for different discretisation schemes at $Pr=1$

For both second order schemes a non-monotonic convergence with a change of sign is observed in particular for $R_{k}$. Figure 3 further demonstrates that the value of $R_{k}$ and $R_{t}$ varies from positive, moderately large values, to negative very small values at the order of 1% on the finest grid. The numerical dissipation of the second order schemes is, as expected, considerably smaller compared to the first order scheme. Furthermore, $R_{t}$ remains negative for all grids for the central differencing applied to the temperature equation. It is worth noting that, for simplicity of notation, the linear interpolation from cell centres to cell faces in OpenFOAM is referred to as central differencing scheme in the current analysis.

The variations of $R_{k}$ and $R_{t}$ for the two-dimensional Rayleigh Bénard convection are shown in Fig. 4a, b for different Pr values and grid resolutions (note that the coarsest grid shown in Fig. 3 is not considered in Fig. 4). It can be seen in Fig. 4 that $R_{k}$ and $R_{t}$ values remain negative (apart from $R_{k}$ value for $Pr=0.1$ for the coarsest grid resolution) for all Pr values considered here. However, these negative values approach to zero with decreasing non-dimensional grid size (i.e. $\Delta x /H$). These results actually confirm that this methodology shows similar trends as Fig. 3 not only for $Pr=1$ but also for $Pr<1$ and $Pr>1$. Positive values of $R_{k}$ and $R_{t}$ indicate a loss of the kinetic energy and thermal variance due to numerical errors, whereas moderate negative values of $R_{k}$ and $R_{t}$ result in a moderate attenuation of these quantities.

Fig. 4

2D case: the variations of a $R_{k}$, b $R_{t}$ values versus non-dimensional grid size (i.e. $\Delta x /H$) for different Pr and central differencing discretisation scheme

In the next section three-dimensional UDNS and DNS results will be presented for mesh sizes of up to $250^3$. The finest grid resolution for the two-dimensional analysis corresponds to a $1000\times 1000$ (i.e. $\Delta x /H=0.001$) uniform mesh. Three dimensional analysis on such grids would be very time consuming and even hardly possible with OpenFOAM due to high memory requirements. Therefore, based on the two-dimensional results, it cannot be expected that $R_{k}$ and $R_{t}$ approach values very close to zero in the three-dimensional case.

5.3 Three-Dimensional Analysis

Finally, three-dimensional DNS and UDNS results will be examined for different grid size and discretisation schemes in this section together with a methodology to assess the quality of UDNS.

The absolute value variations of $R_{k}$, $R_{t}$ are shown in Fig. 5 together with the percentage error in the mean Nusselt numbers (i.e. $e=100\cdot |Nu_{DNS}-Nu_{UDNS}|/ \, Nu_{DNS}$) based on Eqs. 4 and 6 between UDNS and DNS versus non-dimensional grid size (i.e. $\Delta x /H$) for the second order upwind and central differencing schemes at $Pr=1$. It is worth noting that the absolute values of $|R_{k}|$, $|R_{t}|$ are shown here due to their non-monotonic behaviour of sign. Nevertheless, the signed values of $R_{k}$, $R_{t}$ with the mean Nusselt numbers based on Eqs. 4 and 6 are also provided in Tables 1 and 2 for the central differencing and second order upwind schemes, respectively. It has been pointed out earlier that $R_{k}$ and $R_{t}$ values are expected to be proportional to the error in mean Nusselt number (see Eq. 7). However, these equations contain either the numerical dissipation error of kinetic energy or thermal variance, but as the temperature and momentum equation are coupled, the prediction depends on both quantities simultaneously. For this reason Fig. 5 shows the combined contribution (i.e. $0.5 (|R_{k}| + |R_{t}|)$) beside the errors estimated on the individual quantities. It should be noted that evaluating the combined contribution of $|R_{k}|$ and $|R_{t}|$ is an acceptable approach for $Pr=1$ since numerical dissipation from velocity and temperature fields are expected to be equivalent due to the unity ratio of momentum and thermal diffusivity for $Pr=1$. It can be clearly seen from Fig. 5 that the error estimate remains of the same order of magnitude than the percentage error between UDNS and DNS of the mean Nusselt numbers for both the central differencing (Fig. 5a) and second order upwind (Fig. 5b) schemes. Figure  5 further demonstrates that the value of $0.5 (|R_{k}| + |R_{t}|)$ decreases with decreasing non-dimensional grid size while the percentage error between UDNS and DNS of the mean Nusselt numbers also decreases in a similar manner and has a similar magnitude for both the second order central differencing and second order upwind schemes.

Fig. 5

3D case: the variations of absolute values of $|R_{k}|$, $|R_{t}|$ and percentage error in the mean Nusselt numbers (i.e. $e=100\cdot |Nu_{DNS}-Nu_{UDNS}|/ \, Nu_{DNS}$) based on Eqs. 4 and 6 between UDNS and DNS versus non-dimensional grid size (i.e. $\Delta x /H$) for the a central differencing scheme, b second order upwind scheme at $Pr=1$. The values of e are shown with circle, triangle and square pointers and the values of R are shown by the dashed lines

Table 1 The three-dimensional results for a central differencing schemes at $Pr=1$

Table 2 The three-dimensional results for the second order upwind differencing schemes at $Pr=1$

Accordingly, the variations of $R_{k}$ and $R_{t}$ for the three-dimensional Rayleigh Bénard convection in a cubical enclosure are shown in Fig. 6a, b for different Pr values and grid resolutions for the second order upwind discretisation scheme. Fig  6 shows that both $R_{k}$ and $R_{t}$ values converge to zero with decreasing non-dimensional grid size (i.e. $\Delta x /H$). Here, it is also important to emphasize that turbulent flow strengthens with decreasing value of Pr for given value of Ra in the Rayleigh Bénard convection (Petschel et al. 2015). This actually can be explained by the increasing Grashof number (i.e. Gr) with decreasing Pr values for given value of Ra (i.e. $Gr=Ra/Pr$). It is also worth noting that Gr is a dimensionless number which indicates the ratio of buoyancy to viscous forces acting on a fluid. Thus, higher values of Gr signify the relative augmentation of the buoyancy forces in the fluid domain. Correspondingly, the numerical dissipation clearly is higher in the velocity field (i.e. $R_{k}$) for $Pr<1$ than the one for $Pr\ge 1$ as it is noticable from Fig. 6a for the same grid resolution and discretisation scheme. Additionally, Fig. 6b further shows that the magnitude of $R_{t}$ is significantly lower for $Pr>1$ values than the one for $Pr\le 1$. The values of $R_{k}$ and $R_{t}$ can be converted to an alternative quality assessment criterion by estimating effective viscosity and thermal diffusivity here. The effective viscosity ($\nu _{eff}=\nu +\nu _{n}$) of numerical simulations can be determined as the sum of the physical viscosity ($\nu$ is the nominal value used in the calculation) and the numerical viscosity ($\nu _{n}$) due to numerical dissipation. Similarly, the effective thermal diffusivity ($\alpha _{eff}=\alpha +\alpha _{n}$) can also be determined as the sum of physical and numerical thermal diffusivity. This allows to define an effective Rayleigh ($Ra_{eff}$) and Prandtl ($Pr_{eff}$) number in terms of their nominal values (Ra, Pr) as:

Fig. 6

3D case: the variations of a $R_{k}$, b) $R_{t}$ values versus non-dimensional grid size (i.e. $\Delta x /H$) for different Pr and the second order upwind discretisation scheme

$$\begin{aligned} Ra_{eff}&={g\beta \Delta T H^3\over \nu _{eff} \, \alpha _{eff}} = {Ra\over \left( 1+{\nu _n\over \nu }\right) \left( 1+{\alpha _n\over \alpha }\right) } \, , \end{aligned}$$

(8)

$$\begin{aligned} Pr_{eff}&= {\nu _{eff} \over \alpha _{eff}} = Pr {\left( 1+{\nu _n\over \nu }\right) \over \left( 1+{\alpha _n\over \alpha }\right) }. \end{aligned}$$

(9)

The numerical viscosity and thermal diffusivity can be calculated as $\nu _{n}= \nu \,{{{\langle \varepsilon _n^k\rangle }_{(V,t)}}/ {{\langle \varepsilon _\nu \rangle }_{(V,t)}}}= \nu \, {R_{k}}$ and $\alpha _{n}= \alpha \, {{{\langle \varepsilon _n^{t}\rangle }_{(V,t)}}/{{\langle \varepsilon _\alpha \rangle }_{(V,t)}}}=\alpha \, {R_{t}}$ respectively. Accordingly, $Ra_{eff}$ and $Pr_{eff}$ can be written in terms of $R_{k}$ and $R_{t}$ in the following manner:

$$\begin{aligned} Ra_{eff}&=\frac{Ra}{(1+R_{k}) (1+R_{t})} \, , \end{aligned}$$

(10)

$$\begin{aligned} Pr_{eff}&=\frac{Pr(1+R_{k})}{(1+R_{t})}. \end{aligned}$$

(11)

The variations of $R_{k}$, $R_{t}$ values for the three dimensional case are shown in Fig. 7a for different values of non-dimensional grid size (i.e. $\Delta x /H$) and maximum CFL numbers for a central differencing scheme. According to Fig 7a, $R_{t}$ values are insensitive to changes in CFL number whereas $R_{k}$ values decrease with decreasing maximum CFL number, in particular for the coarse grid configurations. Furthermore, Fig. 7b shows the variations of $Ra_{eff}$ and $Pr_{eff}$ values calculated from Eqs. 10 and 11 using the values of $R_{k}$, $R_{t}$ (Fig. 7a for max. CFL=0.5) for different non-dimensional grid sizes. It can be clearly seen from Fig. 7b that $Ra_{eff}$ values are considerably smaller than their nominal value (i.e. Ra) for the coarsest grid configuration where numerical dissipation values are also high (i.e. $R_{k}=20.63$ % and $R_{t}=-5.35$ %). Also, $Ra_{eff}$ increases and then decreases towards its asymptotic value with decreasing non-dimensional grid size and approaches the nominal value where numerical dissipation remains much smaller (i.e. $R_{k}=-3.126$ % and $R_{t}=-3.24$ %). By contrast, $Pr_{eff}$ decreases monotonically with decreasing grid size and approaches its nominal value (i.e. Pr).

Fig. 7

3D case: a the variations of $R_{k}$, $R_{t}$ values for different maximum CFL numbers b The variations of $Ra_{eff}$ and $Pr_{eff}$ values calculated from Eqs. 10 and 11 versus non-dimensional grid size (i.e. $\Delta x /H$) for a central differencing scheme at $Pr=1$

It is well known that the mean Nusselt number scales with a power law of Rayleigh and Prandtl numbers (i.e. $Nu \sim {Ra^a} {Pr^b}$). Dropkin and Somerscales (1965) experimentally analysed turbulent Rayleigh–Bénard convection of Newtonian fluids in a cubic enclosure for $5\times 10^4 \leqslant Ra \leqslant 7.17\times 10^8$ and $0.02 \leqslant Pr \leqslant 11560$. They proposed a mean Nusselt number correlation as $Nu=0.069 \, Ra^{1/3} \, Pr^{0.074}$. Accordingly, $Nu_{eff}$ can be defined based on $Ra_{eff}$ and $Pr_{eff}$ and an error estimate by using the aforementioned quantities can be formulated by comparing ${Nu}_{eff}$ with its nominal value (i.e. Nu) using the mean Nusselt number correlation by Dropkin and Somerscales (1965) as follows:

$$\begin{aligned} Nu={Nu}_{eff}\, [(1+R_{k})^{\frac{1}{3}-0.074} \, (1+R_{t})^{\frac{1}{3}+0.074}] \,. \end{aligned}$$

(12)

Equation 12 can be reformulated to obtain an estimate for the relative error between DNS and UDNS in the following manner where the effective Nusselt number ($Nu_{eff}$) corresponds to the value obtained from UDNS (i.e. $Nu_{eff}=Nu_{\,UDNS}$):

$$\begin{aligned} e_{\,Nu}=\bigg |1-\frac{Nu_{\,UDNS}}{Nu_{\,DNS}}\bigg | \approx \, |1 - [(1+R_{k})^{-\frac{1}{3}+0.074} \, (1+R_{t})^{-\frac{1}{3}-0.074}]| \,. \end{aligned}$$

(13)

Here, $Nu_{UDNS}$ can be taken as one of the values measured directly from UDNS (see Tables 1, 2) i.e. $Nu_{\varepsilon _{\nu }}$ or $Nu_{\varepsilon _{\alpha }}$.

Fig. 8

3D case: the Nusselt numbers obtained from UDNS by using the values of $Nu_{\,UDNS}=Nu_{\varepsilon _{\nu }}$ (blue triangles) and $Nu_{\,UDNS}=Nu_{\varepsilon _{\alpha }}$ are (red circles) versus non-dimensional grid size (i.e. $\Delta x /H$) for different Pr values a $Pr=0.1$, b $Pr=1.0$, c $Pr=10$ for the second order upwind discretisation scheme. The estimated errors based on Eq. 13 are shown with the error bars and the Nusselt number from DNS is shown by dashed lines. The values are slightly shifted along ${\Delta x /H}$ for better visibility

Figure 8 shows the Nusselt numbers obtained from UDNS together with the error estimate obtained from Eq. 13 for different non-dimensional grid size (i.e. $\Delta x /H$) and Pr values. It can be seen from Fig. 8 that Eq. 13 estimates the error quite promisingly based on $Nu_{\,UDNS}=Nu_{\varepsilon _{\alpha }}$ for $Pr<1$, whereas the estimation of Eq. 13 performs better based on $Nu_{\,UDNS}=Nu_{\varepsilon _{\nu }}$ for $Pr>1$. In both cases the statement holds true for the determination of $Nu_{\,UDNS}$ itself not only for the error. For $Pr=1$, the estimation based on both approaches works reasonably well apart from the coarsest grid size (i.e. ${\Delta x /H}=0.0125$), and the deviation significantly decreases with decreasing grid size based on both $Nu_{\,UDNS}=Nu_{\varepsilon _{\alpha }}$ and $Nu_{\,UDNS}=Nu_{\varepsilon _{\nu }}$.

The previous findings (Petschel et al. 2013, 2015; Zhang et al. 2017) on the analysis of kinetic and thermal energy dissipation rates in turbulent Rayleigh–Bénard convection agree that the dissipation rate contributions (i.e. ${\langle \varepsilon _\nu \rangle }_{(V,t)}$ and ${\langle \varepsilon _\alpha \rangle }_{(V,t)}$) that come from the bulk and boundary layer regions are dominated by the boundary layer regions. Therefore, the estimation of ${\langle \varepsilon _\nu \rangle }_{(V,t)}$ (${\langle \varepsilon _\alpha \rangle }_{(V,t)}$) on the boundary layer regions is expected to be more critical for $Pr<1$ ($Pr>1$) due to the inadequate resolution for UDNS analysis. Consequently, the mean Nusselt number can be evaluated more conveniently for UDNS analysis based on $Nu_{\,UDNS}=Nu_{\varepsilon _{\alpha }}$ for $Pr<1$ and based on $Nu_{\,UDNS}=Nu_{\varepsilon _{\nu }}$ for $Pr>1$ as shown in the results of Fig. 8. While error / uncertainty analysis is a very intricate topic and results have always to be treated with care, the findings shown in Fig. 8 are encouraging for future LES quality assesment of turbulent non-isothermal flows.

6 Conclusions

In this study, determination of numerical errors in UDNS of Rayleigh–Bénard convection in enclosed spaces has been analysed by extending the method used by Cadieux et al. (2017). It has been shown that the balance of kinetic energy and thermal variance methodology promisingly works for the quality assessment of UDNS of Rayleigh–Bénard convection. Accordingly, major findings of this study are summarised as follows:

• (i) The numerical dissipation for both kinetic energy ($R_{k}$) and thermal variance ($R_{t}$) can assume positive or negative values depending on the combined effects of temporal and spatial discretisation schemes and grid spacing. Positive values of $R_{k}$ and $R_{t}$ indicate a loss of the kinetic energy and thermal variance due to numerical errors, whereas the negative values of $R_{k}$ and $R_{t}$ result in a moderate attenuation of these quantities.

• (ii) The values of $R_{k}$ and $R_{t}$ remain of the same order of magnitude as the percentage error of the mean Nusselt numbers between UDNS and DNS for both second order discretization schemes used in this work. Both values decrease with decreasing non-dimensional grid size together with the percentage error of the mean Nusselt number between UDNS and DNS.

• (iii) The values of $R_{k}$ and $R_{t}$ can be utilised for evaluating the effective values of the Rayleigh and Prandtl number (i.e. $Ra_{eff}$, $Pr_{eff}$) in UDNS. Therefore, using the scaling of the mean Nusselt number ($Nu \sim {Ra^a} {Pr^b}$), the error between the exact mean Nusselt number of the DNS ($Nu_{\,DNS}$) and the effective mean Nusselt number of the UDNS ($Nu_{\,UDNS}$) can be estimated conveniently for different values of Pr.

• (iv) It has been observed that the Nusselt number definition based on thermal dissipation (viscous dissipation) provides more accurate results for predicting the Nusselt number and its relative error in the case of $Pr<1$ ($Pr>1$).

Finally, the same methodology framework can be extended to explicit LES in the future, to account for additional dissipation from an explicit SGS model beside numerical errors due to grid resolution and discretisation schemes. Further, the methodology will be applied to reactive non-isothermal flows.