An Explanation to the Nusselt–Rayleigh Discrepancy in Naturally Convected Porous Media

Abstract

We model hydrothermal convection using a partial differential equation formed by Darcy velocity and temperature—the velocity formulation. Using the Elder problem as a benchmark, we found that the velocity formulation is a valid model of hydrothermal convection. By performing simulations with Rayleigh numbers in the non-oscillatory regime, we show that multiple quasi-steady-state solutions can be one of the reasons that caused the Nusselt–Rayleigh discrepancy found in previous experiments. The results reveal more understandings about the nature of uncertainty of convection modes in porous media.

1 Introduction

Natural convection in porous media occurs in various geological and industrial settings, such as groundwater, geothermal reservoirs, heat sinks and thermal energy storage. We use the Rayleigh number (Ra) to characterize the material properties of porous media, such as permeability, fluid density and the temperature difference between the top and bottom boundaries. We measure the quality of convective heat transfer using the Nusselt number (Nu). From an engineering perspective, it is beneficial to correlate the Rayleigh number and the Nusselt number.

Cheng (1979) compiled the experimental, analytical and numerical results of the Nusselt number and the Rayleigh number for convection heat transfer in a porous layer heated from below. The compilation showed widespread Nusselt numbers for a particular Rayleigh number. We refer to this phenomenon as the Nu–Ra discrepancy. Experimental and analytical approaches have addressed this problem. Lister (1990) performed experiments using a large porous slab and showed that lateral thermal dispersion is one of the reasons that caused the Nu–Ra discrepancy. Vadasz (2010) and Vadasz and Braester (1992) used analytical techniques to show that boundary and domain imperfections are the causes of the Nu–Ra discrepancy. Karani and Huber (2017b) conducted pore-scale lattice-Boltzmann simulations and concluded that thermal disequilibrium between solid and fluid phases could cause the Nu–Ra discrepancy. Furthermore, Karani et al. (2017) proposed a fractional-order thermal convection model that captures not only the Nu–Ra discrepancy but also the advance and delay of the onset of convection. The authors explain this phenomenon by examining different effects that can be introduced by experimental setups. Therefore, we present numerical simulations, which can be set up as experiments without the imperfections mentioned above.

In a three-dimensional box of saturated porous media, convection can happen in a two or three-dimensional setting (Beck 1972). Straus and Schubert (1978) found out that two-dimensional flows have a larger Nusselt number compared to three-dimensional flows when Ra $\le 97$ using numerical simulations. The numerical simulations from Holst and Aziz (1972) and Horne (1979) also produced this effect. Therefore, if both two and three-dimensional convection can exist in the same box, then the system can at least have two Nusselt numbers. More promisingly, the simulations of Straus and Schubert (1979) indicate that it is always possible to force either steady two-dimensional or steady three-dimensional convection by proper choice of initial conditions. Govorukhin and Shevchenko (2017) also showed that the selection scenarios of a convection mode strongly depend on the initial temperature distribution of the porous media. The cosymmetry theory developed by Yudovich (1991) showed that there exists an infinite number of stable stationary flow for a fixed Rayleigh number.

Therefore, we hypothesize that the non-uniqueness of convection modes in a box of saturated porous media is one of the reasons that caused the Nu–Ra discrepancy. To prove this hypothesis, we perform several numerical simulations using different box sizes and initial conditions of the temperature field. We use the finite element solver—MOOSE Framework (Gaston et al. 2009) to simulate the physics of naturally convected porous media.

Acknowledging that a specific natural or artificial porous media can have multiple heat transfer rates is not too useful for decisions in reservoir engineering. Karani and Huber (2017a) used the basin stability analysis (Menck et al. 2013) to quantify the probability of the occurrence of different convection modes in a two-dimensional setting. We employ the same method for a three-dimensional setting. This probability of convection modes can provide us more information and obtain the expectation value of the Nusselt number.

2 Materials and Methods

2.1 Governing Equations

We present the conservation laws that model natural convection in porous media under the following assumptions mentioned by Horne (1979):

• The Boussinesq approximation.

• Inertial effects are small, or low Reynolds number.

• Viscosity of the fluid is constant.

• Thermal dispersion is negligible.

• Saturating fluid and porous solid are in thermal equilibrium.

$$\begin{aligned} \nabla \cdot \mathbf {q}= & {} 0, \nonumber \\ \nabla ^2 u + \sqrt{\text {Ra}}\frac{\partial ^2 T}{\partial x \partial y}= & {} 0, \nonumber \\ \nabla ^2 v - \sqrt{\text {Ra}}\left( \frac{ \partial ^2 T}{\partial x^2} + \frac{\partial ^2 T}{\partial z^2} \right)= & {} 0, \nonumber \\ \nabla ^2 w + \sqrt{\text {Ra}}\frac{\partial ^2 T}{\partial y \partial z}= & {} 0, \nonumber \\ T_t + \sqrt{\text {Ra}} \left( \mathbf {q}\cdot \nabla T \right) - \nabla ^2 T= & {} 0, \end{aligned}$$

(1)

where $\mathbf {q} = (u,v,w)$, the non-dimensional Darcy fluxes in x, y and z directions. T is non-dimensional temperature, the subscript t is the partial derivatives with respect to time, and Ra is the Rayleigh number. Gravity effects act on the vertical direction y. The Rayleigh number is defined as

$$\begin{aligned} \text {Ra} = \frac{\rho _0^2 g \alpha K_p d \varDelta T c_f}{\mu k_m}, \end{aligned}$$

where $\rho _0$ is the reference density of the fluid, g is the gravity, $\alpha$ is the thermal expansion rate of fluid, $K_p$ is the permeability, d is the vertical height of the porous media, $\varDelta T$ is the temperature difference between the top and bottom boundaries, $\mu$ is the dynamic viscosity and $k_m$ is the overall thermal conductivity. This formulation of velocity and temperature is also used by Florio (2014).

The non-dimensional variables are

$$\begin{aligned} T= & {} T_0 + (\varDelta T)T^*, \nonumber \\ \rho _f= & {} \rho _0(1-\alpha (\varDelta T) T^*), \nonumber \\ \mathbf {q}= & {} [q]\mathbf {q}^*, \nonumber \\ (x, y, z)= & {} d(x^*, y^*, z^*), \nonumber \\ t= & {} [t] t^*, \end{aligned}$$

(2)

where $T_0$ is the reference temperature at the bottom boundary, $\rho _f$ is the fluid density. The non-dimensional variables are defined with asterisks, and we drop them for convenience throughout this paper. The scalings are chosen such that the square root of the Rayleigh number exists in both momentum and energy conservation equations

$$\begin{aligned}&\left[ q\right] = \sqrt{\frac{K B g k_m}{\mu c_f d}}, \nonumber \\&\left[ t\right] = \frac{\phi \rho _0 c_f d^2}{k_m}. \end{aligned}$$

(3)

Since we are using a finite element solver, we require the weak form:

$$\begin{aligned}&\int _{\varOmega } \left( \nabla \omega \cdot \nabla u + \sqrt{\text {Ra}}\frac{\partial \omega }{\partial y} \frac{\partial T}{\partial x}\right) d\varOmega = 0, \nonumber \\&\int _{\varOmega } \left( \nabla \omega \cdot \nabla v - \sqrt{\text {Ra}} \left( \frac{\partial \omega }{\partial x} \frac{\partial T}{\partial x} + \frac{\partial \omega }{\partial z}\frac{\partial T}{\partial z} \right) \right) d\varOmega = 0, \nonumber \\&\int _{\varOmega } \left( \nabla \omega \cdot \nabla w + \sqrt{\text {Ra}}\frac{\partial \omega }{\partial y} \frac{\partial T}{\partial z} \right) d\varOmega = 0, \nonumber \\&\int _{\varOmega } \left( \omega T_t + \omega \sqrt{\text {Ra}} \mathbf {q} \cdot \nabla T + \nabla \omega \cdot \nabla T \right) d\varOmega = 0. \end{aligned}$$

(4)

The weak form is derived by multiplying the strong form by a test function $\omega$ and integrate through the simulation domain. We choose a first order Lagrangian basis function for velocity $\mathbf {q}$ and second order Lagrangian basis function for temperature T. The time integration method of the transient solver is the Crank–Nicolson method. We approximate the CFL number on an element as

$$\begin{aligned} C_{elem} = \varDelta t \frac{|u| + |v| + |w|}{h_{min}}, \end{aligned}$$

where $h_{min}$ is the minimum length of the element w.r.t x, y and z axes. In general, we pick the maximum CFL number of all elements as the representative CFL number, and we make sure during each timestep, the representative CFL number is around 0.5. When the entropy production is constant with respect to time, we claim that the simulation has reached a steady state (Börsing et al. 2017).

2.2 Benchmark Problems and Results

The benchmark results of our code implementation are presented in this subsection. They are essential in this work, but do not serve as our main focus.

2.2.1 The Elder Problem

Elder studied convection caused by localized heating with the Hele-Shaw cell experiment and numerical solutions (Elder 1967). It is one of the benchmark problems in software such as FEFLOW (Trefry and Muffels 2007; Diersch 2014), SUTRA (Voss 1984) and HydroGeoSphere (Brunner and Simmons 2012; Simmons and Elder 2017). Therefore, we benchmark our finite element implementation of the velocity formulation using the Elder problem.

The velocity formulation in 2D is

$$\begin{aligned} \nabla ^2 u + \sqrt{\text {Ra}}\frac{\partial ^2 T}{\partial x \partial y}= & {} 0, \\ \nabla ^2 v - \sqrt{\text {Ra}}\frac{\partial ^2 T}{\partial x^2}= & {} 0, \\ T_t + \sqrt{\text {Ra}} \mathbf {q} \cdot \nabla T - \nabla ^2 T= & {} 0. \end{aligned}$$

The boundary and initial conditions are defined in Fig. 1. Using the physical parameters of the Elder problem defined by Graf and Boufadel (2011), we have a non-dimensionalized problem with Ra $= 521.3$. We use a structured quadrilateral mesh, with the 120 elements in the x direction and 32 elements in the y direction. We choose a uniform timestep $\varDelta t=0.0001$. We compare our solution with the 2, 5 and 10 years simulation results of Graf and Boufadel (2011), which corresponds to 0.01, 0.025 and 0.05 non-dimensionalized time. See Fig. 2 for the benchmark results (Table 1).

Fig. 1

The boundary and initial conditions of the non-dimensionalized Elder problem

Fig. 2

Comparison of numerical solutions of the Elder problem using the velocity formulation. Left: reference solution of Graf and Boufadel (2011). Right: our solution

Table 1 The physical properties of the porous media. Modified from Graf and Boufadel (2011). The thermal expansion rate is calculated using the density of water as a function of temperature from Thiesen et al. (1900), assuming a linear growth of volume

2.2.2 Beck’s Box

For a 3D box of saturated porous media of a given dimension, Beck (1972) derived the preferred cellular mode during the onset of convection. We benchmark our code implementations using different boxes with lengths $h_1$ and widths $h_2$. We use the notation of [$h_1$, $h_2$] to represent the scaled box dimensions normalized by the vertical height, as illustrated in Fig. 4. The boxes we choose are [1.2, 1.2], [2.0, 0.5] and [3.0, 1.0]. The Rayleigh number is set as 42.25, slightly above the critical Rayleigh number, such that the system starts convecting. The initial condition of temperature is the conductive solution with a $\pm 1 \%$ perturbation. We use the notation of (m, n) to represent the cellular modes. The results are in Table 2, and they agree with Beck’s analytical cellular modes.

Table 2 Comparison of Beck’s cellular mode and the simulated cellular mode.

2.3 Quality Measures of Convective Heat Transfer

The Nusselt number and entropy production are used to measure the quality of convective heat transfer in our simulations. The Nusselt number is defined as the ratio of total heat transfer and the stationary conductive heat transfer. Due to simplicity and practical reasons, the Nusselt number is widely used in experiments as a quality measure of convection. However, in numerical simulations, we have access to values of temperature and Darcy fluxes in space. Therefore, we can use another quality measurement—entropy production. Entropy production is generally a better assessment of convection, due to its thermodynamic considerations (Herwig 2016). The Nusselt number is a combination of the quality, given by the entropy production, and the quantity of heat transfer, given by the heat flux involved (Herwig 2016). Though, in order to revisit the Nusselt–Rayleigh relation and gain new insights from it, our simulations still contain the calculation of the Nusselt number

$$\begin{aligned} \text {Nu} = \frac{1}{h_1h_2} \int _0^{h_2}\int _0^{h_1} -\frac{\partial T}{\partial y} \Bigr |_{y=1} dxdz, \end{aligned}$$

(5)

where $h_1$ and $h_2$ are the length of the box with respect to x and z-direction. This is analogous to the calculations of Hewitt et al. (2014), except that we utilize heat flux on the top boundary.

Bejan (2013) formulates the volumetric rate of the total entropy production in a fluid-saturated porous medium as a result of both irreversible heat transfer (subscript therm) and fluid flow friction (subscript visc)

$$\begin{aligned} \frac{dS^{''}_{total}}{dt} = \dot{S}^{''}_{total} = \dot{S}^{''}_{therm} + \dot{S}^{''}_{visc}, \end{aligned}$$

(6)

with $\dot{S}^{''}$ ($\text {J}\,\text {s}^{-1}\,\text {m}^{-3}\,\text {K}^{-1}$). The first term on the right-hand side of Eq. (6) represents the entropy production due to heat transfer irreversibility and invokes the rate of heat flow per unit area and unit time, i.e., Fourier’s law

$$\begin{aligned} \dot{S}^{''}_{therm} = \frac{k_m}{\bar{T}^2}(\nabla T)^2, \text {with } \bar{T} = \frac{T_{top}+T_{bot}}{2}. \end{aligned}$$

(7)

$T_{top}$ and $T_{bot}$ are the temperature on top and bottom boundaries, respectively. Recall non-dimensional variables Eq. (2)

$$\begin{aligned}&T = T_0 + (\varDelta T)T^*, \\&\nabla = \frac{1}{d}\nabla ^*, \\&(\nabla T)^2 = \frac{\varDelta T^2}{d^2} (\nabla ^* T^*)^2. \end{aligned}$$

Substitute $(\nabla T)^2$ in Eq. (7)

$$\begin{aligned} \dot{S}^{''}_{therm} = \frac{k_m \varDelta T^2}{\bar{T}^2 d^2}(\nabla ^* T^*)^2. \end{aligned}$$

(8)

The second term on the right-hand side of Eq. (6) accounts for viscous dissipation effects of the fluid

$$\begin{aligned} \dot{S}^{''}_{visc} = \frac{\mu _f}{K\bar{T}}\mathbf {q}^2 + \frac{\mu _f}{\bar{T}} \varPhi , \end{aligned}$$

(9)

where $\mathbf {q}$ and $\varPhi$ are Darcy flux and the viscous dissipation function, respectively. The second term on the right-hand side of Eq. (9) is only important when the flow tends to behave like non-Darcy flow (Costa 2006) and can be neglected in our study. In general, the viscous dissipation effects become increasingly important for a heterogeneous medium as preferred fluid flow paths lead to a local increase in the fluid velocity. Take the scaling factor of $\mathbf {q}$ in Eq. (3) and plug in Eq. (9)

$$\begin{aligned} \dot{S}^{''}_{visc} = \frac{\alpha \varDelta T g k_m}{\bar{T} c_f d}(\mathbf {q}^*)^2. \end{aligned}$$

We integrate entropy production over the computational domain, assuming homogeneous material properties

$$\begin{aligned} \dot{S}_{total} = \int _{V} \dot{S}^{''}_{total} dV = \frac{k_m (\varDelta T)^2}{\bar{T}^2 d^2}\int _V (\nabla ^* T^*)^2dV + \frac{\alpha \varDelta Tgk_m}{\bar{T} c_f d} \int _V (\mathbf {q}^*)^2 dV, \end{aligned}$$

(10)

where V is the volume of the computational domain. Dropping the asterisks, the dimensionless entropy production, or the entropy generation number is

$$\begin{aligned} N_{\dot{S}} = \frac{d^2 \bar{T}^2}{k_m (\varDelta T)^2} \frac{\dot{S}}{V} = \frac{1}{V}\left( \int _V (\nabla T)^2dV + \frac{\alpha d \bar{T} g}{\varDelta T c_f} \int _V (\mathbf {q})^2 dV \right) . \end{aligned}$$

(11)

We would like to emphasize that the thermal entropy production measures the norm of the linear and nonlinear parts of the temperature gradient, and the Nusselt number measures the flux from nonlinear parts of the temperature from the top boundary. Since the entropy production and the Nusselt number have different meanings in quality measure in convective heat transfer, we calculate both values for all of the simulations.

2.4 Problem Formulation

Börsing et al. (2017) investigated the entropy production of a naturally convecting porous media of various aspect ratios and Rayleigh numbers in a 2D setting. We extended the analysis to a 3D setting and designed numerical experiments over boxes of different dimensions. Figure 3 shows a scatter plot of the box dimensions $[h_1, h_2]$ categorized into the triangular test and the line test. The triangular test aims to inform the entropy production over a wide range of cellular modes. The line test focuses on how entropy production relates to both the total number of convection cells and the cellular modes’ dimensions. Throughout the simulations, the Rayleigh number is set to 42.25, slightly above the critical Rayleigh number. The characteristic length of the mesh is 0.05.

The initial condition of the transient problem is the conductive solution with a $\pm 1 \%$ perturbation. For each point in Fig. 3, we only realize one transient simulation. This is certainly not ideal, as the steady-state solution of the transient problem depends strongly on the initial condition. To compensate, we further analyze the line test by forcing steady-state solutions using the initial conditions

$$\begin{aligned} T= & {} A \sin {(\pi y)} \cos {\left( \frac{m \pi x}{h_1}\right) } \cos {\left( \frac{n \pi z}{h_2}\right) } + 1-y, \nonumber \\ u= & {} -A\frac{m \sqrt{\text {Ra}}}{h_1}\left( \frac{m^2}{h_1^2} + \frac{n^2}{h_2^2} + 1\right) ^{-1} \cos {\left( \pi y\right) }\sin {\left( \frac{m \pi x}{h_1}\right) }\cos {\left( \frac{n \pi z}{h_2}\right) }, \nonumber \\ v= & {} A\sqrt{\text {Ra}}\left( \frac{m^2}{h_1^2} + \frac{n^2}{h_2^2}\right) \left( \frac{m^2}{h_1^2} + \frac{n^2}{h_2^2} + 1\right) ^{-1} \sin {(\pi y)} \cos {\left( \frac{m \pi x}{h_1}\right) } \cos {\left( \frac{n \pi z}{h_2}\right) }, \nonumber \\ w= & {} -A\frac{n \sqrt{\text {Ra}}}{h_2}\left( \frac{m^2}{h_1^2} + \frac{n^2}{h_2^2} + 1\right) ^{-1} \cos {(\pi y)} \cos {\left( \frac{m \pi x}{h_1}\right) } \sin {\left( \frac{n \pi z}{h_2}\right) }, \end{aligned}$$

(12)

where A is the amplitude (Florio 2014). We set the amplitude such that the entropy production does not exceed $\frac{9}{256}$Ra, which is the analytical upper bound of the Nu–Ra relationship (Doering and Constantin 1998).

The observations of varying Nusselt number with respect to 2D or 3D cellular modes (Holst and Aziz 1972; Straus and Schubert 1978; Horne 1979) can be a reason for the Nu–Ra discrepancy. We test this idea by simulating three boxes of different sizes that have 2D or 3D preferred cellular modes during the onset of convection. The long boxes [1.5, 1.0] and [2.3, 0.9] convect with a 2D cellular mode, and the wide boxes [2.5, 1.5] convect with a 3D one. The three boxes are tested in the regions of $4\pi ^2<$ Ra $\le 196$, and the type of convection cells, the Nusselt number and the entropy production are reported. The characteristic length of the mesh is set to 0.1. The goal is to see how the transition of cellular modes due to increasing Rayleigh number influences the Nusselt number, and whether it can explain the Nu–Ra discrepancy.

Fig. 3

Preferred cellular mode, (m, n), as a function of $h_1$, $h_2$. The figure is symmetric with respect to the line $h_1 = h_2$. The change between the modes (m, 0) and $(m+1, 0)$ occurs at $h_1 = \sqrt{m(m+1)}$ (Beck 1972). We examine the convection pattern and entropy production of the plotted points using the proposed finite element solver

Fig. 4

Temperature boundary conditions and side length $(h_1, h_2)$ of Beck’s box. The box is closed with impermeable boundaries

2.5 Basin Stability Analysis Using the Equivalent Entropy Production Initialization

Florio (2014) investigated the probability of convection modes in “critical boxes” using perturbation methods. The critical boxes are boxes with the size that lies on the transition between several convection modes. We investigate this probability using the basin stability analysis similar to Karani and Huber (2017a). Instead of exploring the amplitude space from -1 to 1, we consider the symmetry of amplitudes, and we initialize the amplitudes in the space that is bounded by a certain entropy production.

We pick one of the boxes Florio (2014) used in his analysis $[2^{1/4}, 2^{1/4}]$, approximated as [1.189, 1.189]. The possible convection modes of this box when the Rayleigh number is critical are (0, 1), (1, 0) and (1, 1). The critical Rayleigh number is $\pi ^2(2+3/\sqrt{2}) \approx 40.68$. We set the Rayleigh number to 42.25 to be slightly above the critical Rayleigh number. We define the initial condition and the amplitude of the modes (0, 1) and (1, 1)

$$\begin{aligned}&T_{01} = A_{01} \sin {(\pi y)} \cos {\left( \frac{\pi z}{2^{1/4}}\right) }, \\&T_{11} = A_{11} \sin {(\pi y)} \cos {\left( \frac{\pi x}{2^{1/4}}\right) } \cos {\left( \frac{\pi z}{2^{1/4}}\right) }. \end{aligned}$$

Apply the thermal part of nondimensional entropy production Eq. (11) to $T_{01}$ and $T_{11}$, and sum the results

$$\begin{aligned} N_{\dot{S}\text {sum}}= N_{\dot{S}01} + N_{\dot{S}11} = \frac{\pi ^2(\sqrt{2}+1)}{4\sqrt{2}} (A_{01})^2 + \frac{\pi ^2(\sqrt{2}+2)}{8\sqrt{2}} (A_{11})^2. \end{aligned}$$

Consider the amplitudes $A_{01}$ and $A_{11}$ as axes in the Cartesian coordinates, the combinations of amplitudes that have the same amount of entropy production $N_{\dot{S}\text {sum}}$ are an ellipse

$$\begin{aligned} \frac{(A_{01})^2}{a^2} + \frac{(A_{11})^2}{b^2} = 1, a^2 = \frac{4 \sqrt{2} N_{\dot{S}\text {sum}}}{\pi ^2(\sqrt{2}+1)}, b^2 = \frac{8 \sqrt{2} N_{\dot{S}\text {sum}}}{\pi ^2(\sqrt{2}+2)}, \end{aligned}$$

where a is the semi-major axis and b is the semi-minor axis. We apply a change of variables

$$\begin{aligned}&A_{01} = a \cos {\xi }, A_{11} = b \sin {\xi }, \\&\xi = \tan ^{-1}{\left( \frac{a}{b} \tan {\theta }\right) }, \end{aligned}$$

such that we can characterize the combination of amplitudes with the same amount of entropy production using the angle $\theta$, and the initial entropy production $N_{\dot{S}\text {sum}}$. The initial conditions of velocities u, v, w are the same as Eq. (12). The probability of the (0, 1) mode to occur is defined as the area that is characterized by $\theta =0$ to the separation angle $\theta = \theta ^*$ divided by quarter of the elliptical area

$$\begin{aligned} P(\theta ) = \frac{2}{\pi } \tan ^{-1}{\left( \frac{a}{b} \tan {\theta } \right) } = \frac{2 \xi }{\pi }. \end{aligned}$$

(13)

We divide $\xi$ linearly into 17 parts from 0 to 90 degrees. The entropy productions $N_{\dot{S}\text {sum}}$ we choose to initialize are 1.02, 1.04, 1.08, 1.12 and 1.16. Note that we cannot perform such analysis using the Nusselt number. However, the linear solutions of the convective temperature Eq. (12) have zero contribution to the Nusselt number.

3 Results

3.1 Quality of Convective Heat Transfer in Different Settings of Boxes

Figure 5 shows the cellular mode and the entropy production of the line test with Rayleigh number of 42.25. Figure 6 represents the cellular mode and entropy production of the triangular test with Ra $= 42.25$.

Fig. 5

The line test at Ra=42.25 plotted over Beck’s diagram (Beck 1972). The y-axis of the middle and bottom plots represents the entropy production. The middle plot is solved by the transient solver with the randomly perturbed initial conditions. The bottom plot is solved by the steady-state solver with the initial conditions Eq. (12). The bracketed values (m, n) represent the cellular mode

Fig. 6

The triangular test at Ra=42.25 plotted over Beck’s diagram (Beck 1972). The bottom plot represents the contour of entropy production using linear interpolation. The bracketed values (m, n) represent the cellular mode

3.2 The Nu–Ra Relation

We test the Nusselt number of the boxes [1.5, 1.0], [2.3, 0.9] and [2.5, 1.5] with respect to the Rayleigh numbers in the range between 42.25 and 196. Figures 7 and 8 are the results.

3.3 The Probability of Mode Occurrence Using Basin Stability Analysis

Figure 9 shows the result of the basin stability analysis. We found out that the separation angle of the box $[2^{1/4}, 2^{1/4}]$ between modes (0, 1) and (1, 1) using the basin stability analysis. The separation angle lies between ${44.30}^{\circ }$ and ${49.94}^{\circ }$. We average the two angles and claim the separation angle is ${47.12}^{\circ }$. Using Eq. (13), the probability of mode (0, 1) to occur when the initial condition is the superposition of mode (0, 1) and (1, 1) is 46.8 percent. However, consider the symmetry of modes (0, 1) and (1, 0), we can calculate the probability of mode (1, 1) to occur when the initial condition is the superposition of the three modes, which is 36.2 percent.

Fig. 7

The relation between the Rayleigh number and the entropy production of the boxes [1.5, 1.0], [2.3, 0.9] and [2.5, 1.5]. The bracketed values (m, n) represent the cellular mode

Fig. 8

The Nu–Ra relation of the boxes [1.5, 1.0], [2.3, 0.9] and [2.5, 1.5], plotted over Cheng (1979)’s compilation

Fig. 9

The diagram of basin stability analysis. The red line is drawn as $y = \tan {(\theta ^* x)}$, where $\theta ^*$ is the separation angle. The streamlines show how different initial conditions reach the steady-state solution. The bracketed values (m, n) represent the cellular mode at steady state

4 Discussion

4.1 Quality of Convective Heat Transfer in Different Settings of Boxes

Figure 5 shows that the cellular modes of our simulation match Beck’s prediction well. The slight shift of the boundaries of cellular modes change can be attributed to that our simulations are performed on a porous media with Ra $=42.25>\text {Ra}_c$. We color-coded the cellular modes such that shades of blues and greens are for two-dimensional cellular modes, and gradients of reds and yellows are for three-dimensional cellular modes. By only looking at the entropy production of the blues and greens, the patterns are similar to Börsing et al. (2017)’s tests of 2D boxes. It is clear that box sizes are relevant to the quality of convective heat transfer. Focusing on the whole figure including the three-dimensional cellular modes, the observations are:

• 1. The quality of convective heat transfer of three-dimensional cellular modes is generally worse than those of two-dimensional cellular modes when Ra$=42.25$.

• 2. The entropy production gradually increases from three dimensional cellular modes to two-dimensional cellular modes that have the same total sum of cells $m+n$.

• 3. The transient solution does not necessary converge to cellular modes with better quality of convective heat transfer. See Sect. 4.3 for further discussions about the influence of the initial conditions.

Figure 6 further supports the aforementioned observations. On the boundaries from mode (2, 0) to either (2, 1) or (1, 1), we can also observe a descent of the entropy production. Increasing convection cells do not guarantee a better quality of convective heat transfer. These results give us new insights into the relation between cellular modes and entropy production. We can also view these results as more rigorous benchmarks that compare with theoretical predictions.

4.2 The Nu–Ra Relation

Figure 8 compared the Nu–Ra relation of the three boxes with previous experimental and numerical results. The results do not show the wide scattering of Nusselt number during the onset of convection nor in the region of Ra $\le 100$. This can be attributed to the lack of thermal dispersion effects in our model (Karani et al. 2017).

However, the regions of $140<$Ra$<196$ show a wide scattering of the Nusselt number from 3.212 to 4.145, which can be explained by the multiple steady states of convection pattern. Our results agree with Straus and Schubert (1979)’s numerical tests that both two and three-dimensional convection cells could be obtained at $60 \le$Ra$\le 150$ by perturbing the initial condition of temperature. The results also agree with Borkowska-Pawlak and Kordylewski (1985)’s proof, that continuous transition of pattern flows from two-dimensional to three-dimensional structure is possible (and vice-versa) with Rayleigh number variations.

The hypotheses of box sizes and multiple steady states of convection pattern influencing the Nu–Ra relation are proved using the three boxes experiment. Future experiments of the Nu–Ra relation should also be aware of the chosen box size and the convection pattern of fluids.

4.3 The Probability of Mode Occurrence and its Implications

Our results suggest that we cannot infer the steady-state pattern, given the entropy production or the Nusselt number of the initial condition. The steady-state solution is determined by the initial combination of mode amplitudes. Thus, we investigate the probability of mode occurrence, assuming equal chance for amplitude combinations to occur as initial conditions.

Using the basin stability analysis with the equivalent entropy production initialization, we can calculate the probability of mode occurrence. In our example of the $[2^{1/4}, 2^{1/4}]$ box, the probability of mode (1, 1) to occur is 36.2 percent, which is slightly higher than the modes (0, 1) and (1, 0). From the simulation, we also know that the Nusselt number for mode (0, 1) and (1, 1) is 1.0637 and 1.0546, respectively. Combining with their occurring probability, we can, therefore, calculate the expectation value for the Nusselt number, which is 1.0579. We have provided a method for a new Nu–Ra relationship in the non-oscillating region of the Rayleigh number.

The assumption that we can use a straight line to separate the probability space is only for convenience. It is also possible that the line can be of a higher order. As in Fig. 9, the straight line does not separate the modes. Nevertheless, in this example, it already gives a good approximation and demonstrates how we can apply this method to different box sizes and Rayleigh numbers.

5 Conclusion

We show how the influence of different box sizes and multiple steady states of convection pattern leads to the discrepancy of Nu–Ra relation in the region of moderate Rayleigh number. We also demonstrate the method of basin stability analysis using equivalent entropy production initialization to study the probability of mode occurrences in naturally convected porous media. This method can be utilized for further studies of the Nu–Ra relationship.

Appendix: Methods of Cellular Mode Checking

In the appendix section, we discuss the methods we used for checking cellular modes of a certain convection pattern.

1.1 Counting by Visual Inspection

We can count how many cells a convection pattern has by looking at the temperature profile at $y=0.5$. The method is: count how many times of a temperature change and a 0.5 (-) temperature, in either x or z-direction. See Fig. 10 for better comprehension of this concept. This process is programmed and used to post-process the simulation results in Fig. 5, 6, 7 and 8. The problem with this method is that it cannot differentiate between the 3D cellular modes (bottom left of Fig. 10) and the linear combination of the 2D cellular modes (bottom right of Fig. 10).

1.2 Counting by Fourier Analysis

Our finite element code calculates the temperature in space T(x, y, z) for a given time. The numerical solution consists of the linear solutions Eq. (12); therefore, we can apply the concept of Fourier analysis. We define the convection temperature as the full solution without the conduction solution

$$\begin{aligned} T_{convect}= & {} T(x,y,z) - (1-y) \\= & {} \sum _{m,n=0}^{\infty } A_{mn} T_{mn}, \end{aligned}$$

note that m and n cannot be 0 at the same time. If we want to know if the system is convecting by a certain cellular mode $(m=M,n=N)$, we multiply the convection temperature by the linear solution $T_{MN}$ and integrate through the domain $\varOmega$

$$\begin{aligned} \int _\varOmega T_{MN} T_{convect} d\varOmega= & {} \sum _{m,n=0}^{\infty } \int _\varOmega A_{mn} T_{MN} T_{mn} d\varOmega \\= & {} A_{MN} \int _\varOmega (T_{MN})^2 d\varOmega . \end{aligned}$$

The amplitude $A_{MN}$ is therefore

$$\begin{aligned} A_{MN} = \frac{\int _\varOmega T_{MN} T_{convect} d\varOmega }{\int _\varOmega (T_{MN})^2 d\varOmega }. \end{aligned}$$

If $M \ne 0$ and $N \ne 0$,

$$\begin{aligned} A_{MN} = \frac{8}{h_1 h_2} \int _\varOmega T_{MN} T_{convect} d\varOmega . \end{aligned}$$

If either of M or N is 0,

$$\begin{aligned} A_{MN} = \frac{4}{h_1 h_2} \int _\varOmega T_{MN} T_{convect} d\varOmega . \end{aligned}$$

This Fourier analysis not only gives us the information of which cellular mode the system has, but also the value of amplitudes. The value of amplitudes can also be used to differentiate the 3D cellular modes from the linear combination of the 2D cellular modes. We use this method to analyze Fig. 9.

Fig. 10

Examples of automatic cellular mode checking using different configurations of boxes. The title of each plot represents the cellular mode (m, n)